The Disk Substructures at High Angular Resolution Project (DSHARP): VII.   The Planet-Disk Interactions Interpretation

Shangjia Zhang; Zhaohuan Zhu; Jane Huang; Viviana V. Guzm\'an; Sean M.; Andrews; Tilman Birnstiel; Cornelis P. Dullemond; John M. Carpenter; Andrea; Isella; Laura M. P\'erez; Myriam Benisty; David J. Wilner; Cl\'ement; Baruteau; Xue-Ning Bai; Luca Ricci

arXiv:1812.04045·astro-ph.EP·January 9, 2019

The Disk Substructures at High Angular Resolution Project (DSHARP): VII. The Planet-Disk Interactions Interpretation

Shangjia Zhang, Zhaohuan Zhu, Jane Huang, Viviana V. Guzm\'an, Sean M., Andrews, Tilman Birnstiel, Cornelis P. Dullemond, John M. Carpenter, Andrea, Isella, Laura M. P\'erez, Myriam Benisty, David J. Wilner, Cl\'ement, Baruteau, Xue-Ning Bai, Luca Ricci

PDF

TL;DR

This study uses hydrodynamical simulations to link observed disk substructures with potential young planets, enabling constraints on planet masses and distributions in protoplanetary disks.

Contribution

It provides a systematic method to infer planet properties from disk substructure features, advancing understanding of planet formation and disk-planet interactions.

Findings

01

Gas gap properties correlate with planet mass.

02

Dust gap features can constrain planet mass and orbit.

03

Detected gaps suggest a population of young, wide-orbit planets.

Abstract

The Disk Substructures at High Angular Resolution Project (DSHARP) provides a large sample of protoplanetary disks having substructures which 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 2-D hydrodynamical simulations including both gas and dust components. We present the resulting gas structures, including the relationship between the planet mass and 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 and 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,…

Tables5

Table 1. Table 1: The relation between the gap width Δ Δ \Delta and K ′ superscript 𝐾 ′ K^{\prime}

