Kinematic detection of a planet carving a gap in a protoplanetary disc

TL;DR

This paper reports the direct kinematic detection of a Jupiter-mass planet in a protoplanetary disc, linking the observed dust gap and disturbed gas flow to planet-disc interactions, advancing understanding of planet formation.

Contribution

It provides the first direct kinematic evidence of a planet carving a gap in a protoplanetary disc, confirming the planet's role in creating observed disc substructures.

Findings

Detection of a Jupiter-mass planet at 130 au in HD 97048

Correlation between gas flow disturbances and dust gap

Evidence that some disc gaps are caused by planets

Abstract

We still do not understand how planets form, or why extra-solar planetary systems are so different from our own solar system. But the last few years have dramatically changed our view of the discs of gas and dust around young stars. Observations with the Atacama Large Millimeter/submillimeter Array (ALMA) and extreme adaptive-optics systems have revealed that most --- if not all --- discs contain substructure, including rings and gaps, spirals, azimuthal dust concentrations, and shadows cast by misaligned inner discs. These features have been interpreted as signatures of newborn protoplanets, but the exact origin is unknown. Here we report the kinematic detection of a few Jupiter-mass planet located in a gas and dust gap at 130 au in the disc surrounding the young star HD 97048. An embedded planet can explain both the disturbed Keplerian flow of the gas, detected in CO lines, and the gap detected in the dust disc at the same radius. While gaps appear to be a common feature in protoplanetary discs, we present a direct correspondence between a planet and a dust gap, indicating that at least some gaps are the result of planet-disc interactions.

Tables

Table 1. Table 1: Details of the observations.

UTC Date	Time on source	Nant	Baselines	pwv	Calibrators:
	[minutes]		[m]	[mm]	Bandpass	Flux	Phase
2016 Nov 18	26.3	45	15 to 919	0.60	J0538-4405	J1107-4449	J1058-8003
2017 Nov 24	42.7	49	92 to 8548	0.56	J0522-3627	J0904-5735	J1058-8003
2017 Nov 30	42.7	47	79 to 8283	0.56	J0635-7516	J0904-5735	J1058-8003

Full text

{affiliations}

Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy, Monash University, Clayton Vic 3800, Australia,

Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

Atacama Large Millimeter/Submillimeter Array, Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura 763-0355, Santiago, Chile

Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904,1098XH Amsterdam, The Netherlands

Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

Astronomy Department, University of California, Berkeley, CA 94720-3411, USA

Universidad de Santiago de Chile, Av. Libertador Bernardo O’Higgins 3363, Estación Central, Santiago

Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

A variety of mechanisms have been proposed to explain the formation of rings and gaps in discs, e.g. snow-lines, non-ideal MHD effects, zonal flows, and self-induced dust-traps [9, 10, 11, 12, 13]. The most straightforward explanation is that the gaps are the results of forming planets interacting with the disc [14]. Recent ALMA surveys suggest that planets could indeed be responsible for carving out several of the observed gaps [2, 3, 15], but until now definite evidence has remained elusive. Despite much effort, direct imaging of planets in young discs remains difficult. Many of the claimed detections have been refuted, or require confirmation [16, 17, 18]. The most promising detection to date is a companion imaged in the cleared inner disc around PDS 70 [19, 20, 21, 22]. However, the mass estimate from photometry remains uncertain, and it is not yet clear if PDS 70 b falls within the planetary regime.

Our approach is to search for the dynamical effect of a planet on the surrounding gas disc. Disc kinematics are dominated by Keplerian rotation. Embedded planets perturb the gas flow in their vicinity, launching spiral waves at Lindblad resonances both inside and outside their orbit. The disturbed velocity pattern is detectable by high spectral and spatial resolution ALMA line observations [23]. This technique was used to detect embedded planets in the disc surrounding HD 163296 [24, 25].

