PLANET FORMATION SIGNPOSTS: OBSERVABILITY OF CIRCUMPLANETARY DISKS VIA GAS KINEMATICS

Planets are expected to form during the evolution of circumstellar disks of gas and dust (e.g., Armitage 2011). As a protoplanetary core grows to approximately a Saturn mass, it becomes massive enough to open a gap in the disk (Lin & Papaloizou 1986; Lubow et al. 1999). This process is the result of competition between gravitational, viscous, and pressure torques exerted onto the disk by the planet and by the disk itself (Crida et al. 2006). The gap splits the disk into two radially distinct zones. There are multiple detections of dust-depleted gaps and cavities in well-studied protoplanetary disks (Andrews et al. 2011; Perez et al. 2015), but interestingly, no unambiguous detection of a forming planet has yet been confirmed.

Hydrodynamical models of planet–disk interactions show that a single giant protoplanet continues accreting from the outer disk at formidable rates (∼2 × 10⁻⁴ M_Jup yr⁻¹; Gressel et al. 2013; Shabram & Boley 2013; Szulágyi et al. 2014) even after its gap is evacuated (Kley 1999; Papaloizou & Nelson 2005). Accretion streams converge onto the vicinity of the giant developing a circumplanetary disk (CPD), through which angular momentum disperses, thus regulating planetary accretion (Lubow et al. 1999; Ayliffe & Bate 2009).

In two-dimensional simulations, strong shocks appear near the planet's Hill sphere leading to excessive inflow redirected toward the planet, rapidly depleting the CPD material (Lubow et al. 1999; D'Angelo et al. 2002). These shocks appear much weaker in three-dimensional (3D) calculations, leading to more persistent CPD structures (Bate et al. 2003; Ayliffe & Bate 2009). Grid-based simulations show that the inflow of gas onto the protoplanet mostly happens in the vertical direction, allowing material to cross the shock front near the Hill radius (Machida et al. 2010; Szulágyi et al. 2014). Similar results were found for MHD simulations by Gressel et al. (2013).

Most models show that CPD radii truncate to about one-third (Ayliffe & Bate 2009; Shabram & Boley 2013) or even one-half (Gressel et al. 2013; Szulágyi et al. 2014) of the planet's Hill radius. This implies that a Jupiter-mass planet on a 100 AU orbit could bear a CPD with a diameter of 4.5–7 AU. If located at 100 pc away from Earth, the projected diameter translates into ∼45–70 mas, within the range of modern astronomical instrumentation.

Gas-giant protoplanet candidates have been detected embedded in the HD 100546 and HD 169142 disks, two in each system (Quanz et al. 2013a; Biller et al. 2014; Currie et al. 2014; Reggiani et al. 2014). Both systems are young and nearby Herbig Ae/Be stars (<10 Myr old, <145 pc) bearing large circumstellar disks with confirmed gaps (Osorio et al. 2014; Walsh et al. 2014). The two directly imaged detections are L' (3.8 μm) emission blobs: one is at 05 angular separation (∼52 AU at 100 pc) from HD 100546 (also detected in M' emission at 4.8 μm, Quanz et al. 2014) and the other is at 016 (∼23 AU at 145 pc) from HD 169142. The compact source in HD 169142 falls within a symmetric gap imaged in polarized scattered light by Quanz et al. (2013b). Interestingly, there are no near-IR counterparts for either candidate, supporting the idea that these are in fact accreting gas giants with SEDs driven by circumplanetary accretion (Zhu 2015).

These CPDs have eluded unambiguous detection, mainly because of the lack of predictions on CPD tracers. Isella et al. (2014) presented deep continuum observations to detect a CPD around a protoplanet candidate in LkCa 15, with no positive results. CPD continuum emission may in fact be scant since large dust grains (>100 μm), probed by submillimeter continuum observations, are sieved out of the planet-induced gap by the outer disk pressure bump. Interestingly, only small grains (<100 μm) enter the dust-depleted cavity (Zhu et al. 2012) and make it into the CPD. This dust filtration is evidenced by the numerous dust-depleted cavities seen in submillimeter observations (Andrews et al. 2011). CPD observational signposts to unambiguously confirm forming planet candidates are scant. Realistic 3D simulations coupled with radiative transfer are needed to predict the observational signatures of a CPD at various wavelengths.