Parameters	$Δ_{g}$	$Δ_{d, 0 p 1}$	$Δ_{d, 0 p 3}$	$Δ_{d, 1}$	$Δ_{d, 3}$	$Δ_{d, 10}$	$Δ_{d, 30}$	$Δ_{d, 100}$
$S t_{m a x}$ (DSD1)	–	1.57 $\times$ $10^{- 1}$	5.23 $\times$ $10^{- 2}$	1.57 $\times$ $10^{- 2}$	5.23 $\times$ $10^{- 3}$	1.57 $\times$ $10^{- 3}$	5.23 $\times$ $10^{- 4}$	1.57 $\times$ $10^{- 4}$
$A$ ${p = - 3.5, s_{m a x} = 0.1 m m}$	1.05	1.09	1.73	2.00	1.25	1.18	0.98	1.11
$B$ ${p = - 3.5, s_{m a x} = 0.1 m m}$	0.26	0.07	0.24	0.36	0.27	0.29	0.25	0.29
Uncertainty in $l o g_{10} (K^{'})$	${}_{- 0.12}^{+ 0.03}$	${}_{- 1.12}^{+ 0.86}$	${}_{- 0.63}^{+ 0.53}$	${}_{- 0.50}^{+ 0.21}$	${}_{- 0.16}^{+ 0.22}$	${}_{- 0.17}^{+ 0.14}$	${}_{- 0.14}^{+ 0.16}$	${}_{- 0.16}^{+ 0.13}$
$S t_{m a x}$ (DSD2)	–	–	–	1.57	5.23 $\times$ $10^{- 1}$	1.57 $\times$ $10^{- 1}$	5.23 $\times$ $10^{- 2}$	1.57 $\times$ $10^{- 2}$
$A$ ${p = - 2.5, s_{m a x} = 1 c m}$	–	–	–	–	1.10	1.13	1.55	2.00
$B$ ${p = - 2.5, s_{m a x} = 1 c m}$	–	–	–	–	0.05	0.09	0.23	0.36
Uncertainty in $l o g_{10} (K^{'})$	–	–	–	–	${}_{- 1.06}^{+ 0.80}$	${}_{- 0.77}^{+ 0.70}$	${}_{- 0.59}^{+ 0.65}$	${}_{- 0.29}^{+ 0.28}$

Table 2. Table 2: The relation between δ 𝛿 \delta -1 and K 𝐾 K , where δ 𝛿 \delta is the gap depth.

Parameters	$δ_{g}$ - 1	$δ_{d, 0 p 1}$ - 1	$δ_{d, 0 p 3}$ - 1	$δ_{d, 1}$ - 1	$δ_{d, 3}$ - 1	$δ_{d, 10}$ - 1	$δ_{d, 30}$ - 1	$δ_{d, 100}$ - 1
$S t_{m a x}$ (DSD1)	–	1.57 $\times$ $10^{- 1}$	5.23 $\times$ $10^{- 2}$	1.57 $\times$ $10^{- 2}$	5.23 $\times$ $10^{- 3}$	1.57 $\times$ $10^{- 3}$	5.23 $\times$ $10^{- 4}$	1.57 $\times$ $10^{- 4}$
C ${$ 3.5, 0.1 mm $}$	0.002	14.9	1.18	0.178	0.244	0.135	0.0917	0.0478
D ${$ 3.5, 0.1 mm $}$	2.64	0.926	1.36	1.54	1.25	1.21	1.18	1.23
Uncertainty in $l o g_{10} (K)$	${}_{- 0.13}^{+ 0.08}$	${}_{- 0.71}^{+ 0.82}$	${}_{- 0.57}^{+ 0.74}$	${}_{- 0.48}^{+ 0.53}$	${}_{- 0.26}^{+ 0.32}$	${}_{- 0.23}^{+ 0.19}$	${}_{- 0.24}^{+ 1.20}$	${}_{- 0.18}^{+ 0.18}$
$S t_{m a x}$ (DSD2)	–	–	–	1.57	5.23 $\times$ $10^{- 1}$	1.57 $\times$ $10^{- 1}$	5.23 $\times$ $10^{- 2}$	1.57 $\times$ $10^{- 2}$
C ${$ 2.5, 1 cm $}$	–	–	–	271	998	25.5	1.46	0.069
D ${$ 2.5, 1 cm $}$	–	–	–	1.22	0.533	1.17	1.50	1.94
Uncertainty in $l o g_{10} (K)$	–	–	–	${}_{- 1.26}^{+ 1.45}$	${}_{- 3.38}^{+ 2.64}$	${}_{- 1.18}^{+ 1.08}$	${}_{- 0.56}^{+ 0.78}$	${}_{- 0.52}^{+ 0.56}$

Table 3. Table 3: Inferred Planet Mass from 19 Gaps

Name	$M_{*}$	$r_{g a p}$	width	$Σ_{d u s t}$	$h / r$	$Σ_{g, 0}$	$Σ_{G I}$	$Σ_{g, 0}^{u s e d}$	$S t_{m a x}^{u s e d}$	$M_{p, a m 4}$	$M_{p, a m 3}$	$M_{p, a m 2}$	Uncertainty
	( $M_{⊙}$ )	(au)	( $Δ$ )	( $g {cm}^{- 2}$ )		( $g {cm}^{- 2}$ )	( $g {cm}^{- 2}$ )	( $g {cm}^{- 2}$ )	( $5.23 \times 10^{- 4}$ )	( $M_{J u p}$ )	( $M_{J u p}$ )	( $M_{J u p}$ )	( $l o g_{10} (M_{p})$ )
AS 209	0.83	9	0.42	1.23	0.04	$>$ 100, 100, 100	1278.4	100, 100, 100	0.33, 3, 30	1.00, 0.81, 0.37	2.05, 1.66, 0.76	4.18, 3.38, 1.56	${}_{- 0.16}^{+ 0.13}$ , ${}_{- 0.17}^{+ 0.14}$ , ${}_{- 0.29}^{+ 0.28}$
AS 209	0.83	99	0.31	0.17	0.08	30, 10, 3	19.2^$\dagger$^$\dagger$footnotemark:	10, 10, –	3, 30, –	0.32, 0.18, –	0.65, 0.37, –	1.32, 0.75, –	${}_{- 0.17}^{+ 0.14}$ , ${}_{- 0.50}^{+ 0.21}$ , –
Elias 24	0.78	57	0.32	0.52	0.09	100, 30, 10	58.6^$\dagger$^$\dagger$footnotemark:	30, 30, –	1, 10, –	0.41, 0.19 –	0.84, 0.40, –	1.72, 0.81, –	${}_{- 0.14}^{+ 0.16}$ , ${}_{- 0.16}^{+ 0.22}$ , –
Elias 27	0.49	69	0.18	0.48	0.09	100, 30, 10	25.6^$\dagger$^$\dagger$footnotemark:	10, 10, –	3, 30, –	0.03, 0.02, –	0.06, 0.05, –	0.12, 0.10, –	${}_{- 0.14}^{+ 0.16}$ , ${}_{- 0.50}^{+ 0.21}$ , –
GW Lup^^The gap of the GW Lup at 74 au has width $Δ$ $>$ 0.15, while the gap of the DoAr 25 at 98 au has $Δ$ $<$ 0.15 before rounding.	0.46	74	0.15	0.13	0.08	10, 3, 3	19.8	10, –, –	3, –, –	0.01, –, –	0.03, –, –	0.06, –, –	${}_{- 0.17}^{+ 0.14}$ , –, –
HD 142666	1.58	16	0.20	1.63	0.05	$>$ 100, 100, 100	814.0	100, 100, 100	0.33, 3, 30	0.15, 0.12, 0.09	0.30, 0.25, 0.19	0.62, 0.50, 0.38	${}_{- 0.16}^{+ 0.13}$ , ${}_{- 0.17}^{+ 0.14}$ , ${}_{- 0.29}^{+ 0.28}$
HD 143006	1.78	22	0.62	0.20	0.04	30, 10, 3	442.7	30, 10, –	1, 30, –	9.75, 2.35, –	19.91, 4.80, –	40.64, 9.81, –	${}_{- 0.14}^{+ 0.16}$ , ${}_{- 0.50}^{+ 0.21}$ , –
HD 143006	1.78	51	0.22	0.14	0.05	30, 10, 3	101.6	30, 10, –	1, 30, –	0.16, 0.14 –	0.33, 0.28, –	0.67, 0.57, –	${}_{- 0.14}^{+ 0.16}$ , ${}_{- 0.50}^{+ 0.21}$ , –
HD 163296	2.04	10	0.24	1.43	0.04	$>$ 100, 100, 100	2273.0	100, 100, 100	0.33, 3, 30	0.35, 0.28, 0.19	0.71, 0.58, 0.39	1.46, 1.18, 0.79	${}_{- 0.16}^{+ 0.13}$ , ${}_{- 0.17}^{+ 0.14}$ , ${}_{- 0.29}^{+ 0.28}$
HD 163296	2.04	48	0.34	0.41	0.06	30, 10, 10	146.0	30, 10, –	1, 30, –	1.07, 0.54, –	2.18, 1.10, –	4.45, 2.24, –	${}_{- 0.14}^{+ 0.16}$ , ${}_{- 0.50}^{+ 0.21}$ , –
HD 163296	2.04	86	0.17	0.15	0.07	30, 10, 3	52.6	30, 10, –	1, 30, –	0.07, 0.08, –	0.14, 0.16, –	0.29, 0.34, –	${}_{- 0.14}^{+ 0.16}$ , ${}_{- 0.50}^{+ 0.21}$ , –
SR 4	0.68	11	0.45	1.56	0.05	$>$ 100, 100, 100	792.8	100, 100, 100	0.33, 3, 30	1.06, 0.86, 0.38	2.16, 1.75, 0.77	4.41, 3.57, 1.57	${}_{- 0.16}^{+ 0.13}$ , ${}_{- 0.17}^{+ 0.14}$ , ${}_{- 0.29}^{+ 0.28}$
DoAr 25^^The gap of the GW Lup at 74 au has width $Δ$ $>$ 0.15, while the gap of the DoAr 25 at 98 au has $Δ$ $<$ 0.15 before rounding.	0.95	98	0.15	0.48	0.07	100, 30, 10	20.0^$\dagger$^$\dagger$footnotemark:	10, 10, –	3, 30, –	(– , 0.10, –)	(0.10, –, –)	(– , 0.95, –)	–, –, –
DoAr 25	0.95	125	0.08	0.14	0.07	30, 10, 3	13.1^$\dagger$^$\dagger$footnotemark:	10, –, –	3, –, –	(0.03, –, –)	– , –, –	– , –, –	–, –, –
Elias 20	0.48	25	0.13	0.80	0.08	100, 30, 30	171.9	100, 30, 30	0.33, 10, 100	–, –, –	(0.05, 0.05, 0.05)	– , –, –	–, –, –
IM Lup	0.89	117	0.13	0.20	0.09	30, 10, 3	16.0^$\dagger$^$\dagger$footnotemark:	10, –, –	3, –, –	(0.09 , –, –)	(0.09, –, –)	–, – , –	–, –, –
RU Lup	0.63	29	0.14	1.13	0.07	$>$ 100, 100, 100	144.1	100, 100, 100	0.33, 3, 30	(0.07, –, –)	(–, 0.07, 0.07)	– , –	–, –, –
Sz 114	0.17	39	0.12	0.22	0.10	30, 10, 3	35.3	30, 10, –	1, 30, –	(0.02 , 0.02, –)	–, –, –	–, – , –	–, –, –
Sz 129	0.83	41	0.08	0.47	0.06	100, 30, 10	77.7^$\dagger$^$\dagger$footnotemark:	30, 30, –	1, 10, –	(–, 0.03 , –)	(0.03, –, –)	–, – , –	–, –, –

Table 4. Table 4: Gap Widths for the Gas, DSD1 and DSD2

h/r	$α$	$q$	$Δ_{g}$	$Δ_{d, 0 p 1}$	$Δ_{d, 0 p 3}$	$Δ_{d, 1}$	$Δ_{d, 3}$	$Δ_{d, 10}$	$Δ_{d, 30}$	$Δ_{d, 100}$	$Δ_{d, 1}$	$Δ_{d, 3}$	$Δ_{d, 10}$	$Δ_{d, 30}$	$Δ_{d, 100}$
0.05	$10^{- 4}$	3.3 $\times 10^{- 5}$	0.09	0.69	0.46	0.21	0.16	0.13	0.12	0.12	0.81	0.75	0.64	0.32	0.19
0.05	$10^{- 4}$	1 $\times 10^{- 4}$	0.24	0.80	0.61	0.48	0.25	0.24	0.22	0.21	1.00	0.87	0.76	0.57	0.26
0.05	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.32	0.82	0.57	0.98	0.30	0.29	0.29	0.27	1.00	0.92	0.77	0.55	0.48
0.05	$10^{- 4}$	1 $\times 10^{- 3}$	0.42	0.56	0.56	0.55	0.42	0.40	0.38	0.37	1.00	0.56	0.56	0.56	0.54
0.05	$10^{- 4}$	3.3 $\times 10^{- 3}$	0.55	0.75	0.60	0.97	0.53	0.50	0.49	0.48	1.00	1.00	0.72	0.58	0.98
0.05	$10^{- 3}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.13	0.00	0.00	0.00	0.08	0.00	0.00	0.00	0.00	0.13	0.00
0.05	$10^{- 3}$	1 $\times 10^{- 4}$	0.20	0.77	0.52	0.27	0.20	0.18	0.17	0.16	1.00	0.89	0.72	0.47	0.24
0.05	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.27	0.77	0.54	0.38	0.30	0.27	0.26	0.24	1.00	0.96	0.72	0.51	0.34
0.05	$10^{- 3}$	1 $\times 10^{- 3}$	0.37	0.73	0.56	0.44	0.38	0.36	0.34	0.32	1.00	0.84	0.68	0.54	0.41
0.05	$10^{- 3}$	3.3 $\times 10^{- 3}$	0.57	0.94	0.62	0.54	0.50	0.47	0.46	0.45	0.99	0.98	0.94	0.62	0.53
0.05	$10^{- 2}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.05	$10^{- 2}$	1 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.05	$10^{- 2}$	3.3 $\times 10^{- 4}$	0.21	0.39	0.26	0.21	0.00	0.00	0.00	0.00	0.59	0.52	0.37	0.25	0.19
0.05	$10^{- 2}$	1 $\times 10^{- 3}$	0.30	1.00	0.62	0.38	0.33	0.33	0.31	0.30	1.00	1.00	0.98	0.58	0.35
0.05	$10^{- 2}$	3.3 $\times 10^{- 3}$	0.43	1.00	0.74	0.61	0.53	0.46	0.44	0.42	1.00	1.00	1.00	0.75	0.60
0.07	$10^{- 4}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.11	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.10
0.07	$10^{- 4}$	1 $\times 10^{- 4}$	0.14	0.78	0.68	0.36	0.25	0.18	0.15	0.14	0.96	0.82	0.73	0.61	0.32
0.07	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.33	1.00	0.75	0.63	0.36	0.33	0.31	0.30	1.00	1.00	0.81	0.68	0.36
0.07	$10^{- 4}$	1 $\times 10^{- 3}$	0.42	0.79	0.76	0.48	0.44	0.41	0.39	0.37	1.00	0.82	0.77	0.64	0.44
0.07	$10^{- 4}$	3.3 $\times 10^{- 3}$	0.56	1.00	1.00	0.84	0.58	0.52	0.50	0.48	1.00	1.00	1.00	0.63	0.82
0.07	$10^{- 3}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.07	$10^{- 3}$	1 $\times 10^{- 4}$	0.00	0.00	0.18	0.15	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.17	0.14
0.07	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.28	0.97	0.70	0.43	0.32	0.29	0.27	0.26	1.00	1.00	0.78	0.63	0.38
0.07	$10^{- 3}$	1 $\times 10^{- 3}$	0.36	0.82	0.71	0.53	0.43	0.37	0.35	0.33	1.00	1.00	0.78	0.66	0.49
0.07	$10^{- 3}$	3.3 $\times 10^{- 3}$	0.52	0.97	0.90	0.68	0.57	0.52	0.48	0.47	1.00	1.00	0.97	0.88	0.65
0.07	$10^{- 2}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.07	$10^{- 2}$	1 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.07	$10^{- 2}$	3.3 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.07	$10^{- 2}$	1 $\times 10^{- 3}$	0.28	0.52	0.37	0.15	0.00	0.00	0.00	0.00	0.76	0.65	0.48	0.35	0.13
0.07	$10^{- 2}$	3.3 $\times 10^{- 3}$	0.40	1.00	0.97	0.58	0.46	0.43	0.41	0.41	1.00	1.00	1.00	0.99	0.53
0.10	$10^{- 4}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.10	$10^{- 4}$	1 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.12	0.12	0.11	0.11	0.00	0.00	0.00	0.00	0.10
0.10	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.21	0.83	0.79	0.55	0.43	0.32	0.22	0.19	1.00	0.89	0.80	0.71	0.48
0.10	$10^{- 4}$	1 $\times 10^{- 3}$	0.41	1.00	0.81	0.75	0.50	0.44	0.42	0.40	1.00	1.00	0.84	0.78	0.51
0.10	$10^{- 4}$	3.3 $\times 10^{- 3}$	0.53	0.86	0.84	0.72	0.70	0.54	0.48	0.29	0.75	0.88	0.86	0.72	0.71
0.10	$10^{- 3}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.10	$10^{- 3}$	1 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.10	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.00	0.00	0.00	0.18	0.15	0.15	0.00	0.00	0.00	0.00	0.00	0.00	0.15
0.10	$10^{- 3}$	1 $\times 10^{- 3}$	0.38	1.00	0.82	0.65	0.46	0.40	0.39	0.36	1.00	1.00	1.00	0.76	0.56
0.10	$10^{- 3}$	3.3 $\times 10^{- 3}$	0.49	1.00	0.98	0.74	0.59	0.52	0.49	0.47	1.00	1.00	1.00	0.81	0.68
0.10	$10^{- 2}$	3.3 $\times 10^{- 5}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.10	$10^{- 2}$	1 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.31	0.00	0.00	0.00	0.00
0.10	$10^{- 2}$	3.3 $\times 10^{- 4}$	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
0.10	$10^{- 2}$	1 $\times 10^{- 3}$	0.11	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.20	0.00	0.00	0.00	0.00
0.10	$10^{- 2}$	3.3 $\times 10^{- 3}$	0.38	0.72	0.49	0.40	0.15	0.16	0.00	0.00	1.00	0.88	0.69	0.47	0.38
Kernel	$σ$	$= 0.025 r_{p}$
0.05	$10^{- 4}$	3.3 $\times 10^{- 5}$	0.09	–	–	–	–	0.13	0.12	0.12	–	–	–	–	–
0.05	$10^{- 3}$	3.3 $\times 10^{- 5}$	0.00	–	0.13	0.11	0.11	0.09	0.08	–	–	–	–	0.13	0.11
0.05	$10^{- 2}$	3.3 $\times 10^{- 4}$	0.21	–	–	–	0.20	–	0.11	0.13	–	–	–	–	–
0.07	$10^{- 4}$	3.3 $\times 10^{- 5}$	0.00	–	–	0.10	0.09	0.08	0.09	0.09	–	–	–	0.08	0.09
0.07	$10^{- 3}$	1 $\times 10^{- 4}$	0.00	–	–	0.15	0.14	0.12	0.13	–	–	–	–	0.17	0.14
0.07	$10^{- 2}$	1 $\times 10^{- 3}$	0.28	–	–	–	0.09	0.10	0.09	0.08	–	–	–	–	0.12
0.10	$10^{- 4}$	3.3 $\times 10^{- 5}$	0.00	–	–	0.07	0.08	0.08	0.08	0.08	–	–	–	–	0.07
0.10	$10^{- 4}$	1 $\times 10^{- 4}$	0.00	–	–	0.10	0.12	0.12	0.11	0.12	–	–	–	0.09	0.10
0.10	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.00	–	–	–	–	–	0.15	0.11	–	–	–	0.16	–
0.10	$10^{- 2}$	1 $\times 10^{- 3}$	0.11	–	–	–	0.12	0.08	0.11	–	–	–	–	–	–
0.10	$10^{- 2}$	3.3 $\times 10^{- 3}$	0.38	–	–	–	–	–	0.40	0.13	–	–	–	–	–
Common	Gaps	Separated	by	Horseshoe	$\frac{Δ_{1}}{Δ_{2}}$
0.05	$10^{- 4}$	1 $\times 10^{- 4}$	0.24	–	–	$\frac{0.28}{0.10}$	–	–	–	–	–	–	–	–	–
0.05	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.32	–	–	$\frac{0.32}{0.13}$	–	–	–	–	–	–	–	–	$\frac{0.29}{0.07}$
0.05	$10^{- 4}$	1 $\times 10^{- 3}$	0.42	$\frac{0.10}{0.19}$	$\frac{0.11}{0.20}$	$\frac{0.10}{0.19}$	$\frac{0.25}{0.15}$	$\frac{0.26}{0.15}$	–	–	–	$\frac{0.09}{0.17}$	$\frac{0.10}{0.19}$	$\frac{0.11}{0.20}$	$\frac{0.10}{0.18}$
0.05	$10^{- 4}$	3.3 $\times 10^{- 3}$	0.55	–	–	$\frac{0.56}{0.81}$	–	–	–	–	–	–	–	–	$\frac{0.54}{0.86}$
0.05	$10^{- 3}$	3.3 $\times 10^{- 3}$	0.57	$\frac{0.07}{0.93}$	–	–	–	–	–	–	$\frac{0.09}{0.99}$	$\frac{0.09}{0.97}$	$\frac{0.07}{0.93}$	–	–
0.07	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.33	–	–	$\frac{0.41}{0.19}$	–	–	–	–	–	–	–	–	–
0.07	$10^{- 4}$	1 $\times 10^{- 3}$	0.42	–	–	–	$\frac{0.27}{0.19}$	–	–	–	–	–	–	$\frac{0.51}{0.19}$	$\frac{0.25}{0.17}$
0.07	$10^{- 4}$	3.3 $\times 10^{- 3}$	0.56	–	–	$\frac{0.61}{0.39}$	–	–	–	–	–	–	–	–	$\frac{0.59}{0.27}$
0.07	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.28	$\frac{0.88}{0.74}$	–	–	–	–	–	–	–	–	–	–	–
0.07	$10^{- 2}$	3.3 $\times 10^{- 3}$	0.40	–	$\frac{0.82}{0.62}$	–	–	–	–	–	–	–	–	$\frac{0.81}{0.88}$	–
0.10	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.21	–	–	–	–	–	–	–	–	–	–	$\frac{0.62}{0.18}$	–
0.10	$10^{- 4}$	1 $\times 10^{- 3}$	0.41	–	–	$\frac{0.57}{0.18}$	–	–	–	–	–	–	–	$\frac{0.58}{0.30}$	$\frac{0.27}{0.30}$
0.10	$10^{- 4}$	3.3 $\times 10^{- 3}$	0.53	$\frac{0.16}{0.65}$	$\frac{0.19}{0.66}$	$\frac{0.19}{0.63}$	$\frac{0.16}{0.39}$	$\frac{0.35}{0.13}$	$\frac{0.33}{0.02}$	–	$\frac{0.14}{0.67}$	$\frac{0.15}{0.65}$	$\frac{0.16}{0.65}$	$\frac{0.17}{0.64}$	$\frac{0.16}{0.42}$
0.10	$10^{- 3}$	1 $\times 10^{- 3}$	0.38	–	–	–	–	–	$\frac{0.18}{0.18}$	$\frac{0.17}{0.15}$	–	–	–	–	–
0.10	$10^{- 3}$	3.3 $\times 10^{- 3}$	0.49	–	$\frac{0.85}{0.64}$	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 2}$	1 $\times 10^{- 4}$	0.00	–	–	–	–	–	–	–	$\frac{0.12}{0.12}$	–	–	–	–

Table 5. Table 5: Gap Depths ( l o g 10 ( δ − 1 ) ) 𝑙 𝑜 subscript 𝑔 10 𝛿 1 \big{(}log_{10}(\delta-1)\big{)} for the Gas, DSD1 and DSD2

h/r	$α$	$q$	$δ_{g} - 1$	$δ_{d, 0 p 1} - 1$	$δ_{d, 0 p 3} - 1$	$δ_{d, 1} - 1$	$δ_{d, 3} - 1$	$δ_{d, 10} - 1$	$δ_{d, 30} - 1$	$δ_{d, 100} - 1$	$δ_{d, 1} - 1$	$δ_{d, 3} - 1$	$δ_{d, 10} - 1$	$δ_{d, 30} - 1$	$δ_{d, 100} - 1$
		$(M_{p} / M_{*})$	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )	( $l o g_{10}$ )
0.05	$10^{- 4}$	3.3 $\times 10^{- 5}$	-1.03	1.94	1.25	0.62	0.24	-0.03	-0.21	-0.41	3.66	2.93	2.10	1.25	0.50
0.05	$10^{- 4}$	1 $\times 10^{- 4}$	0.34	1.66	1.05	0.57	0.39	0.29	0.23	0.10	3.67	2.62	1.60	0.87	0.32
0.05	$10^{- 4}$	3.3 $\times 10^{- 4}$	1.15	1.56	1.12	0.82	0.64	0.57	0.53	0.53	2.70	2.10	1.54	0.93	0.54
0.05	$10^{- 4}$	1 $\times 10^{- 3}$	3.43	3.12	2.56	2.14	1.90	1.72	1.65	1.56	3.66	3.47	3.13	2.64	2.05
0.05	$10^{- 4}$	3.3 $\times 10^{- 3}$	2.05	4.55	3.90	3.20	2.67	2.26	1.93	1.73	8.75	8.00	6.93	5.75	4.10
0.05	$10^{- 3}$	3.3 $\times 10^{- 5}$	–	–	-0.55	-0.73	–	–	–	–	–	–	-0.77	-0.61	-0.81
0.05	$10^{- 3}$	1 $\times 10^{- 4}$	-0.27	2.24	1.46	0.72	0.24	0.01	-0.10	-0.26	4.94	3.83	2.55	1.46	0.57
0.05	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.81	3.03	2.26	1.59	1.17	0.82	0.51	0.37	6.67	5.69	4.37	3.03	1.84
0.05	$10^{- 3}$	1 $\times 10^{- 3}$	2.35	2.54	3.08	2.86	2.49	2.17	1.89	1.49	3.40	3.08	2.72	3.23	3.78
0.05	$10^{- 3}$	3.3 $\times 10^{- 3}$	1.96	3.52	2.95	2.43	2.10	1.78	1.48	1.34	–	–	5.44	4.22	3.13
0.05	$10^{- 2}$	3.3 $\times 10^{- 5}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.05	$10^{- 2}$	1 $\times 10^{- 4}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.05	$10^{- 2}$	3.3 $\times 10^{- 4}$	-0.32	0.81	0.16	-0.25	-0.72	–	-0.72	–	1.45	1.20	0.75	0.11	-0.41
0.05	$10^{- 2}$	1 $\times 10^{- 3}$	0.64	1.88	1.24	0.59	0.24	0.04	-0.01	-0.10	–	3.37	2.22	1.27	0.44
0.05	$10^{- 2}$	3.3 $\times 10^{- 3}$	2.69	3.05	2.50	2.01	1.69	1.35	1.05	0.74	4.87	4.12	5.00	3.85	2.60
0.07	$10^{- 4}$	3.3 $\times 10^{- 5}$	–	–	–	-0.49	-0.55	-0.58	-0.65	–	–	–	–	–	-0.61
0.07	$10^{- 4}$	1 $\times 10^{- 4}$	-0.69	2.31	1.92	1.08	0.51	0.15	-0.04	-0.28	4.06	3.26	2.41	1.73	0.88
0.07	$10^{- 4}$	3.3 $\times 10^{- 4}$	0.59	2.13	1.88	1.32	0.93	0.60	0.39	0.25	3.09	2.32	1.94	1.54	1.03
0.07	$10^{- 4}$	1 $\times 10^{- 3}$	1.86	3.34	2.93	2.20	1.81	1.51	1.29	1.02	4.57	3.89	3.44	2.92	2.20
0.07	$10^{- 4}$	3.3 $\times 10^{- 3}$	2.60	5.10	4.63	3.79	3.24	2.75	2.42	2.04	9.89	9.11	8.10	7.02	5.69
0.07	$10^{- 3}$	3.3 $\times 10^{- 5}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.07	$10^{- 3}$	1 $\times 10^{- 4}$	–	–	-0.47	-0.39	-0.63	-0.65	–	–	–	–	-0.56	-0.55	-0.54
0.07	$10^{- 3}$	3.3 $\times 10^{- 4}$	0.03	2.85	2.29	1.38	0.83	0.48	0.28	0.09	6.25	5.40	3.86	2.53	1.33
0.07	$10^{- 3}$	1 $\times 10^{- 3}$	1.04	3.50	2.96	2.12	1.64	1.26	0.93	0.72	7.71	6.91	5.78	4.49	2.86
0.07	$10^{- 3}$	3.3 $\times 10^{- 3}$	2.38	4.86	4.36	3.58	3.12	2.78	2.48	2.11	4.69	–	7.61	6.49	4.89
0.07	$10^{- 2}$	3.3 $\times 10^{- 5}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.07	$10^{- 2}$	1 $\times 10^{- 4}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.07	$10^{- 2}$	3.3 $\times 10^{- 4}$	–	–	–	–	–	–	–	–	-0.45	–	–	–	–
0.07	$10^{- 2}$	1 $\times 10^{- 3}$	-0.09	0.84	0.18	-0.37	-0.71	–	–	–	1.84	1.46	0.79	0.10	-0.49
0.07	$10^{- 2}$	3.3 $\times 10^{- 3}$	1.08	2.09	1.63	1.01	0.58	0.28	0.16	0.11	–	–	2.99	1.84	0.87
0.10	$10^{- 4}$	3.3 $\times 10^{- 5}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 4}$	1 $\times 10^{- 4}$	–	–	–	-0.63	-0.29	-0.25	-0.38	-0.57	–	–	–	–	-0.80
0.10	$10^{- 4}$	3.3 $\times 10^{- 4}$	-0.33	2.51	2.23	1.49	0.73	0.26	0.00	-0.27	4.10	3.50	3.00	2.20	1.15
0.10	$10^{- 4}$	1 $\times 10^{- 3}$	0.79	2.75	2.42	1.76	1.18	0.71	0.45	0.21	4.55	3.74	3.13	2.40	1.46
0.10	$10^{- 4}$	3.3 $\times 10^{- 3}$	1.84	3.71	3.08	2.81	2.36	1.87	1.52	1.14	5.84	4.47	3.77	3.37	4.29
0.10	$10^{- 3}$	3.3 $\times 10^{- 5}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 3}$	1 $\times 10^{- 4}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 3}$	3.3 $\times 10^{- 4}$	–	–	–	-0.41	-0.58	-0.57	-0.62	–	–	–	–	-0.67	-0.67
0.10	$10^{- 3}$	1 $\times 10^{- 3}$	0.17	2.54	2.10	1.33	0.59	0.08	-0.18	-0.34	–	5.06	3.52	2.29	1.14
0.10	$10^{- 3}$	3.3 $\times 10^{- 3}$	1.33	3.19	2.74	2.04	1.40	0.97	0.63	0.26	–	6.64	5.63	4.52	2.92
0.10	$10^{- 2}$	3.3 $\times 10^{- 5}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 2}$	1 $\times 10^{- 4}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 2}$	3.3 $\times 10^{- 4}$	–	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 2}$	1 $\times 10^{- 3}$	-1.45	–	–	–	–	–	–	–	–	–	–	–	–
0.10	$10^{- 2}$	3.3 $\times 10^{- 3}$	0.17	1.24	0.55	-0.07	-0.61	-0.62	-0.62	-0.64	2.87	2.14	1.30	0.48	-0.25