Here, we used ALMA to observe the disc surrounding the young ($\approx$3 Myr) intermediate-mass (2.4 M⊙) star HD 97048 in Band 7 continuum (885 $\mu$m) and in the 13CO J=3-2 transition, with a spectral resolution of 220 m s*-1*. Observations were performed with 3 interferometer configurations sampling baselines from 15 to $\approx$ 8,500 m, resulting in final images with a spatial resolution of 0.07”$\times$0.11” (13$\times$20 au).

We report the detection of a localised deviation from Keplerian flow in the disc. The velocity kink is spatially associated with the gap seen in dust continuum emission (Fig. 1 and S1). The most plausible and simultaneous explanation for those two independent features is the presence of an embedded body of a few Jupiter masses which carves a gap in the dust disc and locally perturbs the gas flow.

The continuum emission shows a system of two rings detected up to $\approx$ 1” from the star. The 13CO emission extends further in radius ($\approx$ 4”), displaying the typical butterfly pattern of a disc in Keplerian rotation (Fig. 1 and 2). No significant brightness variation of the 13CO emission is detected at the location of the gap. In a given spectral channel, the emission is distributed along the corresponding isovelocity curve, i.e. the region of the disc where the projected velocity towards the observer is equal to the channel velocity. The observed East-West asymmetry is characteristic of an optically thick emitting layer located above the midplane. The lower, fainter, disc surface is also detected to the West of the upper disc surface.

The CO emission displays a kink in the upper isovelocity curve, highlighted by the dotted circle in Fig. 1. The velocity kink is detected consistently in channels between +0.7 and +1.1 km/s from the systemic velocity. It is also seen in images reconstructed from individual observing nights, i.e. before combining the data sets. The morphology of the emission around the velocity kink is the same with and without continuum subtraction, indicating that the kink is not the result of optical depth effects (see Supplementary Figure 3). The sensitivity of the ALMA observations allows us to detect the continuum in each individual channel, revealing that the velocity kink is located just above the gap seen in continuum emission, at the same radius. This spatial coincidence points to a common origin for both features. The deformation of the emission is localised to a diameter of $\approx$ 0.3”. Notably, the emission on the opposite side the disc (and at opposite velocity) displays a smooth profile, with no kink. This excludes a large scale perturbation of the disc or an azimuthally symmetric mechanism. The perturbation is similar to the one detected in HD 163296 [24]. In both cases, the kink is only detected over a small range in both radial extent and velocity. A corresponding velocity kink in the lower surface of the disc is not seen as the emission is weaker and masked by the continuum and brighter upper CO surface (Fig. 1 and 2). Using the same procedure as in ref.[26], we measured the altitude of the 13CO layer to be 17$\pm$1 au near the velocity kink (at a distance of 130 au). Assuming that the planet is located in the disc midplane and exactly below the center of the velocity kink, it would be at a projected distance of 0.45$\pm$0.1” and PA = $-55\pm$10∘ from the star.

To infer the mass of the putative planet, we performed a series of 3D global gas and multi-grain dust hydrodynamics simulations, where we embedded a planet on a circular orbit at 130 au with a mass of 1, 2, 3 and 5 MJup (gas disc mass of $10^{-2}$ M⊙). Simulations were performed for approximately 800 orbits ($\approx$ 1 Myr), and then post-processed to compute the thermal structure and resulting continuum emission and CO maps.

The presence of a few Jupiter mass planet produces distinct signatures in the gas and dust (Fig. 3). The embedded planet generates a gap and spirals in the gas, resulting in a non-axisymmetric velocity field. The dynamics of the dust depends on the Stokes number, i.e. the ratio of the gas drag stopping time to the orbital time, which depends on the grain size and dust properties. When the Stokes number is close to unity — corresponding roughly to millimetre sized grains at the gas surface densities considered here if grains are compact and spherical — dust grains form axisymmetric rings inside and outside of the planet orbital radius [14].