In this work, we connect 3D smooth particle hydrodynamics (SPH) of CPD dynamics (Section 2) with 3D radiative transfer (Section 3) in the context of interferometric observations of line emission. We aim to study the observability of CPDs as features imprinted in the gas kinematics (Section 4). We compute predictions tailored for the Atacama Large Millimeter/submillimeter Array (ALMA) for common gas tracers (Section 4.1). Implications and conclusions are discussed in Section 5.

We carried out a set of 3D SPH simulations to characterize the CPD morphology and kinematics, which we feed into a radiative transfer code (see Section 3). We address the question of observability signposts of a single snap CPD embedded in its parent disk after the simulations have reached a reasonably steady state. Exploring how planet–disk interactions evolve with time is beyond the scope of this paper.

2.1. 3D SPH Simulations

We used a modified version of the SPH code Gadget-2 (Springel 2005) to perform a suite of 3D simulations of a planet embedded in a protoplanetary disk. The code has been modified to make it more suitable for simulating planets embedded in disks (see Dunhill et al. 2013). We have simulated the full disk azimuth but with a restricted radial range.

We implement radial boundary conditions in a manner similar to Ayliffe & Bate (2010), reducing the range of the simulated disk while still providing high resolution around the planet. Although spiral waves launched by the planet reflect off the boundary, they do not affect the disk dynamics at the radius of the planet. A more complete description of this method will be given in a forthcoming paper (A. Dunhill et al. 2015, in preparation). We summarize used disk parameters in Table 1.

Table 1.  Simulations Parameters

Parameter	Value
${R}_{\mathrm{in}}$	0.35 Rp
${R}_{\mathrm{out}}$	1.85 Rp
Initial Σ profile	${\rm{\Sigma }}\propto {R}^{-1}$ a
Initial $H/R$ profile	0.05 (constant with R)
Circumstellar disk $\mathrm{masses}$ b
Isothermal (SPH1, SPH2)	$5\times {10}^{-4\;}{M}_{\odot }$
Adiabatic (SPH3)	$7\times {10}^{-5\;}{M}_{\odot }$

Notes.

^aNormalized so that R_p = 1 AU would give ${\rm{\Sigma }}=100\;{\rm{g}}\;{\mathrm{cm}}^{-2}$ at R = R_p. ^bFor the nominal disk model presented in Section 4.

Download table as:  ASCIITypeset image

To reduce the effect of transient waves on the planet at the start of the simulation, we include an initial gap in the disk using the prescription of Lubow & D'Angelo (2006). We model the planet as a point mass potential with a sink radius, inside which any gas particles are swallowed and their mass and momentum are added to the planet. The potential is unsoftened outside the sink radius. Initially, R_sink = 0.05 R_p and decays exponentially to R_sink = 0.001 R_p after approximately eight orbits. This corresponds to ∼2 Jupiter radii for R_p = 1 AU.

We performed three different runs, two locally isothermal (SPH1 and SPH2, where T(R) is enforced so that $H/R$ remains constant) and one with an adiabatic equation of state (SPH3, using an adiabatic index $\gamma =5/3$). It has been shown before that using an isothermal equation of state yields slightly larger CPDs than using more realistic radiation physics (Ayliffe & Bate 2009). SPH1 has a planet mass M_p = 1 M_J, while SPH2 and SPH3 have M_p = 5 M_J, all orbiting a 1 ${M}_{\odot }$ star.

For the simulations with ${M}_{{\rm{p}}}=5\;{M}_{\mathrm{Jup}},$ we model the disk using 2 million SPH particles. In order to achieve an equivalent resolution within the CPD, we used 16 million SPH particles in SPH1. We vertically resolve the CPD into ∼4 SPH smoothing lengths h (typically h ∼10⁻³ R_p in the CPD midplane) ensuring that we do not underestimate the midplane density (Nelson 2006). At this resolution, the artificial viscosity in the simulations gives an effective Shakura–Sunyaev alpha parameter of $\alpha \sim 0.005$ within the CPDs. Due to the extreme computational expense, we halted SPH1 after 10 orbits of the circumstellar disk, although lower-resolution tests indicate that it has settled into a steady state by this time. SPH2 and SPH3 were halted after 50 orbits (see Figure 1).