Equations56

Σ_{g} (r) = Σ_{g, 0} (r / r_{0})^{- 1},

Σ_{g} (r) = Σ_{g, 0} (r / r_{0})^{- 1},

S t = t_{s t o p} Ω = \frac{π s ρ _{p}}{2 Σ _{g a s}} = 1.57 \times 1 0^{- 3} \frac{ρ _{p}}{1 g c m ^{- 3}} \frac{s}{1 mm} \frac{100 g c m ^{- 2}}{Σ _{g}} .

S t = t_{s t o p} Ω = \frac{π s ρ _{p}}{2 Σ _{g a s}} = 1.57 \times 1 0^{- 3} \frac{ρ _{p}}{1 g c m ^{- 3}} \frac{s}{1 mm} \frac{100 g c m ^{- 2}}{Σ _{g}} .

f_{sd}=\int_{min\{s_{min},s_{min}^{code}\}}^{min\{s_{max},s_{min}^{code}\}}s^{3+p}ds\bigg{/}\int_{s_{min}}^{s_{max}}s^{3+p}ds\,,

f_{sd}=\int_{min\{s_{min},s_{min}^{code}\}}^{min\{s_{max},s_{min}^{code}\}}s^{3+p}ds\bigg{/}\int_{s_{min}}^{s_{max}}s^{3+p}ds\,,

f_{l d} = 1 - f_{s d} .

f_{l d} = 1 - f_{s d} .

w_{i} (s) = \frac{M _{t o t}}{N _{p a r t}} \frac{s ^{3 + p} / \int _{ma x {s_{min}, s_{min}^{co d e}}}^{s_{ma x}} s ^{3 + p} d s}{s ^{- 1} / \int _{s_{min}^{co d e}}^{s_{ma x}^{co d e}} s ^{- 1} d s},

w_{i} (s) = \frac{M _{t o t}}{N _{p a r t}} \frac{s ^{3 + p} / \int _{ma x {s_{min}, s_{min}^{co d e}}}^{s_{ma x}} s ^{3 + p} d s}{s ^{- 1} / \int _{s_{min}^{co d e}}^{s_{ma x}^{co d e}} s ^{- 1} d s},

τ_{l d} = f_{l d} \frac{\sum _{i} w _{i} ( s ) κ _{ab s, i} ( s )}{A _{ce l l}},

τ_{l d} = f_{l d} \frac{\sum _{i} w _{i} ( s ) κ _{ab s, i} ( s )}{A _{ce l l}},

τ_{s d} = f_{s d} κ_{ma} Σ_{g} /100

τ_{s d} = f_{s d} κ_{ma} Σ_{g} /100

κ_{ma} = \int_{min {s_{min}, s_{min}^{co d e}}}^{min {s_{ma x}, s_{min}^{co d e}}} κ_{ab s} (s) s^{3 + p} d s

κ_{ma} = \int_{min {s_{min}, s_{min}^{co d e}}}^{min {s_{ma x}, s_{min}^{co d e}}} κ_{ab s} (s) s^{3 + p} d s

τ (r, θ) = τ_{s d} (r, θ) + τ_{l d} (r, θ) .

τ (r, θ) = τ_{s d} (r, θ) + τ_{l d} (r, θ) .

T_{b} (r, θ) = T_{d} (r) (1 - e^{- τ (r, θ)}),

T_{b} (r, θ) = T_{d} (r) (1 - e^{- τ (r, θ)}),

T_{d} (r) = T_{d} (r_{0}) (\frac{r}{r _{0}})^{- 0.5} .

T_{d} (r) = T_{d} (r_{0}) (\frac{r}{r _{0}})^{- 0.5} .

T_{d}(r_{0})=\Big{(}\frac{\phi L_{*}}{8\pi r_{0}^{2}\sigma_{SB}}\Big{)}^{1/4}

T_{d}(r_{0})=\Big{(}\frac{\phi L_{*}}{8\pi r_{0}^{2}\sigma_{SB}}\Big{)}^{1/4}

\frac{v _{ϕ}^{2}}{r} = \frac{v _{K}^{2}}{r} + \frac{1}{ρ _{g a s}} \frac{\partial P}{\partial r} .

\frac{v _{ϕ}^{2}}{r} = \frac{v _{K}^{2}}{r} + \frac{1}{ρ _{g a s}} \frac{\partial P}{\partial r} .

\frac{Δ v _{ϕ}}{v _{K}} \sim \frac{r}{2 ρ _{g a s} v _{K}^{2}} \frac{\partial P}{\partial r},

\frac{Δ v _{ϕ}}{v _{K}} \sim \frac{r}{2 ρ _{g a s} v _{K}^{2}} \frac{\partial P}{\partial r},

K_{v_{r}} = q (h / r)^{p h} α^{p a} .

K_{v_{r}} = q (h / r)^{p h} α^{p a} .

Δ δ v_{r o t} = A K_{v_{r}}^{B} .

Δ δ v_{r o t} = A K_{v_{r}}^{B} .

K_{v_{r}} = q (h / r)^{- 1.27} α^{- 0.41}

K_{v_{r}} = q (h / r)^{- 1.27} α^{- 0.41}

Δ δ v_{r o t} = 0.11 K_{v_{r}}^{0.80} .

Δ δ v_{r o t} = 0.11 K_{v_{r}}^{0.80} .

δ_{Σ} = Σ (r_{p e ak}) /Σ (r_{g a p}),

δ_{Σ} = Σ (r_{p e ak}) /Σ (r_{g a p}),

δ_{I} = I_{mm} (r_{p e ak}) / I_{mm} (r_{g a p}),

δ_{I} = I_{mm} (r_{p e ak}) / I_{mm} (r_{g a p}),

Σ_{e d g e} = \frac{Σ ( r _{p e ak} ) + Σ ( r _{g a p} )}{2},

Σ_{e d g e} = \frac{Σ ( r _{p e ak} ) + Σ ( r _{g a p} )}{2},

I_{e d g e} = \frac{I _{mm} ( r _{p e ak} ) + I _{mm} ( r _{g a p} )}{2} .

I_{e d g e} = \frac{I _{mm} ( r _{p e ak} ) + I _{mm} ( r _{g a p} )}{2} .

Δ = (r_{o u t} - r_{in}) / r_{o u t} .

Δ = (r_{o u t} - r_{in}) / r_{o u t} .

\frac{K^{\prime}}{0.014}=\frac{q}{0.001}\Big{(}\frac{h/r}{0.07}\Big{)}^{-0.18}\Big{(}\frac{\alpha}{10^{-3}}\Big{)}^{-0.31}\,.

\frac{K^{\prime}}{0.014}=\frac{q}{0.001}\Big{(}\frac{h/r}{0.07}\Big{)}^{-0.18}\Big{(}\frac{\alpha}{10^{-3}}\Big{)}^{-0.31}\,.

δ - 1 = C K^{D},

δ - 1 = C K^{D},

\frac{K}{24}=\frac{q}{0.001}\Big{(}\frac{h/r}{0.07}\Big{)}^{-2.81}\Big{(}\frac{\alpha}{10^{-3}}\Big{)}^{-0.38}\,.

\frac{K}{24}=\frac{q}{0.001}\Big{(}\frac{h/r}{0.07}\Big{)}^{-2.81}\Big{(}\frac{\alpha}{10^{-3}}\Big{)}^{-0.38}\,.

1 - \frac{r _{I G 2}}{r _{p}} = 2.3 q^{0.02} (h / r)^{0.58} α^{- 0.01} .

1 - \frac{r _{I G 2}}{r _{p}} = 2.3 q^{0.02} (h / r)^{0.58} α^{- 0.01} .

Σ_{g} = 1700 (\frac{r}{a u})^{- 1.5} g c m^{- 2} .

Σ_{g} = 1700 (\frac{r}{a u})^{- 1.5} g c m^{- 2} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

The Disk Substructures at High Angular Resolution Project (DSHARP):

VII. The Planet-Disk Interactions Interpretation

Shangjia Zhang

Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S. Maryland Pkwy, Las Vegas, NV, 89154, USA

Zhaohuan Zhu

Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S. Maryland Pkwy, Las Vegas, NV, 89154, USA

Jane Huang

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

Viviana V. Guzmán

Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile

Sean M. Andrews

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

Tilman Birnstiel

University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 Munich, Germany

Cornelis P. Dullemond

Zentrum für Astronomie, Heidelberg University, Albert Ueberle Str. 2, 69120 Heidelberg, Germany

John M. Carpenter

Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

Andrea Isella

Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA

Laura M. Pérez

Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

Myriam Benisty

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

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

David J. Wilner

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

Clément Baruteau

CNRS / Institut de Recherche en Astrophysique et Planétologie, 14 avenue Edouard Belin, 31400 Toulouse, France

Xue-Ning Bai

Institute for Advanced Study and Tsinghua Center for Astrophysics, Tsinghua University, Beijing 100084, China

Luca Ricci

Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff Street, Northridge, CA 91130, USA

Abstract

The Disk Substructures at High Angular Resolution Project (DSHARP) provides a large sample of protoplanetary disks having substructures which 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 2-D hydrodynamical simulations including both gas and dust components. We present the resulting gas structures, including the relationship between the planet mass and 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 and 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 probe a wide-orbit planet population (Neptune to Jupiter mass planets beyond 10 au) that is not accessible to other planet searching techniques.

hydrodynamics — planetary systems: protoplanetary disks — planet-disk interactions — submillimeter: planetary systems —

††software: Astropy (Astropy Collaboration et al., 2013), CASA (McMullin et al., 2007), Dusty FARGO-ADSG (Baruteau & Masset, 2008a, b; Baruteau & Zhu, 2016), Matplotlib (Hunter, 2007), Numpy (Van Der Walt et al., 2011), Scipy (Jones et al., 2001–)

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. review by Johansen et al. 2014; Raymond et al. 2014; Chabrier 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 Myrs old) mainly due to their stellar variablity and the presence of the protoplanetary disks. Fewer than 10 young planet candidates in systems $<$ 10 Myrs 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-meter class telescopes) and interferometry at radio wavelengths (especially the ALMA and the VLA) can directly probe the protoplanetary disks down to au-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; Baruteau et al., 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; Gonzalez et al., 2012; Pinilla et al., 2012; 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; Zhu et al. 2015a; Fung & Dong 2015), (sub-)mm dust thermal continuum images (Dipierro et al., 2015; Picogna & Kley, 2015; Dong & Fung, 2017; Dong et al., 2018a), and (sub-)mm 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-)mm wavelengths allows us to probe low mass planets, since 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-)mm 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 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 which 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 §2, we describe our methods. The results are presented in §3. The derived young planet properties for the DSHARP disks are given in §4. After a short discussion in §5, we conclude the paper in §6.

2 Method

We carry out 2-D hydrodynamical planet-disk simulations using the modified version of the grid-based code FARGO (Masset, 2000) called Dusty FARGO-ADSG (Baruteau & Masset, 2008a, b; Baruteau & Zhu, 2016). The gas component is simulated using finite difference methods (Stone & Norman, 1992), while the dust component is modelled 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 §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 Yang & Zhu (2018).

2.1 Setup: Gas and Dust

We adopt polar coordinates ( $r$ , $\theta$ ) 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 non-inertial coordinate frame.

We initialize the gas surface density as

[TABLE]

