Observational Signatures of Planets in Protoplanetary Disks: Planet-Induced Line Broadening in Gaps

Protoplanets can produce structures in protoplanetary disks via gravitational disk-planet interactions. Once detected, such structures serve as signposts of planet formation. Here we investigate the kinematic signatures in disks produced by multi-Jupiter mass ($M_{\rm J}$) planets using 3D hydrodynamics and radiative transfer simulations. Such a planet opens a deep gap, and drives transonic vertical motions inside. Such motions include both a bulk motion of the entire half-disk column, and turbulence on scales comparable to and smaller than the scale height. They significantly broaden molecular lines from the gap, producing double-peaked line profiles at certain locations, and a kinematic velocity dispersion comparable to thermal after azimuthal averaging. The same planet does not drive fast vertical motions outside the gap, except at the inner spiral arms and the disk surface. Searching for line broadening induced by multi-$M_{\rm J}$ planets inside gaps requires an angular resolution comparable to the gap width, an assessment of the gap gas temperature to within a factor of 2, and a high sensitivity needed to detect line emission from the gap.


INTRODUCTION
Detecting protoplanets forming in protoplanetary disks provides fundamental constraints on key parameters governing the planet formation processes, such as location, timescale, and local environment. Apart from directly detecting the photosphere, protoplanets can be detected through emissions from their circumplanetary disks (CPD; e.g., Zhu 2015;Eisner 2015;Szulágyi et al. 2017), as well as planet-induced kinematics (e.g., Perez et al. 2015Perez et al. , 2018 and chemical signatures in the circumstellar disk (Cleeves et al. 2015). Massive efforts have been carried out to look for protoplanets via direct imaging; yet, most yielded null results (e.g., Testi et al. 2015;Canovas et al. 2017;Maire et al. 2017;Uyama et al. 2017). In a few cases, candidate point sources have been detected (e.g., Quanz et al. 2013;Currie et al. 2014;Sallum et al. 2015;Reggiani et al. 2018;Guidi et al. 2018;Keppler et al. 2018), while in most of them the nature of the detection is under debate (e.g., bona fide planets vs disk structures; Thalmann et al. 2016;Follette et al. 2017;Hord et al. 2017;Mendigutía et al. 2018). * NASA Sagan Fellow In this paper, we investigate a large scale kinematic signature that a forming Jovian planet produces in a disk -meridional gas motions in planet-opened gaps (Morbidelli et al. 2014;Fung & Chiang 2016). The mechanism of gap opening can be described as the balance between the planetary torque that repels material from its co-orbital region and the viscous torque that attempts to fill it back in (Lin & Papaloizou 1993). Previous work have shown that under typical conditions, a 1 Jupiter mass (M J ) planet can open a deep gap a few orders of magnitude depleted in gas (e.g., Fung et al. 2014;Fung & Chiang 2016). Under hydrostatic equilibrium, the gas density distribution ρ(r, z) in the vertical direction z at a stellocentric radius r is where ρ 0 (r) is the midplane density at r , and h(r) is the disk scale height at r. Close to the planet's orbit r p , h is nearly constant. The planet's torque is strongest near the disk midplane (z = 0), and diminishes with increasing |z|. The viscous torque is on the other hand more evenly distributed. Even though these two torques are overall balanced, they are unable to cancel each other locally. This leads to a planetary torque dominated midplane where gas flows away from the planet, and a vis-cous torque dominated upper atmosphere that flows toward the planet. The two flows connect via meridional circulations, which can reach velocities comparable to local sound speed c s in gaps opened by super-Jupiter planets (Fung & Chiang 2016). These gas motions may significantly broaden molecular emission lines. This is the key signature that we will examine. The main motion of gas in a disk is Keplerian rotation with a velocity v K . Significant deviations from Keplerian rotation may indicate the presence of substantial density structures (e.g., a warp, Rosenfeld et al. 2014;Casassus et al. 2015;Facchini et al. 2018;Mayama et al. 2018). The thermal motion of gas molecules has a typical velocity of v th ∼ c s ∼ (h/r)v K , usually several percent of v K . Instabilities, such as the magnetorotational instability (MRI), may induce sub-sonic turbulent motions (e.g., Simon et al. 2011).
Planet-induced kinematic signatures have been investigated before. In a pioneering work, Pérez et al. (2015) assessed the detectability of the rotation of CPDs. More recently, Pinte et al. (2018) showed that a planet may introduce local kinks in the isovelocity maps of the circumstellar disk, Teague et al. (2018a) quantified the radial modulation in the circumstellar disk rotation due to gap structures, and Huang et al. (2018b) studied the detectability of the anticyclonic gas motion inside planettriggered vortices. Perez et al. (2018) investigated gas kinematics caused by pressure gradients at gaps, spiral wakes and vortices. While these works mainly focused on planet-induced signatures in intensity weighted velocity (moment-1) maps 1 , we study how meridional motions are manifested in velocity dispersion (moment-2) maps.
Today, Keplerian rotation can be probed precisely (e.g., used to infer the mass of the central star(s); Czekala et al. 2015Czekala et al. , 2016Czekala et al. , 2017. Searching for subsonic turbulent motions such as MRI turbulence (e.g., Bai 2015;Simon et al. 2015Simon et al. , 2017 is an ongoing effort. Pioneering efforts have been carried out using various molecule tracers, such as Hughes et al. (2011, CO, TW Hya and HD 163296), Guilloteau et al. (2012, CS, DM Tau), Flaherty et al. (2015  Our goal -to study the observational signatures of meridional motions in planet-opened gaps -is similar in nature to those in looking for turbulence generated by instabilities. We combine 3D hydrodynamics ( §2.1) and radiative transfer simulations ( §2.2) to produce synthetic line emission observations. We illustrate the key observational signature inside the gap ( §3), while commenting on the variations in the key signature under various circumstances ( §4.1) as well as gas motions outside the gap ( §4.2). A generic strategy in searching for such signatures in real observations is outlined ( §5), and advantages of those signatures as signposts of planets are highlighted ( §5.1).

SIMULATIONS
We adopt a general scheme to produce synthetic line observations similar to Perez et al. (2018). First, we carry out 3D global hydrodynamics simulations using the PEnGUIn code  to calculate the density and velocity structures of disks perturbed by a giant planet. Then, we combine dust radiative transfer simulations using HOCHUNK3D (Whitney et al. 2013) and molecular line radiative transfer simulations using the package "Simulation Package for Astronomical Radiative Xfer" (SPARX 2 ) to produce channel maps.

Hydrodynamics Simulations
The setup we use for our hydrodynamics simulations of gas-opening largely borrows from Fung & Chiang (2016). We perform simulations in spherical coordinates, and the grid has dimensions of 320 (radial) by 810 (azimuthal) by 50 (polar) cells. Our domain extends from 0.3r p to 2r p radially in logarithmic spacing, covers the full 2π in azimuth, and includes 0.25 radians from the midplane to the top of the grid in polar angles. We simulate only one half of the disk, taking advantage of its symmetry around the midplane. This grid setup is nearly identical to the one used by Fung & Chiang (2016), except both the radial and polar domain have been extended to improve our ability to isolate the gap and resolve the disk surface. The resolution we choose have been shown to give good convergence on the gap profile (Fung et al. 2014;Fung & Chiang 2016). Since the gap profile is determined by the balance of planetary torque and disk diffusion, this implies our resolution is also sufficient to constrain the amount of turbulent diffusion in the system.
For boundary conditions, in the radial direction we have fixed boundaries with wave killing zones; in the azimuthal direction it is periodic; and in the polar di-rection it is reflective, both at the midplane to enforce symmetry, and at the top to prevent unrealistic mass influx or outflux. Also, we use an updated version of PEnGUIn (Fung & Lee 2018) that includes the fast orbital advection algorithm (Masset 2000), which helps speed up the code and reduce numerical diffusion.
We carry out two simulations each with one planet on a fixed circular orbit at r p = 30 AU for 1000 planetary orbits. The planets each has a mass of M p = 10 −3 M (Model Q1) and 4 × 10 −3 M (Model Q4); i.e., 1 M J and 4 M J if M = 1M . The disk has an initial surface density profile of: and a time-invariant vertically isothermal temperature profile that gives a sound speed of: This choice of using a vertically isothermal disk is inspired by the fact that at tens of AU the disk is expected to be vertically isothermal for bulk of its mass. A more realistic equation of state (EOS) may introduce an order-of-unit correction to our models, but unlikely to alter our results qualitatively. Our simulations do not account for disk self-gravity, thus we are free to tune Σ 0 to set the total disk mass M disk = 0.001M in both models. We choose this M disk so that the disk is globally optically thin to 13 C 18 O emission ( §3). We include a small viscosity using the Shakura & Sunyaev (1973) prescription, ν = αhc s , where ν is the kinematic viscosity, and α is set to 10 −3 . This amount of viscosity ensures that the disk reaches steady state within 1000 planetary orbits. Note that physically this viscosity corresponds to turbulent diffusion, which would manifest itself in line broadening. This "viscous turbulence" is not included in our line emission simulations, as we focus on planet-induced effects only (nonplanet-induced turbulence in disks is likely weak, see §1). In our plots, the disk rotates counterclockwise.
Naturally, simulation results greatly resemble those by Fung & Chiang (2016) (Figure 1(a)(b)(c)). In Appendix A, Figure A.1 can be compared to Figure 4 of Fung & Chiang (2016), and Figure A.2 illustrates the distribution of kinetic energy due to planet-induced vertical motions in the disk. While we are originally motivated by the meridional circulation in and around the gaps, we find that the gap in model Q4 is turbulent, and this turbulent motion can contribute to line broadening even more greatly than the circulation (see §4 for more discussions).

Radiative Transfer Simulations
We feed the exact 3D hydro density structures into HOCHUNK3D to calculate the disk temperature, following the procedures described in Dong et al. (2016a). The central star, which is the only heating source in our models, is assumed to be a 1 M protostar with a temperature of 4350 K and a radius of 2.3 R . A population of interstellar medium (ISM) dust (Kim et al. 1994, sizes up to ∼ 1µm) is employed to process starlight; they are assumed to be well-mixed with the gas with their density set to 0.1% of the gas density. Both dust scattering (anisotropic; scattering phase function approximated using the Henyey-Greenstein function) and absorption/emission are included (dust properties can be found in Figure 2 in Dong et al. 2012).
Figure 1(d) shows azimuthally averaged midplane temperature (T midplane ) for both models, as well as in a control model with no planet (i.e., a full disk). The gap region is on average about 5 K, or 15%, hotter than in the no planet case. Panel (e) shows azimuthally averaged vertical temperature profiles at r p . The midplane is the coldest, and T rises with z before plateauing at the surface. A nearly isothermal midplane layer with T − T midplane < 10 K encompasses ∼ 95% of the mass.
We note that temperature, and therefore the vertical hydrostatic structure, in our models are not self-consistent. Hydrodynamics simulations assume a parametrized scale height profile (Equation 3), which under the hydrostatic equilibrium condition corresponds to a midplane temperature not fully consistent with the one obtained from radiative transfer (RT) simulations (the former is ∼30% lower at r p in the no planet case; note in self-consistent disk models the difference may be higher, Isella & Turner 2018). Achieving selfconsistency requires either coupled hydro-RT simulations, or iterations, which are beyond our scope. This self-inconsistency however does not affect the key observational signature to be defined in the next section (see more discussions in §4.1).
Synthetic molecular line channel maps and moment maps are produced with SPARX, which calculates the spectral line radiative transfer of specific transitions for molecules of interests. Molecular gas densities and velocities extracted from the hydrodynamics simulations and disk temperature derived from the HOCHUNK3D simulations are provided to SPARX. For demonstration purposes, we consider the rotational transitions in J = 3−2 of CO and its isotopologues C 18 O, and 13 C 18 O. Fractional abundances of CO, C 18 O, and 13 C 18 O relative to H 2 are set at 10 −4 , 2 × 10 −7 , and 2 × 10 −9 , respectively. These transitions are mainly for covering different levels of opacity, and the listed abundances roughly follow the Surface Density and Temperature Structures in Our Models Figure 1. Top: Surface density maps of the two models at 1000 planetary orbits ((a) and (b)); the green plus and dot mark the locations of the star and the planet, respectively. Bottom: Azimuthally averaged radial profiles of the surface density (c) and midplane temperature (d), and azimuthally averaged vertical temperature profiles at rp ((e); θ is the polar angle). See §2 for discussions. measurements in ISM dense clouds (e.g., van Dishoeck & Blake 1998;Wilson & Rood 1994).
While SPARX is capable of evaluating iteratively the non-local-thermodynamical-equilibrium (non-LTE) molecular level populations and the corresponding radiation intensity, the high gas density throughout the disk domain under consideration is sufficient to thermalize the J = 3 − 2 transitions of CO and its isotopolouges, whose critical densities are around 10 4 − 10 5 cm −3 . Thus, we applied the LTE level population following the Boltzmann distribution and performed standard ray-tracing for generating the synthetic images.
A pixel dimension of 300×300 and a velocity resolution of 0.03 km s −1 are used for ensuring the molecular emission features are both spatially and spectrally resolved. Dust continuum emission is included in our line radiative transfer calculations; it is insignificant due to the low ∼mm-wavelength opacity of the assumed ISM dust (on the order of 0.1 cm 2 /g). Millimeter opacity in real disks is expected to be 1-2 orders of magnitude larger; our choice helps isolate and reveal the planet-induced signatures in the cleanest way. Simulated channel maps are not corrupted to incorporate realistic noises and angular resolutions. We defer direct model-observation comparisons to future work.
Finally, we note that moment-2 of line emission from optically thin gas at temperature T with no non-thermal motions equals to the isothermal sound speed c s , where m µ is the molecular weight.

THE KEY OBSERVATIONAL SIGNATURE OF PLANET-INDUCED VELOCITY DISPERSION
In this section, we define and study the key observational signature of planet-induced turbulence inside the gap in Model Q4. We employ a tracer, 13 C 18 O (J = 3 − 2), to which the disk is globally optically thin (optical depth τ ∼ 0.1 at v LOS = 0 outside the gap and τ peaks at ∼ 10 −7 inside the gap). We adopt a face-on geometry, so that the vertical direction in the disk corresponds to the line of sight (LOS), and vertical velocity dispersion corresponds to moment-2. The dependence of the key signature on parameters will be discussed in §4.1. We emphasize that experiments in this section (and most of the paper) are designed to demonstrate the basic principle, rather than to mimic real life situations. In particular, we choose 13 C 18 O, a species almost impossible to detect inside a deep gap at the moment, because it gives the most simple picture of the planetary signature. This simplicity comes from the fact that the thermal broadening would be a smooth and almost monotonic function with radius if, and only if, the emission originates from the midplane everywhere. This makes the planet-induced broadening stands out more clearly. As we will see ( §4.1; Figure 8), planet-induced line broadening inside gaps is clearly present (and retrievable; §5) in emission from tenuous gas tracers such as CO, as long as the gap is optically thin to the emission, which is often easy to satisfy due to the large gap depletion. Figure 2 shows the moment-2 map. A ring of high velocity dispersion around the planet's orbit is prominent. This is the key signature of planet-induced turbulence. Figure 3 compares moment-2 in the no planet case (a) with that in Model Q4 (c). In the former, the velocity dispersion is entirely thermal -Eqn. 4 with T = T midplane (dotted line in the bottom panel) closely traces the measured moment-2 (dashed line). The small difference can be attributed to the higher temperatures of the tenuous gas at high latitudes. In the latter, the velocity dispersion outside the gap (roughly 0.8 − 1.5r p ) is consistent with thermal broadening too, suggesting no substantial planet-induced motions there.
Inside the gap, however, there is a ring of high velocity dispersion. After azimuthal averaging, at its peak at ∼0.9r p the total line broadening is about twice of the thermal broadening, indicating a "mean" non-thermal velocity of ∼1.7c s (total ≈ √ thermal 2 + turbulent 2 for Gaussian-like turbulence). To highlight the "planetinduced" nature of this additional broadening, panel (b) shows the moment-2 map simulated with the vertical velocity v z manually set to 0, so that only thermal broadening contributes. In this case, moment-2 inside the gap is only ∼10% higher than that in the no planet case due to the slightly higher temperature in the former (Figure 1). Inside the planet's Hill radius (marked in Figure 2; resolved by ∼15 cells in hydro simulations) we observe a velocity dispersion exceeding 30%v K , almost entirely non-thermal (note that the heating from the CPD accretion is not taken into account). Figure 4 further illustrates the origin of the planetinduced line broadening. We plot the line profiles at three locations inside the gap. The top row shows the moment-0 (a) and moment-2 (c) maps, as well as the density-weighted vertically averaged v z in the upper (closer to observers) half disk ((b); the full disk has zero average v z everywhere due to its symmetry).
Locations A, B, and C represent three situations with different flow patterns. The line profile of optically thin gas at temperature T follows a 1D Maxwell-Boltzmann distribution For gas at T = T midplane ≈ 40K inside the gap, the fullwidth-half-magnitude (FWHM) of a thermally broadened line is FWHM thermal = 2 2 ln 2 kT midplane /m µ = 0.24 km/s. At Location A, the FWHM is 0.25 km/s in both the half-disk and full disk models ((d), suggesting an absence of both velocity dispersion within the half-disk (intra half-disk) and bulk velocity of the entire half-disk column (also note the zero average v z in the half-disk).
The vertical motion at Location C is completely different (f ). First, the half-disk model line profile has a FWHM of 0.30 km/s, substantially larger than FWHM thermal = 0.24 km/s, indicating the presence of intra-half-disk velocity dispersions. Secondly, the fulldisk line profile is double peaked at ±0.3 km/s, with each peak tracing the bulk motion in one half of the disk (the half-disk profile overlaps with the left peak in the full-disk profile). This indicates that the bulk gas in each half disk is moving at a supersonic bulk velocity of 0.30 km/s ∼ 3c s . Both factors contribute to the resulting full-disk moment-2, 0.33 km/s, significantly larger than that at position A, 0.11 km/s. The sign of the peak v LOS in the half-disk profile indicates a bulk motion towards the midplane in the two half disks.
The flow pattern at Location B is in between. The half disk line profile is kinematically broadened to FWHM = 0.32 km/s, comparable to Location A, but peaks at a lower |v LOS | ∼ 0.1 km/s, indicating a smaller global (entire half-disk) bulk motion. As a result, the full disk line profile remains single peaked, with its moment-2 (0.18 km/s) in between Locations A and C. Note that at Location B the bulk gas motion is away from the midplane, indicated by the negative peak v LOS in the half-disk model; this is different from Location C.

Dependences of the Key Signature on Various Parameters
Temporal Fluctuation - Figure 5 shows moment-2 maps in Model Q4 at 800, 900, and 1000 planetary orbits. Vertical velocity dispersion at different regions inside the gap fluctuates in a seemingly uncorrelated manner. The key signature -the ring of high velocity dispersionpersists, while fluctuating with an azimuthally averaged peak-to-peak amplitude of ∼30% (bottom panel). Note that the change is not monotonic -the bump is lowest at 900 orbits. This temporal fluctuation accords with the interpretation that we find planet-induced turbulence.
The key signature drops dramatically, indicating a transition from turbulent to laminar flow when planet mass decreases from Q4 to Q1. The non-thermal velocity dispersion, once azimuthally averaged, is at the level of ∼10% inside the gap. Narrow streamers of high velocity dispersion, mainly associated with the inner gap edge and spiral arms, are visible (c). However, extremely high angular resolutions are likely needed to resolve them in real observations.
Temperature Structure -As mentioned in §2.2, density and temperature in our models are not self-consistent.
Here we carry out an experiment to qualitatively assess its effect on the key signature. Figure 7 compares two moment-2 maps of Model Q4, one with radiative transfer temperature adopted in the line emission simulation (a), and the other with the temperature corresponding to Equation 3 (hydro input temperature; (b)). Different temperatures result in different velocity dispersions in the thermally broadened inner and outer disks. Inside the gap, however, the two converge as kinematic broadening dominates, and the key signature remains robust.
Tenuous Gas Tracers - Figure 8 compares moment-2 in Model Q4 in CO (a) and in a control case produced with v z manually set to 0 (b). While the gap is still optically thin to CO emission (τ peaks at ∼0.01 at v LOS = 0), the inner and outer disks are extremely optically thick (τ 1000). Therefore outside the gap CO traces a hot surface layer above the midplane, leading to larger thermal broadening than in the midplane-tracing 13 C 18 O emission. Nevertheless, The regions with the highest velocity dispersion inside the gap are still visible in the 2D map (a), and the key planet-induced signature -the difference between the thin solid and dashed curves in (c) -is still prominent. Furthermore, the high velocity dis-Model Q4, 13 C 18 O, Face On (reference: c s (40K) = 0.10 km/s)   3(c)), suggesting an absence of local (intra-half-disk) velocity dispersion as well.
On the other hand, at the two inner spiral arms the half-disk line profile peaks at |v LOS | ∼ 0.05 km/s ∼ 0.5c s , indicating the presence of transonic global vertical motions with respect to the midplane. Such motions have previously been seen (e.g., Lyra et al. 2016;Riols & Latter 2018). Note that shock heating at the arms (e.g., Richert et al. 2015;Hord et al. 2017), expected to further modify the flow pattern, is not included. The line width in the half-disk model at both arms are consistent with thermal broadening, indicating an absence of intra-half-disk velocity dispersion. Interestingly, gas move away from the midplane at the primary arm, and towards the midplane at the secondary arm. The kinematic support at the arms helps lift the disk surface in addition to the usual vertical pressure support (Dong & Fung 2017a, Fig . 3), echoing with previous findings from 3D hydro+RT simulations that the inner arms are more prominent than the outer arms in scattered light (e.g., Dong et al. 2015b;Zhu et al. 2015).

COMMENTS ON MODEL APPLICATIONS
Searching for planet-induced turbulence inside gaps is similar in nature to searching for other types of turbulence induced by instabilities, such as MRI ( §1). Naturally, the techniques developed for the latter applies to the former. Here we briefly outline a generic roadmap for real life applications of our models. azimuthally averaged radial profiles (c). The inner and outer disks are extremely optically thick to CO, which traces warm layers above the midplane; consequently, thermal broadening is larger than in the 13 C 18 O case. The gap region remains optically thin. The narrow arcs and rings with high velocity dispersions in the inner disk in (a) are absent in (b). They are induced by the planet in the surface layer. The key planet-induced signature in velocity dispersion inside the gap, which is the difference between the thin solid and dashed curves in (c), stays robust. Note that the thick black curve in the radial profile plot in Figures 3, 5, 7, and here is the same. See §4.1 for discussions.
Model Q4, Moment-2 (reference: c s (40K) = 0.10 km/s) 13 C 18 O CO  When probing a planet-opened gap using a molecular line to which the gap region is optically thin, the line is broadened by four factors: (A) Keplerian broadening, (B) thermal broadening, (C) planet-induced broadening, and (D) other turbulent broadening. To infer the presence and measure the strength of (C), we need to estimate the contributions from (A), (B), and (D), and take them out of the observed line broadening.
Examples of assessing (A) and (B) using line observations can be found in, e.g., Flaherty et al. (2015Flaherty et al. ( , 2017Flaherty et al. ( , 2018 and Teague et al. (2016Teague et al. ( , 2018b. To estimate thermal broadening, an accurate assessment of the temperature at the emitting regions is needed. Techniques have been developed recently to do so using observations of multiple transitions from a single molecule, and temperatures can now be constrained to within a few K in certain cases (e.g., TW Hya, Teague et al. 2018b). Once the contributions from (A) and (B) are taken out, additional broadening, if present, can be attributed to (C) and (D). Separate efforts are then needed to distinguish (C) from (D). Note that so far the searches for (D) have generally returned null results, putting upper limits on v turb,(D) at a few percent of c s ( §1).
In order to minimize the contributions from (A) and (B), disks close to face-on are preferred, so do heavy molecule as they naturally exhibit lower thermal broadening at a given temperature (turbulent broadening is not expected to depend on molecular weight).
We emphasize again that the experiments here are designed to best illustrate the key planet-induced signature, instead of representing best practices in real life. The detectability of gas emission from a real gap mainly depends on its gas surface density (measurable; e.g., van der Marel et al. 2015van der Marel et al. , 2016, and the choice of the tracer.
In the companion scenario, the properties of featuredriving companions, such as location, orbit, and mass, may be constrained using observed feature morphology (e.g., Fung & Dong 2015;Dong & Fung 2017b,a). This provides an extremely valuable channel in using disk observations to test planet formation models. Unfortunately, planetary and non-planetary mechanisms sometimes produce similar observational features (e.g., Dong et al. 2018a). Except in rare cases (e.g., the stellar mass companion in HD 100453, Dong et al. 2016b;Wagner et al. 2018), it is difficult to confirm the planetary origin of observed features, and subsequently to take advantages of disk observations to study planet formation.
Planet-induced kinematics signatures provide another possible avenue to help distinguish planetary and nonplanetary mechanisms ( §1). The detection of high velocity dispersions inside gaps, as studied here, can potentially lend strong support to the planet scenario. In addition, the strength in these non-thermal dispersions may be used to constrain the planet mass.
Quite a few types of planet-induced kinematics signatures have been explored ( §1). Comparing with them, the signature studied here has four characteristics. First, it is a global signature on scales comparable to the gap width, typically 25%-50%r p , easing the requirement on angular resolution. Further, its signal-to-noise ratio can be boosted by azimuthal averaging. Certain types of planet-induced kinematics, such as the CPD rotation , are local features detectable only under high angular resolutions. Figure 11 shows the moment-1 maps of Model Q4 in both 13 C 18 O (a) and CO (b) at i = 45 • . The CPD rotation, clearly visible, is confined to a region comparable to the planet's Hill sphere, which is ∼10 times smaller than r p for Jovian planets.
Secondly, velocity dispersion in the gap driven by the planet is arguably a more numerically robust feature than other planet-induced kinematic signatures such as the CPD and local kinks in the isovelocity maps. The moment-2 signals reported here have length scales larger than local scales such as the planet's Hill radius or disk scale height, making it easier to be resolved and captured in hydrodynamics simulations.
Thirdly, the requirement on the accuracy of gap temperature assessment is lenient. Turbulent motions can be mildly supersonic, and the total line broadening can double that of thermal. Even in disks with modest inclinations and substantial Keplerian broadening, the planet-induced signature is still easily visible (Figure 9). Knowing the local temperature to within a factor of 2 Model Q4, r p = 30 AU, Moment-1, i = 45 • is likely sufficient to reveal velocity dispersions induced by a multi-M J planet. This is in contrast with certain types of planet-induced kinematics. For example, Teague et al. (2018a) showed that the variations in the gas rotations due to the modulations in the radial pressure gradient induced by Jovian planets (the weak radial structures at close to r p in Figure 11) may only be at the level of 1% v K , or perhaps 10% c s . Measuring such signals demand accurate assessments of disk temperatures to a few percent.
Finally, the bottleneck in real-life application is expected to be the high sensitivity required to probe the weak line emission from deeply depleted gaps. With the current Atacama Large Millimeter/submillimeter Array (ALMA), this poses a strong challenge, if possible at all, for gaps at a few tens of AU, while the situation improves dramatically for giant planets at large radii. For example, if we observe Model Q4 at 140 pc with an angular resolution of 0 .1 (0.5r p ) to just resolve the gap, we need up to a few weeks of ALMA integration time to achieve a 10σ detection of a CO emission line with a spectral resolution of 0.1 km/s inside a τ CO = 0.1 gap (i.e., the scenario in Figure 8). The situation will be much more favorable if the gap is at a larger radiusshould we place the planet at 100 AU with everything else being equal, an angular resolution of 0 .3 would be sufficient to resolve the gap, and the integration time for the same 10σ detection would be reduced to a few hours.

SUMMARY
In this paper we use 3D hydrodynamics and radiative transfer simulations to study the signatures of planetinduced turbulence in gaps in gas line observations. A planet with M p = 4 × 10 −3 M (4 M J around a 1 M star) can induce transonic to mildly supersonic turbulent motions inside its gap, on both local (intra-halfdisk) and global (entire half-disk) scales (Figure 4). The resulting velocity dispersion, after azimuthal averaging, can double that of thermal ( Figure 3). The turbulent motion can result in double-peaked line profiles at certain regions in the gap (Figure 4). This basic key signature is robust against temporal variations ( Figure 5), different treatments of disk temperature (Figure 7), finite inclinations (Figure 9), and is visible in tracers with a variety of opacity as long as the gap is optically thin (Figure 8). For a planet with M p = 10 −3 M , azimuthally averaged turbulent velocity dispersion drops to ∼10%c s ( Figure 6). In the dense outer disk, the planet does not introduce significant vertical motions; this is true in the inner disk too, except at the surface ( Figure 8) and in the spiral arms ( Figure 10; a bulk motion away from the midplane in the primary arm, and opposite in the secondary arm).
Detecting planet-induced velocity dispersion can lend support to the planetary origin of gaps, and its strength may potentially be used for dynamical planet mass assessments. Comparing with other types of planetinduced kinematic signatures, velocity dispersion in gaps has its own characteristics. It is on radial scales comparable to the gap width, and visible after azimuthal averaging. Searching for this signature likely requires an assessment of the gap gas temperature to within a factor of 2, meanwhile a high sensitivity is needed to detect the weak line emission from deep gaps.

A. FLOW PATTERNS IN THE HYDRO MODELS
As discussed in Section 2.1, the two simulations used in this work, Models Q1 and Q4, are essentially replicas of the simulations done by Fung & Chiang (2016) with slightly modified domain size. Comparing their Figure 4 with our Figure A.1, we find the overall meridional flow pattern remains the same. Gas is repelled from the planet's gap region near the midplane, rises near the gap edge to about 2 or 3 scale heights, circulates back into the gap, falls and returns to the midplane to repeat the cycle again.
Despite the similarities, one noticeable difference is that meridional motion in the inner and outer disks are slower in our simulations by a factor of a few. The difference is likely due to both the more extended boundaries reducing the amount of wave reflection, and the more stable advection scheme used in the updated version of PEnGUIn. While this is an improvement, since this work focuses on gas motion within the gap, it has little impact on our results.
where c s is taken from Equation 3. The vertically integrated η is analogous to moment-2 detected face-on using an dense gas tracer; therefore η can be thought of as the vertically differentiated moment-2. Then Figure A.2 essentially illustrates which portion of the gap contributes more to the kinematic fraction of moment-2. Evidently, most of the contribution comes from near the midplane. Comparing it with Figure A.1, we see that the locations with the highest speeds do not necessarily contribute most, as those regions also tend to have low densities.