Fig 4 shows the predicted emission for the various planet masses, in the continuum and for the 13CO line. The channel maps are best reproduced with an embedded planet of 2-3 MJup, giving a velocity kink with amplitude matching the observations. For the 1 MJup planet, the kink is too small. The most massive planet, with 5 MJup, creates a kink too large, and which remains detectable over a range of velocity that is too wide ($\pm$ 1 km/s from the 0.96 km/s channel where the deviation is the strongest). Embedded planets have also been predicted to generate vertical bulk motions and turbulence which should result in detectable line broadening when the planet is massive enough ($>$ a few Jupiter masses) [27]. Analysis of the moment-2 map does not reveal significant line broadening at the location of the gap and are consistent with thermal broadening and Keplerian shearing within the beam. This also rules out the upper end of the range tested in our simulations. HD 97048 was observed with SPHERE on the VLT, resulting in a point source detection limit of $\approx$2 MJup [28] at the location where we detect the velocity kink. This upper limit assumes a hot-start model, and an unattenuated planet atmosphere. Our simulations show that the planet is embedded, with an optical depth of $\gtrsim$0.5 towards the observer at 1.6 $\mu$m, i.e. it would appear about twice as faint as an unobscured planet. This is consistent with a 2-3 MJup planet not being detected by SPHERE. Our kinematic mass determination is also consistent with the planet mass range 0.4 and 4 MJup estimated from the width of the scattered light gap, for a viscosity between $10^{-4}$ and $10^{-2}$ [29].

All of the planet masses explored in our models result in azimuthally symmetric gaps in continuum emission at 885 $\mu$m, as detected by ALMA. At this wavelength, the thermal emission is dominated by dust grains a few hundred microns in size, that decouple from the gas, and form axisymmetric rings, even if the gas flow is locally non-axisymmetric.

The width and/or depth of a gap in sub-millimetre thermal emission depends on the planet mass as well as on the Stokes number of the dust grains that contribute most at the observed wavelength [15]. In most cases, the Stokes number is unknown as the local gas density and dust properties are poorly constrained by observations. Continuum gap width may not provide reliable estimate of planet masses [30]. Conversely, when the mass of the planet is known, for instance via the kinematics of the gas as in this work, the dust continuum observations can be used to directly measure the Stokes number, thus constraining the gas density and dust properties in the vicinity of the gap. For HD 97048, dust grains dominating the thermal emission must have a Stokes number around a few $10^{-2}$ to reproduce the observed gap profile. This excludes compact dust grains with the assumed density in our models. One way to reduce the Stokes number is to increase the gas density, but this causes significant accretion on the planet (see supplementary material). Another possibility is that dust grains of a few 100 microns or millimetres in size consist of fluffy aggregates, as suggested by sub-millimetre polarisation studies [31, 32]. Aggregates have a larger projected area, and experience stronger gas drag than equal mass compact grains. They have smaller Stokes number, and can reproduce the observed dust continuum gap width for a 2 MJup planet (Fig. 4 and S4).

The coincident location of the velocity kink and gap demonstrates that protoplanets are responsible for at least some of the observed gaps in discs. Most of the alternative mechanisms for creating dust gaps in discs — including snow lines, non-ideal MHD, zonal flows and self-induced dust traps — rely on the formation of a pressure bump where the dust grains can be trapped and grow further. While those pressure bumps produce deviations from Keplerian velocity, they are axisymmetric. That is, they do not cause a localised, non-axisymmetric velocity deviation as observed. Other mechanisms might be imagined to create a non-azimuthally symmetric velocity pattern in the disc. Gravitational instabilities or outer companion/flyby create spirals, but these are large-scale structures, and would not produce a velocity kink localised to a small region of the disc. Neither will they result in azimuthally symmetric dust gaps. The interaction of a few Jupiter masses planet with its surrounding disc is to our knowledge the only plausible explanation that can explain both a localised velocity kink and an azimuthally symmetric gap. More systematic kinematic mass estimates may allow us to better connect the population of young embedded planets in discs with the known exoplanet population.