where $r_{0}$ is also the position of the planet and we set $r_{0}=r_{p}=1$ . For studying gaps of individual sources in §4, we scale $\Sigma_{g,0}$ to be consistent with the DSHARP observations. We assume locally isothermal equation of state, and the temperature at radius $r$ follows $T(r)=T_{0}(r/r_{0})^{-1/2}$ . $T$ is related to the disk scale height $h$ as $h/r=c_{s}/v_{\phi}$ where $c_{s}^{2}=RT/\mu=P/\Sigma$ and $\mu$ = 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 $\pi$ in the $\theta$ direction. For low viscosity cases ( $\alpha$ = $10^{-4}$ and $10^{-3}$ ), there are 750 grid points in the radial direction and 1024 grid points in the $\theta$ direction. This is equivalent to 16 grid points per scale height at $r_{0}$ if $h/r=0.1$ . For high viscosity cases ( $\alpha$ =0.01), less resolution is needed so there are 375 and 512 grid points in the radial and $\theta$ direction. For simulations to fit AS 209 in §4.1.1, the resolution is 1500 and 2048 grid points in the radial and $\theta$ 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.12 $r_{0}$ and $r>$ 8 $r_{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 10s of au, 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

[TABLE]

where $\rho_{p}$ is the density of the dust particle, $s$ is the radius of the dust particle, and $\Sigma_{g}$ is the gas surface density. We assume $\rho_{p}$ =1 g cm*-3* in our simulations. We use 200,000 and 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 $\times 10^{-5}$ to 1.57, or physical radii ranging from 1 $\mu$ m to 10 cm if $\Sigma_{g,0}$ =10 $\mathrm{g\,cm^{-2}}$ and $\rho_{p}=1\mathrm{g\,cm^{-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 the every particle varies as being inverse 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 $\alpha$ 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 ( $\nu$ ) 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 $(h/r)_{r_{0}}$ $(0.05,0.07,0.1)$ , five values for the planet-star mass ratio (q $\equiv M_{p}/M_{*}$ = 3.3 $\times$ 10*-5*, 10*-4*, 3.3 $\times$ 10*-4*, 10*-3*, 3.3 $\times$ 10*-3* M*∗*, or roughly $M_{p}$ = 11 $M_{\earth}$ , 33 $M_{\earth}$ , 0.35 $M_{J}$ , 1 $M_{J}$ , 3.5 $M_{J}$ if $M_{*}=M_{\odot}$ ), and three values for the disk turbulent viscosity coefficient ( $\alpha=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, $\alpha=10^{-3}$ (m3 in h5am3p1 means minus 3), $M_{p}/M_{*}=3.3\times 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 §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\,\mathrm{au}$ (D’Alessio et al., 1998). While a moderate $\alpha\sim 10^{-2}$ is preferred to explain the disk accretion (Hartmann et al., 1998), recent works suggest that a low turbulence level ( $\alpha<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 $\alpha$ is smaller than $10^{-4}$ , the viscous timescale over the disk scale height at the planet position ( $H_{p}^{2}/\nu$ ) is longer than $10^{4}/\Omega_{p}$ or 1.6 million years at 100 au, so that the viscosity will not affect the disk evolution significantly. In §4.1, we carry out several simulations with different $\alpha$ 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_{\earth}$ , 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 mm Continuum Intensity Maps

For each simulation, we calculate the mm 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 mm 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 mm are very similar to the gaps in the gas. If the mm emission is dominated by mm 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 mm intensity maps, we choose two very different dust size distributions to generate intensity maps. For the distribution referred to as DSD1, we assume $n(s)\propto 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-)mm 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)\propto s^{-2.5}$ with the maximum grain size of 1 cm ( $p=-2.5$ and $s_{max}$ =1 cm). 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 mm/cm wavelengths (Ricci et al., 2010b, a; Pérez et al., 2015). Both cases assume a minimum grain size of 0.005 $\mu$ 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 $\mathrm{cm^{2}\,g^{-1}}$ , while for the $s_{max}=$ 1 cm case it is 0.46 $\mathrm{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 §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 $\Sigma_{g,0}$ = 0.1 $\mathrm{g\,cm^{-2}}$ , 0.3 $\mathrm{g\,cm^{-2}}$ , 1 $\mathrm{g\,cm^{-2}}$ , 3 $\mathrm{g\,cm^{-2}}$ , 10 $\mathrm{g\,cm^{-2}}$ , 30 $\mathrm{g\,cm^{-2}}$ , and 100 $\mathrm{g\,cm^{-2}}$ (seven groups of models). The maximum-size particle in these disks (0.1 mm), which dominates the total dust mass, corresponds to $St$ = 1.57 $\times 10^{-1}$ , 5.23 $\times 10^{-2}$ , 1.57 $\times 10^{-2}$ , 5.23 $\times 10^{-3}$ , 1.57 $\times 10^{-3}$ , 5.23 $\times 10^{-4}$ , and 1.57 $\times 10^{-4}$ at $r=r_{p}$ . For the $s_{max}=$ 1 cm cases (DSD2), we vary $\Sigma_{g,0}$ as 1 $\mathrm{g\,cm^{-2}}$ , 3 $\mathrm{g\,cm^{-2}}$ , 10 $\mathrm{g\,cm^{-2}}$ , 30 $\mathrm{g\,cm^{-2}}$ , and 100 $\mathrm{g\,cm^{-2}}$ (five groups of models), and the corresponding $St$ for 1 cm particles at $r=r_{p}$ is 1.57, 5.23 $\times 10^{-1}$ , 1.57 $\times 10^{-1}$ , 5.23 $\times 10^{-2}$ , and 1.57 $\times 10^{-2}$ . For each given surface density above, we only select particles with Stokes numbers smaller than the corresponding $St$ in our simulations and use the distribution of these particles to calculate the 1.27 mm intensity maps. For the $s_{max}=$ 1 cm dust distribution (DSD2), we do not have $\Sigma_{g,0}$ = 0.1 $\mathrm{g\,cm^{-2}}$ , 0.3 $\mathrm{g\,cm^{-2}}$ cases since 1 cm particles in these disks have Stokes numbers larger than the largest Stokes number (1.57) in our simulations.

Here, we lay out the detailed steps to scale each simulation to the disks that have surface densities of $\Sigma_{g,0}$ listed above, and then calculate the mm intensity maps for these disks.

First, given a $\Sigma_{g,0}$ , we find the relationship between the particle size in this disk and the Stokes number of super-particles in simulations. For each particle in the simulation, we use its Stokes number in the initial condition to calculate the corresponding particle size $s$ (Equation 2 with known $\Sigma_{g}$ ). The Stokes number of test particles at $r=r_{p}$ in the initial condition ranges from $St_{min}$ = 1.57 $\times$ $10^{-5}$ to $St_{max}$ = 1.57, or in terms of grain size, $s^{code}_{min}$ = $St_{min}\times 2\Sigma_{gas}/(\pi\rho_{p})$ and $s^{code}_{max}$ = $St_{max}\times 2\Sigma_{gas}/(\pi\rho_{p})$ from Equation 2. For instance, a 1 $\mu$ m particle in a disk with $\Sigma_{g}$ = $10\,\mathrm{g\,cm^{-2}}$ at the planet position corresponds to the particle with $St$ = 1.57 $\times$ $10^{-5}$ at $r=r_{p}$ in the initial setup of the simulation. For dust grains with $St$ $<$ $St_{min}$ = 1.57 $\times$ $10^{-5}$ , we use the gas surface density $\Sigma_{g}(r,\theta)$ in our simulations to represent the dust, assuming small dust grains are well coupled with the gas.
Then, with a given $\Sigma_{g,0}$ , we use the assumed particle size distributions (DSD1 and DSD2) in the initial condition to calculate the mass weight for each super-particle in the simulation. Note that during the simulation, the resulting dust size distribution at each radius is different from the initial dust size distribution since particles drift in the disk. As mentioned above, we divide the dust component in the disk into two parts: (a) the small dust particles ( $s$ $<$ $s_{min}^{code}$ ) represented by the gas component in the simulation and (b) large dust particles ( $s$ $\geq$ $s_{min}^{code}$ ) represented by the super-particles in the simulation. We calculate the initial mass fractions of the dust contributed by part (a) and (b). The mass fraction of small particles (part a) with respect to the total dust mass is

[TABLE]

and the mass fraction of large particles using dust super-particles (part b) is

[TABLE]

We want to explore two dust size distributions $n(s)$ $\propto$ $s^{-3.5}$ and $s^{-2.5}$ , given the minimum and maximum dust size $s_{min}$ and $s_{max}$ . However, the super-particles in our setup have a different distribution. The number of super-particles $N(s)$ follows a uniform distribution in the $log(s)$ space, $\int N(s)ds$ $\propto$ $dlog(s)$ . Thus for dust in part (b), if $f_{ld}$ $>$ 0, we give each particle (having size $s$ ) a mass weight to scale them into the desired distribution:

[TABLE]

where $M_{tot}$ is the total dust mass in the disk and $N_{part}$ is the total number of super-particles in the simulation.

Next, we assign the opacity for each particle to derive the total optical depth. DSHARP opacities are produced by Birnstiel et al. (2018), which contains a table of absorption and scattering opacities for a given wavelength and grain size, $\kappa(\lambda,s)$ . For part (b) dust component, we assign each particle a DSHARP absorption opacity $\kappa_{abs,i}(s_{i})$ at 1.27 mm based on the particle’s size, where $s_{i}$ is the $s$ value in the table that is the closest to this particle size. If the particle size is smaller than the minimum size in the opacity table, we take the opacity for the minimum sized particle in the table, namely using a constant extrapolation, since the opacity is already independent of the particle size at the lower size end of the opacity table. We bin all super-particles in each numerical grid cell to derive the total optical depth through the disk for particles in part (b):

[TABLE]

where the sum is adding all particles in the cell, and $A_{cell}$ is the surface area of the grid cell. The optical depth contributed by part (a) is simply

[TABLE]

where

[TABLE]

is the mass-averaged opacity of the small dust within the range of dust sizes in part (a). The final optical depth for each grid cell at ( $r,\theta$ ) is the sum of both components,

[TABLE]

Note that we do not consider dust and gas within one Hill radius $r_{H}$ around the planet for our analysis since our simulations are not able to resolve the circumplanetary region. Thus, we impose the optical depth there to be the minimum optical depth within the annulus ( $r_{0}$ - $r_{H}$ ) $<$ r $<$ ( $r_{0}$ + $r_{H}$ ).

Then, we calculate the brightness temperature or intensity for each grid cell as

[TABLE]

and we assume that the midplane dust temperature follows the assumed disk temperature. Thus,

[TABLE]

Because we seek to derive a scale-free intensity for different systems, the Rayleigh-Jeans approximation is made here. For the young solar system and the HR 8799 calculations in §5.1, and the detailed modeling of AS 209 and Elias 24 in §4.1.1 and §4.1.2, we use the full Planck function at $\nu$ = 240 GHz to derive more accurate intensities.

The normalized brightness temperature ( $T_{b}(r,\theta)/T_{d}(r_{0})$ ) is adequate for the gap width and depth calculation in §3.2.2. But for individual sources, we would like to calculate the absolute brightness temperature. Then, we need to multiply the normalized brightness temperature by the disk temperature at $r_{0}$ ( $T_{d}(r_{0})$ ). We estimate $T_{d}(r_{0})$ using

[TABLE]

where $L_{*}$ is the stellar luminosity and $\phi$ is a constant of 0.02 coming from an estimate from Figure 3 in D’Alessio et al. (2001). This disk mid-plane temperature is the same as Equation 5 in Dullemond et al. (2018), and more details can be found there. We calculate $T_{d}(r)$ for each DSHARP source using the stellar properties ( $L_{*}$ ) listed in Andrews et al. (2018). Knowing $T_{d}(r)$ , we can simply derive $h/r$ at the gap position using $h/r$ = $c_{s}$ / $v_{\phi}$ (the $M_{*}$ that is used to calculate $v_{\phi}$ is also given in Andrews et al. (2018).).

Finally, we convolve these intensity maps with two different Gaussian beams. The beam size is $\sigma=0.06r_{p}$ and $\sigma=0.025r_{p}$ respectively. For a protoplanetary disk 140 pc away, this is equivalent to FWHM (Full Width Half Maximum) beam size of 0.1” and 0.043” if $r_{p}=100$ au, or 0.05” and 0.021” if $r_{p}=50$ au.

3 Simulation Results

3.1 Gas

We will first present results for the gas component in the simulations, including gaseous gap profiles (§3.1.1) and the sub/super Keplerian gas motion at the gap edges (§3.1.2).

3.1.1 Density

Figure 1 shows the two-dimensional gas density maps for all the simulations at 1000 planetary orbits. The left, middle, and right panel blocks show simulations with $h/r$ =0.05, 0.07, and 0.1. Within each panel block, $\alpha=10^{-4}$ , $10^{-3}$ , and $10^{-2}$ cases are shown from left to right. Some large-scale azimuthal structures are evident in the figure. First, low $\alpha$ disks exhibit noticeable horseshoe material within the gap. Since the planet is at ( $x$ =1, $y$ =0) and orbiting around the star in the counterclockwise direction, most horseshoe material is trapped behind the planet (around the L5 point). This is consistent with the shape of the horseshoe streamlines around a non-migrating planet in a viscous disk (Masset, 2002). Second, the gap edge becomes more eccentric and off-centered for smaller $h/r$ , smaller $\alpha$ and larger planet mass cases (especially for $M_{p}\geq$ 3 $M_{J}$ ). Such an eccentric gap edge for the $M_{p}\geq$ 3 $M_{J}$ planet is consistent with previous studies (Lubow, 1991a, b; Kley & Dirksen, 2006; Teyssandier & Ogilvie, 2017). Third, large-scale vortices can be seen at the gap edges for some of the $\alpha=10^{-4}$ cases. Although they are not very apparent in the gas surface density maps, they can trap dust particles azimuthally, causing a large azimuthal contrast in the dust continuum images (as shown in §3.2).

The azimuthally averaged gas surface density profiles for all the models are shown in Figure 2. Several noticeable trends in this figure are:

When the planet mass increases, the gap depth normally increases. However, when the gap is very eccentric (e.g. h5am4p5, h5am3p5), the azimuthally averaged gas surface density at the gap is actually higher than the cases with lower mass planets. This is because azimuthal averaging over an elliptical gap smears out the gap density profile.
With the same planet mass, gaps in $h/r=0.1$ cases are shallower but wider than the $h/r=0.05$ cases. This is consistent with previous studies (Fung et al., 2014; Kanagawa et al., 2015, 2016).
For a given planet mass and $h/r$ , the gaps are shallower and smoother with increasing $\alpha$ . With $\alpha=0.01$ and $10^{-3}$ , the gap edge is smooth, and there is only a single gap at $r/r_{p}\sim 1$ . With $\alpha=10^{-4}$ , there are clearly two shoulders at two edges of the gap, and the material in the horseshoe region still remains in some cases. Especially, for low mass planets in $\alpha=10^{-4}$ disks, the gap at $r/r_{p}\sim 1$ appears to split into two adjacent gaps. This is consistent with non-linear wave steepening theory (Goodman & Rafikov, 2001; Muto et al., 2010; Dong et al., 2011; Duffell & MacFadyen, 2012; Zhu et al., 2013) which suggests that the waves launched by a low mass planet in an inviscid disk need to propagate for some distance to shock and open gaps, leaving the horseshoe region untouched.
For $\alpha=10^{-4}$ cases, we see secondary gaps at $r/r_{p}\sim 0.6$ in $h/r$ =0.05 disks, $r/r_{p}\sim 0.5$ in $h/r=0.07$ disks, and $r/r_{p}\sim 0.4$ in $h/r=0.1$ disks. For some cases, we can even see tertiary gaps at smaller radii. These are consistent with simulations by Bae et al. (2017); Dong et al. (2017) and these gaps are due to the formation of shocks from the secondary and tertiary spirals (Bae & Zhu, 2018a, b).

3.1.2 Kinematics Across the Gap

Recent works by Teague et al. (2018) and Pinte et al. (2018) have shown that, using molecular lines, ALMA can detect the velocity deviation from Keplerian rotation in protoplanetary disks. Such deviations are caused by the radial pressure gradient at the gaseous gap edges,

[TABLE]

In our 2-D simulations, $v_{K}$ is simply $\sqrt{GM_{*}/r}$ and $P$ is $\Sigma c_{s}^{2}$ . Equation 13 suggests that the deviation from the Keplerian motion is

[TABLE]

where $\Delta v_{\phi}=v_{\phi}-v_{K}$ . In a smooth disk where $\partial P/\partial r\sim P/r$ , this deviation is very small, on the order of $(h/r)^{2}$ or $\leq$ 1 $\%$ in a typical protoplanetary disk. But if the gaseous disk has a sharp pressure transition (e.g. at gap edges), the deviation from the Keplerian rotation can be significantly larger. In Figure 3, we plot the azimuthally averaged $\delta v_{rot}\equiv(v_{\phi}-v_{K})/v_{K}$ and $\Sigma$ in run h5am4p4. The directly measured $\delta v_{rot}$ is plotted as the orange curve in the upper panel, while the calculated $\delta v_{rot}$ using the disk surface density profile (presented in the lower panel) and Equation 13 is plotted as the blue curve in the upper panel. We can see that Equation 13 reproduces the measured azimuthal velocity very well, confirming that the sub/super Keplerian motion is due to the radial pressure gradient.

Figure 4 shows $\delta v_{rot}$ for all our cases. As expected, when the gap is deeper due to either smaller $\alpha$ , smaller $h/r$ , or a more massive planet, the amplitude of $\delta v_{rot}$ is larger. However, when the gap becomes very eccentric and off centered (e.g. h5am4p5, h5am3p5), the azimuthally averaged $\delta v_{rot}$ shows a much wider outer bump, indicating an eccentric outer disk. We label these cases as dashed curves in Figure 4 and unfilled markers in panel a of Figure 5. Another interesting feature shown in Figure 4 is that the presence of the gap edge vortices in $\alpha=10^{-4}$ cases does not affect the azimuthally averaged $\delta v_{rot}$ very much. They look similar to the larger $\alpha$ cases without vortices. We interpret this as: if the vortex is strong with fast rotation, it has a smaller aspect ratio so that it is physically small (Lyra & Lin, 2013) and contributes little to the azimuthally averaged gas velocity profile; and if the vortex is weak, although it has a wider azimuthal extent its rotation is small compared with the background shear so again it contributes little to the global velocity profile.

The radial distance and amplitude of the sub/super Keplerian peaks are plotted in panel a of Figure 5. $\Delta\delta v_{rot}$ is the difference between the maximum $\delta v_{rot}$ (at $r>r_{p}$ ) and the minimum $\delta v_{rot}$ (at $r<r_{p}$ ) from Figure 4. Note that these velocity peaks are not peaks (or rings) at mm intensity images that will be presented in §3.2. We first notice that the distance between these peaks in $\Delta r/r$ is roughly 4.4 times $h/r$ , which is not sensitive to either the planet mass or $\alpha$ (upper panel a). Thus, we can use the distance of these sub/super-Keplerian peaks to roughly estimate the disk temperature. On the other hand, the amplitude of the sub/super Keplerian peaks depends on all of these parameters (lower panel a). With increasing planet mass, the amplitude increases until the gap edge becomes eccentric. For the same mass planet in the same $h/r$ disk, the amplitude decreases with increasing $\alpha$ . For the same mass planet in the same $\alpha$ disk, the amplitude decreases with increasing $h/r$ .

Thus, using gas kinematics, we can first use the distance between the peaks to estimate $h/r$ , and then we can use the amplitude together with the estimated $h/r$ and assumed $\alpha$ value to derive the planet mass.

Following Kanagawa et al. (2015, 2016), we seek simple power laws to fit various observable quantities throughout the paper so that the fittings can be easily used by the community. Here, we try to find the best fit for $\Delta\delta v_{rot}$ . We define a $K_{v_{r}}$ parameter that is proportional to $q$ and has power law dependence on $h/r$ and $\alpha$ ,

[TABLE]

We try to find the best fitting parameters $ph$ and $pa$ . If $ph$ =0 or $pa$ =0, it means that the fitting does not depend on the disk $h/r$ or $\alpha$ , respectively. First, we assign values to $ph$ and $pa$ , and we can make the log $\Delta\delta v_{rot}$ - log $K_{v_{r}}$ plot for all the data points. Then, we do a linear-regression fitting for these data points using

[TABLE]

The coefficients in the fitting ( $A$ and $B$ ) are thus determined. The sum of the square difference of the vertical distance between the data points and the fitting is $\sigma$ . Finally, we vary $ph$ and $pa$ and follow the same fitting procedure until the minimum $\sigma$ is achieved. The resulting $ph$ and $pa$ are the best degeneracy parameters, and $A$ and $B$ are the best fitting parameters. For $\Delta\delta v_{rot}$ , the fitting formula is:

[TABLE]

with

[TABLE]

Thus, the sub/super Keplerian motion is most sensitive to $h/r$ , followed by $q$ and $\alpha$ . The fitting formula is shown in panel b of Figure 5 together with all measured $\Delta\delta v_{rot}$ . The uncertainty in $K_{v_{r}}$ is estimated by measuring the horizontal offset between each data point and the fitting line. From the distribution of the offset, the left side error is estimated by the 15.9 percentile of the distribution and the right side error is 84.1 percentile of the distribution. The uncertainty in $log_{10}(K_{v_{r}})$ is ${}_{-0.099}^{+0.103}$ , which is about a factor of 1.25 of $K_{v_{r}}$ .

3.2 Dust Thermal Emission

After exploring the gaseous gaps, we study the gaps in mm dust continuum maps in §3.2.1. We detail our method to fit the gap width and depth in §3.2.2.

3.2.1 Axisymmetric and Non-axisymmetric Features

As discussed in §2.3, we have 45 simulations with different $h/r$ , $\alpha$ , and $M_{p}$ . For each simulation, we generate seven continuum maps for seven $\Sigma_{g,0}$ with the DSD1 dust size distribution and five continuum maps for five $\Sigma_{g,0}$ with the DSD2 dust size distribution. Thus, we produce 45 $\times$ 12 mm maps.

The mm intensity maps for a $\Sigma_{g,0}=3\,\mathrm{g\,cm^{-2}}$ disk with DSD1 and DSD2 dust size distributions are presented in Figures 6 and 7, respectively. We want to emphasize that, if the opacity is a constant with the maximum dust size (which roughly stands when the maximum dust size, $s_{max}$ , is not significantly larger than the wavelength of observation), there is a degeneracy in the relative intensity maps between different $\Sigma_{g}$ and $s_{max}$ because only the Stokes number matters for the gas dynamics. For example, the shapes of intensity maps for the $\Sigma_{g,0}=3\,\mathrm{g\,cm^{-2}}$ and $s_{max}=0.1\,\mathrm{mm}$ cases are very similar to the $\Sigma_{g,0}=300\,\mathrm{g\,cm^{-2}}$ and $s_{max}=1\,\mathrm{cm}$ cases, since they have the same Stokes number. Thus, Figure 6 should be regarded as the dust well-coupled limit, while Figure 7 should be regarded as the dust fast-drifting limit.

Regarding the gaps and rings, there are several noticeable trends:

By comparing these two figures, we can see that the rings are more pronounced when particles with larger Stokes numbers are present in the disk. For the well-coupled case (Figure 6), the gap edge is smoothly connecting to the outer disk and the outer disk is extended. However, for the fast-drift particle cases (Figure 7), there is a clear dichotomy: either the disk does not show the gap or the gap edge becomes a narrow ring. This is because the gap edge acts as a dust trap so that a small gaseous feature can cause significant pileup for fast-drifting particles.
The marginal gap opening cases are in panels that are along the diagonal line in Figures 6 and 7, which are similar to the trend for the gaseous gaps in Figure 1.
The narrow gap edge of the the fast-drifting particle cases (Figure 7) becomes wider with a higher $\alpha$ due to turbulent diffusion. Thus, if we know the particles’ Stokes number at the gap edge, we can use the thickness of the ring to constrain the disk turbulence, as shown in Dullemond et al. (2018).

Besides axisymmetric structures, there are also several non-axisymmetric features to notice:

The gaps in the lower left panels (h5am4p5, h5am3p5) are clearly eccentric and off-centered. We may be able to use the ellipticity of the gap edges to infer the planet properties. Thus, for every mm intensity map, we find the local maximum in each azimuthal angle and use linear fitting method to measure the gap eccentricity and the distance between the center of the ellipse and the star. We find that, even in mm images generated from disks having dramatically different Stokes numbers, the gap eccentricity and off-centered distance are quite similar. However, the lower planet mass cases for the DSD1 have mild dust trapped rings thus having lower SNR, while the higher mass cases for the DSD2 have strong asymmetry, thus leading to half of the rings with the low SNR. Thus, we combine the fitting results for both DSD1 and DSD2 at $\Sigma_{g,0}$ = 3 g cm*-2*, and pick up the smaller values for eccentricity and the off-centered distance (Figure 8). We also test several cases with the ring-fitting method described in §3.1 in Huang et al. (2018a) (a MCMC fitting of the offset $\Delta x$ , $\Delta y$ , the semi-major axis, the aspect ratio and the position angle) and find that the derived eccentricity and the distance from the central star are very similar to those derived here. Clearly, both eccentricity and off-centered distance increase with the planet mass, which is consistent with gas only simulations in Kley & Dirksen (2006); Ataiee et al. (2013); Teyssandier & Ogilvie (2017); Ragusa et al. (2018). These quantities do not quite depend on $h/r$ and $\alpha$ except a weak trend that gaps in larger $\alpha$ disks have higher eccentricities. Unfortunately, due to the limited number of super-particles in the simulations, the Poisson noise in the intensity maps prevents us from measuring the eccentricity very accurately. The adopted Gaussian convolution kernel to reduce the Poisson noise has a $\sigma_{c}$ of 0.06 $r_{p}$ . If the major-axis and the minor-axis have an error of $\sigma_{c}$ /2, the uncertainty of the eccentricity is $\Delta e=\big{(}1-(1-0.03/2)^{2}\big{)}^{\frac{1}{2}}=0.17$ . Thus, any measured eccentricity smaller than 0.15 is consistent with zero eccentricity. For the same reason, any off-centered distance smaller than half of the pixel size (0.015) is consistent with zero. We mark these uncertainties as the light grey area in Figure 8. On the other hand, if the eccentricity and the off-centered distance is above these limits, our results suggest that the eccentric gap edge may be a signature of a massive planet in disks. Eccentric and off-centered gap edges have been measured in HL Tau (ALMA Partnership et al., 2015) and HD 163296 (Isella et al., 2016), which may suggest that these gaps are induced by planets.
For the lowest viscosity cases ( $\alpha=10^{-4}$ ), particle concentration within vortices can be seen at the gap edge. Even a 33 M*⊕* planet can induce particle-concentrating vortices. Interestingly, the vortex sometimes is inside the gap edge, e.g. $h/r=0.05$ , $M_{p}=1M_{J}$ case and $h/r=0.1$ , $M_{p}=3M_{J}$ case. This is probably because large particles are trapped at the gap edges, while small particles move in and get trapped into the vortex. For the majority of cases, the vortices that cause significant asymmetry in mm intensity maps are at the gap edge where $dP/dr=0$ . To characterize such large-scale asymmetries, Figure 9 shows the contrast at the gap edge, which is the ratio between the intensity of the brightest part of the ring over the intensity 180 degree opposites on the previously fitted ellipse. The figure shows that the case with a smaller gas surface density tends to show a higher contrast. We note that the contrast is very large in some cases. A 33 $M_{\earth}$ planet can lead to a factor of 100 contrast at the gap edge for a $h/r=0.05$ disk with $St=0.16$ particles. Thus, a low mass planet may also explain some of the extreme asymmetric systems: e.g. IRS 48 (van der Marel et al., 2013) and HD 142527 (Casassus et al., 2013).
The dust concentration at L5 or both L4/L5 is seen in some $\alpha=10^{-4}$ cases, consistent with previous simulations (Lyra et al., 2009). These features are more apparent than those in the gas (Figure 1). As pointed out by Ricci et al. (2018), such features may be observable. On the other hand, we want to emphasize that the dust concentration at Lagrangian points is not in a steady state, and the amount of dust at those points decreases with time. Thus, in this paper, we will not use these feature to constrain the planet properties.

3.2.2 Fitting Gaps/Rings

To derive the relationship between the gap profiles and the planet mass, we azimuthally average the mm intensity maps as shown in Figure 10. The solid curves are for models with $s_{max}=0.1$ mm (DSD1), while the dashed curves are for models with $s_{max}=1$ cm (DSD2).

We try to find the relationship between the planet mass and the gap properties (such as the gap width $\Delta$ and depth $\delta$ ), using the dust intensity profiles in Figure 10. Previous works such as Kanagawa et al. 2016, 2015; Dong & Fung 2017 studied the relationship between the planet mass and the gaseous gap width and depth. However, mm observations are probing dust with sizes up to mm/cm and these dust can drift in the gaseous disk. Thus, studying only the gaseous gap profiles is not sufficient for explaining mm observations and carrying out a similar study but directly for dust continuum maps is needed. We seek to first find a relationship between disk and planet properties ( $\alpha$ , $h/r$ , and $M_{p}$ ) using the fitting of the azimuthally averaged gas surface density profile, and characterize those three parameters using a single parameter $K$ (for the depth- $K$ relation) or $K^{\prime}$ (for the width- $K$ ’ relation). Then, we fit the azimutally averaged dust intensity profile for our grid of models and find their depth- $K$ and width- $K^{\prime}$ relations. Overall, our fitting follows Kanagawa et al. (2016) and Kanagawa et al. (2015) but extend those relationships to dust particles with different sizes.

The detailed steps are the following:

(1) We measure the gap depth ( $\delta$ ) for both gas surface density profiles (Figure 2) and mm intensity profiles (Figure 10). From the outer disk to the inner disk, we first find the outer peak (the first local maximum, which corresponds to where dust piles up due to the dust trapping) and mark this point as $r_{peak}$ , and then find the bottom of the gap (local minimum) inside $r_{peak}$ and mark it as $r_{gap}$ . $r_{gap}$ is not necessarily $r_{p}$ . As demonstrated in Figure 11, the gap can have the deepest point further out than $r_{p}$ . This is because some gaps have significant horseshoe material in between. In some extreme cases with very shallow gaps, only the outer portion of the gap that is outside the horseshoe region is visible (e.g. the top middle panel in Figure 10). We define the gap depth $\delta$ as

[TABLE]

for the gas surface density profiles, and

[TABLE]

for the dust mm intensity profiles.

(2) Measuring the gap width ( $\Delta$ ) for these profiles. To calculate the width, we first define the edge quantities as the average between the peak and gap surface densities (for gas) or the mm intensities (for intensity maps):

[TABLE]

and

[TABLE]

Then, we find one edge $r_{in}$ at the inner disk and the other $r_{out}$ at the outer disk, where $\Sigma$ ( $r_{in}$ ) = $\Sigma$ ( $r_{out}$ ) = $\Sigma_{edge}$ for the gas surface density or $I(r_{in})$ = $I(r_{out})$ = $I_{edge}$ for the dust intensity (Figure 11). Thus, we define the gap width $\Delta$ for either the gas surface density or the dust intensity as

[TABLE]

Figure 12 shows $\Delta$ for all $\Sigma_{g,0}$ cases with DSD1 (panel a) and DSD2 (panel b) dust distributions. If there is some horseshoe material around $r$ = $r_{p}$ separating the main gap into two gaps, the horizontal $\Sigma_{edge}$ or $I_{edge}$ line will cross through the horseshoe and we treat two individual gaps as a single one (i.e., the $r_{in}$ is taken to be the $r_{in}$ of the inner gap and $r_{out}$ is taken to be the $r_{out}$ of the outer gap), but the individual gaps on either side of the horseshoe region are also plotted in Figure 12 as fainter makers and they are connected to the main gap width using dotted lines.

Note that our definition of gap width is more convenient to use than that in Kanagawa et al. (2016), because the width here is normalized by $r_{out}$ instead of $r_{p}$ as in Kanagawa et al. (2016). In actual observations, we do not have the knowledge of the planet position $r_{p}$ within the gap. Another difference between our defined gap width and the one used in Kanagawa et al. (2016) is that we use $(\Sigma(r_{peak})+\Sigma(r_{gap}))/2$ to define $\Sigma_{edge}$ while Kanagawa et al. (2016) use $\Sigma_{0}/2$ to define the gap edge. Our definition enables us to study shallow gaps that are shallower than $\Sigma_{0}/2$ .

(3) Fitting the width ( $\Delta$ )- $K^{\prime}$ relation. We first use the width $\Delta$ measured from the gas surface density profiles to find the optimal degeneracy parameter $K^{\prime}$ following the same procedure as in Equation 14. Similarly, a least squares fitting was done to minimize the sum of the square difference of the vertical distance between the points and the linear-regression line log( $\Delta$ ) vs. log(K’). With this procedure, we derive that the optimal $K^{\prime}$ is

[TABLE]

With this definition of $K^{\prime}$ , the best fitting relationships ( $\Delta=AK^{\prime B}$ ) are found for each initial gas densities with two dust size distributions DSD1 and DSD2. The resulting $A$ and $B$ for these fits are listed in Table 1. Note that our definition of $K^{\prime}$ is equivalent to the square root of $K^{\prime}$ defined in Kanagawa et al. (2016). Compared with the fitting formula for the gas surface density in Kanagawa et al. (2016), our $K^{\prime}$ is less sensitive to $h/r$ and the gaseous gap width is less sensitive to $q$ . We confirm that this is largely due to our different definition of the gap width (compared with their definition, our normalized gap width is smaller for wide gaps and larger for shallow gaps that are normally narrow.).

Figure 12 shows the fits for all the cases with DSD1 (panel a) and DSD2 (panel b) dust size distributions. We can see that uncertainties of these fittings become large when $\Delta\lesssim 0.15$ . Thus, our fitting procedure does not involve widths that are smaller than 0.15. For these narrow gaps whose widths are smaller than 0.15 (labeled as the open symbols with back numbers in them), their gap profiles start to be affected by the smoothing kernel with $\sigma=0.06r_{p}$ . Thus in Figure 12, we also plot the widths measured from the profiles that are convolved with a $\sigma$ = 0.025 $r_{p}$ kernel. These widths are plotted as open symbols with red numbers in them.

(4) Fitting the depth ( $\delta$ )- $K$ relation. We adopt the same procedure to fit the depth- $K$ as the width- $K^{\prime}$ aforementioned. Since no-gap is equivalent to $\delta$ =1, we try to find the optimal degeneracy parameter $K$ by a least squares fitting for log( $\delta$ - 1) vs. log( $K$ ),

[TABLE]

for various $K$ . The optimal $K$ is fitted to be

[TABLE]

After $K$ is fixed, we use Equation 23 to fit the relationship between $\delta$ - 1 and $K$ for the dust intensity profiles from different $\Sigma_{g,0}$ with DSD1 and DSD2. $C$ and $D$ are found using linear regression. The resulting $C$ and $D$ in different $\Sigma_{g,0}$ cases with either DSD1 or DSD2 are listed in Table 2. Figure 13 show $\delta$ - 1 for all $\Sigma_{g,0}$ cases with DSD1 and DSD2. The best fits are also plotted for each panel. Note that open symbols are not involved in the fitting since these gaps are eccentric and their depths do not follow the trend for other gaps. Clearly, with the Stokes number increasing, the fitting becomes worse. This is expected since particles with larger Stokes numbers drift faster and the gap profile becomes more irregular.

(5) The uncertainty of the fittings. We apply the same measure to calculate the uncertainty of the gap width/depth fitting as that of $\Delta\delta v_{rot}$ - $K_{v_{r}}$ relation mentioned in §3.1.2. That is, we measure the horizontal offset (in $log_{10}(K^{\prime})$ or $log_{10}(K)$ ) between each point and the fitting line at each sets of dust configurations and also the gas surface density. From the distribution of the offset, the left side error is estimated by the 15.9 percentile of the distribution and the right side error is 84.1 percentile of the distribution. These uncertainties are summarized in Table 1 and 2 and marked in grey color at the top of each panel in Figure 12 and 13. For widths that are larger than 0.15, the uncertainties for the fittings are less than a factor of two for $K^{\prime}$ (or $q$ ) when $St\lesssim 5\times 10^{-3}$ and around a factor of three for $K^{\prime}$ (or $q$ ) when $5\times 10^{-3}<St\lesssim 5\times 10^{-2}$ . When $St\gtrsim 10^{-1}$ , particles drift to the central star quickly and most of the gaps only have a single ring left at the outer disk so that $\Delta\sim$ 1 and the uncertainties for $K^{\prime}$ at a given $\Delta$ is very large. For these cases, we cannot use the gap width to estimate the planet mass.

Finally, we summarize all the fits for the width and depth in Figure 14. In the Appendix, we provide gap depth $\delta$ and width $\Delta$ of our whole grid of models. In spite of the dramatically different dust size distributions between DSD1 and DSD2, the fits for DSD1 are quite close to fits for DSD2 as long as the Stokes number for the maximum-size particles is the same (e.g. red solid and dot-dashed lines). This is reasonable since only the Stokes number matters for the dust dynamics, and DSD1 have a similar opacity as DSD2. For 1 mm observations, the opacity is roughly a constant when $s_{max}\lesssim$ 1 cm (the opacity is slightly higher when $s_{max}\sim$ 1mm, see Birnstiel et al. 2018). Thus, different disks with different surface densities ( $\Sigma_{g,0}$ ) and different dust size distributions have the same intensity profiles as long as their Stokes numbers for maximum-size particles (where most of the dust mass is) are the same and $s_{max}\lesssim$ 1 cm. Thus, our derived relationships can be used in other disks with different surface densities and dust size distributions as long as the Stokes number of the maximum-size particles is in our simulated range (1.57 $\times$ 10*-4* to 1.57). For disks with Stokes number smaller than 1.57 $\times$ 10*-4*, their gap profiles should be similar to the disks with $St$ =1.57 $\times$ 10*-4* since dust is well coupled to the gas.

3.2.3 Secondary Gaps/Rings

Previous simulations have shown that a planet can introduce many gaps/rings in disks having very low viscosities (Zhu et al., 2014; Dong et al., 2017; Bae et al., 2017). These gaps can be grouped into two categories: 1) two gaps adjacent to the planet that are separated by the horseshoe material (e.g. two troughs at 0.9 $r_{p}$ and 1.1 $r_{p}$ in Figure 11, also mentioned in §3.1.1), and 2) secondary shallower gaps much further away into the inner and outer disks (e.g. the gap at 0.6 $r_{p}$ in Figure 11). The two gaps in the first category form because: a) the spiral waves, especially excited by low mass planets, need to propagate in the radial direction for some distance to steepen into spiral shocks and induce gaps (Goodman & Rafikov, 2001), b) the horseshoe material has a slow relative motion with respect to the spiral shocks thus this material takes a long time to be depleted. Eventually, these two gaps may merge into one single main gap, which is studied in §3.2.2. The gaps in the second category are induced by additional spiral arms from wave interference (Bae & Zhu, 2018a). Instead of disappearing, these gaps will become deeper with time in inviscid disks. Thus, they are useful to constrain the planet and disk properties (Bae & Zhu, 2018b).