Zoom In Zoom Out Reset image size

These simulations are limited in that we neglect complex radiation physics, including only viscous and shock heating but not passive heating. However, the CPD structures we focus on are still present in radiation hydrodynamic runs as shown by Ayliffe & Bate (2009).

The main driver of this investigation is to study under which conditions an accreting protoplanet would be detectable through ALMA observations of line emission. We have chosen bright and commonly observed CO transitions that lie within the submillimeter range in ALMA. Rotational transitions of ${}^{12}\mathrm{CO}$ and isotopologues, ${}^{13}\mathrm{CO}$ and ${{\rm{C}}}^{18}{\rm{O}}$, are known to trace the gas in protoplanetary disk gaps and cavities (Bruderer 2013; Perez et al. 2015). Most importantly, these lines contain important kinematic information, essential for detecting companion objects embedded in the gas inside dust-depleted cavities.

We consider two model disks for our analysis. A nominal disk located at 100 pc, inclined by 20°, and hosting a planet at 100 AU is used to illustrate the observational features revealed through CO kinematics. The second model is tailored for the HD 100546 disk, with an inclination of 42° and a 5 ${M}_{\mathrm{Jup}}$ planet candidate at ${R}_{{\rm{p}}}=52\;\mathrm{AU}$ (Quanz et al. 2014). Our simulations are resampled accordingly in Cartesian coordinates using a linear interpolation scheme via SPLASH (Price 2007). The Cartesian cells are perfectly cubic, with a cell size of 0.014 R_p, which after scaling becomes 1.4 AU and 0.7 AU for R_p = 100 AU and 52 AU, respectively.

For the isothermal runs (SPH1 and SPH2), temperature is an imposed function of orbital radius, while for the adiabatic run (SPH3) we use the SPH internal energies to calculate temperature by assuming a standard mean molecular weight and adiabatic index. We scale the disks using $T(R)\propto {R}^{-1/2},$ consistent with measurements of flaring disks (Kenyon & Hartmann 1987; Andrews et al. 2011). After scaling, temperatures at 100 AU reach ∼60 K for the isothermal disks and ∼2000 K for SPH3, well above CO freeze-out (20–25 K).

In the event of CO freeze-out, a similar RT calculation can be applied to species with enhanced abundances where CO is depleted, such as DCO⁺ or N₂H⁺, or species that are formed by surface reactions with CO ice, such as H₂CO, which have recently been detected in disks at or beyond the CO snow-line (Mathews et al. 2013; Qi et al. 2013, 2013).

We compute synthetic images in CO(2-1) with the radiative transfer code radmc3D (Dullemond et al. 2015). Line radiative transfer is done in LTE, using molecular data from the LAMDA database.⁶ Line widths include thermal broadening and a local (spatially unresolved) microscopic turbulence set to a constant value of 0.1 $\mathrm{km}\;{{\rm{s}}}^{-1}$. We used a fiducial molecular abundance relative of H₂ relative to ${}^{12}\mathrm{CO}$ of 10⁻⁴. We adopted an ISM isotopic abundance ${}^{12}{\rm{C}}{/}^{13}{\rm{C}}$ of 76 (Stahl et al. 2008) and 500 (Wilson & Rood 1994) for ${}^{13}\mathrm{CO}$ and ${{\rm{C}}}^{18}{\rm{O}}$, respectively.

Channel maps are rendered using radmc3d ray-tracing. The results are synthetic data cubes centered on the star, with a total width in velocity of ∼16 $\mathrm{km}\;{{\rm{s}}}^{-1}$ and individual channels of 0.1 $\mathrm{km}\;{{\rm{s}}}^{-1}$. These data cubes represent our sky model, which is subsequently Fourier transformed and resampled to ALMA's visibility plane.

