Constraining Protoplanetary Disk Accretion and Young Planets Using ALMA Kinematic Observations
MMNRAS , 1–16 (2020) Preprint 16 February 2021 Compiled using MNRAS L A TEX style file v3.0
Constraining Protoplanetary Disk Accretion and Young PlanetsUsing ALMA Kinematic Observations
Ian Rabago, ★ Zhaohuan Zhu Department of Physics and Astronomy, University of Nevada, Las Vegas 4505 S. Maryland Parkway Las Vegas, NV 89154, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Recent ALMA molecular line observations have revealed 3-D gas velocity structure in proto-planetary disks, shedding light on mechanisms of disk accretion and structure formation. 1)By carrying out viscous simulations, we confirm that the disk’s velocity structure differs dra-matically using vertical stress profiles from different accretion mechanisms. Thus, kinematicobservations tracing flows at different disk heights can potentially distinguish different accre-tion mechanisms. On the other hand, the disk surface density evolution is mostly determinedby the vertically integrated stress. The sharp disk outer edge constrained by recent kinematicobservations can be caused by a radially varying 𝛼 in the disk. 2) We also study kinematicsignatures of a young planet by carrying out 3-D planet-disk simulations. The relationshipbetween the planet mass and the “kink” velocity is derived, showing a linear relationship withlittle dependence on disk viscosity, but some dependence on disk height when the planet ismassive (e.g. 10 𝑀 𝐽 ). We predict the “kink” velocities for the potential planets in DSHARPdisks. At the gap edge, the azimuthally-averaged velocities at different disk heights deviatefrom the Keplerian velocity at similar amplitudes, and its relationship with the planet mass isconsistent with that in 2-D simulations. After removing the planet, the azimuthally-averagedvelocity barely changes within the viscous timescale, and thus the azimuthally-averaged ve-locity structure at the gap edge is due to the gap itself and not directly caused to the planet.Combining both axisymmetric kinematic observations and the residual “kink” velocity isneeded to probe young planets in protoplanetary disks. Key words: accretion – accretion disks – astroparticle physics - dynamo - magneto-hydrodynamics (MHD) - instabilities - turbulence
Protoplanetary disks are thought to be the sites of planetary for-mation. In recent years, use of ground-based radio telescopes suchas the Atacama Large Millimeter Array (ALMA) has improvedour ability to resolve protoplanetary disks, revealing a wealth ofstructure in many existing disks, including rings, spirals, and othernon-axisymmetric structures (e.g. ALMA Partnership et al. 2015;Andrews et al. 2016, 2018). It is often suggested that the dark an-nular features observed in the dust continuum are gaps created byprotoplanets as they repel material away from the planet (e.g. reviewby Espaillat et al. 2014; Andrews 2020). Despite the high frequencyof gaps observed in high-resolution imagery, only a small num-ber of forming planets have been discovered to date. One notableprotoplanetary system is the disk of PDS 70, which contains twoplanetary-mass point sources detected by direct observation (Kep-pler et al. 2018; Wagner et al. 2018; Haffert et al. 2019; Christiaenset al. 2019; Isella et al. 2019; Wang et al. 2020). Many other disks ★ E-mail: [email protected] have been suspected to contain protoplanets due to the presence ofgaps and spirals within the disk, but there is much debate over theirexistence as these features are not a unique signature of planetaryformation (e.g. Follette et al. 2017; Rameau et al. 2017; Cugno et al.2019; Brittain et al. 2019).Besides planets, there are many other mechanisms that are ableto form gaps in protoplanetary disks. Gaps and other axisymmet-ric structures can be formed by condensation of materials at snowlines (Zhang et al. 2015; Okuzumi et al. 2016), interaction of thedisk material with magnetic fields (Kretke & Lin 2007; Johansenet al. 2009; Suriano et al. 2017; Bai & Stone 2014; Hu et al. 2019,2020), dust-gas secular gravitational instability (Takahashi & Inut-suka 2014; Tominaga et al. 2020), and indirectly from the spiralarm generated by a planet orbiting at a different radius (Bae et al.2017; Dong et al. 2017). Determining which of the observed gapsare hosts to potential planets is an ongoing process and an area ofactive research.Thanks to ALMA’s high sensitivity to molecular lines, sub-thermal kinematic motion in protoplanetary disks can be probedusing molecular lines, and such kinematic information provides © a r X i v : . [ a s t r o - ph . E P ] F e b I. Rabago and Z. Zhu more direct evidence of the planets in disks (Teague et al. 2018;Pinte et al. 2018). To extract the disk kinematic information frommolecular lines, two different methods have been developed. Pinteet al. (2018) uses the high spectral-resolution channel maps to probethe disk region whose velocity structure is significantly disturbedby the planet. Such a non-axisymmetric disturbance is due to thehorseshoe and circumplanetary motion around the planet (Perezet al. 2015; Pérez et al. 2018), forming a “kink” in the channelmaps. The second method, used by Teague et al. (2018), assumesthat the disk flow is mostly axisymmetric, allowing the axisymmet-ric flow velocity to be calculated by averaging the kinematic in-formation along the disk’s azimuthal direction. With this averagingprocess, sub-thermal axisymmetric deviations from the Keplerianmotion can be extracted. Teague et al. (2018) analyzed the rota-tional velocity structure in the gaps of disk HD 163296 and foundthat the super-Keplerian and sub-Keplerian rotational velocity atthe gap edges could be explained by planets of 1 𝑀 𝐽 at 100 AU and1 . 𝑀 𝐽 at 165 AU. A later work (Teague et al. 2019) examines thetwo-dimensional velocity components of HD 163296, and finds ev-idence of meridional flows in the atmosphere above the gaps. Suchmeridional flow also resembles the flow structure in 3-D planet-diskinteraction simulations (Fung & Chiang 2016).The relationship between the planet and the “kink” velocityhas not been derived. Pinte et al. (2018) and Pinte et al. (2019)derived the planet mass by comparing the observed channel mapswith synthetic channel maps from direct numerical simulations.Although such an approach provides a more robust estimate on theplanet mass, it is time consuming and computationally expensivefor a large sample of disks (e.g. Pinte et al. 2020).The relationship between the axisymmetric deviation from theKeplerian motion at the gap edge and the embedded planet mass hasbeen derived in Zhang et al. (2018) and Gyeol Yun et al. (2019) usinga grid of 𝑅 − 𝜙 𝑅 − 𝜙 Σ ∝ 𝑒𝑥 𝑝 (−( 𝑅 / 𝑅 𝑑 ) ) . Dullemond et al. (2020) suggest thatit is due to outside-in photoevaporation or truncation by an unseencompanion. On the other hand, the sharp outer edge could be dueto the disk evolution itself, as we will argue in §4.1.In this paper, we use 3-D viscous simulations to study therelationship between the planet mass and the disk’s 3-D velocitystructure (both the “kink” velocity and the axisymmetric deviationfrom the Keplerain velocity at the gap edge). We also explore howprobing the velocity structure could differentiate different accretionprocesses in protoplanetary disks. The disk could have complicatedflow structure (e.g. meridional flow) even without gaps (Urpin 1984;Kley & Lin 1992; Rozyczka et al. 1994). We will address how suchflow pattern is related to the stress distribution in the disk, and howthe disk surface density evolution depends on the stress distribution.Finally, we study the differences between the velocity structure ofa gap which has a planet in it and that of a gap without a planet. InSection 2 we outline the viscous theory important to this work. Wedescribe our simulation setup in Section 3, and present our resultsin Section 4. We discuss important results and implications of ourfindings in Section 5. Finally, we conclude in Section 6. Before we carry out viscous simulations for planet-disk interaction,we will first study the evolution of viscous disks without planetsfor two reasons. First, there is a well developed analytical theoryon viscous disk structure/evolution, which can be used to test theviscosity module for our numerical simulations. Second, by mod-ifying the traditional viscous disk theory, we can derive the flowstructure of disks undergoing different instabilities, and how kine-matic observations may be able to constrain these different accretionmechanisms.In this section, a few analytical results from the viscous disktheory are reviewed, starting from one-dimensional surface den-sity evolution to two-dimensional meridional circulation. Variousstress profiles have been adopted to mimic the turbulent structurein disks subject to different disk instabilities. Then, motivated bysome recent MHD simulations showing that the disk mostly accretesin the radial direction under the spherical-polar coordinate system,we analytically derive the meridional circulation pattern under thespherical-polar coordinate system. This new analytical derivationcan also be compare with our spherical-polar simulations more di-rectly in Section 3.
If the viscous disk has an internal viscosity of 𝜈 , the disk’s surfacedensity ( Σ ) follows the diffusion equation 𝜕 Σ 𝜕𝑡 = 𝑅 𝜕𝜕𝑅 (cid:20) 𝑅 𝜕𝜕𝑅 (cid:16) 𝑅 𝜈 Σ (cid:17)(cid:21) . (1)The analytical solution of this equation has been derived in Lynden-Bell & Pringle (1974). For a Keplerian disk with 𝜈 ∝ 𝑅 𝛾 , (2)one similarity solution is (we adopt the derivation from Hartmannet al. 1998) Σ ( 𝑟, 𝑡 ) = 𝐶 𝜋𝜈 𝑇 − / − 𝛾 − 𝛾 𝑒𝑥 𝑝 (cid:32) − ( 𝑅 / 𝑅 ) − 𝛾 𝑇 (cid:33) (3) (cid:164) 𝑀 ( 𝑟, 𝑡 ) = 𝐶𝑇 − / − 𝛾 − 𝛾 𝑒𝑥 𝑝 (cid:32) − ( 𝑅 / 𝑅 ) − 𝛾 𝑇 (cid:33) × (cid:34) − ( − 𝛾 ) ( 𝑅 / 𝑅 ) − 𝛾 𝑇 (cid:35) (4)with 𝑇 = 𝑡 / 𝑡 𝑠 +
1, where the viscous timescale is 𝑡 𝑠 = ( − 𝛾 ) 𝑅 𝜈 , (5)where 𝜈 = 𝜈 ( 𝑅 ) . The scaling factors for the surface density,radius, and time in this similarity solution are 𝐶 , 𝑅 , and 𝑡 𝑠 re-spectively. Note that the surface density profile at one specific timefollows Σ ( 𝑅 ) ∝ 𝜈 − 𝑒𝑥 𝑝 (cid:16) − ( 𝑅 / 𝑅 ) − 𝛾 (cid:17) ∝ 𝑅 − 𝛾 𝑒𝑥 𝑝 (cid:16) − ( 𝑅 / 𝑅 ) − 𝛾 (cid:17) (6)where 𝑅 = 𝑅 𝑇 /( − 𝛾 ) . The widely used relationship Σ ∝ 𝑅 − 𝑒𝑥 𝑝 (− 𝑅 / 𝑅 ) (e.g. Hartmann 1998; Andrews et al. 2009) isthe special case with 𝛾 =
1. Considering that 𝜈 = 𝛼𝑐 𝑠 / Ω , 𝛾 = 𝛼 value in a disk having 𝑇 ∝ 𝑅 − / (closeto the temperature of a passively irradiated disk). MNRAS000
1. Considering that 𝜈 = 𝛼𝑐 𝑠 / Ω , 𝛾 = 𝛼 value in a disk having 𝑇 ∝ 𝑅 − / (closeto the temperature of a passively irradiated disk). MNRAS000 , 1–16 (2020) onstraining Disk Structure 𝑅 − 𝑧 Although the disk’s 1-D evolution is straightforward, the disk’s 2-D( 𝑅 − 𝑧 ) accretion structure is much more complicated. The disk’s 𝑅 − 𝑧 𝑅 − 𝑧 plane. The disk is unlikely to accreteat the same speed at different heights if the disk’s vertical structure isconsidered. The resulting flow pattern shows the meridional circu-lation. Assuming that the initial density profile at the disk midplaneis 𝜌 ( 𝑅, 𝑧 = ) = 𝜌 ( 𝑅 , 𝑧 = ) (cid:18) 𝑅𝑅 (cid:19) 𝑝 , (7)and the temperature varies radially (but constant on cylinders) 𝑇 ( 𝑅, 𝑧 ) = 𝑇 ( 𝑅 ) (cid:18) 𝑅𝑅 (cid:19) 𝑞 , (8)the hydrostatic equilibrium in the 𝑅 − 𝑧 plane requires that (e.g.Nelson et al. 2013) 𝜌 ( 𝑅, 𝑧 ) = 𝜌 ( 𝑅, 𝑧 = ) exp (cid:34) 𝐺 𝑀𝑐 𝑠 (cid:32) √︁ 𝑅 + 𝑧 − 𝑅 (cid:33)(cid:35) , (9)andv 𝜙 ( 𝑅, 𝑧 ) = v 𝐾 (cid:34) ( 𝑝 + 𝑞 ) (cid:18) 𝑐 𝑠 v 𝐾 (cid:19) + + 𝑞 − 𝑞𝑅 √︁ 𝑅 + 𝑧 (cid:35) / , (10)where 𝑐 𝑠 = √︁ 𝑝 / 𝜌 is the isothermal sound speed, v 𝐾 = Ω 𝐾 𝑅 = √︁ 𝐺 𝑀 ∗ / 𝑅 , and 𝐻 = 𝑐 𝑠 / Ω 𝐾 .The disk’s radial velocity under the cylindrical coordinate sys-tem can be derived using the angular momentum equation 𝜕𝜌𝛿 v 𝜙 𝜕𝑡 = − 𝑅 𝜕𝑅 𝑇 𝑅𝜙 𝜕𝑅 − 𝜌 v 𝑅 𝑅 𝜕𝑅 Ω 𝐾 𝜕𝑅 − 𝜕𝑇 𝜙𝑧 𝜕𝑧 − 𝜌 v 𝑧 𝜕𝑅 Ω 𝐾 𝜕𝑧 , (11)where 𝛿 v 𝜙 = v 𝜙 − 𝑅 Ω 𝐾 , and 𝑇 can represent the turbulent Reynoldsstress, Maxwell stress, or/and viscous stress. This equation appliesto not only axisymmetric flows but also non-axisymmetric flowssimply by replacing every term with the azimuthally averaged quan-tities.For steady flows (constant 𝜌𝛿 v 𝜙 ) or axisymmetric flows( 𝛿 v 𝜙 =0), we can derive 𝜌 v 𝑅 𝑅 𝜕𝑅 Ω 𝐾 𝜕𝑅 = − 𝑅 𝜕𝑅 𝑇 𝑅𝜙 𝜕𝑅 − 𝜕𝑇 𝜙𝑧 𝜕𝑧 . (12)after neglecting higher order terms ( ( 𝐻 / 𝑅 ) and above). For viscousfluid having 𝑇 𝑅𝜙 = 𝜇 𝑅𝜙 𝑅𝜕 Ω / 𝜕𝑅 and 𝑇 𝜙𝑧 = 𝜇 𝜙𝑧 𝜕 v 𝜙 / 𝜕𝑧 , if theviscosity is isotropic with 𝜇 𝑅𝜙 = 𝜇 𝜙𝑧 = 𝜇 , we have 𝜌 v 𝑅 𝑅 𝜕𝑅 Ω 𝐾 𝜕𝑅 = − 𝑅 𝜕 (cid:16) 𝑅 𝜇 𝜕 Ω 𝜕𝑅 (cid:17) 𝜕𝑅 − 𝜕 (cid:16) 𝜇 𝜕 v 𝜙 𝜕𝑧 (cid:17) 𝜕𝑧 . (13)If we assume 𝜇 = 𝜌𝜈 , 𝜈 = 𝛼𝑐 𝑠 / Ω , and 𝛼 = 𝛼 ( 𝑅 / 𝑅 ) 𝑠 , and thenplug in the density and velocity profiles, we can derivev 𝑅 = − 𝜈𝑅 (cid:20) 𝑠 + 𝑝 + 𝑞 + + 𝑞 + (cid:16) 𝑧𝐻 (cid:17) (cid:21) . (14)as in Urpin (1984); Kley & Lin (1992); Rozyczka et al. (1994);Takeuchi & Lin (2002). Thus, if 3 𝑠 + 𝑝 + 𝑞 + <
0, the flowis outwards at the disk midplane. For a commonly adopted diskstructure with 𝑠 = 𝑝 = − .
25, and 𝑞 = − .
5, the quantity3 𝑠 + 𝑝 + 𝑞 + = − .
75 so that the flow is outwards at the midplane. Since 5 𝑞 + 𝑧 / 𝐻 is large enough.On the other hand, turbulence in protoplanetary disks may behighly anisotropic. For example, local MHD shearing box simu-lations suggest that turbulence induced by the magneto-rotationalinstability (MRI) generates much stronger 𝑅 − 𝜙 stress than the 𝜙 − 𝑧 stress (Hawley et al. 1995). However, hydrodynamical simulationssuggest that the 𝜙 − 𝑧 stress is much stronger than the 𝑅 − 𝜙 stress forthe vertical shear instability (VSI) (Nelson et al. 2013; Stoll et al.2017; Flock et al. 2020).If we can ignore the 𝜙 − 𝑧 stress (e.g. for MRI generatedturbulence), we can derivev 𝑅 = − 𝜈𝑅 (cid:20) 𝑠 + 𝑝 + 𝑞 + + 𝑞 + (cid:16) 𝑧𝐻 (cid:17) (cid:21) . (15)as in Fromang et al. (2011); Jacquet (2013); Philippov & Rafikov(2017). For cases with 𝜇 𝜙𝑧 significantly larger than 𝜇 𝑅𝜙 (e.g. VSI),we can assume 𝜇 𝜙𝑧 = 𝐶𝜇 𝑅𝜙 (which is equivalent to 𝜈 𝜙𝑧 = 𝐶𝜈 𝑅𝜙 or 𝛼 𝜙𝑧 = 𝐶𝛼 𝑅𝜙 ) and we can derivev 𝑅 = − 𝜈 𝑅𝜙 𝑅 (cid:20) 𝑠 + 𝑝 + 𝑞 + + 𝑞 + (cid:16) 𝑧𝐻 (cid:17) − 𝐶𝑞 (cid:18) − (cid:16) 𝑧𝐻 (cid:17) (cid:19)(cid:21) . (16)where 𝜈 𝑅𝜙 = 𝜇 𝑅𝜙 / 𝜌 is the viscosity generated from the 𝑅 − 𝜙 stress.At the disk midplane, this equation recovers Equation 12 inStoll et al. (2017). When 3 𝑠 + 𝑝 + 𝑞 + − 𝐶𝑞 <
0, the flow isinwards at the midplane. With 𝑞 = − . 𝐶 =
650 (Stoll et al.2017) for VSI due to the significant motion in the vertical direction(Lin & Youdin 2015), the disk flows inwards at the midplane andoutwards at 𝑧 (cid:38) 𝐻 .However, all above derivations are based on the 𝑅 − 𝑧 cylin-drical coordinate system. Some instabilities (e.g. VSI) show cleardistinction between 𝑅 − 𝜙 and 𝜙 − 𝑧 stresses, so that a cylindricalcoordinate system is more appropriate to study disk accretion inthese disks (Stoll et al. 2017). On the other hand, some other disks(e.g. MHD disks) have accretion flow which is in a direction morealigned with the radial direction in the spherical-polar coordinatesystem (e.g. Zhu & Stone 2018), and the 𝑟 − 𝜙 stress and 𝜃 − 𝜙 stressplay more distinct roles on disk accretion. This motivates us to de-rive the meridional circulation under the spherical-polar coordinatesystem, which also makes it easier to compare with our spherical-polar numerical simulations in Section 3. The angular momentumequation under the spherical-polar coordinate system is 𝜕𝜌𝛿 v 𝜙 𝜕𝑡 = − 𝑟 𝜕𝑟 𝑇 𝑟 𝜙 𝜕𝑟 − 𝜌 v 𝑟 𝑟 𝜕𝑟 Ω 𝐾 𝜕𝑟 − 𝑟 sin 𝜃 𝜕 sin 𝜃𝑇 𝜃 𝜙 𝜕𝜃 − 𝜌 v 𝜃 𝑟 sin 𝜃 𝜕 sin 𝜃𝑟 Ω 𝐾 𝜕𝜃 , (17)where v 𝐾 = 𝑟 Ω 𝐾 . Throughout the paper, to distinguish the radialdirection between the cylindrical and spherical-polar coordinatesystems, we use 𝑅 to represent the radial direction in the cylindricalcoordinate system, and 𝑟 for the radial direction in the spherical-polar coordinate system. For steady flows with second and higherorder terms removed, we have 𝜌 v 𝑟 𝑟 𝜕𝑟 Ω 𝐾 𝜕𝑟 = − 𝑟 𝜕𝑟 𝑇 𝑟 𝜙 𝜕𝑟 − 𝑟 sin 𝜃 𝜕 sin 𝜃𝑇 𝜃 𝜙 𝜕𝜃 . (18)Although solving the same fluid equations under different co-ordinate systems won’t change the results, 𝑇 𝜃 𝜙 and 𝑇 𝜙𝑧 are verydifferent at the disk surface due to the different curvatures of these MNRAS , 1–16 (2020)
I. Rabago and Z. Zhu two systems, and ignoring 𝑇 𝜃 𝜙 will lead to a different flow structurethan Equation 15 which ignores 𝑇 𝜙𝑧 . If we ignore the 𝜃 - 𝜙 stress anduse 𝑇 𝑟 𝜙 = 𝜇𝑟𝜕 Ω 𝐾 / 𝜕𝑟 , the disk’s radial velocity is 𝜌 v 𝑟 𝑟 𝜕𝑟 Ω 𝐾 𝜕𝑟 = − 𝑟 𝜕 (cid:16) 𝑟 𝜇 𝜕 Ω 𝐾 𝜕𝑟 (cid:17) 𝜕𝑟 . (19)Assuming 𝜇 = 𝜌𝜈 and 𝛼 = 𝛼 ( 𝑟 / 𝑟 ) 𝑠 in the spherical-polar coor-dinate system, we can derivev 𝑟 = − 𝜈𝑟 (cid:20) 𝑠 + 𝑝 + 𝑞 + + 𝑞 + (cid:16) 𝑧𝐻 (cid:17) (cid:21) . (20)By comparing Equation 15 with Equation 20, we can see that ig-noring 𝑇 𝜃 𝜙 or 𝑇 𝜙𝑧 leads to different flow structures. It is importantto choose the appropriate coordinate system so that the differentcomponents of the stress from turbulence can be more clearly sep-arated (e.g. some instabilities may generate turbulence with 𝑇 𝜃 𝜙 =0but non-zero 𝑇 𝜙𝑧 .).One particular stress profile that is motivated by MHD simu-lations is that the stress is vertically uniform until a certain numberof disk scale heights (z= ℎ 𝑐𝑢𝑡 H) (Fromang et al. 2011). Thus, wewrite 𝜇 in Equation 19 as 𝜇 = 𝛼𝜌 𝑐 𝑠 Ω 𝐾 √ 𝜋 ℎ 𝑐𝑢𝑡 | 𝑧 | ≤ ℎ 𝑐𝑢𝑡 𝐻 ℎ 𝑐𝑢𝑡 𝐻 (21)so that the vertically integrated stress is still ∫ 𝜇𝑑𝑧 = 𝛼 Σ 𝑐 𝑠 / Ω . If 𝜇 can be written as 𝜇 ( 𝑟 / 𝑟 ) 𝜂 within ℎ 𝑐𝑢𝑡 𝐻 , the radial velocitywithin ℎ 𝑐𝑢𝑡 𝐻 isv 𝑟 = −( . + 𝜂 ) 𝜇 (cid:18) 𝑟𝑟 (cid:19) 𝜂 𝑟 𝜌 . (22)where 𝜂 = 𝑠 + 𝑝 + 𝑞 + .
5. Comparing Equation 22 with Equation20, we can see that the uniform stress case has a sharp increaseof velocity at the disk surface (exponential increase due to the 1/ 𝜌 dependence).Different assumptions on the vertical stress profiles due to dif-ferent accretion mechanisms lead to dramatically different merid-ional flow structures. The analytical solutions for the meridionalcirculation (Equations 14, 15, 16, 20, 22) will be compared withour direct numerical simulations in Section 3 to mutually verify theanalytical derivation and numerical simulations. We solve the compressible Navier-Stokes equations using the grid-based magnetohydrodynamics (MHD) code ATHENA++ (Stoneet al. 2020). We adopt the spherical-polar coordinate system ( 𝑟, 𝜃, 𝜙 ) for the simulations. We first carry out 𝑟 − 𝜃 𝑟 = [ . , . ] , and 𝜃 = [ 𝜋 / − . , 𝜋 / + . ] . In the 𝑟 direction, we have 112 grid cellsthat are logarithmically uniformly spaced. In the 𝜃 direction, wehave 48 uniformly spaced cells. The disk setup follows Equations7 to 10 with 𝐻 / 𝑅 = . 𝑝 = − .
25 and 𝑞 = − .
5. Wehave adopted the reflecting boundary condition in the 𝜃 directionand fixed the inner boundary to the initial condition throughout the Figure 1.
Comparison of the disk radial velocity profiles at r=1 in directnumerical simulations with 𝛼 = .
01 (solid curves) to the analytic equations(dotted curves). The disk achieves this profile at a time of 30 orbits at r=1.The three dotted curves which have negative v 𝑟 at the disk surface areanalytical solution from Equation 14, 20, and 22 with 𝛼 = .
01. Thedotted curve with the positive v 𝑟 at the surface is from Equation 16 with 𝛼 = − and C=650 representing the stress derived in the vertical shearinstability simulations. The velocity structure is significantly affected bydifferent vertical stress profiles, although all the 𝛼 = .
01 cases have thesame vertically integrated 𝑟 − 𝜙 stress. simulation. The disk relaxes to the initial temperature at a coolingtime of 𝑡 𝑐 = .
01 orbital time (2 𝜋 / Ω , Zhu et al. 2015). Such coolingtime corresponds to the radiative cooling timescale at 100 au andis not short enough to trigger VSI (Lin & Youdin 2015). The 𝛼 parameter is 0.01. We have carried out three separate simulationswith the full viscous stress (Eq. 13), only the 𝑟 − 𝜙 stress, and thevertically uniform stress (Eq. 21 with ℎ 𝑐𝑢𝑡 =4). The radial velocitiesat 𝑟 = 𝑡 = 𝑇 are shown in Figure 1, where 𝑇 = 𝜋 / Ω is the orbital time at 𝑟 =
1. At 30 orbits, the meridional circulationhas been fully established. The three dotted curves which havenegative v 𝑟 at the surface in Figure 1 are the analytical solutionsfrom Equations 14, 20, 22 with the same parameters . We cansee good agreements between the simulations and the analyticalsolutions. For comparison, the dotted curve with the positive v 𝑟 at the surface is from Equation 16 with 𝛼 = − and C=650representing the stress similar to those derived from the verticalshear instability simulations (Stoll et al. 2017).After testing our code against analytical solutions, we carryout a set of simulations to study the disk surface density evolu-tion for disks having these different vertical stress profiles. Thesesimulations have similar setups as the previous tests except for afew modifications. To study the disk’s viscous spread, We setup aGaussian bump at 𝑟 = .
1. We alsoallow 𝛼 to change with radii ( 𝛾 ≠ 𝑟 = [ . , ] . We have 128 uniformly spaced grid cellsin the 𝜃 direction in the domain 𝜃 = [ 𝜋 / − . , 𝜋 / + . ] . Theoutflow boundary condition and reflecting boundary condition areadopted in the radial and 𝜃 direction. Our outflow boundary con-dition, which does not allow the gas to flow from the ghost zones 𝑧 and 𝑅 in these equations are calculated using 𝑟 and 𝜃 of each cell insimulations. MNRAS000
1. We alsoallow 𝛼 to change with radii ( 𝛾 ≠ 𝑟 = [ . , ] . We have 128 uniformly spaced grid cellsin the 𝜃 direction in the domain 𝜃 = [ 𝜋 / − . , 𝜋 / + . ] . Theoutflow boundary condition and reflecting boundary condition areadopted in the radial and 𝜃 direction. Our outflow boundary con-dition, which does not allow the gas to flow from the ghost zones 𝑧 and 𝑅 in these equations are calculated using 𝑟 and 𝜃 of each cell insimulations. MNRAS000 , 1–16 (2020) onstraining Disk Structure into the computational domain, is different from the default outflowboundary condition in Athena++. The simulations are run for 200 𝑇 . After studying the disk surface density evolution under differentstress profiles, we carry out 3-D viscous simulations to study thedisk’s velocity structure influenced by a young planet. The simula-tion domain covers a range of 𝑟 = [ . , . ] , 𝜃 = [ 𝜋 / − . , 𝜋 / + . ] , and 𝜙 = [ , 𝜋 ] . The domain is divided into 112 logarith-mically spaced cells in 𝑟 , 48 uniformly spaced cells in 𝜃 , and 152uniformly spaced cells in 𝜙 . The number of cells in the 𝜙 directionis half the number required to ensure square cells in the 𝑟 − 𝜙 plane;many quantities examined in this paper are azimuthally averaged,therefore this reduction in cells is expected to have a small effect.To initialize the density profile of the disk, we numericallyintegrate the density at each grid cell to establish vertical hydrostaticequilibrium. We initialize the sound speed profile using the power-law profile 𝑐 𝑠 = 𝑐 𝑠, (cid:18) 𝑟𝑟 (cid:19) − 𝑞 . (23)As with the 2-D simulations, we use power-law exponents of 𝑝 = − .
25 and 𝑞 = − . Σ ∝ 𝑟 − ,consistent with observations of older protoplanetary disks (Andrewset al. 2009). The disk again has a cooling time of 𝑡 𝑐 = .
01 orbitaltime and 𝐻 / 𝑅 = . 𝑟 . We use reflecting boundary conditions inthe 𝜃 direction and fix the radial boundary conditions to their initialvalues throughout the simulation.The disk velocity is initialized with only the rotational compo-nent v 𝜙 . To initialize the radial velocity profile, we allow the diskto evolve for a period of roughly 25 𝑇 and allow the gas to settleinto a steady state. With no initial velocity in the 𝑟 − 𝜃 plane, theradial velocity of the disk will settle towards Equation 14. Once theradial velocity of the disk is properly established, we add a planet ata distance of 𝑟 = 𝑟 𝑝 = 𝑀 𝑝 .The gravitational potential of the planet is written as a second orderpotential (e.g. Dong et al. 2011): Φ 𝑝 = − 𝐺 𝑀 𝑝 (cid:16) 𝑟 + 𝑟 𝑠 (cid:17) / , (24)where 𝑟 𝑠 is the smoothing radius. We choose the value of 𝑟 𝑠 to be0.1 Hill radii for each simulation. After reaching its final mass, theplanet continues to open a gap in the disk for another 450 planetaryorbits, giving the simulation a total time of 500 𝑇 .The 3-D simulations have the isotropic disk viscosity 𝛼 andthe final planet mass 𝑀 𝑝 as free parameters. We choose final planetmasses of 10 − and 10 − central star mass (which are 1 𝑀 𝐽 and10 𝑀 𝐽 if the central star is a solar mass star) in combination witha constant 𝛼 of 10 − and 10 − , creating a set of four differentsimulations. We also run a second set of simulations to study howthe choice of the smoothing length affects the unbound velocity flowin the vicinity of the planet, which we discuss in Section 4.2. Thesesimulations use the same disk viscosities of 𝛼 = − and 10 − ,final planet masses of 1 𝑀 𝐽 and 10 𝑀 𝐽 , but a constant smoothinglength of 2 grid cells. In order to examine the effect of MHD diskson meridian circulation, we run an additional simulation with a Figure 2.
Viscous disk evolution for disks with different stress structures.The Gaussian profile at r=r is the initial condition. All other curves are thedisk surface density at 200 𝑇 where 𝑇 is the orbital time at r=r . final planet mass of 1 𝑀 𝐽 and a vertically-varying 𝛼 according toEquation 21, with ℎ 𝑐𝑢𝑡 = Figure 1 shows that the disk’s velocity structure is significantly af-fected by different vertical stress profiles in disks without planets. Attwo disk scale heights ( 𝜃 ∼ . (cid:38) 𝑐 𝑠 , the 𝑟 − 𝜙 stress only case has an inward velocityof ∼ . 𝑐 𝑠 , the uniform stress case has an inward radial veloc-ity (cid:46) 𝑐 𝑠 , and the VSI case has an outward velocity (cid:38) 𝑐 𝑠 .However, at three disk scale heights, the inward radial velocity inthe uniform stress case quickly increases to (cid:38) 𝑐 𝑠 . Note thatall three cases with 𝛼 = .
01 in Figure 1 have the same verticallyintegrated 𝑟 − 𝜙 stress, which means that they have the same 𝛼 valuein the traditional 1-D viscous disk theory. If we also set 𝛼 = . 𝛼 valuein 1-D disk evolution. We caution that this traditional viscous stressmodel has limitations to capture internal stresses due to large scalemagnetic fields (e.g. Zhu & Stone 2018) or external stresses dueto magnetocentrifugal wind (Bai 2016). Nevertheless, future kine-matic observations using various molecular lines tracing differentdisk heights will not only measure the value of 𝛼 but also constrainthe detailed vertical stress profiles and accretion mechanisms.These largely different velocity structures with different verti-cal stress profiles raise the question as to whether the surface den-sity evolution of 3-D viscous disks are also affected by the detailedvertical stress profiles, so that the traditional 1-D disk evolution-ary model (Equation 3) cannot capture the disk evolution properly.Thus, we have carried out axisymmetric 2-D simulations with dif-ferent stress profiles in the disk to study the spread of a Gaussianbump. The results are shown in Figure 2. With the same verticallyintegrated stress, all three simulations with three different verticalstress profiles (black, blue, and red solid curves) show exactly the MNRAS , 1–16 (2020)
I. Rabago and Z. Zhu
Figure 3.
Radial surface density profiles for the 𝛼 =0.01 simulations (left)and 𝛼 = 0.001 simulations (right) at 𝑡 = 𝑇 . Blue curves indicate aplanet of 1 Jupiter mass, and orange curves indicate a planet of 10 Jupitermasses. same surface density structure at 200 𝑇 , following a power law 𝑅 − before the exponential decrease (Equation 6). This demon-strates that the detailed vertical stress profiles do not affect the disksurface density evolution, and the traditional 1-D viscous model issufficient for studying the disk surface density evolution.As shown in Figure 2, we have also carried out full stresssimulations whose vertically integrated stress varies along the radialdirection steeper or flatter than 𝜈 ∝ 𝑅 . If the vertically integratedstress changes slower with radii (e.g. 𝛾 = 𝛼 ∝ 𝑅 − with 𝑇 ∝ 𝑅 − / ), the surface density follows aflatter power law before a steeper exponential decrease, as predictedby 1-D model in Equation 6.Recent kinematic measurements from Dullemond et al. (2020)have suggested that, at the disk outer edge of HD 163296, thedisk surface density falls faster than 𝑒𝑥 𝑝 (− 𝑅 / 𝑅 𝑑 ) and is closer to 𝑒𝑥 𝑝 (−( 𝑅 / 𝑅 𝑑 ) ) . Although Dullemond et al. (2020) suggests thatthis sharper density drop is due to outside-in photoevaporation ortruncation by an unseen companion, we suggest that it can also bedue to the viscosity’s power law having 𝛾 < 𝛼 decreases withlarger radii. For example, when 𝛾 = 𝛼 ∝ 𝑅 − with 𝑇 ∝ 𝑅 − / ,Equation 6 becomes Σ ∝ 𝑅 𝑒𝑥 𝑝 (−( 𝑅 / 𝑅 ) ) (the dotted curve inFigure 2). This sharper drop is then consistent with observations inDullemond et al. (2020). Another implication of the decreasing 𝛼 with the increasing radius is that the density profile at the inner diskwould be a constant with radii ( ∝ 𝑅 ). Some disks indeed show arelatively flat surface density profile (e.g. Carrasco-González et al.2019). On the other hand, measurements of the disk surface densityat the inner disk could be very uncertain due to the optical deptheffects (Liu 2019; Zhu et al. 2019; Carrasco-González et al. 2019). After presenting the results on the disk structure without planets,we will study how the planet can influence the disk structure.Fig. 3 shows the disk surface density profiles at 𝑡 = 𝑇 when the gap is fully opened by the planet. Blue curves represent1 𝑀 𝐽 simulations while orange curves represent 10 𝑀 𝐽 simulations.The gray shaded regions denote the width of the planet’s Hill spherefor 1 and 10 Jupiter masses. Larger gaps are formed in the presenceof higher-mass planets or lower viscosity disks (e.g. Fung et al.2014; Kanagawa et al. 2015).Figures 4 and 5 show the velocity components of the disk forthe 𝛼 = . , 𝑀 𝐽 and 𝛼 = . , 𝑀 𝐽 simulations, respectively.Each row shows the disk velocities for a different scale height,with panels showing the vertical velocity v 𝑧 , the radial velocity v 𝑟 , and the deviation from Keplerian velocity 𝛿 v 𝜙 . The densitypanel (at top left) shows the gap and spiral wake. The streamlineof a horseshoe trajectory (starting at the red point) outlines the gapregion. To examine the velocity components across different diskscale heights, we cut the disk at equally spaced values of constant 𝜃 . All of our poloidal cuts are taken on the upper half of the disk;any cuts on the lower half would be almost identical, save for achange in the sign of v 𝑧 or v 𝜃 . We note that this is a simple way toanalyze quantities at different disk scale heights and may not reflectthe velocity structure at constant ℎ / 𝑟 throughout the disk.The vertical velocity v 𝑧 (left column) shows the presence ofinflows in the vicinity of the planet (Morbidelli et al. 2014; Funget al. 2015). The radial velocity (middle column) clearly followsthe spirals. For both the inner and outer spirals at the midplane,the disk region outside the spiral has negative v 𝑟 while the regionat 𝑟 smaller than the spiral front has positive v 𝑟 . The gap is notapparent in the v 𝑟 panels. On the other hand, regions of sub/super-Keplerian velocity, shown as 𝛿 v 𝜙 (right column), trace the edgesof the gap very well. The disk is super-Keplerian (positive 𝛿 v 𝜙 ) atthe outer gap edge, while it is sub-Keplerian (negative 𝛿 v 𝜙 ) at theinner gap edge. This is much clearer in Figure 5, as the velocitiesin the 1 𝑀 𝐽 simulation are dominated by the overall sub-Keplerianrotation of the gas. Meanwhile, the spirals are not as visible in the 𝛿 v 𝜙 panel. Overall, we conclude that v 𝑟 traces the spirals while 𝛿 v 𝜙 traces the gap. v 𝑧 has a smaller amplitude than either v 𝑟 or 𝛿 v 𝜙 . Since our temperature profile is nearly isothermal, we do notobserve additional spirals driven by vertical buoyancy resonances(Bae et al. 2021). However, when present these spirals will appearin the disk velocities as additional spiral signatures traced by the v 𝑟 and v 𝑧 components.We also observe some height dependence for all three velocitycomponents, especially for the 10 𝑀 𝐽 case which has a deep gap.Around the planet, v 𝑧 is more negative at 1 scale height than 2scale heights, implying a higher inflow velocity closer to the planet,while the amplitude of 𝛿 v 𝜙 is higher at the midplane decreasingtowards the surface. In the 10 𝑀 𝐽 case, v 𝑟 changes quite dramaticallyat different disk heights. The large v 𝑟 around the planet whichtraces the spirals at the midplane disappears at 1 scale height andonly slightly comes back at 2 scale heights, and the v 𝑟 tracing thespirals at the very inner disk becomes stronger at 2 scale heights.Thus, for massive planets in a deep gap, the kinematic signatures ofthe planet may show some differences by using various moleculartracers probing different depths in the disk.The region around the planet is significantly affected by theplanet in both v 𝑟 and 𝛿 v 𝜙 panels, which can be probed by thedistortion in the velocity channel maps (Pinte et al. 2018). Thecircumplanetary region can be divided into two regions: the regionthat is bound to the planet (e.g. the circumplanetary disk) and theunbound region which is circulating between the planet and thestar (Lubow et al. 1999). The bound circumplanetary disk has acircular motion around the planet, and for our simulations it is notcompletely resolved since the circumplanetary disk is quite small( (cid:46) 𝑅 𝐻 where 𝑅 𝐻 = ( 𝑀 𝑝 / 𝑀 ∗ ) / 𝑎 , Martin & Lubow 2011).Thus, we will only focus on the unbound region for studying thevelocity distortion.The large velocity in the unbound region is associated withthe spirals which are strongest close to the planet. We calculate thedeviation from Keplerian velocity as 𝛿 v 𝐾 = √︃ v 𝑟 + 𝛿 v 𝜙 around theplanet. In Figure 7, we plot v 𝑚𝑎𝑥 = max ( 𝛿 v 𝐾 ) for each simulationagainst the final planet mass 𝑀 𝑝 for several different scale heightsin the disk. This maximum velocity deviation should be close to MNRAS000
Radial surface density profiles for the 𝛼 =0.01 simulations (left)and 𝛼 = 0.001 simulations (right) at 𝑡 = 𝑇 . Blue curves indicate aplanet of 1 Jupiter mass, and orange curves indicate a planet of 10 Jupitermasses. same surface density structure at 200 𝑇 , following a power law 𝑅 − before the exponential decrease (Equation 6). This demon-strates that the detailed vertical stress profiles do not affect the disksurface density evolution, and the traditional 1-D viscous model issufficient for studying the disk surface density evolution.As shown in Figure 2, we have also carried out full stresssimulations whose vertically integrated stress varies along the radialdirection steeper or flatter than 𝜈 ∝ 𝑅 . If the vertically integratedstress changes slower with radii (e.g. 𝛾 = 𝛼 ∝ 𝑅 − with 𝑇 ∝ 𝑅 − / ), the surface density follows aflatter power law before a steeper exponential decrease, as predictedby 1-D model in Equation 6.Recent kinematic measurements from Dullemond et al. (2020)have suggested that, at the disk outer edge of HD 163296, thedisk surface density falls faster than 𝑒𝑥 𝑝 (− 𝑅 / 𝑅 𝑑 ) and is closer to 𝑒𝑥 𝑝 (−( 𝑅 / 𝑅 𝑑 ) ) . Although Dullemond et al. (2020) suggests thatthis sharper density drop is due to outside-in photoevaporation ortruncation by an unseen companion, we suggest that it can also bedue to the viscosity’s power law having 𝛾 < 𝛼 decreases withlarger radii. For example, when 𝛾 = 𝛼 ∝ 𝑅 − with 𝑇 ∝ 𝑅 − / ,Equation 6 becomes Σ ∝ 𝑅 𝑒𝑥 𝑝 (−( 𝑅 / 𝑅 ) ) (the dotted curve inFigure 2). This sharper drop is then consistent with observations inDullemond et al. (2020). Another implication of the decreasing 𝛼 with the increasing radius is that the density profile at the inner diskwould be a constant with radii ( ∝ 𝑅 ). Some disks indeed show arelatively flat surface density profile (e.g. Carrasco-González et al.2019). On the other hand, measurements of the disk surface densityat the inner disk could be very uncertain due to the optical deptheffects (Liu 2019; Zhu et al. 2019; Carrasco-González et al. 2019). After presenting the results on the disk structure without planets,we will study how the planet can influence the disk structure.Fig. 3 shows the disk surface density profiles at 𝑡 = 𝑇 when the gap is fully opened by the planet. Blue curves represent1 𝑀 𝐽 simulations while orange curves represent 10 𝑀 𝐽 simulations.The gray shaded regions denote the width of the planet’s Hill spherefor 1 and 10 Jupiter masses. Larger gaps are formed in the presenceof higher-mass planets or lower viscosity disks (e.g. Fung et al.2014; Kanagawa et al. 2015).Figures 4 and 5 show the velocity components of the disk forthe 𝛼 = . , 𝑀 𝐽 and 𝛼 = . , 𝑀 𝐽 simulations, respectively.Each row shows the disk velocities for a different scale height,with panels showing the vertical velocity v 𝑧 , the radial velocity v 𝑟 , and the deviation from Keplerian velocity 𝛿 v 𝜙 . The densitypanel (at top left) shows the gap and spiral wake. The streamlineof a horseshoe trajectory (starting at the red point) outlines the gapregion. To examine the velocity components across different diskscale heights, we cut the disk at equally spaced values of constant 𝜃 . All of our poloidal cuts are taken on the upper half of the disk;any cuts on the lower half would be almost identical, save for achange in the sign of v 𝑧 or v 𝜃 . We note that this is a simple way toanalyze quantities at different disk scale heights and may not reflectthe velocity structure at constant ℎ / 𝑟 throughout the disk.The vertical velocity v 𝑧 (left column) shows the presence ofinflows in the vicinity of the planet (Morbidelli et al. 2014; Funget al. 2015). The radial velocity (middle column) clearly followsthe spirals. For both the inner and outer spirals at the midplane,the disk region outside the spiral has negative v 𝑟 while the regionat 𝑟 smaller than the spiral front has positive v 𝑟 . The gap is notapparent in the v 𝑟 panels. On the other hand, regions of sub/super-Keplerian velocity, shown as 𝛿 v 𝜙 (right column), trace the edgesof the gap very well. The disk is super-Keplerian (positive 𝛿 v 𝜙 ) atthe outer gap edge, while it is sub-Keplerian (negative 𝛿 v 𝜙 ) at theinner gap edge. This is much clearer in Figure 5, as the velocitiesin the 1 𝑀 𝐽 simulation are dominated by the overall sub-Keplerianrotation of the gas. Meanwhile, the spirals are not as visible in the 𝛿 v 𝜙 panel. Overall, we conclude that v 𝑟 traces the spirals while 𝛿 v 𝜙 traces the gap. v 𝑧 has a smaller amplitude than either v 𝑟 or 𝛿 v 𝜙 . Since our temperature profile is nearly isothermal, we do notobserve additional spirals driven by vertical buoyancy resonances(Bae et al. 2021). However, when present these spirals will appearin the disk velocities as additional spiral signatures traced by the v 𝑟 and v 𝑧 components.We also observe some height dependence for all three velocitycomponents, especially for the 10 𝑀 𝐽 case which has a deep gap.Around the planet, v 𝑧 is more negative at 1 scale height than 2scale heights, implying a higher inflow velocity closer to the planet,while the amplitude of 𝛿 v 𝜙 is higher at the midplane decreasingtowards the surface. In the 10 𝑀 𝐽 case, v 𝑟 changes quite dramaticallyat different disk heights. The large v 𝑟 around the planet whichtraces the spirals at the midplane disappears at 1 scale height andonly slightly comes back at 2 scale heights, and the v 𝑟 tracing thespirals at the very inner disk becomes stronger at 2 scale heights.Thus, for massive planets in a deep gap, the kinematic signatures ofthe planet may show some differences by using various moleculartracers probing different depths in the disk.The region around the planet is significantly affected by theplanet in both v 𝑟 and 𝛿 v 𝜙 panels, which can be probed by thedistortion in the velocity channel maps (Pinte et al. 2018). Thecircumplanetary region can be divided into two regions: the regionthat is bound to the planet (e.g. the circumplanetary disk) and theunbound region which is circulating between the planet and thestar (Lubow et al. 1999). The bound circumplanetary disk has acircular motion around the planet, and for our simulations it is notcompletely resolved since the circumplanetary disk is quite small( (cid:46) 𝑅 𝐻 where 𝑅 𝐻 = ( 𝑀 𝑝 / 𝑀 ∗ ) / 𝑎 , Martin & Lubow 2011).Thus, we will only focus on the unbound region for studying thevelocity distortion.The large velocity in the unbound region is associated withthe spirals which are strongest close to the planet. We calculate thedeviation from Keplerian velocity as 𝛿 v 𝐾 = √︃ v 𝑟 + 𝛿 v 𝜙 around theplanet. In Figure 7, we plot v 𝑚𝑎𝑥 = max ( 𝛿 v 𝐾 ) for each simulationagainst the final planet mass 𝑀 𝑝 for several different scale heightsin the disk. This maximum velocity deviation should be close to MNRAS000 , 1–16 (2020) onstraining Disk Structure Figure 4.
Midplane cuts of the 𝛼 = . , 𝑀 = 𝑀 𝐽 simulation at 𝑡 = 𝑇 . Each row shows the velocity components of the disk at different diskheights (midplane, 1 scale height, and 2 scale heights). Top left:
Density plot, with a horseshoe trajectory also plotted within the gap region, beginning at thered marker.
Left column: vertical velocity v 𝑧 . Middle column:
Radial velocity. The inward/outward radial velocity trace the inner/outer spirals.
Right column:
Deviation from Keplerian velocity. Regions of sub-Keplerian and super-Keplerian velocity trace the inner and outer gap edges. All velocities are scaled to thelocal Keplerian velocity at 𝑅 = , 1–16 (2020) I. Rabago and Z. Zhu
Figure 5.
Same as Figure 4, but for the 𝛼 = . , 𝑀 = 𝑀 𝐽 simulation. MNRAS000
Same as Figure 4, but for the 𝛼 = . , 𝑀 = 𝑀 𝐽 simulation. MNRAS000 , 1–16 (2020) onstraining Disk Structure Figure 6.
Vertical slice of the 𝛼 = . , 𝑀 = 𝑀 𝐽 simulation at 𝑡 = 𝑇 . Background contours show gas density. Dashed circle on the right wedgedenotes the planet’s Hill radius. Vectors denote the gas velocity, scaled to the local sound speed. Figure 7.
Maximum deviation from Keplerian velocity v 𝑚𝑎𝑥 versus planetmass 𝑀 𝑝 , across multiple scale heights of the disk. v 𝑚𝑎𝑥 is scaled to thelocal Keplerian velocity. Triangle and pentagon markers correspond to diskviscosities of 𝛼 = .
001 and 𝛼 = .
01, respectively. The gray line representsa linear relationship between v 𝑚𝑎𝑥 and 𝑀 𝑝 , as described in the text. the maximum “kink” velocity measured in observations (e.g. Fig-ure 1 of Pinte et al. 2018). For this Figure, we plot the 1 𝑀 𝐽 and3 𝑀 𝐽 simulations with a smoothing length of 2 grid cells, as wellas the 10 𝑀 𝐽 simulations with the smoothing length scaled to theplanet mass. To ensure we do not measure the circumplanetary disk, we ignore values within the planet’s Hill sphere while calculatingv 𝑚𝑎𝑥 . We tried excluding different sized regions (e.g. 0.4 𝑅 𝐻 ) orusing simulations with different smoothing lengths, and found theresults to be almost identical. Above 3 scale heights, v 𝑚𝑎𝑥 increasesdramatically and no longer follows the relationship found at lowerscale heights. We believe that these large velocities are numerical innature, created by the gas density dropping to the simulation den-sity floor in the disk atmosphere. Thus, we did not plot the valuesbeyond 3 scale heights in Figure 7. For the case with the 10 𝑀 𝐽 planet in an 𝛼 = .
001 disk, we did not plot v 𝑚𝑎𝑥 at 2 scale heightsand above where the density reaches the floor inside the deep gap.Figure 7 shows that, for 1 and 3 𝑀 𝐽 cases, v 𝑚𝑎𝑥 does notdepend on 𝛼 or the disk height where it is measured. But for 10 𝑀 𝐽 cases, v 𝑚𝑎𝑥 decreases from the midplane to 2 scale heights andthen increases beyond 2 scale heights. This is also apparent fromFigure 5 where v 𝑟 becomes small at 1 and 2 scale heights.Our best fit for the values taken at the disk midplane is plottedas the gray line, which isv 𝑚𝑎𝑥 = 𝐾 𝑀 𝑝 𝑀 ∗ , (25)where v 𝐾 is the planet’s Keplerian velocity around the star. If weapply Equation 25 to the observation in Pinte et al. (2018), the 0.15v 𝐾 deviation in their Figure 1 corresponds to a 5.4 𝑀 𝐽 planet arounda 1.9 𝑀 (cid:12) star. This is higher than the 2 𝑀 𝐽 planet that Pinte et al.(2018) derived by comparing channel maps between observationsand direct simulations. Although this difference may be attributedto that we didn’t carry out synthetic channel maps, we notice thatour 𝛿 v 𝐾 around the planet ( 𝑀 𝑝 / 𝑀 ∗ = .
001 cases) is indeed lowerthan that shown in Pinte et al. (2018) (their Figure 4) by a factor of2. On the other hand, our values seem to be more consistent withPérez et al. (2018) (their Figure 1). Thus, more detailed comparisonsamong these simulations with the same parameters are needed infuture. We note that at the disk midplane, our fit in Figure 7 haslittle dependence on the choice of 𝛼 . Even if 𝛼 changes by one MNRAS , 1–16 (2020) I. Rabago and Z. Zhu order of magnitude ( 𝛼 = − and 10 − ), the maximum velocitiesof 𝛿 v 𝐾 for both 1 𝑀 𝐽 and 10 𝑀 𝐽 cases are almost identical. Thus,we can safely use this relationship to derive the embedded planetmass, with a minimum impact from the unknown disk viscosity.From this perspective, using the “kink velocity” serves as a morerobust method to derive the planet mass than using properties of thegap (e.g. Zhang et al. 2018) which depends on both 𝛼 and the gasscale height.The independence of v 𝑚𝑎𝑥 on 𝛼 indicates that the maximumvelocity of 𝛿 v 𝐾 is not related to the disk’s accretion process. In-stead, it is more likely determined by the planet-disk gravitationalinteraction. The linear dependence between v 𝑚𝑎𝑥 and 𝑀 𝑝 seemsto be consistent with the linear theory (Goldreich & Tremaine 1979)which shows that, for small amplitude perturbations, the velocityperturbation is proportional to the planet mass. However, the spiralscan quickly become highly non-linear when the planet mass is large(Goodman & Rafikov 2001; Muto et al. 2010; Dong et al. 2011;Zhu et al. 2013). The spirals steepen to shocks when they propagateaway from the planet at a distance of | Δ 𝑟 𝑠ℎ | ≈ . (cid:32) 𝛾 + / 𝐺 𝑀 𝑝 Ω 𝑝 𝑐 𝑠 (cid:33) − / 𝐻 , (26)where 𝛾 is the adiabatic index. Thus, the shocking distance is ∼ 𝐻 and ∼ . 𝐻 for the 1 𝑀 𝐽 and 10 𝑀 𝐽 planet in our 𝐻 / 𝑅 = . 𝑚𝑎𝑥 in the region beyond the planet’s Hillradius (which is 0.7 𝐻 and 1.5 𝐻 for 1 𝑀 𝐽 and 10 𝑀 𝐽 cases), themeasured spirals are in the non-linear regime (especially for 10 𝑀 𝐽 cases). Furthermore, a deep gap is induced by the massive planet,which clearly indicates that the linear theory cannot be applied to themassive planet case. Thus, we need to sort a different explanationfor the measured v 𝑚𝑎𝑥 with a massive planet. For massive planets(e.g. 10 𝑀 𝐽 ), the flow may be quite dynamic. The circumplanetaryregion is well separated from the background disk, and we expectthat the maximum velocity of disk material in a horseshoe orbitshould scale with the free fall velocity to the circumplanetary diskedge. Since the circumplanetary disk has a size of (cid:46) 𝑅 𝐻 , theunbound flow will have a maximum velocity at the closest approachat ∼ 𝑅 𝐻 , similar to the escape velocity at ∼ 𝑅 𝐻 , whichis v 𝑒𝑠𝑐 = . 𝐾 ( 𝑀 𝑝 / 𝑀 ∗ ) / . For a 10 𝑀 𝐽 planet around a solarmass star, v 𝑒𝑠𝑐 = . 𝐾 , which is close to our measurement of0.5 v 𝐾 . On the other hand, the escape velocity scales with 𝑀 / 𝑝 ,while our fitting shows a linear scaling with 𝑀 𝑝 . For a 𝑀 𝐽 planet,the escape velocity at 0.4 𝑅 𝐻 is higher than our measured v 𝑚𝑎𝑥 .We attribute this to the fact that the gas flow is significantly affectedby the gas pressure within the shallow gap induced by a low massplanet and the linear perturbation theory may apply. Figure 8 shows the velocity structure of the disk, azimuthally aver-aged in the 𝑟 − 𝜃 plane at 𝑡 = 𝑇 . All velocities are scaled relativeto the local sound speed. The poloidal velocity is shown as blackarrows, while the background color denotes the gas rotational speed.Blue and red areas denote regions of sub- and super-Keplerian ro-tation, respectively. Meridional flows are visible in all simulationsaround the planet’s Hill radius (denoted by the orange circle at 𝑅 = 𝜙 = 𝜙 𝑝 + 𝜋 , there is no such meridionalmotion present. Thus, understanding whether the observed merid-ional motion is truly axisymmetric in disks (e.g. similar meridionalpatterns at different values of 𝜙 ) or whether it only occurs in somelocal disk region is important for confirming its planet origin.A detailed examination of the azimuthally-averaged velocity isshown in Fig. 9. We examine each velocity component across theentire (cylindrical) radial extent and vertical extent of the disk, usingthe same polar cuts used in Figures 4 and 5 to examine differentscale heights. The disk scale heights are visible in Figure 8 as graydashed lines in the upper half of the disk. We plot the azimuthallyaveraged velocities at the end of the simulation ( 𝑡 = 𝑇 ), andremove all cells within the Hill sphere of the planet to remove localeffects from the planet’s vicinity. Since the number of 𝜃 cells in oursimulations is even, the midplane velocity is calculated by averagingthe two cells just above and below the midplane.The 𝛿 v 𝜙 panels (top panels) in Figure 9 show a vertical shiftas one moves vertically upward from the disk midplane, a featurewe attribute to the gas pressure gradient. The displacements aremore noticeable in disks with smaller perturbations relative to thedisk Keplerian speed (this is also true relative to the local soundspeed 𝑐 𝑠 , since ℎ / 𝑟 = 𝑐 𝑠 / v 𝜙 is assumed a constant vertically). Themagnitude of the velocity shift induced by the pressure gradient is atmost a few percent of the Keplerian velocity. In spite of the verticalshift, the shape of the curve is unchanged below 3 scale heights.By comparing these curves with the surface density plots (redrawnat the top of the figure), we see that the points of minimum andmaximum deviation correspond to the points of steepest pressuregradient, and the magnitude of the deviation is roughly equal to thevalue of the background pressure gradient at the gap edges, similarto that described in Teague et al. (2018). Above 3 scale heights, theresidual velocities in the inner disk become highly negative, largerthan expected for the given planet and disk parameters.Zhang et al. (2018) and Gyeol Yun et al. (2019) have stud-ied the deviation from Keplerian velocity produced at the gap edgeand how the velocities change with different disk parameters. Onthe other hand, these works are based on 2-D simulations, so thatthey focus on the velocity deviation at the disk midplane, whilemost ALMA molecular line observations actually probe the disksurface (Teague et al. 2018, 2019). These studies have found re-lationships between the parameters of the planet-disk system andthe amplitude of the velocity deviations Δ ( 𝛿 v 𝑟𝑜𝑡 ) (measured fromthe sub-Keplerian minimum interior to 𝑅 𝑝 to the super-Keplerianmaximum exterior to 𝑅 𝑝 ). The dependence of the planet-disk pa-rameters are commonly expressed in terms of a fitting parameter 𝐾 ,where 𝐾 = 𝑞 𝑝𝑎 (cid:16) ℎ𝑟 (cid:17) 𝑝𝑏 𝛼 𝑝𝑐 and the power-law exponents 𝑝𝑎 , 𝑝𝑏 ,and 𝑝𝑐 are to be determined.Gyeol Yun et al. (2019) run a parameter study using the fittingparameter 𝐾 = 𝑞 (cid:16) ℎ𝑟 (cid:17) − 𝛼 − from Kanagawa et al. (2015). Theyfind a fit of Δ ( 𝛿 v 𝑟𝑜𝑡 ) = (cid:18) ℎ𝑟 (cid:19) . 𝐾 . + . 𝐾 . (27) MNRAS000
001 cases) is indeed lowerthan that shown in Pinte et al. (2018) (their Figure 4) by a factor of2. On the other hand, our values seem to be more consistent withPérez et al. (2018) (their Figure 1). Thus, more detailed comparisonsamong these simulations with the same parameters are needed infuture. We note that at the disk midplane, our fit in Figure 7 haslittle dependence on the choice of 𝛼 . Even if 𝛼 changes by one MNRAS , 1–16 (2020) I. Rabago and Z. Zhu order of magnitude ( 𝛼 = − and 10 − ), the maximum velocitiesof 𝛿 v 𝐾 for both 1 𝑀 𝐽 and 10 𝑀 𝐽 cases are almost identical. Thus,we can safely use this relationship to derive the embedded planetmass, with a minimum impact from the unknown disk viscosity.From this perspective, using the “kink velocity” serves as a morerobust method to derive the planet mass than using properties of thegap (e.g. Zhang et al. 2018) which depends on both 𝛼 and the gasscale height.The independence of v 𝑚𝑎𝑥 on 𝛼 indicates that the maximumvelocity of 𝛿 v 𝐾 is not related to the disk’s accretion process. In-stead, it is more likely determined by the planet-disk gravitationalinteraction. The linear dependence between v 𝑚𝑎𝑥 and 𝑀 𝑝 seemsto be consistent with the linear theory (Goldreich & Tremaine 1979)which shows that, for small amplitude perturbations, the velocityperturbation is proportional to the planet mass. However, the spiralscan quickly become highly non-linear when the planet mass is large(Goodman & Rafikov 2001; Muto et al. 2010; Dong et al. 2011;Zhu et al. 2013). The spirals steepen to shocks when they propagateaway from the planet at a distance of | Δ 𝑟 𝑠ℎ | ≈ . (cid:32) 𝛾 + / 𝐺 𝑀 𝑝 Ω 𝑝 𝑐 𝑠 (cid:33) − / 𝐻 , (26)where 𝛾 is the adiabatic index. Thus, the shocking distance is ∼ 𝐻 and ∼ . 𝐻 for the 1 𝑀 𝐽 and 10 𝑀 𝐽 planet in our 𝐻 / 𝑅 = . 𝑚𝑎𝑥 in the region beyond the planet’s Hillradius (which is 0.7 𝐻 and 1.5 𝐻 for 1 𝑀 𝐽 and 10 𝑀 𝐽 cases), themeasured spirals are in the non-linear regime (especially for 10 𝑀 𝐽 cases). Furthermore, a deep gap is induced by the massive planet,which clearly indicates that the linear theory cannot be applied to themassive planet case. Thus, we need to sort a different explanationfor the measured v 𝑚𝑎𝑥 with a massive planet. For massive planets(e.g. 10 𝑀 𝐽 ), the flow may be quite dynamic. The circumplanetaryregion is well separated from the background disk, and we expectthat the maximum velocity of disk material in a horseshoe orbitshould scale with the free fall velocity to the circumplanetary diskedge. Since the circumplanetary disk has a size of (cid:46) 𝑅 𝐻 , theunbound flow will have a maximum velocity at the closest approachat ∼ 𝑅 𝐻 , similar to the escape velocity at ∼ 𝑅 𝐻 , whichis v 𝑒𝑠𝑐 = . 𝐾 ( 𝑀 𝑝 / 𝑀 ∗ ) / . For a 10 𝑀 𝐽 planet around a solarmass star, v 𝑒𝑠𝑐 = . 𝐾 , which is close to our measurement of0.5 v 𝐾 . On the other hand, the escape velocity scales with 𝑀 / 𝑝 ,while our fitting shows a linear scaling with 𝑀 𝑝 . For a 𝑀 𝐽 planet,the escape velocity at 0.4 𝑅 𝐻 is higher than our measured v 𝑚𝑎𝑥 .We attribute this to the fact that the gas flow is significantly affectedby the gas pressure within the shallow gap induced by a low massplanet and the linear perturbation theory may apply. Figure 8 shows the velocity structure of the disk, azimuthally aver-aged in the 𝑟 − 𝜃 plane at 𝑡 = 𝑇 . All velocities are scaled relativeto the local sound speed. The poloidal velocity is shown as blackarrows, while the background color denotes the gas rotational speed.Blue and red areas denote regions of sub- and super-Keplerian ro-tation, respectively. Meridional flows are visible in all simulationsaround the planet’s Hill radius (denoted by the orange circle at 𝑅 = 𝜙 = 𝜙 𝑝 + 𝜋 , there is no such meridionalmotion present. Thus, understanding whether the observed merid-ional motion is truly axisymmetric in disks (e.g. similar meridionalpatterns at different values of 𝜙 ) or whether it only occurs in somelocal disk region is important for confirming its planet origin.A detailed examination of the azimuthally-averaged velocity isshown in Fig. 9. We examine each velocity component across theentire (cylindrical) radial extent and vertical extent of the disk, usingthe same polar cuts used in Figures 4 and 5 to examine differentscale heights. The disk scale heights are visible in Figure 8 as graydashed lines in the upper half of the disk. We plot the azimuthallyaveraged velocities at the end of the simulation ( 𝑡 = 𝑇 ), andremove all cells within the Hill sphere of the planet to remove localeffects from the planet’s vicinity. Since the number of 𝜃 cells in oursimulations is even, the midplane velocity is calculated by averagingthe two cells just above and below the midplane.The 𝛿 v 𝜙 panels (top panels) in Figure 9 show a vertical shiftas one moves vertically upward from the disk midplane, a featurewe attribute to the gas pressure gradient. The displacements aremore noticeable in disks with smaller perturbations relative to thedisk Keplerian speed (this is also true relative to the local soundspeed 𝑐 𝑠 , since ℎ / 𝑟 = 𝑐 𝑠 / v 𝜙 is assumed a constant vertically). Themagnitude of the velocity shift induced by the pressure gradient is atmost a few percent of the Keplerian velocity. In spite of the verticalshift, the shape of the curve is unchanged below 3 scale heights.By comparing these curves with the surface density plots (redrawnat the top of the figure), we see that the points of minimum andmaximum deviation correspond to the points of steepest pressuregradient, and the magnitude of the deviation is roughly equal to thevalue of the background pressure gradient at the gap edges, similarto that described in Teague et al. (2018). Above 3 scale heights, theresidual velocities in the inner disk become highly negative, largerthan expected for the given planet and disk parameters.Zhang et al. (2018) and Gyeol Yun et al. (2019) have stud-ied the deviation from Keplerian velocity produced at the gap edgeand how the velocities change with different disk parameters. Onthe other hand, these works are based on 2-D simulations, so thatthey focus on the velocity deviation at the disk midplane, whilemost ALMA molecular line observations actually probe the disksurface (Teague et al. 2018, 2019). These studies have found re-lationships between the parameters of the planet-disk system andthe amplitude of the velocity deviations Δ ( 𝛿 v 𝑟𝑜𝑡 ) (measured fromthe sub-Keplerian minimum interior to 𝑅 𝑝 to the super-Keplerianmaximum exterior to 𝑅 𝑝 ). The dependence of the planet-disk pa-rameters are commonly expressed in terms of a fitting parameter 𝐾 ,where 𝐾 = 𝑞 𝑝𝑎 (cid:16) ℎ𝑟 (cid:17) 𝑝𝑏 𝛼 𝑝𝑐 and the power-law exponents 𝑝𝑎 , 𝑝𝑏 ,and 𝑝𝑐 are to be determined.Gyeol Yun et al. (2019) run a parameter study using the fittingparameter 𝐾 = 𝑞 (cid:16) ℎ𝑟 (cid:17) − 𝛼 − from Kanagawa et al. (2015). Theyfind a fit of Δ ( 𝛿 v 𝑟𝑜𝑡 ) = (cid:18) ℎ𝑟 (cid:19) . 𝐾 . + . 𝐾 . (27) MNRAS000 , 1–16 (2020) onstraining Disk Structure Figure 8.
Azimuthally averaged velocity structure for each of the 3D simulations. Black arrows denote the velocity flow in the ( 𝑟, 𝜃 ) plane. Background colorsdenote the deviation from normal Keplerian rotational velocity, with red colors indicating super-Keplerian flow and blue colors indicating sub-Keplerian flow.The orange circle denotes the planet’s Hill radius. Gray dashed lines mark the disk scale heights. Some vectors are omitted for clarity; see text for details. In a different parameter study by Zhang et al. (2018), both 𝐾 and Δ ( 𝛿 v 𝑟𝑜𝑡 ) are fit to the simulation data, with: 𝐾 v 𝑟 = 𝑞 (cid:18) ℎ𝑟 (cid:19) − . 𝛼 − . , (28)and Δ ( 𝛿 v 𝑟𝑜𝑡 ) = . 𝐾 . 𝑟 . (29)Both of these parameter studies used 2D hydrodynamic sim-ulations, and thus were unable to study how Δ ( 𝛿 v 𝑟𝑜𝑡 ) changes atdifferent disk scale heights. In Figure 10, we plot Δ ( 𝛿 v 𝑟𝑜𝑡 ) at dif- ferent disk heights with respect to 𝐾 and 𝐾 v 𝑟 in our simulations,together with the two fitting formulas (Equation 27 on the left andEquation 29 on the right). We can see that both formulas fit equallywell at large 𝐾 and 𝐾 v 𝑟 values. But Equation 27 fits slightly betterthan Equation 29 for small 𝐾 and 𝐾 v 𝑟 (which means smaller massplanets). Most importantly, Δ ( 𝛿 v 𝑟𝑜𝑡 ) is almost a constant within 3disk scale heights. For larger values of 𝐾 and 𝐾 v 𝑟 (e.g. disks withlarger planets, smaller scale heights, or lower values of 𝛼 ), the valueof Δ ( 𝛿 v 𝑟𝑜𝑡 ) begins to deviate above 3 scale heights. This indicatesthat the deviation from the Keplerian velocity increases towards thedisk surface faster in a deeper gap. MNRAS , 1–16 (2020) I. Rabago and Z. Zhu
Figure 9.
Azimuthally averaged velocity profiles for the 3D simulations. Each row plots a velocity component versus cylindrical
R across multiple scale heightsof the disk.
Top Row : differential azimuthal velocity. Gray curves are surface density curves shown in Fig. 3.
Middle Row : Radial velocity.
Bottom Row : Poloidalvelocity. All speeds are normalized to the local Keplerian speed. Different colors represent different scale heights in the disk. Shaded regions represent theextent of the planet’s Hill radius.
Figure 10.
The deviation amplitude Δ ( 𝛿 v 𝑟𝑜𝑡 ) plotted for different fitting parameters 𝐾 . Left : Using the fitting parameter from Kanagawa et al. (2015) andbest-fit equation from Gyeol Yun et al. (2019). Right: Using the fitting parameter and best-fit equations from Zhang et al. (2018). Best-fit equations are shownin gray. Colors are used to represent different scale heights, as in Fig. 9. The triangle and pentagon markers correspond to disk 𝛼 -viscosities of 10 − and 10 − ,respectively, following the convention used in Zhang et al. (2018). MNRAS000
The deviation amplitude Δ ( 𝛿 v 𝑟𝑜𝑡 ) plotted for different fitting parameters 𝐾 . Left : Using the fitting parameter from Kanagawa et al. (2015) andbest-fit equation from Gyeol Yun et al. (2019). Right: Using the fitting parameter and best-fit equations from Zhang et al. (2018). Best-fit equations are shownin gray. Colors are used to represent different scale heights, as in Fig. 9. The triangle and pentagon markers correspond to disk 𝛼 -viscosities of 10 − and 10 − ,respectively, following the convention used in Zhang et al. (2018). MNRAS000 , 1–16 (2020) onstraining Disk Structure At 4 scale heights, closer to the disk atmosphere, we find that Δ ( 𝛿 v 𝑟𝑜𝑡 ) departs strongly from the expected behavior and increasessubstantially to roughly order unity. Since 𝛿 v 𝜙 = v 𝜙 − 𝑅 Ω 𝐾 = v 𝜙 − v 𝐾 , the large values of Δ ( 𝛿 v 𝑟𝑜𝑡 ) suggests that 𝛿 v 𝜙 is changingby nearly the local Keplerian velocity in the vicinity of the gap.Examining the azimuthal velocities in Figure 9 shows that the super-Keplerian peak in 𝛿 v 𝜙 is (cid:46)
10% of Keplerian velocity at 4 scaleheights, suggesting that most of the contribution to Δ ( 𝛿 v 𝑟𝑜𝑡 ) is fromsub-Keplerian motion. Thus, a possible interpretation of this featureis that, in the vicinity of the gap, material in the disk atmosphere haslost nearly all of its rotational velocity and is instead falling directlytowards the star.Overall, as long as the molecular tracers (Isella et al. 2018) aretracing the disk region within ∼ ∼ 𝑟 𝑝 , there is adecrease in the radial velocities (even to negative values), whileexterior to 𝑟 𝑝 there is an increase. This is most apparent in the 𝑀 𝑝 = 𝑀 𝐽 simulations. These flows are strongest at the diskmidplane and weaken as one moves vertically in the disk; for mostof our simulations, these flows are no longer visible above 2 scaleheights. Some simulations show a reversed radial flow direction inthe disk atmosphere, with v 𝑟 positive interior to 𝑟 𝑝 and negativeexterior to 𝑟 𝑝 . This reversal may correspond to the material at thesurface flowing into the gap. The appearance of these local extremain the radial velocity channel is a signature unique to the planet andmay be important in determining the existence of a planet in anobserved gap. Section 5 examines this idea further by comparingthe velocity fields of planetary and non-planetary gaps.Gas inflows are also visible in the disk poloidal velocities (Fig.9, bottom row), as velocity peaks at the planet’s orbital radius (exceptfor the velocity at 4 scale heights). Since the data in Figure 9 aretaken from the upper half of the disk, the inflows appear as a positivebump near 𝑟 =
1. The inflow pattern is visible for several scaleheights, though the magnitude of these features is not as large as the 𝛿 v 𝜙 signature. We note that, for the 10 𝑀 𝐽 simulations, additionaldownward motion at 𝑅 ∼ . The velocity disturbances created by the planet in Figures 4 and5 trace features that span across the entire disk, some of whichare present across multiple scale heights. These disturbances arestrongest in the vicinity of the planet, and so good azimuthal res-olution is desirable in observations. For disks with high azimuthalresolution (roughly (cid:46) 𝑅 𝐻 ), deviations from Keplerian velocitycan be localized to specific areas in the disk, which has been used to infer the presence and position of planets using velocity channelmaps (Pinte et al. 2018, 2020).Given a mass for the embedded planet, Equation 25 can beused to predict the velocity deviation created in the disk. In Table 1,we calculate the expected velocity deviations (’kink velocity’) forplanets in selected disks. We use planet masses from Zhang et al.(2018). For the 22 AU gap in HD 143006, our calculated value for 𝛿 v is roughly consistent with the measured value obtained in Pinteet al. (2020). On the other hand, most planet masses constrainedby Zhang et al. (2018) for the DSHARP sample are less than 𝑀 𝐽 ,which corresponds to v 𝑚𝑎𝑥 / v 𝐾 (cid:46) Although Figure 10 links the planet mass with the gas kinematicsat the gap edge, we want to note that the gas kinematics measure-ments in Teague et al. (2018) and Teague et al. (2019) are actuallysignatures of a local decrease in the gas surface density, and theyare not direct probes for the embedded planet. For the gas, there isa radial force due to the pressure gradient 𝑑𝑃𝑑𝑟 , and the change in thepressure gradient due to the gap causes a change in the azimuthalvelocity v 𝜙 . The strength of the meridional flows observed in thedisk atmosphere is also a function of gap depth, as the flows anglemore strongly in the polar direction over the gap. However, neitherof these is a direct signature of an embedded planet; a gap createdby non-planetary means with a similar gap depth would exhibit thesame deviations in azimuthal velocity and gap inflows, despite noplanet existing within the gap.Thus, we want to compare the velocity fields in gaps createdby planetary and non-planetary methods in order to identify anyvelocity signatures that are unique to the planet. To do this, we extendour simulations by removing the planetary mass and simulating anadditional time of 50 𝑇 . Without the gravitational influence ofthe planet, any unique dynamic planetary signatures are dispersedwithin several orbits. This roughly mimics the gaps created by non-planetary methods which may operate more axisymmetrically, suchas snow lines and MRI, even though no additional physics has beenimplemented in our model. We note that the gap slowly closes overtens ( 𝛼 = .
01 cases) or hundreds ( 𝛼 = .
001 cases) of orbits dueto the viscous spreading.A snapshot of the original, “planetary” gap is compared tothe new, “non-planetary” gap in Figure 11 for the simulation with 𝛼 = .
001 and 𝑀 𝑝 = 𝑀 𝐽 . The “non-planetary” gap is chosen at10 orbits after removing the planet from the gap. The snapshot ofthe “planetary” gap is chosen such that the gap depth Σ / Σ in bothsnapshots is nearly the same. Though we only show one simulationin this figure, the behavior is roughly the same across all of oursimulations. The left column shows both 1D and 2D profiles forthe planetary gap, while the right column shows the same for thenon-planetary gap. Comparing the azimuthal velocities, we see thatthe deviations from normal rotational velocity are nearly identicalfor both gaps along the full vertical extent of the disk. Thus, a MNRAS , 1–16 (2020) I. Rabago and Z. Zhu
Disk Gap(AU) 𝑀 ∗ ( 𝑀 (cid:12) ) 𝑀 𝑝,𝑚𝑖𝑛 𝑀 𝑝 𝑀 𝑝,𝑚𝑎𝑥 𝛿 v 𝑚𝑖𝑛 𝛿 v 𝛿 v 𝑚𝑎𝑥 𝛿 v 𝑜𝑏𝑠 AS 209 9 0.83 0.37 3.38 4.18 0.0241 0.22 0.272 -AS 209 99 0.83 0.18 0.75 1.32 0.0117 0.0488 0.0859 -Elias 24 57 0.78 0.19 0.81 1.72 0.0132 0.0561 0.119 -Elias 27 69 0.49 0.02 0.1 0.12 0.0022 0.011 0.0132 ?GW Lup 74 0.46 0.01 0.06 0.06 0.00117 0.00704 0.00704 < 0.3HD 142666 16 1.58 0.09 0.5 0.62 0.00308 0.0171 0.0212 -HD 143006 22 1.78 2.35 9.81 40.6 0.0713 0.298 1.23 ≈ ≈ Table 1.
Estimated velocity deviations for selected disks using Equation 25, in units of the local Keplerian velocity. Distance and stellar mass data fromAndrews et al. (2018) and planet mass estimates from Zhang et al. (2018). Gap distances are measured in AU and used as the planet-star distance. Star massesare given in 𝑀 (cid:12) , and planet masses are given in 𝑀 𝐽 . The “middle” planet mass estimates are taken from the middle values of Column 13 in Zhang et al.(2018). Values in 𝛿 v 𝑜𝑏𝑠 are the maximum velocity deviations from Pinte et al. (2020), where applicable. measurement of 𝛿 v 𝜙 alone is not sufficient evidence to confirm theexistence of an embedded planet.Both the v 𝑟 and v 𝜃 components are also similar between the“planetary” gap and the “non-planetary” gap with the circulationpattern at the gap edge. The only noticeable difference is that, inthe planetary gap, both the v 𝑟 and v 𝜃 components at the midplane(black curves) show peaks in the vicinity of the planet’s orbital ra-dius. The v 𝑟 peak represents the repelling radial motion away fromthe planet which is probably driven by the spirals or the gravitationalinfluence of the planet, and the v 𝜃 peak represents infall onto theplanet. Thus these peaks are signatures of the planet itself. However,these signatures are roughly an order of magnitude smaller than the 𝛿 v 𝜙 signature, making their detection in velocity maps particularlychallenging. Their vertical extent is also quite limited; the outflowspresent in the v 𝑟 component weaken quickly with height (disap-pearing at one scale height), and so are likely not visible with COmeasurements. The inflows that cause a signature in the v 𝜃 com-ponent have a larger vertical extent, of 2-3 scale heights. Anotherplanetary signature which is not captured in our azimuthally aver-aged approach is the broader line width within the gap due to strongturbulence induced by the planet (Dong et al. 2019).We also compare our constant- 𝛼 simulations (with 𝛼 = . 𝛼 simulation, which uses the stress profile describedby Equation 21. Both simulations have 𝑀 𝑝 = 𝑀 𝐽 . The two modelsare very similar, with most of the global features remaining the same.The gap in the variable- 𝛼 model is somewhat shallower than theconstant- 𝛼 model, by a factor of about 1.4. The deviation amplitudein 𝛿 v 𝜙 is slightly smaller than the constant- 𝛼 simulation, and sothe measured value of Δ ( 𝛿 v 𝑟𝑜𝑡 ) would also be smaller. Close to 𝑧 = 𝐻 = ℎ 𝑐𝑢𝑡 𝐻 , the region where the stress profile changes, thegas velocity increases sharply. We have studied how ALMA kinematic observations can help us toconstrain 1) the accretion mechanisms, and 2) the planet properties.We have examined how different disk stress profiles affect itsvelocity structure, and how the velocity structure in turn affectsits subsequent evolution. We find that the radial velocity profileof the disk is very sensitive to the stress profile that is chosen.Thus, future kinematic observations using various molecular linestracing different disk heights will not only measure the value of 𝛼 but also constrain the detailed vertical stress profiles and accretionmechanisms.On the other hand, as long as the vertically integrated stressis the same, the evolution of the disk’s global surface density isunaffected. We also find that steep dropoffs at the outer edge of thedisk, which are normally explained by phenomenon external to thedisk such as truncation or photoevaporation, can also be explainedby a disk with a radially-varying 𝛼 profile.In our study of three-dimensional velocity flows of a planetopening a gap in the disk, we are able to see how different velocitycomponents follow different components of the planet-disk system.The radial velocity v 𝑟 traces out the spiral wake of the planet,while the sub-/super-Keplerian velocity 𝛿 v 𝜙 traces out the edgesof the gap. We observe some dependence of these features on thedisk height, especially for massive planet cases which have deepgaps. The linear relationship between the planet mass and the “kinkvelocity” is derived, and it is independent from the disk viscosityand the disk height (except for very massive planet cases). Usingsuch a relationship, we predict the “kink velocity” for the planetcandidates in the DSHARP sample.The velocity deviation at the gap edge within 3 disk scaleheights is consistent with previous 2D studies at the midplane.We see meridional circulation in the azimuthally averaged velocitymaps, and are able to identify components of the circulation processacross the vertical extent of the disk. However, by comparing thegap carved by the planet to a gap that is non-planetary in nature, wefind that the deviation from Keplerian rotation and the meridionalcirculation is a feature of the gap and not necessarily a signature ofthe planet itself. Examining the velocity field throughout the disk re-veals some velocity signatures that are unique to the planet, but theyare highly non-axisymmetric, limited in their vertical extent, and arean order of magnitude smaller than the deviations caused at the gapedge. Thus, they are not easily detectable in azimuthally averagedvelocity maps and may also be difficult to detect with current ob-servations. Combining both axisymmetric kinematic observationsand the residual “kink” velocity is needed to probe young planets inprotoplanetary disks. More specifically, we can derive the gaseousgap depth and width using the planet mass constrained by the “kink”velocity, and then study if the predicted velocity structure at the gapedge are consistent with axisymmetric kinematic observations. MNRAS000
Estimated velocity deviations for selected disks using Equation 25, in units of the local Keplerian velocity. Distance and stellar mass data fromAndrews et al. (2018) and planet mass estimates from Zhang et al. (2018). Gap distances are measured in AU and used as the planet-star distance. Star massesare given in 𝑀 (cid:12) , and planet masses are given in 𝑀 𝐽 . The “middle” planet mass estimates are taken from the middle values of Column 13 in Zhang et al.(2018). Values in 𝛿 v 𝑜𝑏𝑠 are the maximum velocity deviations from Pinte et al. (2020), where applicable. measurement of 𝛿 v 𝜙 alone is not sufficient evidence to confirm theexistence of an embedded planet.Both the v 𝑟 and v 𝜃 components are also similar between the“planetary” gap and the “non-planetary” gap with the circulationpattern at the gap edge. The only noticeable difference is that, inthe planetary gap, both the v 𝑟 and v 𝜃 components at the midplane(black curves) show peaks in the vicinity of the planet’s orbital ra-dius. The v 𝑟 peak represents the repelling radial motion away fromthe planet which is probably driven by the spirals or the gravitationalinfluence of the planet, and the v 𝜃 peak represents infall onto theplanet. Thus these peaks are signatures of the planet itself. However,these signatures are roughly an order of magnitude smaller than the 𝛿 v 𝜙 signature, making their detection in velocity maps particularlychallenging. Their vertical extent is also quite limited; the outflowspresent in the v 𝑟 component weaken quickly with height (disap-pearing at one scale height), and so are likely not visible with COmeasurements. The inflows that cause a signature in the v 𝜃 com-ponent have a larger vertical extent, of 2-3 scale heights. Anotherplanetary signature which is not captured in our azimuthally aver-aged approach is the broader line width within the gap due to strongturbulence induced by the planet (Dong et al. 2019).We also compare our constant- 𝛼 simulations (with 𝛼 = . 𝛼 simulation, which uses the stress profile describedby Equation 21. Both simulations have 𝑀 𝑝 = 𝑀 𝐽 . The two modelsare very similar, with most of the global features remaining the same.The gap in the variable- 𝛼 model is somewhat shallower than theconstant- 𝛼 model, by a factor of about 1.4. The deviation amplitudein 𝛿 v 𝜙 is slightly smaller than the constant- 𝛼 simulation, and sothe measured value of Δ ( 𝛿 v 𝑟𝑜𝑡 ) would also be smaller. Close to 𝑧 = 𝐻 = ℎ 𝑐𝑢𝑡 𝐻 , the region where the stress profile changes, thegas velocity increases sharply. We have studied how ALMA kinematic observations can help us toconstrain 1) the accretion mechanisms, and 2) the planet properties.We have examined how different disk stress profiles affect itsvelocity structure, and how the velocity structure in turn affectsits subsequent evolution. We find that the radial velocity profileof the disk is very sensitive to the stress profile that is chosen.Thus, future kinematic observations using various molecular linestracing different disk heights will not only measure the value of 𝛼 but also constrain the detailed vertical stress profiles and accretionmechanisms.On the other hand, as long as the vertically integrated stressis the same, the evolution of the disk’s global surface density isunaffected. We also find that steep dropoffs at the outer edge of thedisk, which are normally explained by phenomenon external to thedisk such as truncation or photoevaporation, can also be explainedby a disk with a radially-varying 𝛼 profile.In our study of three-dimensional velocity flows of a planetopening a gap in the disk, we are able to see how different velocitycomponents follow different components of the planet-disk system.The radial velocity v 𝑟 traces out the spiral wake of the planet,while the sub-/super-Keplerian velocity 𝛿 v 𝜙 traces out the edgesof the gap. We observe some dependence of these features on thedisk height, especially for massive planet cases which have deepgaps. The linear relationship between the planet mass and the “kinkvelocity” is derived, and it is independent from the disk viscosityand the disk height (except for very massive planet cases). Usingsuch a relationship, we predict the “kink velocity” for the planetcandidates in the DSHARP sample.The velocity deviation at the gap edge within 3 disk scaleheights is consistent with previous 2D studies at the midplane.We see meridional circulation in the azimuthally averaged velocitymaps, and are able to identify components of the circulation processacross the vertical extent of the disk. However, by comparing thegap carved by the planet to a gap that is non-planetary in nature, wefind that the deviation from Keplerian rotation and the meridionalcirculation is a feature of the gap and not necessarily a signature ofthe planet itself. Examining the velocity field throughout the disk re-veals some velocity signatures that are unique to the planet, but theyare highly non-axisymmetric, limited in their vertical extent, and arean order of magnitude smaller than the deviations caused at the gapedge. Thus, they are not easily detectable in azimuthally averagedvelocity maps and may also be difficult to detect with current ob-servations. Combining both axisymmetric kinematic observationsand the residual “kink” velocity is needed to probe young planets inprotoplanetary disks. More specifically, we can derive the gaseousgap depth and width using the planet mass constrained by the “kink”velocity, and then study if the predicted velocity structure at the gapedge are consistent with axisymmetric kinematic observations. MNRAS000 , 1–16 (2020) onstraining Disk Structure Figure 11.
Comparison of the disk velocity for "planetary" and "non-planetary" gaps in the 𝛼 = . 𝑀 = 𝑀 𝐽 simulation. Left column: snapshot with a planet opening the gap.
Right column: snapshot 10 orbitsafter removing the planet from the gap.
From top to bottom:
Azimuthallyaveraged velocity profile in the 𝑟 − 𝜃 plane, zoomed in on 𝑟 =
1; Azimuthallyaveraged azimuthal, radial, and poloidal velocity components, as presentedin Fig. 9; surface density profiles. The snapshots are chosen so that the gapdepths for the disks with and without the planet are nearly identical.
ACKNOWLEDGEMENTS
All simulations are carried out using computers supported by theTexas Advanced Computing Center (TACC) at the University ofTexas at Austin through XSEDE grant TG-AST130002 and fromthe NASA High-End Computing (HEC) program through the NASAAdvanced Supercomputing (NAS) Division at Ames Research Cen-ter. Z. Z. acknowledges support from the National Science Founda-tion under CAREER grant AST-1753168.
DATA AVAILABILITY
The data used in this paper is available upon request to the corre-sponding author.
REFERENCES
ALMA Partnership et al., 2015, ApJ, 808, L3Andrews S. M., 2020, arXiv e-prints, p. arXiv:2001.05007Andrews S. M., Wilner D. J., Hughes A. M., Qi C., Dullemond C. P., 2009,ApJ, 700, 1502Andrews S. M., et al., 2016, ApJ, 820, L40Andrews S. M., et al., 2018, ApJ, 869, L41Bae J., Zhu Z., Hartmann L., 2017, ApJ, 850, 201Bae J., Teague R., Zhu Z., 2021, arXiv e-prints, p. arXiv:2102.03899Bai X.-N., 2016, ApJ, 821, 80Bai X.-N., Stone J. M., 2014, ApJ, 796, 31Brittain S. D., Najita J. R., Carr J. S., 2019, ApJ, 883, 37Carrasco-González C., et al., 2019, ApJ, 883, 71Christiaens V., et al., 2019, MNRAS, 486, 5819Cugno G., et al., 2019, A&A, 622, A156Dong R., Rafikov R. R., Stone J. M., Petrovich C., 2011, ApJ, 741, 56Dong R., Li S., Chiang E., Li H., 2017, ApJ, 843, 127Dong R., Liu S.-Y., Fung J., 2019, ApJ, 870, 72Dullemond C. P., Isella A., Andrews S. M., Skobleva I., Dzyurkevich N.,2020, A&A, 633, A137Espaillat C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P.,Henning T., eds, Protostars and Planets VI. p. 497 ( arXiv:1402.7103 ),doi:10.2458/azu_uapress_9780816531240-ch022Flaherty K. M., et al., 2017, ApJ, 843, 150Flock M., Turner N. J., Nelson R. P., Lyra W., Manger N., Klahr H., 2020,ApJ, 897, 155Follette K. B., et al., 2017, AJ, 153, 264Fromang S., Lyra W., Masset F., 2011, A&A, 534, A107Fung J., Chiang E., 2016, ApJ, 832, 105Fung J., Shi J.-M., Chiang E., 2014, ApJ, 782, 88Fung J., Artymowicz P., Wu Y., 2015, ApJ, 811, 101Goldreich P., Tremaine S., 1979, ApJ, 233, 857Goodman J., Rafikov R. R., 2001, ApJ, 552, 793Gyeol Yun H., Kim W.-T., Bae J., Han C., 2019, ApJ, 884, 142Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., GirardJ. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749Hartmann L., 1998, Cambridge Astrophysics Series, 32Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998, ApJ, 495, 385Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742Hu X., Zhu Z., Okuzumi S., Bai X.-N., Wang L., Tomida K., Stone J. M.,2019, ApJ, 885, 36Hu X., Wang L., Okuzumi S., Zhu Z., 2020, arXiv e-prints, p.arXiv:2002.01583Isella A., et al., 2018, ApJ, 869, L49Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L.,2019, ApJ, 879, L25Jacquet E., 2013, A&A, 551, A75Johansen A., Youdin A., Klahr H., 2009, ApJ, 697, 1269Kanagawa K. D., Muto T., Tanaka H., Tanigawa T., Takeuchi T., TsukagoshiT., Momose M., 2015, ApJ, 806, L15Keppler M., et al., 2018, A&A, 617, A44Kley W., Lin D. N. C., 1992, ApJ, 397, 600Kretke K. A., Lin D. N. C., 2007, ApJ, 664, L55Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17Liu H. B., 2019, ApJ, 877, L22Lubow S. H., Seibert M., Artymowicz P., 1999, ApJ, 526, 1001Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603Martin R. G., Lubow S. H., 2011, MNRAS, 413, 1447Morbidelli A., Szulágyi J., Crida A., Lega E., Bitsch B., Tanigawa T., Kana-gawa K., 2014, Icarus, 232, 266Muto T., Suzuki T. K., Inutsuka S.-i., 2010, ApJ, 724, 448Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610Okuzumi S., Momose M., Sirono S.-i., Kobayashi H., Tanaka H., 2016, ApJ,821, 82Perez S., Dunhill A., Casassus S., Roman P., Szulágyi J., Flores C., MarinoS., Montesinos M., 2015, ApJ, 811, L5Pérez S., Casassus S., Benítez-Llambay P., 2018, MNRAS, 480, L12MNRAS , 1–16 (2020) I. Rabago and Z. Zhu
Philippov A. A., Rafikov R. R., 2017, ApJ, 837, 101Pinte C., et al., 2018, ApJ, 860, L13Pinte C., et al., 2019, Nature Astronomy, 3, 1109Pinte C., et al., 2020, ApJ, 890, L9Rameau J., et al., 2017, AJ, 153, 244Rosotti G. P., Teague R., Dullemond C., Booth R. A., Clarke C. J., 2020,MNRAS, 495, 173Rozyczka M., Bodenheimer P., Bell K. R., 1994, ApJ, 423, 736Stoll M. H. R., Kley W., Picogna G., 2017, A&A, 599, L6Stone J. M., Tomida K., White C. J., Felker K. G., 2020, arXiv e-prints, p.arXiv:2005.06651Suriano S. S., Li Z.-Y., Krasnopolsky R., Shang H., 2017, MNRAS, 468,3850Szulágyi J., Masset F., Lega E., Crida A., Morbidelli A., Guillot T., 2016,MNRAS, 460, 2853Takahashi S. Z., Inutsuka S.-i., 2014, ApJ, 794, 55Takeuchi T., Lin D. N. C., 2002, ApJ, 581, 1344Tanigawa T., Ohtsuki K., Machida M. N., 2012, ApJ, 747, 47Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018,ApJ, 860, L12Teague R., Bae J., Bergin E. A., 2019, Nature, 574, 378Tominaga R. T., Takahashi S. Z., Inutsuka S.-i., 2020, arXiv e-prints, p.arXiv:2008.02564Urpin V. A., 1984, Soviet Ast., 28, 50Wagner K., et al., 2018, ApJ, 863, L8Wang J. J., et al., 2020, AJ, 159, 263Zhang K., Blake G. A., Bergin E. A., 2015, ApJ, 806, L7Zhang S., et al., 2018, ApJ, 869, L47Zhu Z., Stone J. M., 2018, ApJ, 857, 34Zhu Z., Stone J. M., Rafikov R. R., 2013, ApJ, 768, 143Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88Zhu Z., et al., 2019, ApJ, 877, L18This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000