We label the positions of all these additional gaps and rings in Figure 15. We find that the positions of these rings and gaps in dust intensity radial profiles are similar to those in gas surface density profiles. Thus, we plot the positions based on the gas density profiles. It turns out that only disks with $\alpha\leq 10^{-4}$ can form noticeable multiple gaps. Thus, if we find a system with multiple gaps induced by a single planet (e.g. AS 209 in the next section), the disk viscosity has to be small. From Figure 15, we can see that distance between the secondary gap and the main gap mainly depends on the disk scale height ( $h$ ).

For the secondary gap at $\sim 0.5-0.7$ , following our fitting procedure before, we find that the position of the secondary gap ( $r_{IG2}$ ) and $r_{p}$ is best fitted with

[TABLE]

This clearly shows that the position of the secondary gap is almost solely determined by the disk scale height. Thus, if the secondary gap is present, we can use its position to estimate the disk scale height ( $h/r$ ). The fitting is given in Figure 16. The $\alpha=10^{-5}$ cases are the AS 209 cases which will be discussed in the next section. We caution that the fitting has some scatter. Within each $h/r$ group in Figure 16, the $r_{IG2}/r_{p}$ depends on the planet mass. But this dependence seems to be different for different $h/r$ groups, so that the fitting using all $h/r$ suggests a weak dependence on the planet mass. We also note that our fit is different from the recent fit by Dong et al. (2018b) which has a $q^{-0.2}(h/r)^{1.3}$ dependence (note that their planet mass is normalized by the thermal mass). The difference may be due to: 1) The disks in Dong et al. (2018b) are thinner, where their main set of simulations uses $h/r$ =0.03, 2) Dong et al. (2018b) fit the gap positions at different times for different simulations while we fit the gap positions at the same time in the simulations.