For completeness, we calculate the continuum assuming a simple dust distribution consisting of 30% amorphous carbon grains (Li & Greenberg 1997) and 70% astronomical silicates (Draine & Lee 1984), following the gas density. Grain size distribution follow a power law with exponent 3.5 from 0.05 to 1000 μm. We compute dust opacities with Mie theory. CPD dust continuum predictions are addressed elsewhere in the literature (see Wolf & D'Angelo 2005; Isella et al. 2014), but are strongly affected by gas pressure bumps (D'Angelo & Podolak 2015). Continuum emission in submillimeter observations arises from dust thermal radiation and does not contain kinematic information, thereby it is not within the scope of this paper.

Figure 2 presents our ${}^{13}\mathrm{CO}$ channel map predictions at velocities −1.0, 0.0, +1.0, +2.0, and +3.0 $\mathrm{km}\;{{\rm{s}}}^{-1}$ for the nominal disk model. The top and middle panels are isothermal simulations for 1 and 5 ${M}_{\mathrm{Jup}}$ planets (SPH1 and SPH2, respectively), while the bottom panels show adiabatic results (SPH3).

Zoom In Zoom Out Reset image size

We find that the presence of a CPD produces deviations from circumstellar Keplerian kinematics. The CPD is in itself a miniature Keplerian accretion disk embedded in the gap, and it can be separated in velocity from the overall Keplerian pattern of the circumstellar disk. This velocity separation can be seen in Figure 2 as compact emission at the CPD position. This compact emission is persistent over a velocity range given by the CPD kinematics.

Spectra extracted from the vicinity of the CPDs produce a broader profile when compared with a spectrum extracted from the point-symmetric location at the opposite side of the disk (see Figure 3). The profile extracted from the vicinity of the 1 ${M}_{\mathrm{Jup}}$ planet shows a broad core on top of even broader wings (top spectrum in Figure 3). The middle panel in Figure 3 shows that the isothermal CPD around the 5 ${M}_{\mathrm{Jup}}$ planet exhibits a distinct double-peaked profile, with a peak separation of ∼1.5 $\mathrm{km}\;{{\rm{s}}}^{-1}$. Full CPD spectra span over ${\rm{\Delta }}v\;\sim$ 2.5 $\mathrm{km}\;{{\rm{s}}}^{-1}$ for the isothermal 1 ${M}_{\mathrm{Jup}}$ CPD, while the 5 ${M}_{\mathrm{Jup}}$ covers more than ${\rm{\Delta }}v\;\sim$ 4 $\mathrm{km}\;{{\rm{s}}}^{-1}$. The double-peaked profile for the 1 ${M}_{\mathrm{Jup}}$ case is unresolved at 0.1 $\mathrm{km}\;{{\rm{s}}}^{-1}$ resolution. The width of these spectral features may inform on the size of their respective CPDs and, ultimately, on the mass of the accreting planet via their Hill radii. Figure 3 also shows that the CPD line wings for SPH1 and SPH2 end in an abrupt shoulder. This is likely due to our planet accretion model, a point mass with a sink radius inside which particles are swallowed and their kinematics cannot be sampled.

Zoom In Zoom Out Reset image size

The adiabatic CPD spectrum does not reveal distinctive features (bottom spectra in Figure 3). The compact CPD emission is also much less distinctive in the adiabatic disk (bottom panels in Figure 2). In opposition to the isothermal disk where all compressive work is radiated away immediately, the adiabatic disk cannot cool, reaching temperatures of ∼2000 K around the CPD. This causes the hot gas to rapidly fill the gap preventing a CPD to fully develop (see the second moment map in Figure 3), which hinders clear identification of spectral features in an adiabatic flow. Indeed, adiabatic and isothermal disks represent two extrema of the phenomenon that we are modeling.

A second signpost of planet formation arises when the Doppler-shifted line emission of the circumstellar disk's Keplerian pattern sweeps across the CPD location. The butterfly pattern becomes strongly bent and twisted, while the point-symmetric location at the opposite side of the disk remains undisturbed. This can be seen in the central panels of Figure 2. The insets show the same twisted pattern, but for ${}^{12}\mathrm{CO}$. Optically thicker than ${}^{13}\mathrm{CO}$, the ${}^{12}\mathrm{CO}$ maps still reveal the kinematic bend, even for a 1 ${M}_{\mathrm{Jup}}$ planet whose gap is shallower and its ${}^{12}\mathrm{CO}$ appears optically thick around the CPD. For SPH3, the CPD vicinity is much hotter than in the isothermal cases, producing enhanced ${}^{12}\mathrm{CO}$ emission (see the bottom inset in Figure 2).

