#5E?60??@><;<:E.0A9@E#A/960.@6<;? #5E?60?.;1?@><;<:E
&526?8%A/?@>A0@A>2?.@645;4A9.>$2?<9A@6<;#><720@&526?8%A/?@>A0@A>2?.@645;4A9.>$2?<9A@6<;#><720@
%$#(&52#9.;2@G6?8;@2>.0@6<;?;@2>=>2@.@6<;%$#(&52#9.;2@G6?8;@2>.0@6<;?;@2>=>2@.@6<;
%5.;476.,5.;4
';6B2>?6@E<3!2B.1..?(24.?
,5.<5A.;,5A
';6B2>?6@E<3!2B.1..?(24.?
F5.<5A.;F5AA;9B21A
.;2A.;4
';6B2>?61.1125692
(6B6.;.(AF:I;
<6;@ "/?2>B.@<>E#<;@6H06.';6B2>?61.1.@K960.125692
%2.; ;1>2C?
.>B.>1%:6@5?<;6.;2;@2>3<>?@><=5E?60?
%22;2D@=.423<>.116@6<;.9.A@5<>?
<99<C@56?.;1.116@6<;.9C<>8?.@5@@=?1646@.9?05<9.>?56=A;9B21A=5E?.?@>-3.0-.>@6092?
#.>@<3@52?@><=5E?60?.;1?@><;<:E<::<;?
$2=<?6@<>E6@.@6<;$2=<?6@<>E6@.@6<;
,5.;4%,5A,A.;4AF:I;((;1>2C?% 6>;?@629&A992:<;1#.>=2;@2>
?299.#J>2F 2;6?@E )69;2>.>A@2.A.6*$6006&526?8
%A/?@>A0@A>2?.@645;4A9.>$2?<9A@6<;#><720@%$#(&52#9.;2@G6?8;@2>.0@6<;?
;@2>=>2@.@6<;
?@><=5E?60.9<A>;.92@@2>?
5@@=1D1<6<>4..3
&56?>@60926?=><@20@21/E0<=E>645@.;1<>>29.@21>645@?@5.?/22;/><A45@@<E<A/E646@.9%05<9.>?56='!(
C6@5=2>:6??6<;3><:@52>645@?5<912>?+<A.>23>22@<A?2@56?>@60926;.;EC.E@5.@6?=2>:6@@21/E@52
0<=E>645@.;1>29.@21>645@?9246?9.@6<;@5.@.==962?@<E<A>A?2<><@52>A?2?E<A;221@<</@.6;=2>:6??6<;3><:
@52>645@?5<912>?16>20@9EA;92??.116@6<;.9>645@?.>26;160.@21/E.>2.@6B2<::<;?9602;?26;@52>20<>1.;1
<><;@52C<>86@?293
&56?>@60925.?/22;.002=@213<>6;09A?6<;6;#5E?60??@><;<:E.0A9@E#A/960.@6<;?/E.;.A@5<>6F21
.1:6;6?@>.@<><3646@.9%05<9.>?56='!(<>:<>26;3<>:.@6<;=92.?20<;@.0@1646@.9?05<9.>?56=A;9B21A
A@5<>?A@5<>?
%5.;476.,5.;4,5.<5A.;,5A.;2A.;4(6B6.;.(AF:I;%2.; ;1>2C?&69:.;6>;?@629
<>;296?#A992:<;1<5; .>=2;@2>;1>2.?299..A>. #J>2F E>6.:2;6?@E.B61)69;2>
9J:2;@.>A@2.A*A2!6;4.6.;1A0.$6006
&56?.>@60926?.B.69./92.@646@.9%05<9.>?56='!(5@@=?1646@.9?05<9.>?56=A;9B21A=5E?.?@>-3.0-.>@6092?
The Disk Substructures at High Angular Resolution Project (DSHARP).
VII. The Planet–Disk Interactions Interpretation
Shangjia Zhang
1
, Zhaohuan Zhu
1
, Jane Huang
2
, Viviana V. Guzmán
3,4
, Sean M. Andrews
2
, Tilman Birnstiel
5
,
Cornelis P. Dullemond
6
, John M. Carpenter
3
, Andrea Isella
7
, Laura M. Pérez
8
, Myriam Benisty
9,10
, David J. Wilner
2
,
Clément Baruteau
11
, Xue-Ning Bai
12
, and Luca Ricci
13
1
Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S. Maryland Parkway, Las Vegas, NV 89154, USA; zhaohuan.zhu@unlv.edu
2
Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
3
Joint ALMA Observatory, Avenida Alonso de Crdova 3107, Vitacura, Santiago, Chile
4
Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile
5
University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 Munich, Germany
6
Zentrum für Astronomie, Heidelberg University, Albert Ueberle Str. 2, D-69120 Heidelberg, Germany
7
Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA
8
Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile
9
Unidad Mixta Internacional Franco-Chilena de Astronomía (CNRS, UMI 3386), Departamento de Astronomía,
Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile
10
Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France
11
CNRS/Institut de Recherche en Astrophysique et Planétologie, 14 avenue Edouard Belin, F-31400 Toulouse, France
12
Institute for Advanced Study and Tsinghua Center for Astrophysics, Tsinghua University, Beijing 100084, People’s Republic of China
13
Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff Street, Northridge, CA 91130, USA
Received 2018 October 8; revised 2018 November 26; accepted 2018 November 26; published 2018 December 26
Abstract
The Disk Substructures at High Angular Resolution Project (DSHARP) provides a large sample of protoplanetary
disks with substructures that could be induced by young forming planets. To explore the properties of planets that
may be responsible for these substructures, we systematically carry out a grid of 2D hydrodynamical simulations,
including both gas and dust components. We present the resulting gas structures, including the relationship
between the planet mass, as well as (1) the gaseous gap depth/width and (2) the sub/super-Keplerian motion
across the gap. We then compute dust continuum intensity maps at the frequency of the DSHARP observations.
We provide the relationship between the planet mass, as well as (1) the depth/width of the gaps at millimeter
intensity maps, (2) the gap edge ellipticity and asymmetry, and (3) the position of secondary gaps induced by the
planet. With these relationships, we lay out the procedure to constrain the planet mass using gap properties, and
study the potential planets in the DSHARP disks. We highlight the excellent agreement between observations and
simulations for AS 209 and the detectability of the young solar system analog. Finally, under the assumption that
the detected gaps are induced by young planets, we characterize the young planet population in the planet mass–
semimajor axis diagram. We find that the occurrence rate for >5 M
J
planets beyond 5–10 au is consistent with
direct imaging constraints. Disk substructures allow us to probe a wide-orbit planet population (Neptune to Jupiter
mass planets beyond 10 au) that is not accessible to other planet searching techniques.
Key words: hydrodynamics – planet–disk interactions – planets and satellites: detection – planets and satellites:
formation – protoplanetary disks – submillimeter: planetary systems
1. Introduction
Discoveries over the past few decades show that planets are
common. The demographics of exoplanets have put constraints
on planet formation theory (e.g., see the review by Chabrier
et al. 2014; Johansen et al. 2014; Raymond et al. 2014).
Unfortunately, most discovered exoplanets are billions of years
old and have therefore been subject to significant orbital
dynamical alteration after their formation (e.g., review by
Davies et al. 2014). To test planet formation theory, it is crucial
to constrain the young planet population right after they are
born in protoplanetary disks. However, the planet search
techniques that have discovered thousands of exoplanets
around mature stars are not efficient at finding planets around
young stars (<10 Myr old) mainly due to their stellar variablity
and the presence of the protoplanetary disks. Fewer than 10
young planet candidates in systems <10 Myr have been
detected so far (e.g., CI Tau b, Johns-Krull et al. 2016;V
830 Tau b, Donati et al. 2016; Tap 26 b, Yu et al. 2017; PDS 70
b, Keppler et al. 2018; LkCa 15 b, Sallum et al. 2015).
On the other hand, recent high resolution imaging at near-IR
wavelengths (with the new adaptive optics systems on 10 m
class telescopes) and interferometry at radio wavelengths
(especially the ALMA and the VLA) can directly probe the
protoplanetary disks down to astronomical unit scales, and a
variety of disk features ( such as gaps, rings, spirals, and large-
scale asymmetries) have been revealed (e.g., Casassus et al.
2013; van der Marel et al. 2013
; ALMA Partnership et al. 2015;
Andrews et al. 2016; Garufi et al. 2017). Despite that there are
other possibilities for producing these features, they may be
induced by young planets in these disks, and we can use these
features to probe the unseen young planet population.
Planet–disk interactions have been studied over the past three
decades with both analytical approaches (Goldreich & Tremaine
1980; Tanaka et al. 2002) and numerical simulations (Kley &
Nelson 2012;Baruteauetal.2014). While the earlier work focused
on planet migration and gap opening, more recently efforts have
been dedicated to studying observable disk features induced by
planets (Wolf & D’Angelo 2005; Dodson-Robinson & Salyk 2011;
Zhu et al. 2011;Gonzalezetal.2012; Pinilla et al. 2012;
The Astrophysical Journal Letters, 869:L47 (32pp), 2018 December 20 https://doi.org/10.3847/2041-8213/aaf744
© 2018. The American Astronomical Society. All rights reserved.
1
Ataiee et al. 2013; Bae et al. 2016; Kanagawa et al. 2016; Rosotti
et al. 2016; Isella & Turner 2018), including the observational
signatures in near-IR scattered light images (e.g., Dong et al. 2015;
Fung & Dong 2015;Zhuetal.2015a), (sub)millimeter dust
thermal continuum images (Dipierro et al. 2015; Picogna & Kley
2015;Dong&Fung2017;Dongetal.2018),and(su b)millimeter
molecular line channel maps that trace the gas kinematics at the
gap edges or around the planet (Perez et al. 2015; Pinte et al. 2018;
Teague et al. 2018).
Among all these indirect methods for probing young planets
at various wavelengths, only dust thermal emission at ( sub)
millimeter wavelengths allows us to probe low-mass planets,
because a small change in the gas surface density due to the
low-mass planet can cause dramatic changes in the dust surface
density (Paardekooper & Mellema 2006; Zhu et al. 2014).
However, this also means that hydrodynamical simulations
with both gas and dust components are needed to study the
expected disk features at (sub)millimeter wavelengths. Such
simulations are more complicated due to the uncertainties about
the dust size distribution in protoplanetary disks. Previously,
hydrodynamical simulations have been carried out to explain
features in individual sources (e.g., Jin et al. 2016; Dipierro
et al. 2018; Fedele et al. 2018). With many disk features
revealed by the Disk Substructures at High Angular Resolution
Project (DSHARP; Andrews et al. 2018), a systematic study of
how the dust features relate to the planet properties is desirable.
By conducting an extensive series of disk models spanning a
substantial range in disk and planet properties, we can enable a
broad exploration of parameter space that can then be used to
rapidly infer young planet populations from the observations,
and we will also be more confident that we are not missing
possible parameter space for each potential planet.
In this work, we carry out a grid of hydrodynamical
simulations including both gas and dust components. Then,
assuming different dust size distributions, we generate intensity
maps at the observation wavelength of DSHARP. In Section 2,
we describe our methods. The results are presented in
Section 3. The derived young planet properties for the
DSHARP disks are given in Section 4. After a short discussion
in Section 5, we conclude the Letter in Section 6.
2. Method
We carry out 2D hydrodynamical planet–disk simulations
using the modified version of the grid-based code FARGO
(Masset 2000) called Dusty FARGO-ADSG (Baruteau &
Masset 2008a, 2008b; Baruteau & Zhu 2016). The gas
component is simulated using finite difference methods (Stone
& Norman 1992), while the dust component is modeled as
Lagrangian particles. To allow our simulations to be as scale-
free as possible, we do not include disk self-gravity, radiative
cooling, or dust feedback. These simplifications are suitable for
most disks observed in DSHARP. Most of the features in these
disks lie beyond 10 au where the irradiation from the central
star dominates the disk heating such that the disk is nearly
vertically isothermal close to the midplane ( D’Alessio et al.
1998). Although the dust dynamical feedback to the gas is
important when a significant amount of dust accumulates at gap
edges or within vortices (Fu et al. 2014; Crnkovic-Rubsamen
et al. 2015), simulations that have dust particles but do not
include dust feedback to the gas (so-called “passive dust”
models) serve as reference models and allow us to scale our
simulations freely to disks with different dust-to-gas mass
ratios and dust size distributions. As shown in Section 4,
passive dust models are also adequate in most of our cases
(especially when the dust couples with the gas relatively well).
Simulations with dust feedback will be presented in C. Yang &
Z. Zhu
(2019, in prepration).
2.1. Setup: Gas and Dust
We adopt polar coordinates (r, θ) centered on the star and fix
the planet on a circular orbit at r=1. Since the star is
wobbling around the center of mass due to the perturbation by
the planet, indirect forces are applied to this noninertial
coordinate frame.
We initialize the gas surface density as
rrr,1
gg,0 0
1
S=S
-
() ( ) ()
where r
0
is also the position of the planet and we set
r
r 1
p0
==
.
For studying gaps of individual sources in Section 4,wescaleΣ
g,0
to be consistent with the DSHARP observations. We assume
locally isothermal equation of state, and the temperature at radius r
follows
Tr T rr
00
12
=
-
() ( )
. T is related to the disk scale height
h as h/r=c
s
/v
f
where
c
RT P
s
2
m==S
and μ=2.35. With
our setup, h/r changes as r
1/4
. In the rest of the text, when we
give a value of h/r, we are referring to h/r at r
0
.
Our numerical grid extends from 0.1 r
0
to 10 r
0
in the radial
direction and 0 to 2π in the θ direction. For low viscosity cases
(α=10
−4
and 10
−3
), there are 750 grid points in the radial
direction and 1024 grid points in the θ direction. This is
equivalent to 16 grid points per scale height at r
0
if h/r=0.1.
For high viscosity cases (α=0.01), less resolution is needed
so there are 375 and 512 grid points in the radial and θ
direction. For simulations to fit AS 209 in Section 4.1.1, the
resolution is 1500 and 2048 grid points in the radial and θ
direction to capture additional gaps at the inner disk. We use
the evanescent boundary condition, which relaxes the fluid
variables to the initial state at r<0.12r
0
and r>8r
0
.A
smoothing length of 0.6 disk scale height at r
0
is used to
smooth the planet’s potential (Müller et al. 2012).
We assume that the dust surface density is 1/100 of the gas
surface density initially. The open boundary condition is
applied for dust particles, so that the dust-to-gas mass ratio for
the whole disk can change with time.
The dust particles experience both gravitational forces and
aerodynamic drag forces. The particles are pushed at every
timestep with the orbital integrator. When the particle’s
stopping time is smaller than the numerical timestep, we use
the short friction time approximation to push the particle. Since
we are interested in disk regions beyond tens of astronomical
units, the disk density is low enough that the molecular mean-
free path is larger than the size of dust particles. In this case, the
drag force experienced by the particles is in the Epstein regime.
The Stokes number St for particles (also called particles’
dimensionless stopping time) is
t
s
s
St
2
1.57 10
1gcm 1mm
100gcm
,
2
pp
g
stop
gas
3
3
2
pr r
=W
=
S
=´
S
-
-
-
()
where ρ
p
is the density of the dust particle, s is the radius of
the dust particle, and Σ
g
is the gas surface density. We
assume ρ
p
=1gcm
−3
in our simulations. We use 200,000 and
2
The Astrophysical Journal Letters, 869:L47 (32pp), 2018 December 20 Zhang et al.
100,000 particles for high and low resolution runs, respectively.
Each particle is a super particle representing a group of real dust
particles having the same size. The super particles in our
simulations have Stokes numbers ranging from 1.57×10
−5
to
1.57, or physical radii ranging from 1 μmto10cmifΣ
g,0
=
10 g cm
−2
and ρ
p
=1gcm
−3
. We distribute super particles
uniformly in log(s) space, which means that we have the same
number of super particles per decade in size. Since dust-to-gas
back reaction is not included, we can scale the dust size
distribution in our simulations to any desired distribution.
During the simulation, we keep the size of the super-particle
the same no matter where it drifts to. Thus, the super-particle’s
Stokes number changes when this particle drifts in the disk,
because the particle’s Stokes number also depends on the local
disk surface density (Equation (2)). More specifically, during
the simulation, the Stokes number of every particle varies so as
to be inversely proportional to the local gas surface density.
Turbulent diffusion for dust particles is included as random
kicks to the particles (Charnoz et al. 2011; Fuente et al. 2017).
The diffusion coefficient is related to the α parameter as in
Youdin & Lithwick (2007) through the so-called Schmidt
number Sc. In this work, Sc is defined as the ratio between the
angular momentum transport coefficient (ν) and the gas
diffusion coefficient (D
g
). We set Sc=1, which serves as a
good first-order approximation, although that Sc can take on
different values and its value can differ between the radial and
vertical directions (Zhu et al. 2015b; Yang et al. 2018),
2.2. Grid of Models
To explore the full parameter space, we choose three values
for
hr 0.05, 0.07, 0.1
r
0
(
)( )
, five values for the planet–star
mass ratio (q≡M
p
/M
*
=3.3×10
−5
,10
−4
, 3.3×10
−4
,
10
−3
, 3.3×10
−3
M
*
, or roughly M
p
=11 M
⊕
,33M
⊕
,
0.35 M
J
,1M
J
, 3.5 M
J
if M
*
= M
e
), and three values for the
disk turbulent viscosity coefficient (α = 0.01, 0.001, 0.0001).
Thus, we have 45 simulations in total. We label
each simulation in the following manner: h5am3p1 means
h/r=0.05, α=10
−3
(m3 in h5am3p1 means minus 3),
M
p
/M
*
=3.3×10
−5
M
*
(p1 refers to the lowest planet mass
case). We also run some additional simulations for individual
sources (e.g., AS 209, Elias 24) which will be presented in
Section 4.1 and Guzmán et al. (2018).
This parameter space represents typical disk conditions.
Protoplanetary disks normally have h/r between 0.05 and 0.1
at
r
10 au>
(D’Alessio et al. 1998). While a moderate
α∼10
−2
is preferred to explain the disk accretion (Hartmann
et al. 1998), recent works suggest that a low turbulence level
(α < 10
−2
) is needed to explain molecular line widths in TW
Hya (Flaherty et al. 2018) and dust settling in HL Tau (Pinte
et al. 2016). When α is smaller than 10
−4
, the viscous timescale
over the disk scale height at the planet position (
H
p
2
n
) is
longer than 10
4
/Ω
p
or 1.6 million years at 100 au, so that the
viscosity will not affect the disk evolution significantly. In
Section 4.1, we carry out several simulations with different α
values to extend the parameter space for some sources in the
DSHARP sample. As shown below, when the planet mass is
less than 11 M
⊕
, the disk features are not detectable with
ALMA. When the planet mass is larger than 3.5 M
J
, the disk
features have strong asymmetries, and we should be able to
detect the planet directly though direct imaging techniques.
We run the simulations for 1000 planetary orbits (1000 T
p
),
which is equivalent to 1 Myr for a planet at 100 au or 0.1 Myr
for a planet at 20 au. These timescales are comparable to the
disk ages of the DSHARP sources.
2.3. Calculating Millimeter Continuum Intensity Maps
For each simulation, we calculate the millimeter continuum
intensity maps assuming different disk surface densities and
dust size distributions. Since dust-to-gas feedback is neglected,
we can freely scale the initial disk surface density and dust size
distribution in simulations to match realistic disks.
Both the disk surface density and dust size distribution have
large impacts on the millimeter intensity maps. If the dust
thermal continuum is mainly from micron sized particles and
the disk surface density is high, these dust particles have small
Stokes numbers (Equation (2)). Consequently, they couple to
the gas almost perfectly and the gaps revealed in millimeter are
very similar to the gaps in the gas. If the millimeter emission is
dominated by millimeter sized particles and the disk surface
density is low, the dust particles can have Stokes numbers close
to 1 and they drift very fast in the disk. In this case, they can be
trapped at the gap edges, producing deep and wide gaps. To
explore how different dust size distributions can affect the
millimeter intensity maps, we choose two very different dust
size distributions to generate intensity maps. For the distribu-
tion referred to as DSD1, we assume n(s)∝s
−3.5
with a
maximum grain size of 0.1 mm in the initial condition
(
p 3.5=-
and s
max
=0.1 mm. This is motivated by recent
(sub)millimeter polarization measurements (Kataoka et al.
2017; Hull et al. 2018), which indicate that the maximum
grain size in a variety of disks is around 0.1 mm. In the other
case referred to as DSD2, we assume n(s)∝s
−2.5
with the
maximum grain size of 1 cm (p = −2.5 and s
max
=1cm). This
shallower dust size distribution is expected from dust growth
models (Birnstiel et al. 2012) and consistent with SED
constraints (D’Alessio et al. 2001) and the spectral index at
millimeter/centimeter wavelengths (Ricci et al. 2010a, 2010b;
Pérez et al. 2015). Both cases assume a minimum grain size of
0.005 μm. We find that the minimum grain size has no effect
on the dust intensity maps since most dust mass is in larger
particles. Coincidentally, these two size distributions lead to the
same opacity at 1.27 mm (1.27 mm is the closest wavelength to
1.25 mm in the table of Birnstiel et al. 2018) in the initial
condition (the absorption opacity for the s
max
= 0.1 mm case is
0.43 cm
2
g
−1
, while for the s
max
= 1 cm case it is 0.46 cm
2
g
−1
based on Birnstiel et al. 2018). More discussion on how to
generalize our results to disks with other dust size distributions
can be found in Section 3.2.2.
For each simulation, we scale the simulation to different disk
surface densities. Then for each surface density, we calculate
the 1.27 mm intensity maps using DSD1 or DSD2 dust size
distributions. For the s
max
=0.1 mm dust size distribution
(DSD1), we calculate the 1.27 mm intensity maps for disks
with Σ
g,0
=0.1, 0.3, 1, 3, 10, 30, and 100 g cm
−2
(seven
groups of models). The maximum-size particle in these disks
(0.1 mm) that dominates the total dust mass corresponds to
St=1.57×10
−1
, 5.23×10
−2
, 1.57×10
−2
, 5.23×10
−3
,
1.57×10
−3
, 5.23×10
−4
, and 1.57×10
−4
at r=r
p
. For the
s
max
=1 cm cases (DSD2), we vary Σ
g,0
as 1, 3, 10, 30, and
100 g cm
−2
(five groups of models), and the corresponding St
for 1 cm particles at r=r
p
is 1.57, 5.23×10
−1
, 1.57×10
−1
,
5.23×10
−2
, and 1.57×10
−2
. For each given surface density
3
The Astrophysical Journal Letters, 869:L47 (32pp), 2018 December 20 Zhang et al.