4 Planet Properties

With all the relationships derived in previous sections regarding the planet mass and gap profiles, we can now put them together to constrain the mass of potential planets in the DSHARP disks. We use the measured radial intensity profiles from Figure 2 in Huang et al. (2018a). These profiles are derived by deprojecting the observed images to the face-on view and then averaging the intensity in the azimuthal direction. Details regarding generating the radial intensity profiles are given in Huang et al. (2018a). By using these intensity profiles, we can derive the planet mass following the flowchart given in Figure 17.

First, for each source, we plot the observed radial intensity profile and identify gaps that have $\Delta$ $\geq$ 0.15. As shown in Figure 12, $\Delta$ $\lesssim$ 0.15 have large scatter and are sensitive to the size of the convolution beam. By examining the surface density profiles in detail, we find that such narrow gaps are also very shallow and they are actually the outer one of the double gaps around the horseshoe region. Since these gaps are very shallow, the inner one does not cause enough disk surface density change to be identified as a gap. Thus, for narrow gaps with $\Delta$ $\lesssim$ 0.15, we do not use the fitting formula to derive the planet mass. Instead, we try to directly match the gap $\Delta$ with data points in Figure 12 by eye to get a rough planet mass estimate. For these narrow gaps, the size of the convolution beam matters. Thus, if the gap is at 10s of au, we use the widths derived in images with the $\sigma=0.06r_{p}$ beam, and if the gap is at $\sim$ 100 au we use the widths derived in images with the $\sigma=0.025r_{p}$ beam.

Second, we estimate the gas surface density, using the observed mm flux at the outer disk and/or some other constraints. We integrate the observed intensity from 1.1 $r_{gap}$ to 2 $r_{gap}$ where $r_{gap}$ is the gap center. Using $T_{d}$ derived by Equation 12 and the dust opacity of 0.43 $\mathrm{cm^{2}\,g^{-1}}$ (§2.3), we calculate the averaged dust surface density ( $\Sigma_{d}$ ) from 1.1 to 2 $r_{gap}$ . We have done the same exercise for all our simulations, and Figure 18 shows the relationship between $\Sigma_{g,0}$ and the averaged $\Sigma_{d}$ at the outer disk for the simulations. Figure 18 indicates that, with a smaller gas surface density or larger particles (higher Stokes numbers), the ratio between $\Sigma_{d}$ and $\Sigma_{g,0}$ increases because particles with larger Stokes numbers are more easily trapped at the gap edges. We can then use Figure 18 to estimate $\Sigma_{g,0}$ based on the derived $\Sigma_{d}$ from the observation, the estimated $h/r$ , and the assumed $\alpha$ and planet mass. After we derive the planet mass, we will go back to this step to see if the derived planet mass is consistent with our assumed mass. Otherwise, we iterate these processes again with the new assumed planet mass. On the other hand, this estimate is prone to large errors. If we have more ways to estimate the gas surface density, such as using molecular tracers or constraints from the gravitational instability, we should adopt these constraints.

Third, with known $\Sigma_{g,0}$ and the assumed dust size distribution, we can calculate $St_{max}$ and use the $\Delta$ - $K^{\prime}$ relationship (§3.2.2 and Table 1) to derive the $K^{\prime}$ parameter. Given the sensitivity limits of ALMA, we decide not to use the gap depth ( $\delta$ ) to estimate the $K$ parameter. For example, two gaps with different depths, one being a factor of 105 deep and the other being a factor of 103 deep, can look similar if the S/N of the observation is 100.

Next, we need to constrain the disk scale height and the disk $\alpha$ parameter to break the degeneracy of $K^{\prime}$ in order to derive $q$ . For each major gap, if there is a shallower gap at $r/r_{p}\sim$ 0.5-0.7, the shallower gap may be the secondary gap induced by the planet. The distance between the secondary gap and $r_{p}$ is very sensitive to $h$ (§3.2.3 and Equation 25). Thus, the presence of the secondary gap at the right radii not only makes the planet gap-opening scenario more plausible but also gives constraints on the disk scale height. If there is no secondary gap, we may need to use radiative transfer calculations or Equation 12 to estimate the disk temperature. The existence of the secondary gap also implies that the disk viscosity parameter $\alpha\lesssim 10^{-4}$ . Without the presence of the secondary gap, the $\alpha$ parameter can then be constrained by the symmetry of the disk structures. If the rings/gaps are highly axisymmetric, $\alpha$ is likely to be larger than $10^{-4}$ .

Finally, we can use Equation 22 to calculate $q$ and thus the planet mass. With $M_{p}$ derived, we can go back to Step 2 to estimate a more accurate gas surface density. We can also do a consistency check with the derived $M_{p}$ . For example, we can check if the sub/super-Keplerian motion at the gap edge could be detected (§3.1.2, Equation 16), if the planet should produce large-scale asymmetries (e.g. eccentricity, vortices §3.2.1, Figure 8), and if the gap depth is consistent with observations (Table 2).

Following this procedure (Figure 17), we identify potential planets in the DSHARP disks (as summarized in Table 3) using the intensity profiles from Huang et al. (2018a). All the gaps with $\Delta$ $\geq$ 0.15 in the DSHARP sample have been carefully measured for their widths and then we use the fitting formula to estimate the planet mass based on their widths. These are shown in the upper part of Table 3. Since each fitting line with a Stokes number comes with an uncertainty in $K^{\prime}$ (See §3.2.2, and Table 1), the uncertainties of the planet mass with the given $\alpha$ and $h/r$ are also included in the table. For shallow gaps with $\Delta$ $\leq$ 0.15, our fitting formulae fail to fit the gap widths from the simulations and the gap width is also sensitive to the convolution beam size (Figure 12). Thus, we only choose those that look similar to shallow gaps in our grid of numerical simulations and compare them directly with simulations. Thus, only a subset of the shallow gaps in DSHARP sample have been fitted. They are shown in the lower part of Table 3. Since we compare these shallow gaps with the simulations by eye, proper error estimate can not be provided. Thus, they are considered not robust and complete, and will not be included in the statistical study later. This also means that our statistical study may miss low mass planets. In the next section, we will comment on each case in detail.

Table 3 gives the gap positions, measured gap widths, outer disk dust surface densities and estimated $h/r$ . Using the dust-to-gas mass ratio (Figure 18) in simulations with different dust size distributions (DSD1 and DSD2), the gas surface densities are also provided. If the gas surface density is above the gravitational instability (GI) limit with $Q=1$ , we use the GI limit as the gas surface density. Then with $St_{max}$ calculated for DSD1 and DSD2, we derive $K^{\prime}$ for DSD1 and DSD2 using $\Delta-K^{\prime}$ relationships. To break the degeneracy in $K^{\prime}$ to derive $q$ , we need to know the disk viscosity. Thus, for either DSD1 or DSD2, we provide three possible planet masses with the disk $\alpha$ =10*-2*, 10*-3*, and 10*-4*. These three masses are labeled as $M_{p,am2}$ , $M_{p,am3}$ , and $M_{p,am4}$ , which are listed in Table 3. The inferred planet mass is roughly twice as high if $\alpha$ is 10 times larger. This is because $K^{\prime}=q(h/r)^{-0.18}\alpha^{-0.31}$ , so that $q\propto\alpha^{0.31}$ with a given $K^{\prime}$ and $h/r$ . As shown in Table 3, many gaps (especially having low $\Sigma_{g,0}$ ) cannot be fit using DSD2 dust size distribution. This is because the Stokes number for dust in DSD2 is very large, so that particles in the inner disk quickly drift to the central star forming a cavity with a single ring at the gap edge. This is consistent with the conclusion in Dullemond et al. (2018) that large particles (cm-sized) are not preferred in the DSHARP disks.

As can be seen from Equation 2 and Table 3, the Stokes number estimated from DSD1 and DSD2 can differ by three orders of magnitude. DSD1 with $s_{max}$ = 0.1 mm and DSD2 with $s_{max}$ = 1 cm can be seen as two extreme cases. Dust with $s_{max}<$ 0.1 mm should have similar profiles as DSD1 since 0.1 mm particles already couple with the gas well in the sample. Dust with $s_{max}$ = 1 cm already drifts very fast and we can hardly find a mass solution for most of our disks. To cover a more comprehensive parameter space, we add a new set of planet masses estimated assuming $s_{max}$ = 1 mm (”1 mm” hereafter). The estimated initial gas density $\Sigma_{g,0}$ are used between the values of DSD1 and DSD2. Holding $\Sigma_{g,0}$ constant, $St_{max}$ for ”1 mm” is 10 times larger than that of the DSD1 or 10 times smaller for DSD2. Thus, the Stokes number of the ”1 mm” models are in between those two extremes. The gap width-K’ relation of the ”1 mm” models are taken from the corresponding $St_{max}$ fits in DSD1. The justification is that only the Stokes number matters regarding the gap width, as discussed at the end of §3.2.2 and demonstrated in Figure 14. The estimated $\Sigma_{g,0}$ , $St_{max}$ , three planet masses given $\alpha=10^{-4},10^{-3},10^{-2}$ and their uncertainties are all given in Table 3 in the order of DSD1, ”1 mm” and DSD2 (ascending $s_{max}$ ). Among the nine planet masses estimated for each source, we prefer $M_{p,am3}$ with DSD1 size distribution. The main reason that $\alpha=10^{-3}$ is preferred is that most rings of the DSHARP sample do not show significant asymmetry, indicating that $\alpha\gtrsim 10^{-3}$ . On the other hand, if the gaps are shallow, low mass planets in $\alpha=10^{-4}$ disks can also produce axisymmetric gaps/rings.

4.1 Comments on Individual Sources

4.1.1 AS 209

AS 209 is a system with many gaps. Fedele et al. (2018) found two gaps at 62 au and 103 au and they proposed that a 0.7 $M_{Saturn}$ planet at $\sim$ 103 au can explain both gaps. Huang et al. (2018a) and Guzmán et al. (2018) identified many gaps in this system including dark annuli at 9, 24, 35, 61, 90, 105 and 137 au. Following our procedure (Figure 17), we first derive the $K^{\prime}$ parameter for the main gap at $\sim$ 100 au. The narrow width of the gap suggests that it is a sub-Jupiter mass planet. Then we find that the gap at r = 61 au is shallower than the main gap, and it is at 0.5-0.7 $r_{p}$ . Thus, we treat it as a secondary gap induced by the planet. The distance between the secondary and primary gaps suggests that $h/r\sim 0.05-0.06$ (Equation 25 and Figure 16). This $h/r$ is slightly smaller than the simple estimate with Equation 12, but the faint emission at the near-IR scattered light image (Avenhaus et al., 2018) may support that the disk is indeed thin (another possibility is that the disk is significantly less flared.). With this $h/r$ and $K^{\prime}$ , we derive that the 100 au planet has a mass of $q=3\times 10^{-4}$ in a $\alpha=10^{-4}$ disk or $q=10^{-4}$ in a $\alpha=10^{-5}$ disk. Motivated by the smaller gaps at 24 and 35 au from the DSHARP data (Guzmán et al., 2018), we carry out several additional simulations extending the range of $\alpha$ to $10^{-5}$ . Since a smaller $\alpha$ is used, we double the numerical resolution for all simulations that are constructed for AS 209. Surprisingly, the $q=10^{-4}$ planet in a $\alpha=10^{-5}$ and $h/r=0.05$ disk can explain all 5 gaps at 24, 35, 62, 90 and 105 au (Figure 19). Although we assume that there is another planet at 9 au to explain the 9 au gap, it is possible that the 9 au gap is also produced by the main planet at 99 au, considering that our simulation domain does not extend to 9 au. We want to emphasize that our simulation with one planet at 99 au not only matches the primary gap around 100 au, but also matches the position and amplitude of secondary (61 au), tertiary (35 au) and even the fourth (24 au) inner gaps. This makes AS 209 the most plausible case that there is indeed a planet within the 100 au gap.