As noted in Section 3, our choice of temperature profile in SPH1 and SPH2 is inconsistent with the scaling used for the RT calculations. Self-consistency should result in a CPD approximately twice as thick as in the SPH, reducing the midplane densities and emitted fluxes by a similar factor. It is possible that the signatures highlighted in Figures 2 and 3 would thus be at a lower contrast to the background, but still present at detectable levels, as the velocities are largely unaffected. CO freeze-out at the CPD's location would be prevented by adding a background temperature to account for accretion radiation feedback (Montesinos et al. 2015) and incident radiation from their environment (see Shabram & Boley 2013). Including self-consistent thermal physics in future simulations will settle this discrepancy.

4.1. HD 100546 through a 15 km Baseline Submillimeter Observation

To assess the observability of these kinematic CPD signposts, we performed a second calculation tailored for the protoplanet candidate in HD 100546, based on the SPH2 run. We tied the fluxes of our model to match approximately previous CO observations of this source (Pineda et al. 2014; Walsh et al. 2014). We filtered our sky model using the uv coverage from the long baseline (∼15 km) Science Verification observations of HL Tau (Partnership et al. 2015). We corrupted our model with thermal and phase noises extracted directly from the HL Tau data set. The simulated observation was then self-calibrated and CLEANed using routines in the CASA package.

Figure 4 shows the HD 100546 simulated ALMA observation. Left and right panels show selected channels for ${}^{12}\mathrm{CO}$ and ${}^{13}\mathrm{CO}$ emission. The upper panels illustrate the recovery of the bent Keplerian locus, while the lower panels show how the CPD compact emission can also be detected. Both emission features are recovered at the 5σ level (the rms noise is ∼1 mJy beam⁻¹).

Zoom In Zoom Out Reset image size

The presence of a CPD produces distinct signposts in simulated CO maps: a striking compact emission separated in velocity from the overall Keplerian pattern of the circumstellar disk, a strong influence on the velocity pattern of the gas when the Doppler-shifted line emission sweeps across the CPD location, and a local increase in the velocity dispersion. Moreover, for the locally isothermal simulation with 5 ${M}_{\mathrm{Jup}},$ even the CPD spectra exhibit a double-peaked profile. These distinctive features rely on kinematics and can reveal the presence of an embedded CPD perturber, even in optically thick tracers like ${}^{12}\mathrm{CO}$.

The feasibility of an ALMA observation of HD 100546 was assessed by corrupting our synthetic visibilities with realistic phase and thermal noises extracted from the HL Tau long baseline campaign. We found that these CPDs are detectable in all three signposts with ALMA Cycle 3 capabilities for both 1 and 5 ${M}_{\mathrm{Jup}}$ protoplanets, when embedded in a locally isothermal disk. On the other hand, in the pessimistic case of an adiabatic disk, the CPD formation was hampered, hence detectability is scant. Previous radiation hydrodynamic calculations show that CPDs are better described by isothermal disks and that the adiabatic simulation is in fact a rather pessimistic case (see Ayliffe & Bate 2009, their Figure 12).

The immediate vicinity of the planet's Hill sphere, including the CPD itself, offers an environment for gas-phase physics that produces distinctive kinematic observational features. Future ALMA long baseline observations of gas tracers could detect these signposts of planet formation and could provide not only confirmation of forming planets but also valuable kinematic information on CPD physics.

We thank the referee for a careful review. S.P. acknowledges financial support by FONDECYT grant 3140601. Financial support was provided by Millennium Nucleus RC130007 (Chilean Ministry of Economy), and additionally by FONDECYT grants 1130949, 1141175, 3140634. A.D. and P.R. acknowledge ALMA-CONICYT grants 31120007 and 31120006. S.M. acknowledges CONICYT-PCHA 2014-22140628. M.M. acknowledges CONICYT-Gemini grant 32130007. The Geryon2 cluster housed at Centro de Astro-Ingenieria UC was used for the SPH calculations. The BASAL PFB-06 CATA, Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL-130008 provided funding for improvements to the Geryon2 cluster.