0.1 HD 97048.

HD 97048 is located in the Chameleon I cloud, at a distance of 185 pc [33]. It has a spectral type Be9.5/AO, an effective temperature of 10,000 K, and a luminosity of 40 solar luminosities. The star is surrounded by a large disc [34, 35], extending up to at least 850 au in 12CO J=2-1 emission which is seen at an inclination $\approx 40^{\circ}$. Previous ALMA observations of the dust sub-disc (probing dust grains a few hundred microns in size) revealed a central cavity and two bright rings separated by a gap centered at $\approx$130 au from the star [35]. The two bright rings are also detected in scattered light with VLT/SPHERE [28] (Supplementary figure 1), where the observations probes the distribution of sub-micron sized dust grains. Such small grains experience a high gas drag and tend to follow closely the gas spatial distribution, suggesting that the gap is also present in the gas disc. Two extra rings extending up to 2.2” from the star are also detected in scattered light. The disc is detected in PAH emission up to $\approx$ 650 au, revealing a flaring surface [34]. Scattered light images confirmed that the disc surface shows a significant flaring [28].

0.2 Observations and data reduction.

We observed HD 97048 with ALMA in band 7 in the C40-4 (1 execution) and C40-7 (2 executions) configurations reaching a total time on source of 112 minutes (ALMA program #2016.1.00826.S, PI: G. van der Plas). The details of the observations can be found in Table 1. One of the spectral windows for each observation was centered at the 13CO J=3-2 rest frequency with an individual channel width of 122 kHz, resulting in a 244 kHz spectral resolution (220 m s*-1*) after Hanning smoothing. The three other spectral windows were used for continuum with a bandwidth of 1.875 GHz each.

We performed one round of phase self calibration on the continuum data set observed November 24th, and two rounds on the other two data sets, and applied the self calibration solutions to the line data. We imaged the visibilities at a 120 m s*-1* velocity spacing using Briggs weighing, resulting in a beam size of 0.11” $\times$ 0.07” arcsec at a PA of -38∘. Typical RMS values obtained are 3.6 mJy/beam per channel for the combined data. Channel maps are presented in Supplementary Figure 2.

The kink in CO emission discussed in this manuscript is present in all 3 individual executions and in both the continuum subtracted and non-subtracted images (Supplementary Figure 3).

0.3 3D modelling procedure.

We performed a series of 3D global simulations using the phantom Smoothed Particle Hydrodynamics (SPH) code [36]. We performed a series of multi-grain gas+dust simulations using the algorithm described in [37, 38], using 2 million SPH particles and following the dust fraction of particles of sizes ranging from 1.5625 to 1600 $\mu$m (11 bins total, each bin doubling the grain size). Each dust species experiences a different gas drag depending on the grain size, resulting in differential vertical setting and radial migration. All 11 populations of dust grains were evolved simultaneously with the gas. We thus self-consistently took account of the cumulative backreaction on the gas.

We assume a central mass of 2.4 M⊙ and a distance of 185 pc [33]. We set the initial disc inner and outer radii to 40 au and 700 au, respectively. We set the gas mass to $10^{-2}$ M⊙, and use an exponentially tapered power-law surface density profile with a critical radius of 500 au, power-law index of $-0.5$. The disc aspect ratio was set to 0.06 at 40 au (consistent with the scattered light images [28]), with a vertically isothermal equation of state, and sound speed power-law index of $-0.25$. We set the artificial viscosity in the code to obtain an average Shakura-Sunyaev[39] viscosity of $10^{-3}$ [40].

We embedded two planets in the disc orbiting at 30 au, with a mass of 1 MJup, and at 130 au with a mass of either 1, 2, 3, or 5 MJup. The inner planet is used to carve a central cavity in the disc, as seen in the ALMA observations, but we did not vary the mass of the inner planet in this work. The presence of an inner planet is also suggested by the torqued intra-cavity HCO+ velocity field [35]. We used sink particles [41] to represent the star and planets. We set the accretion radius of the planets to 0.125 times the Hill radius, with an accretion radius of 10 au for the central star. The model surface density is plotted in the top left panel of figure 3 for the 2 MJup planet, along with the radial velocity (top centre) and predicted deviation from Keplerian flow (top right panel) and dust column density for 1.56, 6.25 and 200$\mu$m grains (bottom row). We evolved the models for 800 orbits of the outer planet ($\approx 1$ million years). The flow pattern around the planet establishes itself over a much shorter timescale, but the gap carving is set by the viscous timescale and establishment of dust gaps at lower Stokes numbers requires a large number of orbits.

The planets accrete a moderate amount of gas from the disc. The final masses after 800 orbits are 1.38, 2.52, 3.63 and 5.77 MJup for the 1, 2, 3, 5 MJup planets respectively. Migration is negligible with all planets migrating by less than 1 au by the end of the calculations.

Additional simulations were also performed with a disc gas mass of $10^{-1}$ M⊙. As they result in significant accretion on the outer planet, we explore a range of planet masses. Planets with initial masses of 0.1, 0.15, 0.2, 0.25, and 0.5 MJup reach a mass of 0.11, 0.18, 0.30, 2.0, 3.6, and 5.9 MJup, respectively, after 800 orbits. Migration remains limited to less than 3 au for all planet masses, except for the most massive planet which migrated by 9 au. We can produce a velocity kink matching the observations using a planet with a final mass of 2 MJup, giving us a similar planet mass estimate as with the lower disc gas mass models. However, in the high gas mass models the accretion rate on the planet increases with time, leading to a runaway accretion process and a high final planet mass.

To compute the disc thermal structure, continuum images and synthetic line maps, we used the mcfost Monte Carlo radiative transfer code [42, 43], assuming $T_{\mathrm{gas}}=T_{\mathrm{dust}}$, and local thermodynamic equilibrium as we are looking at low-$J$ CO lines. The central star was represented by a sphere of radius 2.25 R⊙, radiating isotropically with a Kurucz spectrum at 10,000 K [44]. To avoid interpolating the density structure between the SPH and radiative transfer code, we used a Voronoi tesselation where each mcfost cell corresponds to a SPH particle. We set the 13CO abundance to 7 10*-7* and follow the prescription described in Appendix B of [26] to account for freeze-out where T $<20$ K, photo-dissociation and photo-desorption in locations where the UV radiation is high. We adopted a turbulent velocity of 50 m s*-1*.

We adopted a fixed dust mixture composed of 60 % silicate, 15 % amorphous carbon [45]. Each grain size is represented by a distribution of hollow spheres with a maximum void fraction of $0.8$. We use a grain population with 100 logarithmic bins in size ranging from 0.03 to 3000$\mu$m. At each point in the model, the density of a given grain size is obtained by interpolating between the SPH dust grain sizes, assuming grains smaller than half the smallest SPH grain size, ie  0.78$\mu$m, follow the gas distribution, and grains larger than 1.6 mm follow the distribution of the 1.6 mm grains. The total grain size distribution is normalised by integrating over all grain sizes, assuming a power-law $\mathrm{d}n(a)\propto a^{-3.5}\mathrm{d}a$, and over all cells, where we set the total dust mass to 1/100 of the total SPH gas mass. We computed the dust optical properties using the Mie theory.

CO maps were generated at a spectral resolution of 50 m/s, binned at the observed resolution and Hanning smoothed to match the observed spectral resolution. Continuum and CO maps were then convolved with a Gaussian beam matching the ALMA CLEAN beam.

The spatial distribution of the dust grains is set by their Stokes number. In phantom, we assume compact spheres, but grains with different shape and mass will follow the same spatial distribution if they have the same Stokes number. To study the impact of the Stokes number on the thermal emission maps, we can mimic fluffyness by shifting the SPH grain sizes before interpolating the grain size distribution in mcfost. Thermal emission at 885 $\mu$m is dominated by dust grains a few hundred microns in size. For our model with a 10*-2* M⊙ gas mass, compact grains of this size would have a Stokes number close to 1. The corresponding gas width appears larger than in the observations. Supplementary Figure 4 shows that if the emitting grains have a Stokes number around $10^{-2}$ instead of $1$, for instance if they are fluffy aggregates, the sub-millimetre gap and ring widths are in much better agreement with the observations for the same mass planet.

0.4 Radiative-equilibrium hydrodynamics calculations.

The 13CO emission originates from between 1 and 2 hydrostatic scale heights[26], at an altitude where the temperature is higher than in the midplane. To test the validity of the vertically isothermal structure used in the SPH calculations, we also performed a set of phantom simulations where the temperature is regularly updated by mcfost. The two codes have been interfaced to run simultaneously. Thanks to the fast mapping between the distribution of SPH particles and radiative transfer Voronoi mesh, we can perform frequent radiative transfer calculations within the SPH simulation. At a specified time interval, phantom passes the local density, grain distribution, and sink particles properties to mcfost which returns the 3D disc temperature (assuming the gas temperature is equal to the dust temperature). Between two mcfost calls, phantom evolves the disc with the temperature of each particle held constant. This main advantage of this method is that we include the full frequency dependence, as well as light scattering, which are critical for accurate temperature calculations in protoplanetary discs. We assume radiative equilibrium at each call of mcfost, which is only a valid approximation if the radiative timescale is much smaller than the dynamical timescale. Due to the limited optical depths of the models of HD 97048, these conditions are satisfied here. The temperature structure is updated every $1/10^{th}$ of the outer planet orbit. Simulations were also performed with calls to mcfost once per orbit, producing almost indistinguishable results.

Supplementary figures 5, 6 and 7 compares the velocity fields for the 2 MJup model between the vertically isothermal and vertically stratified (ie with regular mcfost temperature updates) cases. The isothermal simulation was designed to have a similar midplane temperature as the mcfost+phantom model, i.e. with $h/r=0.06$ at $r=40$ au. This is the same simulation as presented in the rest of the paper. As expected, differences increase with altitude, with larger deviations in the vertically stratified case. In this model, they remain limited to $\approx$0.1 and 0.05 km/s for the azimuthal and radial velocities, respectively.

In the case of the 13CO emission of HD 97048, the vertically isothermal structure appears as a reasonable approximation, as long as the $h/r$ is chosen sensibly. Hence vertical temperature stratification does not significantly affect our planet mass estimate.

0.5 Impact of observational noise and $uv$ plane sampling.

As we do not aim to perform a detailed fitting of the data, all models so far were presented with a simple Gaussian convolution to compare with observations, i.e. showing noise-free and with fully-sampled $uv$-plane synthetic maps.

To assess if observational artefacts could affect the images and in particular the detection of the kink, we also post-processed the 2 MJup model through a modified version of the CASA ALMA simulator. Synthetic visibilities were computed at the same ($u$,$v$) coordinates as the data. A precipitable water vapor of 0.6 mm was used to set the thermal noise. The resulting synthetic visibilities were CLEANed using the same parameters as observed visibilities.

A comparison of the Gaussian-convolved and CLEANed synthetic images is showed in supplementary figure 8. The shape of the velocity kink is not affected by observational effects, indicating that a simple convolution is a good approximation for qualitatively comparing models to data.

Data availability

Raw data is publicly available via the ALMA archive under project id 2016.1.00826.S. Final reduced and calibrated data cubes are available with the DOI 10.6084/m9.figshare.8266988.

Code availability

Phantom is publicly available at https://bitbucket.org/danielprice/phantom. mcfost is currently available under request and will be made open-source soon. Figures were generated with splash[46] (http://users.monash.edu.au/~dprice/splash/) and pymcfost (https://github.com/cpinte/pymcfost), which are both open-source.