Although the above model reproduces the positions and intensities of gaps and rings very well, its synthetic image (the upper middle panel in Figure 19) shows a noticeable horseshoe region and some degree of asymmetry in the rings. Such asymmetry disappears when $\alpha\gtrsim 10^{-3}$ . On the other hand, the presence of the tertiary and the forth inner gaps requires a small $\alpha$ . Thus, we carry out a simulation with a radially varying $\alpha$ ( $\alpha=3\times 10^{-4}(r/r_{p})^{2}$ ). This model reproduces the 2-D intensity maps better, as shown in the right panels of Figure 19 and also presented in Guzmán et al. (2018). Such a radially varying $\alpha$ disk has also been suggested to explain HD 163296 (Liu et al., 2018). If these models are correct, they suggest that $\alpha$ in protoplanetary disks is not a constant throughout, supporting the idea that different accretion mechanisms are operating at different disk regions (Turner et al., 2014).

Dullemond et al. (2018) constrained that the $\alpha/St$ for the ring at 74 au has a range roughly between 0.03 and 0.7 from the limits of pressure bump width argument (See Table 3 therein). Such constraint is derived using the particle trapping model and does not depend on the origin of the ring. In our $\alpha=10^{-5}$ model, $\alpha/St_{max}\approx 0.003$ and in our $\alpha$ varying model, $\alpha/St_{max}\approx 0.02$ . The actual characteristic $St$ can be smaller, considering that the $St_{max}$ here is the maximum Stokes number at the position of the planet in the initial condition ( $t_{0}$ ). Since for both models $n(s)\propto s^{-3.5}$ , 50% of the dust mass in $t_{0}$ at $r_{p}$ have $St\leq 0.25\ St_{max}$ . Adopting these values, their $\alpha/St\approx$ 0.012 and 0.08, respectively. Thus, the $\alpha=10^{-5}$ model is off the lower limit of $\alpha/St$ by a factor of 3, whereas the $\alpha$ varying model is safely above the lower limit. Considering that the turbulent diffusion with the small $\alpha$ ( $\alpha=10^{-5}$ ) in our simulations may have not reached to a steady state, we conclude that these models are consistent with Dullemond et al. (2018).

4.1.2 Elias 24

Elias 24 (Cieza et al., 2017) is another system that looks very similar to our planet-disk interaction simulations. It has a deep gap at 57 au, a narrow ring at 77 au, and an extended outer disk (Huang et al., 2018a). The narrowness of the ring is suggestive of particle trapping at the gap edge. Dipierro et al. (2018) estimated that there is a 0.7 $M_{J}$ mass planet at 57 au, while Cieza et al. (2017) suggested that the mass of the 57 au planet is 1-8 $M_{J}$ . Our estimate is roughly consistent with these previous estimates. The planet mass is $\sim$ 0.8 $M_{J}$ with $\alpha=10^{-3}$ and DSD1. On the other hand, the clear signature of dust pile-up at the outer gap edge may indicate that dust is larger than 0.1 mm as used in DSD1. If dust particles in Elias 24 are larger than 0.1 mm, the planet mass can be lower than our estimates. Based on our grid of simulations, we run an additional simulation with $\alpha=5\times 10^{-4}$ , $h/r=0.07$ and $M_{p}=0.16\ M_{J}$ ( $q=0.2\ M_{J}/M_{*}$ ). We put the single planet at the 57 au gap and the result is shown in Figure 20. The dust distribution is $n(s)\propto s^{-3.5}$ , $s_{max}$ = 2 mm, and initial gas surface density $\Sigma_{g,0}=15$ $\mathrm{g\,cm^{-2}}$ , hence $St_{max}=2.09\times 10^{-2}$ . Dullemond et al. (2018) estimated that the $\alpha/St$ is between 0.077 to 0.66 at the 77 au bright ring. Our estimated $\alpha/St_{max}=2.39\times 10^{-2}$ is roughly consistent with their lower limit considering that 50% of the dust mass has $\alpha/St>0.096$ under the dust size distribution $p=-3.5$ .

4.1.3 Elias 27

The spiral arms detected in Elias 27 (Pérez et al., 2016) suggest that the disk may be undergoing gravitational instability or there is a massive companion at the outer disk (Meru et al., 2017). Besides the spirals, there is a shallow annular gap at 70 au (Huang et al., 2018b). If we follow our procedure to fit this gap, the planet mass is 0.06 $M_{J}$ using $\alpha=10^{-3}$ and DSD1. Such a low mass planet can not induce the large-scale spirals as observed (Zhu et al., 2015a). On the other hand, detecting this shallow gap means that if there are massive companions in the system within 200 au (e.g. with masses larger than 0.06 $M_{J}$ ), we should be able to see the induced gaps at the mm continuum images. The lack of deep gaps suggests that there are no massive companions in this disk within 200 au. The spirals must be induced by a massive companion outside 200 au or by some other mechanisms (e.g. GI).

4.1.4 GW Lup

GW Lup has two narrow gaps at 74 and 103 au. The former gap is barely above $\Delta=$ 0.15 and the latter is extremely narrow with $\Delta\lesssim$ 0.15. We decide to only fit the 74 au gap since the 103 au gap is too shallow to fit with any of our models. To produce the 74 au gap, the planet mass must be very small ( $\sim$ 0.03 $M_{J}$ or 10 $M_{\earth}$ ). If both 74 and 103 au gaps are part of a wide gap separated by the horseshoe region, the planet will be at $\sim$ 85 au with $M_{p,am3}$ = 0.36 $M_{J}$ or $M_{p,am4}$ = 0.18 $M_{J}$ . The $K$ parameter (Equation 14) is thus $\sim$ 11 and the gaseous gap depth $\delta$ is $\sim$ 2, which is roughly consistent with the observations (Huang et al., 2018a). Thus, this more massive planet solution remains a possibility.

4.1.5 HD 142666

HD 142666 has several shallow dark annuli at 16, 36, and 55 au (Huang et al., 2018a). The outer two dark annuli (36 and 55 au) as identified in Huang et al. (2018a) have widths of 0.05 and 0.04 by our definition, less than the minimum width measured in our models. Thus, we do not fit those two gaps either. We only fit the 16 au gap, and it suggests that $M_{p,am3}$ is 0.3 $M_{J}$ with DSD1 and 0.2 $M_{J}$ with DSD2.

4.1.6 HD 143006

HD 143006 has two wide gaps at r = 22 au and r = 51 au (Pérez et al., 2018). The gap at r = 22 au has the widest relative width ( $\Delta$ ) in all DSHARP disks, which also leads to the highest inferred planet mass with $M_{p,am4}$ = 10 $M_{J}$ and $M_{p,am3}$ = 20 $M_{J}$ . Both submm continuum observations (Pérez et al., 2018) and the near-IR scattered light observations (Benisty et al., 2018) have suggested that the inner disk inside 10 au is misaligned with the outer disk. If such misalignment is caused by a planet on an inclined orbit, the planet mass needs to be larger than 2 $M_{J}$ in an $\alpha=10^{-3}$ disk (Zhu 2018), which is consistent with the high planet mass derived from fitting the gap profile here. With such a massive planet predicted, HD 143006 is a prime target to look for exoplanets with direct imaging techniques.

The outer gap at 51 au can be explained by a sub-Jovian planet in the disk. The 51 au gap also has an interesting arc feature at the outer edge, which implies that the disk viscosity may be low ( $\alpha\lesssim 10^{-4}$ ) and $M_{p,am4}$ are preferred in this system.

Note that such high inferred planet-stellar mass ratio at 22 au exceeds the largest $q$ (3 $M_{J}/M_{*}$ ) in our grid of simulations. This brings more uncertainties to the estimated planet mass. Nevertheless, we believe that our extrapolation of Equation 22 to $q=0.01$ is justifiable since the dust is well coupled to the gas due to the small Stokes number under DSD1, and the previous study with a grid of much higher $q$ (Fung et al., 2014) showed that the relation between gaseous gap properties and the planet mass can extend to $q$ =0.01.

4.1.7 HD 163296

HD 163296 is another system with multiple gaps. The DSHARP observations (Huang et al. 2018a; Isella et al. 2018) reveal 4 gaps at 10 au, 48 au, 86 au and 145 au. Based on the gap widths, we estimate that the planets at 10 au, 48 au, and 86 au have masses of 0.71, 2.18, 0.14 $M_{J}$ in an $\alpha=10^{-3}$ disk with DSD1 dust. If the disk $\alpha=10^{-4}$ , the planet masses are 0.35, 1.07, 0.07 $M_{J}$ with DSD1 dust. Except the 10 au gap, the rest gaps have been revealed by previous ALMA observations (Isella et al., 2016). Isella et al. (2016) estimated that the 48 au planet has a mass between 0.5 and 2 $M_{J}$ and the 86 au planet has a mass between 0.05 and 0.3 $M_{J}$ , which are roughly consistent with our estimate. Our derived gas surface density ( $\Sigma_{g,0}$ ) of 3-30 g cm*-2* at 48 au and 86 au is also consistent with $\sim$ 10 g cm*-2* derived in Isella et al. (2016). Teague et al. (2018) studied the deviation from the Keplerian velocity profile as measured from CO line emission and inferred that the planet at 86 au has a mass around $M_{J}$ , which is larger than our derived $M_{p,am2}$ by a factor of 3. However, the planet mass assuming $\alpha=10^{-2}$ and 1 mm sized particles including $1\sigma$ error can reach to $\sim$ 0.6 $M_{J}$ . Considering that the uncertainty is a factor of two in Teague et al. and also the uncertainties in our adopted gas density, dust size distribution and disk viscosity, these results are still consistent. Liu et al. (2018) has adopted a disk with an increasing $\alpha$ from 10*-4* at 48 au to $10^{-2}$ at 86 au, and estimated that planets at 48 au and 86 au have masses of 0.46 and 0.46 $M_{J}$ (their same values were purely a coincidence). This is consistent with our estimate if we adopt the same $\alpha$ values.

An asymmetric structure is discovered at the outer edge of the 48 au gap (Isella et al., 2018), implying that the disk viscosity $\alpha\lesssim 10^{-4}$ . Thus, the Mp,am4 may be more representative for the 48 au gap.

4.1.8 SR 4

SR 4 has a wide single gap at 11 au. We estimate its mass $M_{p,am3}$ = 2.16 $M_{J}$ with DSD1 and 0.77 $M_{J}$ with DSD2. The gap is also quite deep, consistent with the presence of a Jovian mass planet. Thus, SR 4 may be an interesting source to follow up to study its gas kinematics or detect the potential planet with direct imaging observations.

4.1.9 DoAr 25, Elias 20, IM Lup, RU Lup, Sz 114 and Sz 129

These six systems have shallow gaps with $\Delta<0.15$ . Thus, we compare the observed gap widths directly with those derived in numerical simulations (Figure 12). The inferred planet mass is less than 0.1 $M_{J}$ for all these gaps. The smallest planet is 0.02 $M_{J}$ or 6.4 $M_{\earth}$ . Note also that IM Lup features intricate spiral arms inside the gap fit at 117 au (Huang et al., 2018b).

On the other hand, DoAr 25, Elias 20, and RU Lup have adjacent double gaps, similar to GW Lup. If we treat these double gaps as one main gap which is separated by the horseshoe material, we can derive the planet mass under this scenario. To explain both the 98 and 125 au gaps in DoAr 25 using a single planet, the planet is at 111 au with $M_{p,am3}=0.73\,M_{J}$ or $M_{p,am4}=0.36\,M_{J}$ . To explain the 25 and 33 au gaps in Elias 20, the planet is at 29 au with $M_{p,am3}=0.57\,M_{J}$ or $M_{p,am4}=0.28\,M_{J}$ . To explain the 21 and 29 au gaps in RU Lup, the planet is at 24 au with $M_{p,am3}=1.18\,M_{J}$ or $M_{p,am4}=0.58\,M_{J}$ . To make the gaps as shallow as possible, we assume DSD1 dust distribution here. Even so, the corresponding gap depth $\delta$ is larger than 2 with these planet masses. By comparing with the intensity profiles in Huang et al. (2018a), DoAr 25 has gaps that could be deep enough, while the gaps in both Elias 20 and RU Lup are too shallow and this scenario seems unlikely.

4.2 Young Planet Population

Now, we can put these potential young planets in the exoplanet mass-semimajor axis diagram (Figure 21). Considering most of these systems do not show asymmetric structures, we pick the planet mass that is derived using $\alpha=10^{-3}$ and DSD1. The mass errorbar is chosen as the minimum and maximum planet mass among all the nine masses that have constrained values in Table 3 (columns 11 to 13), adding up the additional uncertainty due to the fitting from the column 14 of the table. Thus, this is a comprehensive estimate of the error covering different disk $\alpha$ (from 10*-4* to $10^{-2}$ ), particle sizes ( $s_{max}$ from 0.1 mm to 1 cm), and the errors of the fitting. The planet masses that are from very narrow gaps in the lower part of Table 3 (the ones with brackets) are labeled with light circles, and we do not count them in the statistical study below since the narrowness of the gaps leads to large uncertainties in the mass estimate. Bae et al. (2018) has collected young planets from previous disk observations in the literature (most are Herbig Ae/Be stars). Here, we only consider the DSHARP sample (Andrews et al., 2018). Although this sample is more homogeneous with similar observation requirements, it is still slightly biased towards bright disks and thus high accretion rate disks around more massive stars.

Since the DSHARP observations have resolutions of $\sim$ 3-5 au and most disks only extend to 200 au in the dust continuum images, the planet population we can probe lies between 5 and 200 au. The probed mass limit is around the Neptune mass in the outer disk and a little bit higher (a factor of $\sim$ 2) in the inner disk ( $<$ 10 au, with a larger beam size). If there are planet-induced gaps in the disk, we should always detect them at almost all the viewing angles unless the disk is very edge on. Thus, the probability that we are missing gap-induced planets due to the observational bias is small. Under this circumstance, we can simply estimate the planet occurrence rate through dividing the number of planets by the total number of disks observed. Although DSHARP observes 20 disks, 2 are certainly in multiple star systems (Kurtovic et al., 2018). Since we only focus on single star systems here, the total number of disks is 18.

Since the gaps in protoplanetary disks may not be due to young planets, our derived planet occurrence rates should be considered as the upper limits. On the other hand, we may miss planets at the mass detection limit ( $\sim$ Neptune mass), as evidenced by that we do not include those planets that are fitted by eye and have no error estimates. Thus, the planet occurrence rates for Neptune mass planets may be higher than our estimates.

By comparing with exoplanets discovered with other methods, we find that:

First, we only have one planet that is more massive than 5 $M_{J}$ . Thus, the occurrence rate for $>$ 5 $M_{J}$ planets beyond 5-10 au is 1/18 or 6%. Wide-orbit giant planets are very rare. This is consistent with the direct imaging constraints that the occurrence rate for 5-20 $M_{J}$ planets at $>$ 5-10 au is 1-10% (Meshkat et al., 2017; Vigan et al., 2017; Bowler & Nielsen, 2018).

Second, using disk features, we may be probing a planet population that is not accessible by other planet searching techniques. These are Neptune to Jupiter mass planets beyond 10 au. Young planetary systems may harbor Uranus and Neptune mass planets beyond 10 au similar to our Solar System. The occurrence rate for 0.2 $M_{J}\lesssim M_{p}\lesssim$ 5 $M_{J}$ planets beyond 5-10 au is 8/18 or 44%, and the occurrence rate for all the planets more massive than Neptune and less than 5 $M_{J}$ beyond 5-10 au is 10/18 or 56%. These rates are comparable to the 31% giant planet ( $>$ 0.1 $M_{J}$ ) occurrence rates (Clanton & Gaudi, 2014) within $10^{4}$ days ( $<$ 9 au for solar mass stars). If we consider that our derived planets spread from 5 au to 200 au, the occurrence rate per decade of semi-major axis is 27% and 35%, respectively. This rate is comparable to the occurrence rate (20%) for giant planets ( $>$ 0.1 $M_{J}$ ) with period between $10^{3}$ and $10^{4}$ days. Thus, giant planet distribution may be flat beyond several au to $\sim$ 100 au.

