Observational Signatures of Tightly Wound Spirals Driven by Buoyancy Resonances in Protoplanetary Disks
DDraft version February 9, 2021
Typeset using L A TEX twocolumn style in AASTeX62
OBSERVATIONAL SIGNATURES OF TIGHTLY WOUND SPIRALS DRIVEN BY BUOYANCY RESONANCESIN PROTOPLANETARY DISKS
Jaehan Bae, ∗ Richard Teague, and Zhaohuan Zhu Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154, USA
ABSTRACTBesides the spirals induced by the Lindblad resonances, planets can generate a family of tightlywound spirals through buoyancy resonances. The excitation of buoyancy resonances depends on thethermal relaxation timescale of the gas. By computing timescales of various processes associated withthermal relaxation, namely, radiation, diffusion, and gas-dust collision, we show that the thermalrelaxation in protoplanetary disks’ surface layers (
Z/R (cid:38) . ) and outer disks ( R (cid:38) au) is limitedby infrequent gas-dust collisions. The use of isothermal equation of state or rapid cooling, commonin protoplanetary disk simulations, is therefore not justified. Using three-dimensional hydrodynamicsimulations, we show that the collision-limited slow thermal relaxation provides favorable conditions forbuoyancy resonances to develop. Buoyancy resonances produce predominantly vertical motions, whosemagnitude at the CO emission surface is of order of
100 m s − for Jovian-mass planets, sufficientlylarge to detect using molecular line observations with ALMA. We generate synthetic observations anddescribe characteristic features of buoyancy resonances in Keplerian-subtracted moment maps andvelocity channel maps. Based on the morphology and magnitude of the perturbation, we propose thatthe tightly wound spirals observed in TW Hya could be driven by a (sub-)Jovian-mass planet at 90 au.We discuss how non-Keplerian motions driven by buoyancy resonances can be distinguished from thosedriven by other origins. We argue that observations of multiple lines tracing different heights, withsufficiently high spatial/spectral resolution and sensitivity to separate the emission arising from thenear and far sides of the disk, will help constrain the origin of non-Keplerian motions. Keywords:
Protoplanetary disks (1300), Spiral arms (1559), Hydrodynamical simulations (767), Sub-millimeter astronomy (1647) INTRODUCTIONRecent observations have revealed a plethora of sub-structures in protoplanetary disks (e.g., ALMA Partner-ship et al. 2015; Andrews et al. 2018; Avenhaus et al.2018). Spiral arms, along with concentric rings andgaps, are found to be one of the most common typesof substructures. They are often interpreted as an out-come of the gravitational interaction between the diskand embedded planets therein (e.g., Bae et al. 2016c),although we cannot rule out other possibilities, such asgravitational instabilities (e.g., Meru et al. 2017; Hallet al. 2018), stellar flybys (e.g., Cuello et al. 2019, 2020),
Corresponding author: Jaehan [email protected] ∗ NASA Hubble Fellowship Program Sagan Fellow or infalling materials (e.g., Lesur et al. 2015) until we di-rectly detect companions.So far, most of the spiral arms are detected inoptical/near-infrared scattered light or (sub-)millimetercontinuum observations observations (e.g., Hashimotoet al. 2011; Muto et al. 2012; Grady et al. 2013; Benistyet al. 2015; Pérez et al. 2016; Benisty et al. 2017; Krauset al. 2017; Andrews et al. 2018; Canovas et al. 2018;Huang et al. 2018a; Reggiani et al. 2018; Uyama et al.2018; Gratton et al. 2019; Monnier et al. 2019; Keppleret al. 2020; Muro-Arena et al. 2020). Spirals in these ob-servations could be results of the increase in the densityof emitting materials (i.e., dust grains), the increase ofthe temperature at the shock front, and/or the increasein the height of the scattering surface.With the unprecedented high spatial/spectral res-olution and sensitivity the Atacama Large Millime- a r X i v : . [ a s t r o - ph . E P ] F e b Bae et al. ter/submillimeter Array (ALMA) offers, it is now possi-ble to use molecular line observations to probe kinemat-ics associated with spiral arms (e.g., Christiaens et al.2014; Tang et al. 2017; Teague et al. 2019; Huang et al.2020; Phuong et al. 2020), which can help better under-stand the origin of the spirals. It is also worth men-tioning so-called velocity kinks (Pinte et al. 2018, 2019,2020) and Doppler flips (Pérez et al. 2018b; Casassus& Pérez 2019; Pérez et al. 2020) seen in ALMA molec-ular line observations. While these are localized fea-tures rather than large-scale spirals, they are generallyinterpreted as the velocity perturbations associated withplanet-induced spirals.Using ALMA CO line observations, Teague et al.(2019) reported three spiral arms in the velocity andtemperature space in the TW Hya disk. One of theinteresting features about the spirals is that the pitchangle is very small, decreasing from 9 to 3 degrees be-tween 70 and 200 au. This tightly wound morphologydistinguishes themselves from a large number of spiralshaving pitch angles of 10 to 30 degrees in other proto-planetary disks (e.g., Huang et al. 2018a; Reggiani et al.2018; Uyama et al. 2018; Monnier et al. 2019; Yu et al.2019), bringing into question their origin. Because Lind-blad spirals’ pitch angle decreases as a function of thedistance from the planet (e.g. Zhu et al. 2015; Bae &Zhu 2018a), it is not completely impossible to explainsuch a tightly wound morphology with the traditionalview of Lindblad resonance-driven spirals (Goldreich &Tremaine 1979, 1980; Ogilvie & Lubow 2002; Bae & Zhu2018a,b). One may argue that a small pitch angle couldbe reconciled if the planet is located sufficiently far in-ward of the observed spirals. In this case, however, it isunclear why we do not observe spirals near the planetbut only far from it because we expect spirals generatestronger perturbations closer to the planet.1.1.
Theoretical Background: Buoyancy Resonances
Here we examine an alternative: buoyancy resonances(Zhu et al. 2012; Lubow & Zhu 2014; McNally et al.2020). Let us consider a gas parcel that is verticallydisplaced from the equilibrium position. Its responseto the perturbation can be described with the verticalbuoyancy frequency (a.k.a. Brunt-Bäisälä frequency) N Z = gγ ∂∂Z (cid:20) ln (cid:18) Pρ γ (cid:19)(cid:21) , (1)or N Z = gγ ∂∂Z (cid:20) ln (cid:18) Tρ γ − (cid:19)(cid:21) (2)adopting the ideal gas law. In the above equations, g is the gravitational acceleration, γ is the adiabatic in-dex, P is the gas pressure, ρ is the gas density, and T is the gas temperature. When the disk temperature isdominated by the stellar irradiation, the disk is hotternear the surface and colder near the midplane (Chiang& Goldreich 1997; D’Alessio et al. 1998; Dartois et al.2003; Rosenfeld et al. 2013). With a positive verticaltemperature gradient and γ ≥ , Equation (2) yields N Z > . We thus expect the gas parcel to vertically os-cillate around its equilibrium position with a frequency N Z .When the buoyancy frequency matches with the forc-ing frequency which, in this case, is planet’s orbital fre-quency, buoyancy resonances can develop (Zhu et al.2012; Lubow & Zhu 2014). Buoyancy resonances givea rise to the density and velocity perturbations along afamily of trailing spirals. As we will show below, one ofthe main characteristics of buoyancy spirals is that theyare very tightly wound compared with Lindblad spirals,in particular in the vicinity of the planet.It is worth pointing out that the thermodynamic prop-erties of the disk is important in the development ofbuoyancy resonances. In order for buoyancy resonancesto fully develop, the timescale for the gas to respondto thermal perturbations (hereafter relaxation timescale t relax ) has to be longer than the timescale associatedwith the buoyancy: t relax (cid:38) N − Z . When t relax (cid:28) N − Z ,the gas behaves isothermally and buoyancy resonancesare expected to be weak or absent.In this paper, we show planets can excite spiralsthrough buoyancy resonances, which can be detectableusing molecular line observations with ALMA. The pa-per is organized as follows. In Section 2, we describethree processes that determine the relaxation timescaleof the gas – radiation, diffusion, and gas-dust collision– and compute the corresponding timescales using aTW Hya disk model. We show that thermal relaxationin the surface layers ( Z/R (cid:38) . ) and outer regions( R (cid:38) au) of protoplanetary disks can be limited byinsufficient gas-dust collision. In Section 3, we presentthree-dimensional hydrodynamic simulations and showthat the collision-limited slow thermal relaxation pro-vides favorable conditions for buoyancy resonances. InSection 4, we generate synthetic observations and showthat the spirals driven by buoyancy resonances are ob-servable with molecular line observations using ALMA.Based on the tightly wound morphology and the magni-tude of velocity perturbations, we propose that the ve-locity spiral seen in TW Hya could be driven by a (sub-)Jovian-mass planet at 90 au. In Section 5, we show thatthe relaxation timescale is comparable to or longer thanthe dynamical timescale under a broad range of condi-tions, and discuss its implications for the development ofbuoyancy resonances, planet-induced gap profiles, and bservational Signatures of Buoyancy Resonances THERMAL RELAXATION OF THE DISK GAS2.1.
Disk Model
As one of the motivations of this work is to explainthe tightly wound spirals observed in the TW Hya disk(Teague et al. 2019), we adopt the disk density and tem-perature profiles similar to the one constrained for thedisk in Huang et al. (2018b).The gas surface density follows Σ g ( R ) = Σ p ( R/R p ) − p . (3)where Σ p is the surface density at R p = 90 au and p = 0 . . We choose Σ p = 3 . − such that thetotal disk gas mass within 210 au is 0.05 M (cid:12) , broadlyconsistent with observational estimates for the TW Hyadisk (Bergin et al. 2013). Throughout this paper we use R = r sin θ for the cylindrical radius and Z = r cos θ forthe height, where r and θ are the radius and the polarangle in the spherical coordinates, respectively.Following the prescription in Dartois et al. (2003), thegas temperature is parameterized as follows. T g ( Z ) = T atm + ( T mid − T atm ) cos (cid:16) π ZZ q (cid:17) if Z < Z q T atm if Z ≥ Z q (4)Here, the midplane and atmosphere temperatures are afunction of the cylindrical radius R following T atm ( R ) = T atm , p ( R/R p ) − q (5)and T mid ( R ) = T mid , p ( R/R p ) − q , (6)where T atm , p = 44 . K, T mid , p = 14 . K, and q = 0 . (Zhang et al. 2017; Huang et al. 2018b). In Equation(4), Z q ( R ) = 4 H mid ( R ) where H mid is the gas scaleheight determined with the disk midplane temperature.In terms of the disk aspect ratio, the above tempera-ture profile at the midplane corresponds to H mid /R =0 . × ( R/R p ) . , assuming . M (cid:12) for the stellarmass (Andrews et al. 2012; Huang et al. 2018b) and 2.4for the mean molecular weight of the gas.Using the above surface density and temperature pro-files, we construct the three-dimensional gas density dis-tribution that satisfies the vertical hydrostatic equilib-rium: − GM ∗ Z ( R + Z ) / − ρ g ∂P∂Z = 0 . (7) Solving the above equation results in the vertical densitydistribution that follows ρ g ( Z ) = ρ g , mid c s, mid c s ( Z ) exp (cid:34) − (cid:90) Z c s ( Z (cid:48) ) GM ∗ Z (cid:48) ( R + Z (cid:48) ) / d Z (cid:48) (cid:35) , (8)where c s, mid and c s ( Z ) denote the sound speed at thedisk midplane and at height Z . We then compute theangular velocity Ω that satisfies the radial force balance,taking into account the gas pressure gradient: Ω = Ω K sin θ + 1 ρ g r sin θ ∂P∂r . (9)Here, Ω K = (cid:112) GM ∗ /R is the Keplerian angular veloc-ity. The initial radial and meridional velocities are setto zero. 2.2. Thermal Relaxation Timescale
Protoplanetary disks are a mixture of gas and dust.Hydrogen molecules dominate the total mass and ther-mal energy, but they are inefficient at emitting radiation.Let us consider a situation where a gas parcel has beenperturbed from its equilibrium state to have a highertemperature. In order for the gas to lose its thermal en-ergy, hydrogen molecules first have to transfer their ki-netic energy to the surrounding dust grains. Dust grainsthen radiate away the excess thermal energy. This pro-cess can thus be understood as a sequential, two-stepprocess (see Figure 1 for a schematic diagram).When the collision between hydrogen molecules anddust grains is sufficiently frequent, the gas cools at therate dust grains radiate away their thermal energy; inthis case, gas and dust are thermally coupled. In theoptically thin regime, thermal photons emitted by dustgrains can freely escape from the disk. The relaxationtimescale of the gas can be described by the radiativetimescale of dust grains t rad = c V κ P σ SB T , (10)where c V denotes the specific heat capacity of the gas, κ P is the Planck mean opacity of the dust, and σ SB isthe Stefan-Boltzmann constant.In the optically thick regime, thermal photons emit-ted from dust grains are absorbed and emitted by othergrains multiple times before they eventually escape thedisk. In this case, the cooling timescale can be character-ized by the diffusion of photons. The diffusion timescaleassociated with the length scale λ diff can be written as t diff = λ (cid:104) D (cid:105) , (11) Bae et al. thinthick ⓵ t coll H H ⓵ t coll H H ⓶ t diff, ∥ ⓶ t diff, ⟂ ⓶ t rad Figure 1.
A schematic diagram showing the cooling processof the gas in protoplanetary disks. The cooling of the gas isa sequential, two-step process. Hydrogen molecules – themost abundant but poorly emitting species – first have tocollide with dust grains to lose their energy. Dust grainssubsequently radiate away the excess thermal energy. Thethermal photons from dust grains escape the disk freely whenthe disk is optically thin (the upper half of the diagram),while they do so after multiple re-absorption/emission (i.e.,diffusion) when the disk is optically thick (the lower half ofthe diagram). where D is the diffusion coefficient defined as D = 16 σ SB T c V κ R ρ (12)with κ R being the Rosseland mean opacity of the dustand the brackets (cid:104) (cid:105) denoting an average over λ diff (Ma-lygin et al. 2017). While the density, temperature, andopacity are expected to vary moderately along the radialand azimuthal directions, their variation can be largeralong the vertical direction due to the vertical strati-fication of the disk. We thus define in-plane diffusiontimescale t diff , (cid:107) and vertical diffusion timescale t diff , ⊥ ,for which the diffusion coefficient is averaged along theradial and vertical directions, respectively. The overalldiffusion timescale considering both in-plane and verti-cal diffusion can be estimated as t diff = (cid:18) t diff , (cid:107) + 1 t diff , ⊥ (cid:19) − . (13)For a given length scale λ diff , we find that t diff , (cid:107) and t diff , ⊥ are comparable when λ diff (cid:28) H g , not surpris-ingly, where H g is the gas scale height. However, one can be larger than the other by a factor of a few when λ diff (cid:38) H g . It is also worth mentioning that it is pos-sible that the length scale of the temperature perturba-tion along the radial/azimuthal direction and that alongthe vertical direction differ (see e.g., Miranda & Rafikov2020a).We note that, with the assumption that the diskhas a constant temperature vertically, along with τ R ≡ (cid:82) κ R ρ g d Z and Σ g (cid:39) ρ g H g , combining Equations (11)and (12) approximates to t diff , ⊥ (cid:39) c V Σ g σ SB T τ R , (14)a formula that is often adopted to account for opti-cally thick cooling in vertically-integrated, one-/two-dimensional framework.Let us now move on to the situation where the gas-dust collision is not sufficiently frequent. In this case,cooling of the gas is limited by the gas-dust collisionrate (assuming lack of other cooling mechanisms, suchas molecular/atomic line cooling; see below). Since theassumption of thermal equilibrium is no longer valid, wedefine the dust temperature T d . The gas-dust collisionaltimescale t coll can be written as t coll = ρ g c V Λ coll | T g − T d | , (15)where Λ coll is the cooling rate per unit volume via gas-dust collisions for which we follow Burke & Hollenbach(1983): Λ coll = (cid:90) a max a min k B | T g − T d | ¯ α T n g n d ( a ) σ d ( a ) v th d a. (16)In the above equation, a min and a max denote the min-imum and maximum dust grain sizes, k B is the Boltz-mann constant, ¯ α T = 0 . is the thermal accommoda-tion coefficient that characterizes the efficiency of theheat transfer between gas molecules and dust grains, n g is the number density of gas molecules, n d ( a ) is thenumber density of dust particles with size a , σ d = πa is the geometrical cross section of dust particles, and v th = (8 k B T g /πm g ) / is the thermal velocity of thegas. Using second and third moment of the dust grainsize, which are defined as (cid:104) a (cid:105) = (cid:90) a max a min a n d ( a )da (17)and (cid:104) a (cid:105) = (cid:90) a max a min a n d ( a )da , (18) bservational Signatures of Buoyancy Resonances Λ coll = 2 k B | T g − T d | ¯ α T (cid:18) ρ s (cid:19) (cid:18) (cid:104) a (cid:105)(cid:104) a (cid:105) (cid:19) (cid:32) ρ v th m g (cid:33) (cid:18) ρ d ρ g (cid:19) , (19)where ρ s is the bulk density of dust grains while ρ d is themass density of dust grains. Then, combining Equations(15) and (19), the collisional timescale can be written as t coll = c V k B ¯ α T (cid:18) ρ s (cid:19) (cid:18) (cid:104) a (cid:105)(cid:104) a (cid:105) (cid:19) (cid:18) m g ρ g v th (cid:19) (cid:18) ρ d ρ g (cid:19) − . (20)Note that the collisional timescale depends on the meandust grain size and the local dust-to-gas mass ratio,which are dependent upon the grain size distributionand the level of turbulence of the disk among many oth-ers.Then, taking into account radiation, diffusion, andgas-dust collisional timescales, the relaxation timescaleof the gas can be written as t relax = t coll + max( t rad , t diff ) . (21)Again, this equation shows that the collision between gasand dust has to precede thermal emission of dust grains.When the gas-dust collision is not sufficiently frequent,the collisional energy exchange can be the bottleneck ofthe cooling process and thus the overall cooling timescaleis determined by the collisional timescale (i.e., t relax (cid:39) t coll ). When the gas-dust collision is frequent, the gasand dust are in thermal equilibrium and the gas coolsover the timescale t rad or t diff depending on the opticallythickness of the disk (i.e., t relax (cid:39) max( t rad , t diff ) ). Thetransition between the optically thick and thin regimesoccur when t diff (cid:39) t rad (Malygin et al. 2017).In order to make a quantitative comparison between t rad , t coll , and t diff , we compute the timescales usingthe disk model described in Section 2.1. To do so,we first need to define the spatial/size distribution ofdust grains as it dictates the gas-dust collision rate butalso the radiative/diffusion timescales through the opac-ity. We adopt a maximum grain size that is decreas-ing over radius: a max ( R ) = 1 mm ( R/
30 au) − . Thischoice is motivated by the fact that the (sub-)millimetercontinuum emission of the TW Hya disk is confinedwithin about 60 au (e.g., Andrews et al. 2016; Tsuk-agoshi et al. 2016). With the vertically-integrated totaldust-to-gas mass ratio fixed to 0.01 at each radius (i.e., Σ d / Σ g = 0 . ), we distribute the dust mass between . µ m and a max adopting a power-law dust size dis-tribution with a power-law index − . : n d ( a ) ∝ a − . .We then determine the vertical scale height of each dust species assuming that vertical settling is balanced byturbulence mixing characterized by α = 10 − . This re-sults in a dust scale height of H d ( R, a ) = H g ( R ) × min (cid:18) , (cid:114) α min(St( a ) , / a ) ) (cid:19) (22)(Dullemond & Dominik 2004; Birnstiel et al. 2010),where St ( a ) ≡ ( ρ s a/ρ g v th )Ω K is the Stokes number ofparticle having size a and we assume a grain internaldensity ρ s = 1 .
67 g cm − (see below). Then, the verti-cal density distribution of each dust species is obtainedfollowing ρ d ( R, Z, a ) = Σ d ( R, a ) √ πH d ( R, a ) exp (cid:18) − Z H d ( R, a ) (cid:19) . (23)The resulting mean dust particle size (cid:104) a (cid:105) / (cid:104) a (cid:105) anddust-to-gas mass ratio ρ d /ρ g are shown in Figure 2a andb. Larger grains settle near the midplane due to shortersettling times. Since larger grains contain a larger frac-tion of the total dust mass, the dust-to-gas mass ratiodecreases over height. Dust grains with sizes (cid:38) µ m,which are believed to dominate the (sub-)millimeter con-tinuum emission, are confined in radius within the inner ∼ au, while µ m-sized grains extend much further out,consistent with both (sub-)millimeter and optical/near-infrared observations of TW Hya (e.g., Andrews et al.2016; Tsukagoshi et al. 2016; Debes et al. 2013, 2017;van Boekel et al. 2017).Next, we calculate the dust opacity in each grid cellbased on the dust distribution obtained as above, adopt-ing the opacity model from the DSHARP collaboration(Birnstiel et al. 2018). In this opacity model, grainsare assumed as a mixture of water ice, astronomical sil-icates, troilite, and refractory organic material, havinga bulk density of .
67 g cm − . The optical constantsto compute the DSHARP opacity are originally fromHenning & Stognienko (1996); Draine (2003); Warren &Brandt (2008). The absorption and scattering opacitiesare calculated following κ abs , sca ν ( R, Z ) = (cid:80) i κ abs , sca ν ( a i ) ρ d ( R, Z, a i ) (cid:80) i ρ d ( R, Z, a i ) , (24)where the subscription ν shows the frequency depen-dency of the opacity. The Rosseland and Planck meanopacities are calculated as κ R ( R, Z ) = (cid:34) (cid:82) ∞ κ ext ν ( R,Z ) dB ν dT dν (cid:82) ∞ dB ν dT dν (cid:35) − (25)and κ P ( R, Z ) = (cid:82) ∞ κ abs ν ( R, Z ) B ν dν (cid:82) ∞ B ν dν , (26) Bae et al. log (/)/1 µ m -1 0 1 2 30 30 60 90 120 150 180radius [au]0.00.10.20.30.4 Z / R (a) log ρ d / ρ g -8 -7 -6 -5 -4 -3 -20 30 60 90 120 150 180radius [au]0.00.10.20.30.4 Z / R - - - (b) log t diff /t orb -5 -4 -3 -2 -1 0 1 2 3 4 50 30 60 90 120 150 180radius [au]0.00.10.20.30.4 Z / R - - (d) log t rad /t orb -5 -4 -3 -2 -1 0 1 2 3 4 50 30 60 90 120 150 180radius [au]0.00.10.20.30.4 Z / R - (c) log t coll /t orb -5 -4 -3 -2 -1 0 1 2 3 4 50 30 60 90 120 150 180radius [au]0.00.10.20.30.4 Z / R (e) log t relax /t orb -5 -4 -3 -2 -1 0 1 2 3 4 50 30 60 90 120 150 180radius [au]0.00.10.20.30.4 Z / R (f) Figure 2. (a) The mean dust grain size ¯ a ≡ (cid:104) a (cid:105) / (cid:104) a (cid:105) . (b) The dust-to-gas mass ratio ρ d /ρ g . The (c) radiative, (d) diffusion,and (e) collisional timescales in units of the local orbital time. (f) The relaxation timescale t relax , defined as in Equation (21).In panel (f), the shaded region in gray shows where t diff ≥ t coll . We emphasize that the thermal relaxation of the gas is limitedby infrequent gas-dust collision in the most region of the disk ( Z/R (cid:38) . ) and t relax is comparable to or longer than the orbitaltimescale. The white curves in panel (f) show the emission surfaces of (upper) CO J = 3 − and (lower) CO J = 3 − lines from the synthetic observation presented in Section 4, based on the standard model with a . M Jup planet. bservational Signatures of Buoyancy Resonances Z / R with gas-dust collision -5-4-3-2-1012345 l og t r e l a x N Z Z / R without gas-dust collision - - - Figure 3. t relax N Z in a logarithmic scale (left) considering gas-dust collision and (right) without considering gas-dust collision.Note that when gas-dust collision is considered t relax (cid:38) N − Z in the entire disk except near the midplane, providing favorableconditions for buoyancy resonances to develop. The shaded region in gray shows where t diff ≥ t coll . where κ ext ν = κ abs ν + κ sca ν .The radiation timescale t rad calculated as in Equation(10) is shown in Figure 2c. The radiation timescale is or-ders of magnitude shorter than the dynamical timescaleeverywhere in the disk. At a given radius, the radia-tion timescale decreases over height because of its steeptemperature dependency ( t rad ∝ T − ).Figure 2d presents the diffusion timescale t diff . Diffu-sion is slower when the optical depth is larger. So t diff is greater at smaller radii and near the midplane, andis a decreasing function of R and Z . In particular, notethat t diff drops exponentially over height because of thedensity dependency ( t diff ∝ ρ ). We also note that thediffusion timescale depends on the temperature pertur-bation length scale as t diff ∝ λ . Here, we opt to use λ diff = H g . In reality, the length scale of any perturba-tions can range from (cid:28) H g to the thickness of the disk,which is a few scale heights. However, as we will showbelow (see also Section 5.1), the thermal relaxation inthe surface layers ( Z/R (cid:38) . ) is limited by infrequentgas-dust collision, insensitive to the choice of λ diff .Figure 2e shows the gas-dust collisional timescale t coll .For given gas density and temperature structures, thecollisional timescale is set by the mean grain size and thedust-to-gas mass ratio ( t coll ∝ (cid:104) a (cid:105) / (cid:104) a (cid:105) , ( ρ d /ρ g ) − ).As the dust-to-gas mass ratio drops exponentially overheight, gas molecules have significantly less frequentcollisions with dust grains. This makes the collisionaltimescale orders of magnitude longer than the dynami-cal timescale in the surface layers.The relaxation timescale t relax , calculated as in Equa-tion (21), is shown in Figure 2f. The plot clearly showsthat the assumption of thermal equilibrium between gasand dust is not necessarily valid in the outer and sur-face regions of the disk because of the long collisional timescale there. Note that this picture is in a goodagreement with previous studies of thermal relaxationin protoplanetary disks (e.g., Malygin et al. 2017; Bar-ranco et al. 2018; Pfeil & Klahr 2019). In Section 5.1,we explore how different assumptions on the diffusionlength scale, level of disk turbulence, grain size distribu-tion, disk mass, and the existence of a gap in the diskcan affect the cooling timescales. As we will show, thethermal relaxation of the gas in surface layers of proto-planetary disks is limited by infrequent gas-dust collisionover a broad range parameter space.As we mentioned in Section 1, buoyancy resonancesrequire adiabatic responses to thermal perturbations tofully develop (i.e., t relax (cid:38) N − Z ). To see how t relax com-pares with N − Z , we present t relax N Z in Figure 3. Asshown, when gas-dust collision is taken into account, t relax (cid:38) N − Z in the entire disk except at the midplanewhere N Z is zero due to the symmetry across the mid-plane. This suggests that most part of the disk, in-cluding the surface layers CO lines probe, has favorableconditions for buoyancy resonances to develop. In con-trast, if gas-dust collision is neglected, the surface layershave t relax (cid:28) N − Z suggesting that buoyancy resonancesare weak or unlikely to develop there.In addition to the cooling processes we consideredabove, atomic and molecular line cooling plays a rolein the surface layers of protoplanetary disks (e.g., Gortiet al. 2011; Du & Bergin 2014; Kama et al. 2016; Fac-chini et al. 2018). The exact height beyond whichline cooling dominates depends on the underlying ther-mal/chemical properties of the disk as well as the stel-lar/external irradiation. While including comprehen-sive thermo/photo-chemistry requires full thermochem-ical radiative transfer calculations, which is beyond thescope of the paper, here we test the potential effect line Bae et al. cooling would have on buoyancy resonances by assum-ing that the disk gas cools efficiently beyond a certainheight of the disk Z = Z line . Taking this into account,we adopt the following form for the relaxation timescaleconsidering diffusion, radiation, gas-dust collision, andline cooling t relax = [ t coll + max( t rad , t diff )] exp (cid:34) − (cid:18) ZZ line (cid:19) − (cid:35) . (27)The exponential term on the right-hand-side of the equa-tion is added to mimic the effect of line cooling .In our fiducial models, we adopt Z line = 4 H g whichcorresponds to Z/R = 0 . at the radial location of theplanet in our simulations, 90 au. This choice is mo-tivated by the thermochemical model of TW Hya pre-sented in Kama et al. (2016), where the major atomicline emission, including [C I], [C II], and [O I] lines,originates from Z/R ∼ . − . (see their Figure D.2.).However, it is important to note that the exact line cool-ing rate is dependent upon various factors, includingthe gas phase abundance of the coolants and electronnumber density which is determined by UV flux as wellas full chemical chains. To test the effect of line cool-ing, we ran additional simulations adopting Z line = 2 H g ( Z/R = 0 . at 90 au; Section 5.2). In this model, the CO molecular line probes the layers where cooling israpid, dominated by line cooling.For the sake of reproducing the hydrodynamic simu-lations we will present in the following sections, we pro-vide parameterized fits to the diffusion and collisionaltimescales: t diff t orb = 2 . (cid:18) R
90 au (cid:19) − . exp (cid:34) − (cid:18) ZZ diff (cid:19) . (cid:35) , (28)where Z diff = 2 . R/
90 au) . , and t coll t orb = 0 . (cid:18) R
90 au (cid:19) − . exp (cid:34)(cid:18) ZZ coll (cid:19) . (cid:35) , (29)where Z coll = 11 . R/
90 au) . . HYDRODYNAMIC SIMULATIONS3.1.
Hydrodynamic Equations Solved
We solve the hydrodynamic equations for mass,momentum, and energy conservation in the three-dimensional spherical coordinates ( r, θ, φ ) using FARGO The exponent − is chosen such that t relax falls sufficientlyrapidly over height so t relax (cid:28) t orb at the upper boundary of thesimulation domain.
3D (Benítez-Llambay & Masset 2016; Masset 2000): ∂ρ g ∂t + ∇ · ( ρ g v ) = 0 , (30) ρ g (cid:18) ∂v∂t + v · ∇ v (cid:19) = −∇ P − ρ g ∇ (Φ ∗ +Φ p )+ ∇· Π , (31) ∂e∂t + ∇ · ( ev ) = − P ∇ · v + Q relax . (32)In the above equations, ρ g is the gas density, v is thevelocity vector, P is the gas pressure, Φ ∗ = − GM ∗ /r is the gravitational potential of the central star havingmass M ∗ , Φ p is the gravitational potential of the planet, Π is the viscous stress tensor, e is the internal energyper unit volume, and Q relax is the rate at which the gasthermally relaxes to the initial state (see Section 2.2).Note that Q relax can be either a positive or a negativevalue depending on whether the gas is colder or hotterthan the initial equilibrium temperature.The gravitational potential of the planet is computedas Φ p ( r, θ, φ ) = − GM p ( | r − r p | + s ) / , (33)where M p is the mass of the planet, r and r p are three-dimensional radius vectors of the center of the grid cellin question and of the planet, and s is the smoothinglength. Since the smoothing length in three-dimensionalcalculations is used only to avoid the singularity in thepotential on the grid scale, we adopt the cell diago-nal size at the position of the planet for the smoothinglength. We insert the planet at R p = 90 au with a fixed,circular orbit. We use three planet masses: M p = 0 . , ,and M Jup . We run simulations for 500 planetary or-bits. The planet mass is linearly increased over the first5 planetary orbits.For the thermal evolution, we adopt an adiabaticequation of state with an adiabatic index γ = 1 . .The gas pressure and the internal energy are relatedas P = ( γ − e . In addition to the thermal energyevolution via P d V work, which is accounted for by thefirst term of the right-hand-side of Equation (32), cool-ing/heating of the gas is realized through the relaxationof the temperature T g towards the initial temperature T g , init (described in Equations 4 - 6) over the thermalrelaxation timescale t relax computed in Section 2. Thethermal relaxation rate can be written as Q relax = − ρ g c V T g − T g , init t relax . (34)In practice, we use max( t relax , ∆ t hydro ) in the denomina-tor of Equation (34) where ∆ t hydro is the hydrodynamic bservational Signatures of Buoyancy Resonances Figure 4.
Results from the standard model with M p = 0 . M Jup taken at t = 5 t orb . From left to right, the perturbed density δρ/ (cid:104) ρ (cid:105) φ , perturbed temperature δT / (cid:104) T (cid:105) φ , and vertical velocity v Z in units of local sound speed at Z/R = 0 . ( (cid:39) H at 90 au)in a φ − r plane. Here, (cid:104)(cid:105) φ denotes average over the azimuthal direction and δρ ≡ ρ − (cid:104) ρ (cid:105) φ and δT ≡ T − (cid:104) T (cid:105) φ . The planet islocated in the midplane, at ( φ, r ) = (0 radian ,
90 au) . The vertical velocity v Z is computed from radial and meridional velocitiesin the spherical coordinates. For a more straightforward comparison with observations later in Section 4, motions toward thedisk midplane are shown with red colors and positive numbers while motions toward the disk surface are shown with blue colorsand negative numbers, such that they match with the red- and blue-shift convention. timestep, in order to avoid over-relaxation that can hap-pen when t relax < ∆ t hydro .For our standard model, we take into account ra-diation, diffusion, gas-dust collision, and line cooling(Section 3.3.1). In practice, this is done by adopt-ing a prescribed relaxation timescale using Equation(27), along with the fits in Equations (28) and (29).As a comparison to the standard model, we addition-ally carry out simulations without gas-dust collision:i.e., t relax = max( t rad , t diff ) exp[ − ( Z/Z line ) − ] (Section3.3.2; hereafter gas-dust thermal equilibrium model).This model is similar to what is often used in three-dimensional protoplanetary disk simulations where it isimplicitly assumed that the gas and dust have instanta-neous energy balance.3.2. Simulation Setup
The simulation domain extends from 30 au (= . R p )to 210 au (= . R p ) in r , from π/ − . to π/ in θ which covers 5.6 scale heights at the radial locationof the planet, and from 0 to π in φ . We adopt 460logarithmically-spaced grid cells in the radial direction,96 uniformly-spaced grid cells in the meridional direc-tion, and 1482 uniformly-spaced grid cells in the az-imuthal direction. With this choice, one gas scale height at the location of the planet is resolved with about 18grid cells in all direction.At the radial boundaries, we adopt a wave-dampingzone to suppress wave reflection (de Val-Borro et al.2006). At the lower meridional boundary, which isthe disk midplane, we adopt the symmetric boundarycondition for all variables but the meridional velocityfor which we apply the reflecting boundary condition.At the upper meridional boundary we adopt the zero-gradient boundary condition.We adopt a kinematic viscosity characterized by α =10 − , a value broadly consistent with the level of turbu-lence observationally constrained for the TW Hya disk(Teague et al. 2016; Flaherty et al. 2018).In order to ensure that no other hydrodynamic insta-bilities operate in the disk, we ran both standard modeland gas-dust thermal equilibrium model in the absence aplanet. We found that no hydrodynamic instabilities de-velop in these runs. The vertical shear instability (Urpin& Brandenburg 1998; Nelson et al. 2013) is suppresseddue to the non-zero viscosity and a long thermal relax-ation time, in agreement with previous studies (Nelsonet al. 2013). 3.3. Simulation Results Bae et al.
Standard Model
We start by discussing the results from our standardmodel where we consider radiation, diffusion, gas-dustcollision, and line cooling for the thermal relaxation ofthe disk gas. Figure 4 shows the perturbed density, per-turbed temperature, and vertical velocity in a φ − r planet at Z/R = 0 . for the M p = 0 . M Jup model.This height corresponds to about three scale heightsabove the midplane at the radial location of the planetand is close to the CO emission surface in the syntheticobservation we will present in Section 4.As most clearly shown in the perturbed density distri-bution, the planet excites a pair of primary Lindblad spi-rals, one in the inner disk (the arc crossing the r = 40 auboundary at φ ∼ π/ ) and one in the outer disk (thearc crossing the r = 140 au boundary at φ ∼ − π/ ).The Lindblad spirals are nearly perpendicular to the az-imuth axis near the planet, suggesting that they havea large pitch angle close to 90 ◦ there. In addition tothe primary Lindblad spirals, the planet excites a sec-ondary Lindblad spiral in the inner disk, which emergesat r ∼ au and φ ∼ π . At Z/R = 0 . , density, tem-perature, vertical velocity perturbations created by theLindblad spirals of the . M Jup planet are a few to
10 % of the background values or the local sound speed.Along with the Lindblad spirals, the planet excites afamily of spirals via buoyancy resonances, which can bemost clearly seen in the vertical velocity plot. Two mainfeatures that distinguish buoyancy spirals from Lindbladspirals are (1) the tightly wound morphology and (2) thelarge vertical motions, which we will explain in detailone by one.In Figure 5, we show the measured pitch angle of Lind-blad and buoyancy spirals with a 5 au interval in radius.For Lindblad spirals, we measure the pitch angle usingthe peak in the density perturbation. For buoyancy spi-rals, we measure the pitch angle using the peak in posi-tive vertical velocities. We opt to use the vertical veloc-ity instead of the density because (1) buoyancy spiralsare more clearly identifiable with the vertical velocityand (2) we can make a more direct comparison to theTW Hya spiral detected in the velocity space (Section4.1.1).We also present the pitch angle derived with lineartheory. For Lindblad spirals, the phase angle φ L (i.e.,the azimuthal angle from the spiral to the planet) as afunction of the cylindrical disk radius R is φ L ( R ) = − sgn( R − R p ) π m − (cid:90) RR ± m Ω( R (cid:48) ) c s ( R (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:32) − R (cid:48) / R / p (cid:33) − m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) / d R (cid:48) (35)
30 60 90 120 150radius [au]010203040 p it c h a ng l e [ d e g ]
30 60 90 120 150radius [au]010203040 p it c h a ng l e [ d e g ] Lindblad spiralbuoyancy spiralTW Hya spiralz=1H g z=3H g Figure 5.
The pitch angles of (black curves) Lindblad and(red curves) buoyancy spirals calculated based on the lineartheory. For each spirals, the lower and upper curves show thelinear theory prediction at Z = 1 H g and H g , respectively.Square and diamond symbols present measured pitch anglesfrom the hydrodynamic simulation with M p = 0 . M Jup at(square) Z = 1 H g and (diamond) H g . Blue dots presentthe pitch angle of the spirals observed in the TW Hya disk(Teague et al. 2019). The gray shaded regions show ± ± ± following Bae & Zhu (2018a), where m = (1 / H/R ) − p and R ± m = (1 ± /m ) / R p .For buoyancy spirals, the phase angle φ B is φ B ( R ) = ± nπ | Ω p − Ω K | Ω K H g Z (cid:18) Z R (cid:19) / (36) × (cid:20) γ − γ + 2 c s g ∂c s ∂Z (cid:21) − / , following Zhu et al. (2015), where n is a positive inte-ger and g is the gravity from the star. Compared withthe analytic form given in Zhu et al. (2015), note thatEquation (36) has an additional term in the parenthesisof the right-hand-side, ( c s /g ) ∂c s /∂Z , which takes thevertical temperature stratification into account. Whenthere is no vertical temperature stratification Equation(36) reduces to the phase angle derived in Zhu et al.(2015). The pitch angles of Lindblad and buoyancy spi-rals are computed as β ≡ − tan − (d R/R d φ ) , using thephase equations in Equations (35) and (36).As shown in Figure 5, buoyancy spirals’ pitch angle issmaller than Lindblad spirals’ pitch angle over a broadrange of distance from the planet. In particular, withina few scale heights from the planet Lindblad spirals’pitch angle increases to ◦ toward the planet, whereasbuoyancy spirals’ pitch angle remains (cid:46) ◦ . We notethat Figure 5 shows pitch angles of the first-order buoy-ancy spirals only (i.e., n = 1 in Equation 36) for bothlinear theory prediction and measurement from the sim- bservational Signatures of Buoyancy Resonances Figure 6.
The two-dimensional distributions of (first row) perturbed density δρ/ (cid:104) ρ (cid:105) φ , (second row) vertical velocity v Z , (thirdrow) radial velocity v R , and (fourth row) perturbed azimuthal velocity v φ − (cid:104) v φ (cid:105) φ , taken at the end of the simulation with M p = 0 . M Jup . From left to right in each row, three panels present the distributions at
Z/R = 0 (the midplane), at
Z/R = 0 . ( (cid:39) H at 90 au), and Z/R = 0 . ( (cid:39) H at 90 au), respectively. The two arrows in the upper left panel point to the inner andouter Lindblad spirals, while the four arrows in the upper middle panel point to the first- and second-order, inner and outerbuoyancy spirals. Dotted circles show the planet’s orbit. The buoyancy spirals can be best distinguished from Lindblad spiralsusing vertical velocities. Bae et al. δρ /< ρ > φ - π πφ [radian]-0.20.00.20.40.60.8 midplane Z/R=0.08Z/R=0.23 δ T/
The azimuthal profiles of (from left to right) perturbed density δρ/ (cid:104) ρ (cid:105) φ , perturbed temperature δT / (cid:104) T (cid:105) φ , verticalvelocity v Z , radial velocity v R , and perturbed azimuthal velocity v φ − (cid:104) v φ (cid:105) φ from the standard model with M p = 0 . M Jup . Theupper panels present the profiles at two scale heights beyond the planet ( R = R = Z/R = 0 . ( (cid:39) H g at 90 au), and Z/R = 0 . ( (cid:39) H g at 90 au), respectively. The red dashed lines showthe azimuthal locations where the vertical velocity peaks. bservational Signatures of Buoyancy Resonances n > ),the pitch angle is smaller than that of first-order buoy-ancy spirals by a factor of (cid:39) /n , meaning that they aremore tightly wound.It is also worth pointing out that buoyancy spirals’pitch angle varies continuously as a function of radiuswithout a singularity at the radial location of the planet,unlike Lindblad spirals. From the observational point ofview, this implies that the inner and outer buoyancyspirals can appear connected to each other as a singlespiral, especially when the spatial resolution is insuffi-cient.Another characteristic that distinguishes buoyancyspirals from Lindblad spirals is the large vertical mo-tions. In order to visualize the density and velocityperturbations at different heights, we present the per-turbed density and vertical, radial, azimuthal velocitiesat the midplane, Z/R = 0 . , and Z/R = 0 . in Figure6. The azimuthal profiles of the perturbed density, per-turbed temperature, and vertical, radial, and azimuthalvelocities at R = R p ± H g are presented in Figure 7.Focusing on the Lindblad spirals first, both densityand velocity perturbations driven by Lindblad spiralsis the strongest at the midplane and decreases overheight. The only exception is the near-zero vertical ve-locity at the midplane, which is because of the symme-try across the midplane. For the buoyancy spirals, onthe other hand, perturbations are the smallest at themidplane because N Z = 0 there. Between Z/R = 0 . and 0.23, perturbations remain comparable or becomestronger than their Lindblad counterparts. In particu-lar, we note that the vertical velocity perturbation asso-ciated with buoyancy spirals become stronger and moreextended in azimuth over height. Compared to Lind-blad spirals, buoyancy spirals generally produce smallerperturbations, but the vertical motions associated withbuoyancy spirals can be stronger especially as we moveto the surface layers. These characteristics of buoyancyspirals suggest that the best strategy to observe buoy-ancy spirals is to look for vertical velocity perturbationsin the surface layers of face-on disks.Finally, we note that there is a π/ phase shift be-tween the vertical velocity and density/temperature per-turbations driven by buoyancy resonances. At R =103 . au, for example, the vertical velocity has peaksat φ = − . , − . , − . , and − . radians and thedensity and temperature perturbation is close to zero atthose azimuthal locations (see the red dashed lines inFigure 7).3.3.2. Gas-dust Thermal Equilibrium Model
We now turn our discussion to the gas-dust thermalequilibrium model in which the gas is assumed to bethermally coupled with the dust. With this assumptionthe relaxation timescale in the surface layers is muchshorter than the dynamical time (Figure 2c) and buoy-ancy resonances are expected to be weak or absent inthe surface layers (Figure 3 right).The perturbed density, perturbed temperature, andvertical velocity at
Z/R = 0 . is presented in Figure8. In Figure 9 we present the azimuthal profiles of theperturbed density, perturbed temperature, and vertical,radial, and azimuthal velocities at R = R p ± H g fromthe planet in the radial direction. As apparent from thefigures, buoyancy spirals do not generate clear densityand temperature perturbations. Vertical velocity per-turbations arising from buoyancy resonances are seen,but they are much weaker compared with the standardmodel. We note that the resonance is not completelysuppressed in this model, as suggested by the verticalvelocity perturbation, because the buoyancy frequencyis not strictly zero with the imposed stratified disk tem-perature (see Equation 2).Speaking of Lindblad spirals, we find that they pro-duce smaller perturbations as we move toward the sur-face layers, in agreement with what is seen in the stan-dard model. However, we note that the level of densityperturbations at Z/R = 0 . and . are larger thanthe standard model by about a factor of two. Note alsothat the Lindblad spirals are more tightly wound – com-pare the azimuthal angles where Lindblad spirals meet r = 40 and 140 au in Figure 4 and 8. This is becausethe gas behaves (nearly) isothermally and thus the soundspeed is smaller than the standard model where the gasbehaves adiabatically. SIMULATED OBSERVATIONSIn order to examine the observability of buoyancy spi-rals, we carry out synthetic observations of the CO J = 3 − and CO J = 3 − lines using the 3Dradiative transfer code RADMC-3D (Dullemond et al.2012). To do so, we first add an inner disk inward of thecomputational domain, from 1 to 30 au, using the gasdensity profile described with Equation (3). We assumeCO is photodissociated at the surface layers, where thesum of the vertically integrated gas column density (toconsider external irradiation) and the radially integratedgas column density from the central star (to consider thecentral star’s irradiation) is less than cm − (Visseret al. 2009). We assume CO is frozen onto grains in the In practice, this is done along the meridional direction in thespherical coordinates. Bae et al.
Figure 8.
Same as Figure 4, but for the gas-dust thermal equilibrium model. Due to the short thermal relaxation timescalein the surface layers, buoyancy resonances do not fully develop. δρ /< ρ > φ - π πφ [radian]-0.20.00.20.40.60.8 midplane Z/R=0.08Z/R=0.23 δ T/
Same as Figure 7, but for the gas-dust thermal equilibrium model. regions where the temperature is below 21 K (Schwarzet al. 2016). We adopt a CO-to-H ratio of − , whichis smaller than the canonical interstellar medium valueof − , motivated by the fact that CO in the TW Hya isknown to be largely depleted in the gas phase (Schwarzet al. 2016; Zhang et al. 2019). For CO, we adopt CO/ CO ratio of 1/60. As opposed to running Monte Carlo calculations tocompute the dust temperature and assuming the gasand dust are thermally coupled, an approach often takenin the post-processing of hydrodynamic simulations, weadopt the gas temperature from the hydrodynamic sim-ulations. This is because the gas and dust do not nec-essarily have the same temperature as we discussed inSection 2.2. bservational Signatures of Buoyancy Resonances Figure 10.
The vertical velocity v Z in a X − Y plane in units of m s − at (upper panels) Z/R = 0 . ( (cid:39) H g at 90 au) and(lower panels) Z/R = 0 . ( (cid:39) H g at 90 au), after 500 orbits. From left to right, results from standard 0.5 M Jup , 1 M Jup ,and 2 M Jup models, respectively. The planet is located at ( X, Y ) = ( −
90 au , . Note that the vertical motions driven byLindblad spirals are smaller than that those driven by buoyancy spirals, especially at Z/R = 0 . . The black and white dashedcurves in the lower left panel trace the buoyancy and Lindblad spirals, respectively. For the fiducial disk geometry, we use a disk positionangle (PA) of ◦ , defined as the angle from the northto the redshifted major axis in the counter-clockwise di-rection, and an inclination of ◦ , comparable to thoseof the TW Hya disk. We test the influence of varyingdisk inclination on the kinematic signatures of buoyancyspirals in Section 4.3. The planet is placed at PA= ◦ (i.e., East) in all cases. The disk rotates clockwise onthe sky.We create image cubes at 10 m s − velocity resolutionand average down to the desired velocity resolution of100 m s − . This is because most radiative transfer codesincluding RADMC-3D return emission at the centralfrequency of each channel, rather than the integratedemission across the channel (see Rosenfeld et al. 2013for a demonstration of this effect). We then convolvethe image cubes with a circular Gaussian beam with aFWHM of . (cid:48)(cid:48) . Observations at this high spatial res- olution will be the product of multiple executions withdiffering array configurations and observing conditionsin order to fill in the uv plane. As we are not tryingto recreate specific observations, where this would be anecessary step, but rather provide a quantitative pre-diction, we opt to simply convolve each channel with acircular Gaussian beam. We add correlated (both spa-tially due to the beam and spectrally due to the Hanningsmoothing) noise to each channel with specified RMS of1 K for CO and 0.3 K for CO, which correspond to0.94 and 0.28 mJy beam − , respectively.We provide simulated CO and CO cubes (averageddown to 100 m s − velocity resolution but without beamconvolution and correlated noise) at https://doi.org/10.5281/zenodo.4361639.4.1. Buoyancy Spirals in Keplerian-subtracted MomentMaps Bae et al.
Figure 11.
Keplerian-subtracted moment maps with simulated (upper panels) CO and (lower panels) CO observationsadopting 0.1” spatial resolution, 100 m s − spectral resolution, and 1.0 K and 0.3 K rms noise level, respectively. From left toright, results with standard 0.5 M Jup , 1 M Jup , and 2 M Jup models. The beam is shown in the lower left corner of each panel.
Before we present Keplerian-subtracted momentmaps, we show the vertical velocity distribution fromstandard . M Jup , M Jup , and M Jup models in Fig-ure 10. Buoyancy resonances produce vertical motionsof order of 100 m s − , with a larger magnitude for moremassive planets. As discussed earlier, Lindblad spiralsproduce much smaller vertical perturbations comparedwith buoyancy spirals, especially in the surface layers.We thus do not expect to detect them in CO in face-ondisks.We generate Keplerian-subtracted moment mapsusing the quadratic method described in Teague &Foreman-Mackey (2018). Using the Python package eddy (Teague 2019), we fit a Keplerian rotation profileto the rotation maps, allowing the source center, thedisk inclination and position angle and the systemic velocity to vary. Given the low inclination of TW Hya( i ∼ ◦ ) we do not include any terms describing anelevated emission surface.The Keplerian-subtracted moment maps are shown inFigure 11. We note that the residual maps clearly showa coherent, tightly wound spiral structure. It is alsoworth pointing out the excellent agreement between theinput velocity field from hydrodynamic simulations andthe retrieved velocities.4.1.1. Case Study: TW Hya
Interestingly enough, the tightly wound morphologyof buoyancy spirals resemble to the spiral seen in theKeplerian-subtracted moment map of TW Hya (Teagueet al. 2019). As presented in Figure 5 the buoyancy spi-rals driven by a planet at 90 au can explain the small bservational Signatures of Buoyancy Resonances − − O ff s e t [ a r c s ec ] − − O ff s e t [ a r c s ec ] − − − − O ff s e t [ a r c s ec ] − − − − − − − − I n t e n s i t y [ m J y b e a m − ] Figure 12.
Channel maps from a simulated CO (3-2) line observation based on the hydrodynamic simulation with 2 M Jup planet. The channel velocity, relative to the central channel, is presented in the upper left corner of each panel. Buoyancyspirals are pointed out with white/black arrows in relevant panels. The beam is presented in the lower left corner of each panel. pitch angle as well as the monotonically decreasing pat-tern over radius. The magnitude of observed velocityperturbations ( ∼
40 m s − ) is also in a good agreementwith our simulations.We note that it is not impossible to explain the smallpitch angle of the TW Hya spirals with Lindblad res-onance; one may place a planet far inward of the ob-served spirals. However, if this has to be the case, itis unclear why we do not see any large perturbationsin the inner disk from neither line nor continuum ob-servations because we expect the largest perturbationswould arise at the vicinity of the planet (Figure 6). Inaddition, given the low inclination of the TW Hya diskwe are likely seeing vertical motions. Our simulationsshow that the outer Lindblad spiral produces little ver-tical motions ( (cid:28)
50 m s − ) and it is unlikely that thevelocity spiral in the TW Hya disk is associated withLindblad spirals.We thus propose that, if the spirals in TW Hya aredriven by an embedded planet, the spirals could be as-sociated with buoyancy resonances driven by a (sub-)Jovian-mass planet at around 90 au.4.2. Buoyancy Spirals in Velocity Channel Maps
When the perturbations driven by buoyancy spiralsare sufficiently large, the spirals can be seen in velocitychannel maps. We present channel maps of the synthetic CO line observation from the standard M Jup modelin Figure 12. Channel maps from standard . M Jup and M Jup models are presented in Figures 22 and 23in Appendix B.The main characteristic of buoyancy spirals in chan-nel maps is wedge-like features standing out of the so-called butterfly pattern of the Keplerian disk. In Fig-ure 12 these non-Keplerian features are most clearlyseen close to the planet, on the North-East side ofthe disk at v = 0 . − . − . Along the major-axis of the disk, buoyancy spirals can appear as an arcbridging the Keplerian wings (e.g., South-East side at v = − . − , North-West side at v = 0 . − )or as an arc beyond the inner disk (e.g., South-East sideat v = 0 . − . − ).4.3. Effect of Disk Inclination
In order to examine the observational appearance ofbuoyancy spirals in more inclined disks, we generate ad-ditional image cubes adopting 15, 30, 45, and 60 de-grees inclinations. Keplerian-subtracted moment mapsfor standard M Jup model are shown in Figure 13.8
Bae et al.
100 50 0 50 100
Residual (ms ) 100 50 0 50 100 Residual (ms ) 100 50 0 50 100 Residual (ms )202 Offset (arcsec)202 O ff se t ( a r cs e c ) (a) i =5 ◦ Offset (arcsec) (b) i =15 ◦ Offset (arcsec) (c) i =30 ◦ Figure 13.
Keplerian-subtracted moment maps with (a) 5, (b) 15, and (c) 30 degrees inclination for the M Jup model. Thedotted ellipse in each panel shows planet’s radial location r = 1 . (cid:48)(cid:48) on sky. The feather-like feature in panel (c) is due to thechannelization effect owing to the channel spacing of the data. Higher spectral resolution data can suppress this effect (seeChristiaens et al. (2014) for an example). − − O ff s e t [ a r c s ec ] i = 5 ◦ i = 15 ◦ i = 30 ◦ i = 45 ◦ i = 60 ◦ − − O ff s e t [ a r c s ec ] − − − − O ff s e t [ a r c s ec ] − − − − − − − − I n t e n s i t y [ m J y b e a m − ] Figure 14.
Channel maps from a simulated CO (3-2) line observation based on the hydrodynamic simulation with 2 M Jup planet, using various disk inclination: (from left to right) 5, 15, 30, 45, and 60 degrees. The channel velocity, relative to thecentral channel, is presented in the upper left corner of each panel in units of km s − . Buoyancy spirals are pointed out withwhite/black arrows in relevant panels. The beam is shown in the lower left corner of each panel. Although non-Keplerian velocity components arisingfrom buoyancy resonances are still present, they appear more as an ellipse for more inclined geometry, making it bservational Signatures of Buoyancy Resonances i = 15 ◦ and i = 30 ◦ cases, the red-shifted buoyancy spirals appearstronger than i = 5 ◦ case. This is because we are moresensitive to in-plane velocity perturbations for inclineddisks. As the planet opens a gap around its orbit it altersthe gas pressure profile from the unperturbed one, suchthat the disk gas rotates at sub-Keplerian speed insideof the planet’s orbit and at super-Keplerian outside ofthe planet’s orbit (Teague et al. 2018). For an inclineddisk, sub-Keplerian rotation would appear as a blue-/red-shifted semi-ellipse on the red-/blue-shifted side ofthe disk, while super-Keplerian rotation would appear asa red-/blue-shifted semi-ellipse on the red-/blue-shiftedside of the disk (see Figure 5 of Teague et al. 2019). Thepattern we see across the planet’s orbit in the Keplerian-subtracted moment maps in Figure 13 exactly matcheswith the aforementioned expectation. If we first look atthe i = 15 ◦ model, there is an elongated blue-shiftedarc on the red-shifted side (South-East side) of the disk,right inward of the planet’s orbit. On the blue-shiftedside of the disk, we see an elongated red-shifted arc in-ward of the planet’s orbit, along with the red-shiftedbuoyancy spiral. Right outside of the planet’s orbit onthe red-shifted side, the super-Keplerian rotation addsred-shifted residuals to the buoyancy spiral, enhancingthe overall magnitude of the residual. The same patternis consistently observed in the i = 30 ◦ model and therotation modulation is stronger in this case.It is thus reasonable to conclude that disentangling thevertical motions induced by buoyancy resonances andthe modulation in rotational motions associated withthe gap is generally more challenging for inclined disks.However, if a buoyancy spiral is sufficiently extended inazimuth such that the spiral crosses the minor axis ofthe disk, this offers a possibility to disentangle verticalmotions from rotational motions. This is because ver-tical motions do not change their sign across the minoraxis while rotational motions do change their sign.We now turn our attention to velocity channel maps.Representative velocity channels for different inclinationare shown in Figure 14. Channel maps for standard . M Jup and M Jup models are presented in AppendixB. While the characteristic features of buoyancy spi-rals are still present, they are clearly weaker for inclineddisks due to the sin i components of the projection, inan agreement with what we see in Keplerian-subtractedmoment maps. This suggests that disks with small in-clination of (cid:46) ◦ offer the best opportunities to searchfor buoyancy spirals. DISCUSSION5.1.
Thermal Relaxation Timescale
As we have shown with hydrodynamic simulations inSection 3.3, the development of buoyancy resonances de-pends on the local thermal relaxation timescale. To ex-amine the development and observability of buoyancyresonances under various disk conditions, we explore abroad range of parameter space and compute the re-laxation timescale following the approach described inSection 2. Specifically, we vary (1) the diffusion lengthscale λ diff , (2) the dust scale height, (3) the maximumdust grain size, (4) the power-law slope in the dust sizedistribution, and (5) the disk mass. The resulting ther-mal relaxation timescales are shown in Figure 15. Wecompare t relax with N − Z in Appendix A and Figure 21.In what follows, we focus our discussion on the surfacelayers of the disk ( Z/R (cid:38) . ) as we are mostly con-cerned about the observability of buoyancy resonancesusing optically thick CO lines which will trace elevatedregions of the disk.Due to the steep dependency ( t diff ∝ λ ), havinga short diffusion length scale of λ = 0 . H g results ina diffusion timescale that is much shorter than the dy-namical timescale almost everywhere in the disk. Evenin such a case, however, the relaxation of the gas is lim-ited by long collisional timescales between gas moleculesand dust grains. In particular, the relaxation timescalein the surface layers where CO lines probe is deter-mined by the collisional timescale and is insensitive tothe choice of λ diff .The gas-dust collision rate is dependent upon the de-tailed spatial and size distribution of grains. We firstvary the scale height of dust grains by changing α inEquation (22). This has two opposite effects on the colli-sional timescale. For an increased scale height, the meandust size increases which would lengthen the collisionaltimescale (Equation 20). At the same time, the dust-to-gas mass ratio also increases which would shorten thecollisional timescale. The two effects effectively cancelout and the collisional timescale in the surface layersremains sufficiently long for buoyancy resonances to de-velop, as can be seen in Figures 15 and 21.When a max is decreased, the mean grain size de-creases. The dust-to-gas mass ratio increases in the sur-face layers because, with a smaller a max , a larger frac-tion of the total dust mass is in the grains that can belofted sufficiently high. Together, this shortens the col-lisional timescale, resulting in t relax < N − Z in a largerregion of the disk as shown in Figure 21. Nevertheless, t relax (cid:38) N − Z near the CO emission surface, suggestingthat buoyancy resonances likely develop there.Changes in the power-law index of the dust size distri-bution work in a similar way. When the dust size distri-bution follows a steeper power-law distribution, a larger0
Bae et al. standard model Z / R λ diff = 0.1H g λ diff = 3H g α = 10 -2 Z / R α = 10 -4 × a max × a max Z / R n d (a) ∝ a -2.5 n d (a) ∝ a -4.5 M disk = 0.017 M Sun Z / R M disk = 0.15 M Sun with a gap, M p = 1M Jup -5 -4 -3 -2 -1 0 1 2 3 4 5 log t relax /t orb Figure 15.
Same as Figure 2f, which is shown again in the upper left panel, but the relaxation timescale calculated with variousparameters. Unless otherwise stated in the title of each panel, the fiducial parameters adopted to compute the diffusion andcollisional timescales are λ diff = H g , α = 10 − , a max ( R ) = 1 mm ( R/
30 au) − , and n d ( a ) ∝ a − . . The shaded region in grayshows where t diff ≥ t coll . In the bottom right panel, the two dashed lines show the radial locations at which the gap depth is thehalf of the maximum (75 and 105 au). We note that, for the broad range of parameter space explored, the thermal relaxationof the gas in the surface layers of protoplanetary disks (i.e., Z/R (cid:38) . ) is limited due to infrequent gas-grain collisions. bservational Signatures of Buoyancy Resonances t relax (cid:38) N − Z in the majority of the disk regions.Lastly, we adopt the azimuthally-averaged gas densityfrom the standard M Jup model to examine the influ-ence the gap has on the relaxation timescale. Becausedust grains are depleted within the gap the collisionaltimescale becomes longer, facilitating the developmentof buoyancy resonances there.In summary, we argue that the relaxation timescaleof the gas in the surface layers (
Z/R (cid:38) . ) is limitedby infrequent gas-dust collisions under a broad rangeof conditions applicable to protoplanetary disks (Figure15). This is because the dust-to-gas mass ratio is small( ρ d /ρ g (cid:46) − ) in the surface layers due to the verticalsettling of dust grains. The resulting thermal relaxationtimescale is comparable to or longer than the timescaleassociated with buoyancy oscillations N − Z (Figure 21),suggesting that the surface layers have a favorable con-dition for buoyancy resonances to develop.5.1.1. Implications for planet-induced gaps
While we focused our analysis mainly on the surfacelayers so far, it is worth pointing out the finite relax-ation timescale near the disk midplane. Even at largeradii for which it is typically thought that less absorbingmaterials would result in more efficient cooling, we findthat the infrequent gas-dust collision could prevent thegas from cooling instantaneously.The finite relaxation timescale in the main body ofa disk can have important implications for the forma-tion of gaps by planets. Adopting a constant dimension-less relaxation timescale β ≡ πt relax /t orb in vertically-integrated two-dimensional disks, Miranda & Rafikov(2020b) and Zhang & Zhu (2020) independently showedthat the isothermal assumption does not provide a goodapproximation when β (cid:38) − – note again that we con-sistently find β (cid:38) − under various disk conditions(see Figure 15). In particular, when β ∼ the radia-tive dissipation of Lindblad spirals becomes importantand the gap around the planet becomes narrower thanotherwise.In Figure 16, we present the radial profiles of the gassurface density from the standard model and the gas-dust equilibrium model. We additionally ran a simula-tion adopting an isothermal equation of state and in-
30 60 90 120 150 180radius [au]0.00.20.40.60.81.01.2 Σ g / Σ g , i n it standardgas-dustequilibriumisothermal30 60 90 120 150 180radius [au]0.00.20.40.60.81.01.2 Σ g / Σ g , i n it Figure 16.
The azimuthally-averaged gas surface densityas a function of radius. The hill sphere is excluded whenazimuthally averaging the density. Red, blue, black curvesshow the profiles from the standard model, the gas-dust ther-mal equilibrium model, and the isothermal model, respec-tively. The vertical dashed lines present the gap width atthe half maximum. The gray shaded regions show ± ± ± cluded the density profile from the isothermal simula-tion in the same figure. The gap widths measured atthe half maximum are ∆ gap = Implications for hydrodynamic instabilities
Hydrodynamic instabilities are an important sourceof turbulence in protoplanetary disks. Vertical shear in-stability can be largely suppressed with finite thermalrelaxation timescales (Nelson et al. 2013; Lin & Youdin2015; Malygin et al. 2017; Pfeil & Klahr 2019, 2020).Spiral wave instability is known to have a forbidden re-gion near the disk surface where the buoyancy frequencyis larger than a half of the Doppler-shifted forcing fre-quency (Bae et al. 2016a). The forbidden region extendstoward the midplane with more adiabatic gas response2
Bae et al.
Figure 17.
Same as Figure 10 but for models with rapid cooling near the surface using Z line = 2 H g ( Z/R (cid:39) . . (Bae et al. 2016a), so we can infer that spiral wave in-stability operates in a more confined region with finitethermal relaxation timescales. On the other hand, finitethermal relaxation timescales can help other instabil-ities that operate with slower thermal relaxation, suchas convective over stability (Klahr & Hubbard 2014) andzombie vortex instability (Marcus et al. 2015; Barrancoet al. 2018).5.2. Effect of Efficient Line Cooling Near the Surface
As discussed in Section 2.2, exactly at which heightatomic/molecular line cooling becomes the dominantcooling mechanism depends upon the balance betweenvarious heating and cooling mechanisms, which in turnis determined by various factors including the amountof dust grains, gas temperature, abundance of coolantatoms/molecules, and electron number density. Here,we test the effect of efficient line cooling near the sur-face, by adopting Z line = 2 H g in Equation (27) (c.f., Z line = 4 H g in the standard model).Figure 17 presents the vertical velocity distribution.Not surprisingly, buoyancy resonances are weaker at Z = 3 H g ( Z/R = 0 . due to the rapid cooling inthe surface layers. However, we note that buoyancyresonances are not completely suppressed. This is be-cause even when cooling is rapid (effectively γ (cid:39) ),the buoyancy frequency is not strictly zero due to thestratified temperature profile (see Equation 2). At Z = 1 H g ( Z/R = 0 . , buoyancy resonances developat the level they develop in the fiducial model adopting Z line = 4 H g (see Figure 10). This suggests that the de-velopment of buoyancy resonances depends on the local thermodynamic properties and that, as far as line cool-ing is not efficient all the way to the midplane, there willbe regions where buoyancy resonances would develop.From the observational point of view, this suggeststhat we can choose optically thiner lines that probe theadequate heights. To support this argument, we gener-ate Keplerian-subtracted moment maps from the mod-els with efficient line cooling, which are shown in Figure18. The residual velocities are smaller than the fiducialmodel in CO because the line probes the surface layerswhere cooling is rapid. On the other hand, the morphol- bservational Signatures of Buoyancy Resonances Figure 18.
Same as Figure 11 but for models with rapid cooling near the surface using Z line = 2 H g . ogy and magnitude of the buoyancy spirals in the COmaps are nearly identical to the standard model (seeFigure 11).5.3.
Can we observationally distinguish the origin ofnon-Keplerian motions?
Here we discuss potential ways to discriminate thenon-Keplerian motions driven by buoyancy resonancesfrom those driven by other origins, Lindblad resonance,corrugated vertical flows, and gas pressure changes.
Lindblad spirals:
The vertical dependency of the per-turbation driven by buoyancy and Lindblad spirals isthe key to discriminate the two. Because buoyancy fre-quency is strictly zero at the disk midplane, we expectno or weak buoyancy resonances there. On the otherhand, perturbations driven by Lindblad spirals are thestrongest at the disk midplane and decrease over height(Figures 6 and 7). Observations of multiple lines tracing different heights in the disk, for instance CO vs. COor C O, will thus help discriminate between buoyancyspirals and Lindblad spirals.
Corrugated vertical flows:
Various (magneto-)hydrodynamicprocesses are known to create radially-alternating cor-rugated vertical flow patterns, including vertical shearinstability (Nelson et al. 2013), spiral wave instabil-ity (Bae et al. 2016a,b), and magnetically-driven zonalflows (Johansen et al. 2009; Flock et al. 2015; Riolset al. 2020). Vertical velocity perturbations driven bybuoyancy resonances (and Lindblad resonance) are sym-metric against the midplane, so at a given radius the gasmotion will be either toward the midplane or toward thesurface. In contrast, vertical velocity perturbations as-sociated with corrugated vertical flows are not symmet-ric against the midplane. Rather, instabilities developinto corrugation modes and the entire column oscillatesvertically (Nelson et al. 2013; Bae et al. 2016a,b). This4
Bae et al. − O ff s e t [ a r c s ec ] .
64 km s − .
96 km s − .
28 km s − .
60 km s − .
92 km s − − O ff s e t [ a r c s ec ] .
24 km s − .
56 km s − .
88 km s − .
20 km s − .
52 km s − I n t e n s i t y [ m J y b e a m − ] −
101 Offset [arcsec] − O ff s e t [ a r c s ec ] .
84 km s − −
101 Offset [arcsec]9 .
16 km s − −
101 Offset [arcsec]9 .
48 km s − −
101 Offset [arcsec]9 .
80 km s − −
101 Offset [arcsec]continuum 0 . . . . . . I n t e n s i t y [ m J y b e a m − ] Figure 19. CO channel maps of HD 143006, originally presented in Pérez et al. (2018a). Contours are drawn at 2.4, 4.8, and . − . The rms noise is . − . The bottom right panel shows the continuum emission of the disk. Notethe CO arcs connecting the Keplerian wings at ∼ . (cid:48)(cid:48) and . (cid:48)(cid:48) from the center, most clearly seen on the North side of the diskat .
24 km s − and on the South side of the disk at .
20 km s − . These non-Keplerian velocity components resemble what weexpect from buoyancy spirals. difference suggests that we can distinguish the two sce-narios if we probe velocity perturbations at upper andlower surfaces of the disk separately. Optically thicktracers (e.g., CO), a sufficiently high spatial resolu-tion, and a moderately inclined disk geometry wouldprovide the best chance to separate the upper and lowersurface emission.
Gas pressure changes:
Gas pressure changes acrossradius can lead sub-/super-Keplerian rotation of thegas to maintain the radial force balance (Teague et al.2018). As we discussed in Section 4.3, distinguishingnon-Keplerian motions associated with buoyancy spiralsfrom rotation modulation can become a challenge in in-clined disks as we become sensitive to both vertical andazimuthal velocities. As an example, we present chan-nel maps of CO line emission of HD 143006 in Figure19. Morphologically, the arcs connecting the Keplerianwings (most clearly seen on the North side of the diskat v = 7 .
24 km s − and on the South side of the diskat v = 8 .
20 km s − ) show a good resemblance to thosefeatures expected from buoyancy resonance. However,these arcs can instead be interpreted as rotation mod-ulation. The non-Keplerian features are most promi- nently seen as red-shifted arcs in blue-shifted channels(e.g., v = 7 .
24 km s − ) and blue-shifted arcs in red-shifted channels (e.g., v = 8 .
20 km s − ), which are con-sistent with the expected modulation associated withsub-Keplerian rotation (Teague et al. 2019). In fact, theinner arcs at ∼ . (cid:48)(cid:48) coincide with the outermost con-tinuum ring, suggesting that the arcs could arise fromthe rotation modulation around the pressure peak. Theouter arcs at ∼ . (cid:48)(cid:48) is close to where CO emissionfades, suggesting they could be due to a rapid drop inthe gas density at that radius. Because of the degener-acy between vertical and azimuthal velocities in channelmaps, we recommend to use Keplerian-subtracted mo-ment maps rather than velocity channel maps when itcomes to distinguishing buoyancy spirals and rotationmodulation arising from gas pressure changes (see Sec-tion 4.3).5.4.
Buoyancy Resonances and Dust
How would buoyancy resonances appear in scatteredlight and (sub-)millimeter continuum observations? In-terestingly enough, van Boekel et al. (2017) reported atightly wound spiral in the near-infrared polarized in-tensity map of the TW Hya disk. The spiral is located bservational Signatures of Buoyancy Resonances Figure 20.
Cartoon channel maps summarizing main features expected from buoyancy spirals: (left) spurs around the centralvelocity channel, standing out of the Keplerian wings; (middle) an arc connecting Keplerian wings across the semi-major axisof the disk; and (right) an arc beyond the inner disk across the semi-major axis of the disk. at the outer edge of an annular gap centered at about93 au , and extends about 90 degrees in azimuth on theSouth-West side of the disk. While we opt out of makingsimulated scattered light observations from our modelsbecause our simulations do not include dust grains, itis interesting to point out that a buoyancy spiral pro-duces positive density perturbations at the exact loca-tion where the scattered light spiral is revealed (see thelower-right quadrant of δρ/ (cid:104) ρ (cid:105) φ at Z/R = 0 . and . in Figure 6). It is interesting to speculate that the veloc-ity spiral in the CO observation and the density spiralin the scattered light observation are probing buoyancyresonances driven by a planet embedded within the gapat ∼ au. Future simulations including dust particleswill help further investigate this possibility.On the other hand, given that buoyancy resonancesare weak or absent near the midplane and that buoyancyspirals are confined in the corotating region where large(sub-)millimeter-sized grains are expected to depleteddue to radial drift, we believe observing buoyancy res-onances in (sub-)millimeter continuum observations isless likely. SUMMARY AND CONCLUSIONAlong with the spirals driven by the well-known Lind-blad resonance, we showed that planets can excite spi-rals via buoyancy resonances, which we can detect usingmolecular line observations. We summarize our findingsbelow.(1) Under a broad range of conditions applicable toprotoplanetary disks, we showed that infrequent gas- This number is updated from van Boekel et al. (2017) in ac-cordance with the
Gaia distance of 60.1 pc (Bailer-Jones et al.2018). dust collision can be the bottleneck in the energy ex-change between the gas and dust in the surface layers(
Z/R (cid:38) . ; Figures 2 and 15). Although this hasbeen previously suggested (Malygin et al. 2017; Bar-ranco et al. 2018; Pfeil & Klahr 2019), to our knowledge,it is the first time that this effect is taken into accountin planet-disk interaction simulations.(2) The collision-limited slow thermal relaxation pro-vides favorable conditions for buoyancy resonances todevelop (Figures 3 and 21). Adopting the thermal re-laxation timescale estimated by considering radiation,diffusion, and gas-dust collision, we showed that planetscan excite a family of tightly wound spirals via buoyancyresonances, in addition to those excited by Lindblad res-onance (Figures 4 and 6).(3) Two main characteristics of buoyancy spirals aretheir small pitch angles and large vertical motions.Buoyancy spirals have a pitch angle of a few to 10 de-grees in the corotating region of the planet (Figure 5).The vertical motions associated with buoyancy reso-nances is of order of
100 m s − for Jovian-mass plan-ets, corresponding to about
20 % of the sound speedor a few % of the Keplerian speed. This is comparableto or larger than the velocity perturbations driven byLindblad resonance (Figures 6 and 7).(4) By generating synthetic ALMA observations,we showed that the non-Keplerian motions associatedwith buoyancy resonances is detectable. Buoyancy spi-rals would appear as tightly-wound arcs in Keplerian-subtracted moment maps (Figure 11). In velocity chan-nel maps, buoyancy spirals appear as spurs aroundthe central velocity channel, arcs connecting Keplerianwings or an arc beyond the inner disk across the semi-major axis of the disk. We summarize these featureswith a cartoon in Figure 20.6 Bae et al. (5) Because buoyancy resonances predominantly pro-duce vertical velocity perturbations, face-on disks pro-vide the best opportunities to search for their signatures.Based on the morphology and the magnitude of velocityperturbations, we propose that the tightly wound spi-rals seen in the near-face-on TW Hya disk (Teague et al.2019) could be driven by a (sub-)Jovian-mass planet at ∼ au.(6) Along with the implications for buoyancy reso-nances, the finite relaxation timescale has importantimplications for planet-induced gaps and. As shown inMiranda & Rafikov (2020b) and Zhang & Zhu (2020),when the cooling of the disk gas is moderate (i.e., β ≡ πt relax /t orb ∼ ) the gap around the planet be-comes narrower than that in fully isothermal ( β (cid:28) )or fully adiabatic ( β (cid:29) ) simulations. The finite re-laxation timescale implies that the mass of planets re-sponsible for gaps seen in continuum observations canbe underestimated if inferred based on isothermal sim-ulations.(7) We discussed potential ways to distinguish non-Keplerian motions driven by buoyancy resonances andthose driven by other mechanisms: Lindblad resonance,corrugated vertical flows, and gas pressure changes. Werecommend the community to observe multiple linestracing different heights in the disk. It is also crucialto have sufficiently high spatial/spectral resolution andsensitivity to separate the emission arising from the nearand far sides of the disk.We conclude by emphasizing that numerical simula-tions have to include more realistic and complete treat-ments for thermodynamics to fully capture planet-diskinteraction. Planet-disk interaction simulations often(but not always) adopt a vertically isothermal tempera-ture structure and/or an isothermal equation of state.When it comes to buoyancy resonances, such simpli- fied models can completely suppress the resonance. Weshould point out that our simulations have caveats. Weadopted a prescribed, fixed thermal relaxation model.In reality, the spatial and size distribution of dust grainswould evolve over time, and this is neglected in currentsimulations. It will be also interesting to implement athermo-chemistry model that evolves over time, coupledwith the hydro evolution. Future simulations with morecomplete treatments for dust- and thermo-dynamics willhelp better interpret state-of-the-art observations.We thank the anonymous referee for providing us ahelpful report that improved the initial manuscript. JBis thankful to Myriam Benisty, Stefano Facchini, andSteve Lubow for helpful conversations. JB acknowl-edges support by NASA through the NASA HubbleFellowship grant Software:
FARGO3D (Benítez-Llambay & Masset2016),
RADMC-3D (Dullemond et al. 2012), eddy (Teague2019)APPENDIX A. THE RELAXATION TIMESCALEIn Figure 21 we present t relax N Z for the disk models discussed in Section 5.1. For the broad range of parameterspace we explored, t relax (cid:38) N − Z in Z/R (cid:38) . , suggesting that buoyancy resonances likely develop in the surface layersof protoplanetary disks. B. ADDITIONAL CHANNEL MAPSFigures 22 and 23 present channel maps from synthetic CO observations of models with . M Jup and M Jup planets. Simulated cubes are publicly available at https://doi.org/10.5281/zenodo.4361639.REFERENCES
ALMA Partnership, Brogan, C. L., Pérez, L. M., et al.2015, ApJL, 808, L3, doi: 10.1088/2041-8205/808/1/L3 Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012,ApJ, 744, 162, doi: 10.1088/0004-637X/744/2/162 bservational Signatures of Buoyancy Resonances standard model Z / R λ diff = 0.1H g λ diff = 3H g α = 10 -2 Z / R α = 10 -4 × a max × a max Z / R n d (a) ∝ a -2.5 n d (a) ∝ a -4.5 M disk = 0.017 M Sun Z / R M disk = 0.15 M Sun with a gap, M p = 1M Jup -2 -1 0 1 2 3 4 5 6 log t relax N Z Figure 21.
Similar to Figure 15, but showing t relax N Z .Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL,820, L40, doi: 10.3847/2041-8205/820/2/L40Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL,869, L41, doi: 10.3847/2041-8213/aaf741 Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ,863, 44, doi: 10.3847/1538-4357/aab846Bae, J., Nelson, R. P., & Hartmann, L. 2016a, ApJ, 833,126, doi: 10.3847/1538-4357/833/2/126 Bae et al. − − O ff s e t [ a r c s ec ] − − O ff s e t [ a r c s ec ] − − − − O ff s e t [ a r c s ec ] − − − − − − − − I n t e n s i t y [ m J y b e a m − ] Figure 22.
Same as Figure 14, but with M p = 0 . M Jup . − − O ff s e t [ a r c s ec ] − − O ff s e t [ a r c s ec ] − − − − O ff s e t [ a r c s ec ] − − − − − − − − I n t e n s i t y [ m J y b e a m − ] Figure 23.