Finally, the planet’s mass distribution is almost flat from Neptune to Jupiter mass. We have $\sim$ 5 planets with 0.03 $M_{J}\lesssim M_{p}\lesssim$ 0.3 $M_{J}$ , and 6 planets with 0.3 $M_{J}\lesssim M_{p}\lesssim$ 3 $M_{J}$ .

We bin the planet masses in decade in part due to the number of sources available and in part because of the uncertainties of the mass range for each planet (see Figure 21). The uncertainties for most of the planet masses are around a factor of 10. We want to emphasize that the derived planet mass has larger uncertainties due to the unknown disk $\alpha$ and dust size distribution. On the other hand, as long as all these disks have similar $\alpha$ values among each other, the derived planet mass will systematically shift up and down with the same fraction (e.g. decreasing the $\alpha$ value by a factor of 10 will decreasing the planet mass by a factor of two for all the planets).

5 Discussion

5.1 Our Solar System and HR 8799 Analogs in Taurus

Exoplanetary systems are very diverse with systems having multiple low-mass planets within 1 au (as probed by the $Kepler$ spacecraft) or systems having multiple giant planets beyond 10s of au (e.g. HR 8799). Our solar system has both terrestrial and giant planets. Are any of the DSHARP sources analogous to our Solar System when it was young? Is DSHARP capable of detecting young Solar System analog or HR 8799 analog?

To answer these questions, we embed planets in our Solar System and HR 8799 into a protoplanetary disk having a minimum mass solar nebulae surface density

[TABLE]

To maximize our chances to detect disk features, we use DSD2 dust size distribution ( $s_{max}$ =1 cm). The initial dust-to-gas mass ratio is 1/100. We run simulations with both $\alpha=10^{-2}$ and $10^{-4}$ to explore the parameter space slightly. The mass of the HR 8799 central star is 1.47 M*⊙*, and the four giant planets in HR 8799 are chosen as 7 $M_{J}$ at 14.5 au, 7 $M_{J}$ at 24 au, 7 $M_{J}$ at 38 au, and 5 $M_{J}$ at 68 au Marois et al. (2010). The inner and outer boundary of these simulations are 0.1 $r_{0}$ and 10 $r_{0}$ , where $r_{0}$ = 10 au for two young solar system runs and $r_{0}$ = 20 au for two HR 8799 runs. The $\alpha=10^{-4}$ run for the solar system has 1500 and 2048 grid points in the radial and $\theta$ direction, whereas the three other models have 750 and 1024 grids in the radial and $\theta$ direction. The Solar System simulation runs for $\sim$ 500 orbits at 10 au (due to the higher resolution and computational cost) and the HR 8799 simulation runs for $\sim$ 1000 orbits at 20 au. The mm intensity images are calculated using the temperature structure from Equation 12 with luminosities at 1 Myr found from D’Antona & Mazzitelli (1994) given current masses. Before making the ALMA synthetic images, the dust emission for the young solar system and HR 8799 runs are convolved with a 2-D Gaussian FWHM 1.4 au and 2.8 au, respectively.

Then, we use the CASA simobserve task to generate synthetic observations with sensitivities and angular resolutions comparable to those of the DSHARP observations, which are shown in Figure 22. The angular resolutions in FWHM are equivalent to $\sim$ 5 au in distance and are marked in the lower left corners in the figure. Each set of synthetic observations consist of 12 minutes of on-source integration time with the Cycle 5 C43-5 antenna configuration, 35 minutes on source in the C43-8 configuration, and 35 minutes on-source in the C43-9 configuration. A precipitable water vapor level of 1.0 mm is adopted throughout. The resulting synthetic visibilities are imaged in the same manner as the DSHARP sources, as described in Andrews et al. (2018). Clearly the DSHARP observational setup is capable of detecting both our Solar System analogs and HR 8799 analogs at a distance of 140 pc away.

The four giant planets induce a wide gap in the HR 8799 analog. When the disk viscosity is high ( $\alpha=10^{-2}$ ), the disk has an annular ring with an inner cavity, similar to transitional disks (Espaillat et al., 2014). When the disk viscosity is low ( $\alpha=10^{-4}$ ), we see bright arcs. We also see bright sources at the inner disk, which are vortices at the gap edge between the adjacent pair of planets and the horseshoe region of the planets. In actual observations, we may misinterpret them as planets or circumplanetary disks. One way to distinguish these possibilities is studying if the bright sources are spatially resolved (Zhu et al., 2018). Either the planet or circumplanetary disks should be smaller than the planet’s Hill radius. If the structures within the gap are spatially resolved, it is likely that they are not from the planets or the circumplanetary disks.

For the Solar System analog, when the disk viscosity is high ( $\alpha=0.01$ ), we can only observe the gap induced by Jupiter. When the viscosity is low ( $\alpha=10^{-4}$ ), the common gap induced by Jupiter and Saturn can be seen. Gap edge vortices and horseshoe regions can also be seen in this case. From the synthetic observations, we can barely see the disk features induced by Uranus and Neptune. Even by examining the radial intensity profiles, we can only see an extremely shallow dimple at the Neptune position. Thus, Uranus and Neptune in our Solar System analogs are not detectable with DSHARP. The reason we have Neptune mass planet candidates in Table 3 and Figure 21 is because either the planet is further away or the central stellar mass is lower (so that $q$ is larger and gaps are deeper).

5.2 Caveats

Although we seek to explain gaps with young planets, we want to point out that there are many other possible mechanisms to produce gaps and rings, such as ice lines (Zhang et al., 2015; Okuzumi et al., 2016), the dead zone transition (Pinilla et al., 2016), MHD zonal flows (Flock et al., 2015; Ruge et al., 2016), the secular gravitational instability (Takahashi & Inutsuka, 2014), disk winds (Bai, 2017; Suriano et al., 2018) and so on. On the other hand, quantitative predictions from these mechanisms are desired for the future so that we can test various ideas and understand the nature of these gaps and rings.

Another major caveat in this work is that we fit the gap profiles at 1000 planetary orbits. The gap depth and width $do$ change with time (Rosotti et al., 2016). To get a rigorous comparison between simulations and observations, we need to know when planets formed in the disk and how planets grew in time (Hammer et al., 2017), which we have little knowledge about. We can only assume that the gap opening timescale is similar to the disk lifetime. Although 1000 orbits at $\sim$ 100 au is close to the disk lifetime, it is only 10% of the disk lifetime for a planet at 20 au. A study similar to this work but also including the gap’s change with time is needed in future. On the other hand, we can do some analytical estimates on the relationship between the gap width and time. First, we do not expect that the gap profile can change dramatically over several thousand orbits if the disk has a large $\alpha$ (e.g. $\alpha>10^{-3}$ ) and small particles (e.g. $St<10^{-3}$ ). This is because, in these disks, the viscous timescale over the gap width is much shorter than 1000 planetary orbits and the gas disk has already reached the steady state. Small particles couple with the gas relatively well and their drift timescale is much longer than several thousand orbits. Dust turbulent diffusion with the large $\alpha$ can further smooth out dust features (Zhu et al., 2012). Second, for particles which are marginally coupled to the gas ( $St\gtrsim 10^{-2}$ ), they drift fast in the disk and we expect that the gap width will increase with time. As long as the gas profile is fixed (e.g. $\alpha\sim 10^{-3}$ ), particles will drift twice further away from the planet over twice amount of time. On the other hand, particles with twice $St$ will drift twice further way from the planet over the same amount of time. Thus, we expect that the gap width is proportional to $St\times t$ for fast drifting particles. We have done a test for disks at different orbits and with fast drifting dusts having different Stokes numbers using the Elias 24 simulation above. We find that if the gas profile is about the same, the time $t$ and the Stokes number $St$ indeed play the same role in widening the gap: the gap width at $2t$ is similar to the gap width at $t$ from particles with $2St$ . However, we have not explored the full parameter space, and the results may change with some other disk parameters. Especially, if $\alpha$ is small, the dramatic change in the gas profile with time will complicate the issue and break the degeneracy between $St$ and $t$ . A detailed study requires adding the time dimension in the parameter space and is beyond the scope of this paper.

Dust evolution and feedback to the gas is ignored in our study so that we can scale the simulations. In reality, particles are trapped at the gap edges which will promote its growth. When a significant amount of dust is trapped at the gap edge, the dust-to-gas feedback can affect the gap depth and width (Yang & Zhu in prep.) or even trigger streaming instability (Youdin & Goodman, 2005). A proper study with all these effects considered is difficult for 2-D numerical simulations. But it can be incorporated into 1-D dust evolutionary models.

We want to emphasize that, as shown in §4, it is straightforward to derive the planet mass assuming other dust size distributions besides DSD1 and DSD2. As shown in Figure 14, only the maximum Stokes number affects the gap profiles. Thus, we can calculate the Stokes number for any given dust size distribution, and then use the fits to derive the planet mass.

6 Conclusion

DSHARP provides a homogeneous sample of young protoplanetary disks showing a variety of substructures, e.g. rings, gaps, spirals, and small scale asymmetry (Andrews et al., 2018). If these substructures are induced by forming young planets, they are revealing a hidden young planet population which has not been probed by direct planet searching techniques.

To explore the potential planet population that is responsible to observed features in the DSHARP disks, we carry out two-dimensional hydrodynamical simulations including dust particles to study the relationships between the gap properties and the planet mass. We systematically study a grid of 45 gas models (as in §2.2), with three values of $\alpha$ ( $10^{-4}$ , $10^{-3}$ , $10^{-2}$ ), three values of $h/r$ (0.05, 0.07, 0.10), and five values of planet mass (from 10 $M_{\earth}$ to 3 $M_{J}$ ). For each model, we scale the dust distribution in the simulation to disks with different surface densities and different dust size distributions. Two different dust size distributions motivated by (sub-)mm polarization measurements (DSD1: $s_{max}$ =0.1 mm, $p$ =-3.5) and (sub-mm) dust thermal continuum observations (DSD2: $s_{max}$ =1 cm, $p$ =-2.5) are considered. Overall, for each model, we generate 12 millimeter images including 7 images using the DSD1 dust size distribution and 5 images using the DSD2 dust size distribution.

•

First, we study the gas structure in these 45 simulations. Overall, the gap becomes deeper with higher $q$ , smaller $h/r$ , and lower $\alpha$ . But when $q\gtrsim 3M_{J}$ in a low $\alpha$ disk, the gap edge becomes eccentric and the gap depth starts to decrease. These are all consistent with previous studies.

•

We study the sub/super-Keplerian motion at the gap edges. We confirm that the deviation from the Keplerian motion is due to the gas radial pressure gradient. The distance between the sub/super-Keplerian motion peaks is roughly 4.4 times $h$ , with a weak dependence on $\alpha$ and $q$ . The amplitude of the sub/super-Keplerian motion peaks is fitted with Equation 16, which shows a strong dependence on $h/r$ .

•

Then, we study the mm intensity maps for all our simulations. The gap edge becomes more eccentric and off-centered with the increasing planet mass. The eccentricity and off-centered distance are provided (Figure 8). Large eccentricity and off-centered distance may be indications of planets in disks.

•

Particle trapping in gap edge vortices and the horseshoe region are apparent in mm intensity maps for disks with $\alpha=10^{-4}$ , leading to large-scale asymmetries in the images. For some parameters, even a $33M_{\earth}$ planet can lead to a factor of 100 contrast between different azimuthal parts of the disk. In some cases, the vortex shows up at smaller radii than the gap edge (similar to the arc structure in HD 163296).

•

We derive several empirical relationships between the width/depth of the gaps in mm intensity maps and the planet/disk properties. All the fits for the width and depth are given in Table 1 and 2 and shown in Figure 14. We show that different disks with different surface densities and different dust size distributions have the same gap shape as long as their Stokes numbers for the maximum-size particles (where most of the dust mass is) are the same. Thus, our derived relationships can be used in other disks with different surface densities and dust size distributions.

•

A single planet can open multiple gaps. The position of the secondary gap is fitted with Equation 25. We find that the position of the secondary gap is almost solely determined by the disk scale height. Thus, if the secondary gap is present, we can use its position to estimate the disk scale height ( $h/r$ ).

•

With all these relationships, we lay out the procedure to constrain the planet mass using gap properties (the flowchart is presented in Figure 17).

•

Applying these steps, we identify potential planets in the DSHARP disks. We provide planet masses that are derived using three different values of $\alpha$ and three dust size distributions.

•

We comment on the potential planets in each disk. Particularly, for AS 209, we point out that our simulation matches not only the primary gap, but also the position and amplitude of the secondary (61 au), tertiary (35 au) and even the fourth (24 au) inner gaps. This makes AS 209 the most plausible case that there is indeed a planet within the 100 au gap (also in Guzmán et al. (2018)). The best fit model also suggests that the disk $\alpha$ increases with radii in AS 209, which may have implications for studying disk accretion theory.

•

We make synthetic observations for HR 8799 and Solar System analogs to show that DSHARP is capable of detecting giant planets in these systems.

•

We plot these potential young planets in the exoplanet mass-semimajor axis diagram (Figure 21). We find that the occurrence rate for $>$ 5 $M_{J}$ planets beyond 5-10 au is $\sim$ 6%, consistent with direction imaging constraints. Using disk features, we can probe a planet population which is not accessible by other planet searching techniques. These are Neptune to Jupiter mass planets beyond 10 au. The occurrence rate is $\sim$ 50%, suggesting a flat distribution beyond several au and planets with Neptune mass and above are common. On the other hand, we caution that there are large uncertainties for both the origin of these gaps and the inferred planet mass.

S. Z. and Z. Z. thank Lee Hartmann for very useful discussions. S. Z and Z. Z. thank the referee for prompt and constructive comments. Z. Z. acknowledges support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with Grant No. NNX17AK40G and Sloan Research Fellowship. Simulations are carried out with the support from the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG- AST130002. J.H. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144152. 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. 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. 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. 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. 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. L.P. acknowledges support from CONICYT project Basal AFB-170002 and from FCFM/U. de Chile Fondo de Instalación Académica. M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming disks). 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. This paper makes use of ALMA data ADS/JAO.ALMA #2016.1.00484.L.

Appendix

The fitted gap widths and depths for all the models are listed in Table 4 and Table 5. Column (4) shows the gap widths/depths of the gas; Column (5-11) show the gap widths/depths of the dust emission with increasing initial gas surface density $\Sigma_{g,0}$ (decreasing Stokes number $St_{max}$ ) under dust size distribution DSD1; similarly Column (12-16) show the gap widths/depths of the dust under DSD2. All widths/depths shown in Table 4 and 5 are derived from the images with a Gaussian convolution $\sigma$ = 0.06 $r_{p}$ (the larger kernel), except for the bottom of Table 4 (below the horizontal line and above the double horizontal lines) where widths are derived using $\sigma$ = 0.025 $r_{p}$ (the smaller kernel). These widths with a smaller beam are listed only if the gap widths $\Delta$ $<$ 0.15 using the larger kernel ( $\sigma$ = 0.06 $r_{p}$ ). Rows below the double lines show the individual widths of the gaps whose common gap is separated into two due to the horseshoe. The value on top the bar shows the width of the inner gap ( $\Delta_{1}$ ), whereas the value under the bar shows the width of the outer gap ( $\Delta_{2}$ ).

\startlongtable

Bibliography130

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, Ap J, 808, L 3, doi: 10.1088/2041-8205/808/1/L 3 · doi ↗
2Andrews et al. (2018) Andrews, S. M., Huang, J., & et al. 2018, Ap J
3Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, Ap J, 820, L 40, doi: 10.3847/2041-8205/820/2/L 40 · doi ↗
4Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A 33, doi: 10.1051/0004-6361/201322068 · doi ↗
5Ataiee et al. (2013) Ataiee, S., Pinilla, P., Zsom, A., et al. 2013, A&A, 553, L 3, doi: 10.1051/0004-6361/201321125 · doi ↗
6Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, Ap J, 863, 44, doi: 10.3847/1538-4357/aab 846 · doi ↗
7Bae et al. (2018) Bae, J., Pinilla, P., & Birnstiel, T. 2018, Ap J, 864, L 26, doi: 10.3847/2041-8213/aadd 51 · doi ↗
8Bae & Zhu (2018 a) Bae, J., & Zhu, Z. 2018 a, Ap J, 859, 118, doi: 10.3847/1538-4357/aabf 8c · doi ↗