Observational Signature of Tightly Wound Spirals Driven by Buoyancy Resonances in Protoplanetary Disks

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 vertically displaced from the equilibrium position. Its response to the perturbation can be described with the vertical buoyancy frequency (a.k.a. Brunt–Väisälä frequency):

$$\begin{eqnarray}&&{N}_{Z}^{2}=\displaystyle \frac{g}{\gamma }\displaystyle \frac{\partial }{\partial Z}\left[\mathrm{ln}\left(\displaystyle \frac{P}{{\rho }^{\gamma }}\right)\right]\end{eqnarray} \tag{ 1 }$$

or

$$\begin{eqnarray}&&{N}_{Z}^{2}=\displaystyle \frac{g}{\gamma }\displaystyle \frac{\partial }{\partial Z}\left[\mathrm{ln}\left(\displaystyle \frac{T}{{\rho }^{\gamma -1}}\right)\right]\end{eqnarray} \tag{ 2 }$$

adopting the ideal gas law. In the above equations, g is the gravitational acceleration, γ is the adiabatic index, P is the gas pressure, ρ is the gas density, and T is the gas temperature. When the disk temperature is dominated by the stellar irradiation, the disk is hotter near 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 vertical temperature gradient and γ ≥ 1, Equation (2) yields ${N}_{Z}^{2}\gt 0$. We thus expect the gas parcel to vertically oscillate around its equilibrium position with a frequency N_Z .

When the buoyancy frequency matches with the forcing frequency which, in this case, is the planet's orbital frequency, buoyancy resonances can develop (Zhu et al. 2012; Lubow & Zhu 2014). Buoyancy resonances give rise to the density and velocity perturbations along a family of trailing spirals. As we will show below, one of the main characteristics of buoyancy spirals is that they are very tightly wound compared with Lindblad spirals, in particular in the vicinity of the planet.

It is worth pointing out that the thermodynamic properties of the disk are important in the development of buoyancy resonances. In order for buoyancy resonances to fully develop, the timescale for the gas to respond to thermal perturbations (hereafter relaxation timescale t_relax) has to be longer than the timescale associated with the buoyancy: ${t}_{\mathrm{relax}}\gtrsim {N}_{Z}^{-1}$. When ${t}_{\mathrm{relax}}\ll {N}_{Z}^{-1}$, the gas behaves isothermally, and buoyancy resonances are expected to be weak or absent.

In this paper, we show that planets can excite spirals through buoyancy resonances, which can be detectable using molecular line observations with ALMA. The paper is organized as follows. In Section 2, we describe three processes that determine the relaxation timescale of the gas—radiation, diffusion, and gas–dust collision—and compute the corresponding timescales using a TW Hya disk model. We show that thermal relaxation in the surface layers (Z/R ≳ 0.1) and outer regions (R ≳ 100 au) of protoplanetary disks can be limited by insufficient gas–dust collision. In Section 3, we present three-dimensional hydrodynamic simulations and show that the collision-limited slow thermal relaxation provides favorable conditions for buoyancy resonances. In Section 4, we generate synthetic observations and show that the spirals driven by buoyancy resonances are observable with molecular line observations using ALMA. Based on the tightly wound morphology and the magnitude of velocity perturbations, we propose that the velocity spiral seen in TW Hya could be driven by a (sub-)Jovian-mass planet at 90 au. In Section 5, we show that the relaxation timescale is comparable to or longer than the dynamical timescale under a broad range of conditions and discuss its implications for the development of buoyancy resonances, planet-induced gap profiles, and hydrodynamic instabilities. We also discuss potential ways to discriminate non-Keplerian motions driven by buoyancy resonances from those driven by other mechanisms, including Lindblad resonance, corrugated vertical flows, and gas pressure changes. We summarize our findings and conclude in Section 6.

As one of the motivations of this work is to explain the tightly wound spirals observed in the TW Hya disk (Teague et al. 2019), we adopt the disk density and temperature profiles similar to the one constrained for the disk in Huang et al. (2018a).

The gas surface density follows

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}(R)={{\rm{\Sigma }}}_{p}{\left(R/{R}_{p}\right)}^{-p},\end{eqnarray} \tag{ 3 }$$

where Σ_p is the surface density at R_p = 90 au and p = 0.9. We choose Σ_p = 3.8 g cm⁻² such that the total disk gas mass within 210 au is 0.05 M_⊙, broadly consistent with observational estimates for the TW Hya disk (Bergin et al. 2013). Throughout this paper, we use $R=r\sin \theta$ for the cylindrical radius and $Z=r\cos \theta$ for the height, where r and θ are the radius and the polar angle in the spherical coordinates, respectively.

Following the prescription in Dartois et al. (2003), the gas temperature is parameterized as follows:

$$\begin{eqnarray}{T}_{{\rm{g}}}(Z)=\left\{\begin{array}{ll}{T}_{\mathrm{atm}}+({T}_{\mathrm{mid}}-{T}_{\mathrm{atm}}){\cos }^{4}\left(\displaystyle \frac{\pi }{2}\displaystyle \frac{Z}{{Z}_{q}}\right) & \mathrm{if}\,Z\lt {Z}_{q}\\ {T}_{\mathrm{atm}} & \mathrm{if}\,Z\geqslant {Z}_{q}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

Here, the midplane and atmosphere temperatures are a function of the cylindrical radius R following

$$\begin{eqnarray}&&{T}_{\mathrm{atm}}(R)={T}_{\mathrm{atm},{\rm{p}}}{\left(R/{R}_{p}\right)}^{-q}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}(R)={T}_{\mathrm{mid},{\rm{p}}}{\left(R/{R}_{p}\right)}^{-q},\end{eqnarray} \tag{ 6 }$$

where T_atm,p = 44.5 K, T_mid,p = 14.2 K, and q = 0.47 (Zhang et al. 2017; Huang et al. 2018a). In Equation (4), Z_q (R) = 4H_mid(R), where H_mid is the gas scale height determined with the disk midplane temperature. In terms of the disk aspect ratio, the above temperature profile at the midplane corresponds to ${H}_{\mathrm{mid}}/R=0.075\times {(R/{R}_{p})}^{0.265}$, assuming 0.88 M_⊙ for the stellar mass (Andrews et al. 2012; Huang et al. 2018a) and 2.4 for the mean molecular weight of the gas.

Using the above surface density and temperature profiles, we construct the three-dimensional gas density distribution that satisfies the vertical hydrostatic equilibrium:

$$\begin{eqnarray}&&-\displaystyle \frac{{{GM}}_{* }Z}{{\left({R}^{2}+{Z}^{2}\right)}^{3/2}}-\displaystyle \frac{1}{{\rho }_{{\rm{g}}}}\displaystyle \frac{\partial P}{\partial Z}=0.\end{eqnarray} \tag{ 7 }$$

Solving the above equation results in the vertical density distribution that follows

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}(Z)={\rho }_{{\rm{g}},\mathrm{mid}}\displaystyle \frac{{c}_{s,\mathrm{mid}}^{2}}{{c}_{s}^{2}(Z)}\exp \left[-{\int }_{0}^{Z}\displaystyle \frac{1}{{c}_{s}^{2}(Z^{\prime} )}\displaystyle \frac{{{GM}}_{* }Z^{\prime} }{{\left({R}^{2}+Z{{\prime} }^{2}\right)}^{3/2}}{dZ}^{\prime} \right],\end{eqnarray} \tag{ 8 }$$

where c_s,mid and c_s (Z) denote the sound speed at the disk midplane and at height Z. We then compute the angular velocity Ω that satisfies the radial force balance, taking into account the gas pressure gradient:

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{2}={{\rm{\Omega }}}_{K}^{2}\sin \theta +\displaystyle \frac{1}{{\rho }_{{\rm{g}}}r{\sin }^{2}\theta }\displaystyle \frac{\partial P}{\partial r}.\end{eqnarray} \tag{ 9 }$$

Here, ${{\rm{\Omega }}}_{K}=\sqrt{{{GM}}_{* }/{R}^{3}}$ is the Keplerian angular velocity. The initial radial and meridional velocities are set to zero.

Protoplanetary disks are a mixture of gas and dust. Hydrogen molecules dominate the total mass and thermal energy, but they are inefficient at emitting radiation. Let us consider a situation where a gas parcel has been perturbed from its equilibrium state to have a higher temperature. In order for the gas to lose its thermal energy, hydrogen molecules first have to transfer their kinetic energy to the surrounding dust grains. Dust grains then radiate away the excess thermal energy. This process can thus be understood as a sequential, two-step process (see Figure 1 for a schematic diagram).

Zoom In Zoom Out Reset image size

When the collision between hydrogen molecules and dust grains is sufficiently frequent, the gas cools at the rate dust grains radiate away their thermal energy; in this case, gas and dust are thermally coupled. In the optically thin regime, thermal photons emitted by dust grains can freely escape from the disk. The relaxation timescale of the gas can be described by the radiative timescale of dust grains,

$$\begin{eqnarray}&&{t}_{\mathrm{rad}}=\displaystyle \frac{{c}_{V}}{16{\kappa }_{{\rm{P}}}{\sigma }_{\mathrm{SB}}{T}_{{\rm{g}}}^{3}},\end{eqnarray} \tag{ 10 }$$

where c_V denotes the specific heat capacity of the gas, κ_P is the Planck mean opacity of the dust, and σ_SB is the Stefan–Boltzmann constant.

In the optically thick regime, thermal photons emitted from dust grains are absorbed and emitted by other grains multiple times before they eventually escape the disk. In this case, the cooling timescale can be characterized by the diffusion of photons. The diffusion timescale associated with the length scale λ_diff can be written as

$$\begin{eqnarray}&&{t}_{\mathrm{diff}}=\displaystyle \frac{{\lambda }_{\mathrm{diff}}^{2}}{\langle D\rangle },\end{eqnarray} \tag{ 11 }$$

where D is the diffusion coefficient defined as

$$\begin{eqnarray}&&D=\displaystyle \frac{16{\sigma }_{\mathrm{SB}}{T}_{{\rm{g}}}^{3}}{3{c}_{V}{\kappa }_{{\rm{R}}}{\rho }_{{\rm{g}}}^{2}},\end{eqnarray} \tag{ 12 }$$

with κ_R being the Rosseland mean opacity of the dust and the brackets 〈〉 denoting an average over λ_diff (Malygin et al. 2017). While the density, temperature, and opacity are expected to vary moderately along the radial and azimuthal directions, their variation can be larger along the vertical direction due to the vertical stratification of the disk. We thus define the in-plane diffusion timescale t_diff,∥ and vertical diffusion timescale t_diff,⊥, for which the diffusion coefficient is averaged along the radial and vertical directions, respectively. The overall diffusion timescale considering both in-plane and vertical diffusion can be estimated as

$$\begin{eqnarray}&&{t}_{\mathrm{diff}}={\left(\displaystyle \frac{1}{{t}_{\mathrm{diff},\parallel }}+\displaystyle \frac{1}{{t}_{\mathrm{diff},\perp }}\right)}^{-1}.\end{eqnarray} \tag{ 13 }$$

For a given length scale λ_diff, we find that t_diff,∥ and t_diff,⊥ are comparable when λ_diff ≪ H_g, not surprisingly, where H_g is the gas scale height. However, one can be larger than the other by a factor of a few when λ_diff ≳ H_g. It is also worth mentioning that it is possible that the length scale of the temperature perturbation along the radial/azimuthal direction and that along the vertical direction differ (see e.g., Miranda & Rafikov 2020a).

We note that, with the assumption that the disk has a constant temperature vertically, along with τ_R ≡ ∫κ_R ρ_g dZ and Σ_g ≃ ρ_g H_g, combining Equations (11) and (12) approximates to

$$\begin{eqnarray}&&{t}_{\mathrm{diff},\perp }\simeq \displaystyle \frac{3{c}_{V}{{\rm{\Sigma }}}_{{\rm{g}}}}{16{\sigma }_{\mathrm{SB}}{T}_{{\rm{g}}}^{3}}{\tau }_{{\rm{R}}},\end{eqnarray} \tag{ 14 }$$

a formula that is often adopted to account for optically thick cooling in a 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 collision rate (assuming lack of other cooling mechanisms, such as molecular/atomic line cooling; see below). Because the assumption of thermal equilibrium is no longer valid, we define the dust temperature T_d. The gas–dust collisional timescale t_coll can be written as

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}=\displaystyle \frac{{\rho }_{{\rm{g}}}{c}_{V}}{{{\rm{\Lambda }}}_{\mathrm{coll}}}| {T}_{{\rm{g}}}-{T}_{{\rm{d}}}| ,\end{eqnarray} \tag{ 15 }$$

where Λ_coll is the cooling rate per unit volume via gas–dust collisions for which we follow Burke & Hollenbach (1983):

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{coll}}={\int }_{{a}_{\min }}^{{a}_{\max }}2{k}_{B}| {T}_{{\rm{g}}}-{T}_{{\rm{d}}}| \bar{{\alpha }_{T}}{n}_{{\rm{g}}}{n}_{{\rm{d}}}(a){\sigma }_{{\rm{d}}}(a){v}_{\mathrm{th}}{da}.\end{eqnarray} \tag{ 16 }$$

In the above equation, ${a}_{\min }$ and ${a}_{\max }$ denote the minimum and maximum dust grain sizes, k_B is the Boltzmann constant, $\bar{{\alpha }_{T}}=0.5$ is the thermal accommodation coefficient that characterizes the efficiency of the heat transfer between gas molecules and dust grains, n_g is the number density of gas molecules, n_d(a) is the number density of dust particles with size a, σ_d = π a² is the geometrical cross section of dust particles, and ${v}_{\mathrm{th}}={(8{k}_{B}{T}_{{\rm{g}}}/\pi {m}_{g})}^{1/2}$ is the thermal velocity of the gas. Using the second and third moments of the dust grain size, which are defined as

$$\begin{eqnarray}&&\langle {a}^{2}\rangle ={\int }_{{a}_{\min }}^{{a}_{\max }}{a}^{2}{n}_{{\rm{d}}}(a){da}\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&\langle {a}^{3}\rangle ={\int }_{{a}_{\min }}^{{a}_{\max }}{a}^{3}{n}_{{\rm{d}}}(a){da},\end{eqnarray} \tag{ 18 }$$

Equation (16) can be written as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{coll}}=2{k}_{B}| {T}_{{\rm{g}}}-{T}_{{\rm{d}}}| \bar{{\alpha }_{T}}\left(\displaystyle \frac{3}{4{\rho }_{s}}\right)\left(\displaystyle \frac{\langle {a}^{2}\rangle }{\langle {a}^{3}\rangle }\right)\left(\displaystyle \frac{{\rho }_{{\rm{g}}}^{2}{v}_{\mathrm{th}}}{{m}_{{\rm{g}}}}\right)\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right),\end{eqnarray} \tag{ 19 }$$

where ρ_s is the bulk density of dust grains while ρ_d is the mass density of dust grains. Then, combining Equations (15) and (19), the collisional timescale can be written as

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}=\displaystyle \frac{{c}_{V}}{2{k}_{B}\bar{{\alpha }_{T}}}\left(\displaystyle \frac{4}{3}{\rho }_{s}\right)\left(\displaystyle \frac{\langle {a}^{3}\rangle }{\langle {a}^{2}\rangle }\right)\left(\displaystyle \frac{{m}_{{\rm{g}}}}{{\rho }_{{\rm{g}}}{v}_{\mathrm{th}}}\right){\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right)}^{-1}.\end{eqnarray} \tag{ 20 }$$

Note that the collisional timescale depends on the mean dust grain size and the local dust-to-gas mass ratio, which are dependent upon the grain size distribution and the level of turbulence of the disk among many others.

Then, taking into account radiation, diffusion, and gas–dust collisional timescales, the relaxation timescale of the gas can be written as

$$\begin{eqnarray}&&{t}_{\mathrm{relax}}={t}_{\mathrm{coll}}+\max ({t}_{\mathrm{rad}},{t}_{\mathrm{diff}}).\end{eqnarray} \tag{ 21 }$$

Again, this equation shows that the collision between gas and 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 of the cooling process, and thus, the overall cooling timescale is determined by the collisional timescale (i.e., t_relax ≃ t_coll). When the gas–dust collision is frequent, the gas and dust are in thermal equilibrium and the gas cools over the timescale t_rad or t_diff depending on the optical thickness of the disk (i.e., ${t}_{\mathrm{relax}}\simeq \max ({t}_{\mathrm{rad}},{t}_{\mathrm{diff}})$). The transition between the optically thick and thin regimes occur when t_diff ≃ 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 using the disk model described in Section 2.1. To do so, we first need to define the spatial/size distribution of dust grains as it dictates the gas–dust collision rate but also the radiative/diffusion timescales through the opacity. We adopt a maximum grain size that is decreasing over radius: ${a}_{\max }{(R)=1\mathrm{mm}(R/30\mathrm{au})}^{-2}$. This choice is motivated by the fact that the (sub)millimeter continuum emission of the TW Hya disk is confined within about 60 au (e.g., Andrews et al. 2016; Tsukagoshi et al. 2016). With the vertically integrated total dust-to-gas mass ratio fixed to 0.01 at each radius (i.e., Σ_d/Σ_g = 0.01), we distribute the dust mass between 0.1 μm and ${a}_{\max }$ adopting a power-law dust size distribution with a power-law index −3.5: n_d(a) ∝ a^−3.5. We then determine the vertical scale height of each dust species assuming that vertical settling is balanced by turbulence mixing characterized by α = 10⁻³. This results in a dust scale height of

$$\begin{eqnarray}&&\begin{array}{l}{H}_{{\rm{d}}}(R,a)={H}_{{\rm{g}}}(R)\\ \ \times \ \min \left(1,\sqrt{\displaystyle \frac{\alpha }{\min (\mathrm{St}(a),1/2)(1+\mathrm{St}{\left(a\right)}^{2})}}\right)\end{array}\end{eqnarray} \tag{ 22 }$$

(Dullemond & Dominik 2004; Birnstiel et al. 2010), where St (a) ≡ (ρ_s a/ρ_g v_th)Ω_K is the Stokes number of particle having size a, and we assume a grain internal density ρ_s = 1.67 g cm⁻³ (see below). Then, the vertical density distribution of each dust species is obtained following

$$\begin{eqnarray}&&{\rho }_{d}(R,Z,a)=\displaystyle \frac{{{\rm{\Sigma }}}_{d}(R,a)}{\sqrt{2\pi }{H}_{d}(R,a)}\exp \left(-\displaystyle \frac{{Z}^{2}}{2{H}_{d}{\left(R,a\right)}^{2}}\right).\end{eqnarray} \tag{ 23 }$$

The resulting mean dust particle size 〈a³〉/〈a²〉 and dust-to-gas mass ratio ρ_d /ρ_g are shown in Figures 2(a) and (b). Larger grains settle near the midplane due to shorter settling times. Because larger grains contain a larger fraction of the total dust mass, the dust-to-gas mass ratio decreases over height. Dust grains with sizes ≳100 μm, which are believed to dominate the (sub)millimeter continuum emission, are confined in radius within the inner ∼60 au, while micron-sized grains extend much farther out, consistent with both (sub)millimeter and optical/near-infrared observations of TW Hya (e.g., Debes et al. 2013, 2017; Andrews et al. 2016; Tsukagoshi et al. 2016; van Boekel et al. 2017).

Zoom In Zoom Out Reset image size

Next, we calculate the dust opacity in each grid cell based on the dust distribution obtained as above, adopting the opacity model from the DSHARP collaboration (Birnstiel et al. 2018). In this opacity model, grains are assumed to be a mixture of water ice, astronomical silicates, troilite, and refractory organic material, having a bulk density of 1.67 g cm⁻³. The optical constants to compute the DSHARP opacity are originally from Henning & Stognienko (1996), Draine (2003), and Warren & Brandt (2008). The absorption and scattering opacities are calculated following

$$\begin{eqnarray}&&{\kappa }_{\nu }^{\mathrm{abs},\mathrm{sca}}(R,Z)=\displaystyle \frac{{\displaystyle \sum }_{i}{\kappa }_{\nu }^{\mathrm{abs},\mathrm{sca}}({a}_{i}){\rho }_{{\rm{d}}}(R,Z,{a}_{i})}{{\displaystyle \sum }_{i}{\rho }_{{\rm{d}}}(R,Z,{a}_{i})},\end{eqnarray} \tag{ 24 }$$

where the subscript ν shows the frequency dependence of the opacity. The Rosseland and Planck mean opacities are calculated as

$$\begin{eqnarray}&&{\kappa }_{{\rm{R}}}(R,Z)={\left[\displaystyle \frac{{\int }_{0}^{\infty }\tfrac{1}{{\kappa }_{\nu }^{\mathrm{ext}}(R,Z)}\tfrac{{{dB}}_{\nu }}{{dT}}d\nu }{{\int }_{0}^{\infty }\tfrac{{{dB}}_{\nu }}{{dT}}d\nu }\right]}^{-1}\end{eqnarray} \tag{ 25 }$$

and

$$\begin{eqnarray}&&{\kappa }_{{\rm{P}}}(R,Z)=\displaystyle \frac{{\int }_{0}^{\infty }{\kappa }_{\nu }^{\mathrm{abs}}(R,Z){B}_{\nu }d\nu }{{\int }_{0}^{\infty }{B}_{\nu }d\nu },\end{eqnarray} \tag{ 26 }$$

where ${\kappa }_{\nu }^{\mathrm{ext}}={\kappa }_{\nu }^{\mathrm{abs}}+{\kappa }_{\nu }^{\mathrm{sca}}$.

The radiation timescale t_rad calculated as in Equation (10) is shown in Figure 2(c). The radiation timescale is orders of magnitude shorter than the dynamical timescale everywhere in the disk. At a given radius, the radiation timescale decreases over height because of its steep temperature dependency (${t}_{\mathrm{rad}}\propto {T}_{{\rm{g}}}^{-3}$).

Figure 2(d) presents the diffusion timescale t_diff. Diffusion is slower when the optical depth is larger. So, t_diff is greater at smaller radii and near the midplane, and is a decreasing function of R and Z. In particular, note that t_diff drops exponentially over height because of the density dependency (${t}_{\mathrm{diff}}\propto {\rho }_{{\rm{g}}}^{2}$). We also note that the diffusion timescale depends on the temperature perturbation length scale as ${t}_{\mathrm{diff}}\propto {\lambda }_{\mathrm{diff}}^{2}$. Here, we opt to use λ_diff = H_g. In reality, the length scale of any perturbations can range from ≪ H_g to the thickness of the disk, which is a few scale heights. However, as we will show below (see also Section 5.1), the thermal relaxation in the surface layers (Z/R ≳ 0.1) is limited by infrequent gas–dust collision, insensitive to the choice of λ_diff.

Figure 2(e) shows the gas–dust collisional timescale t_coll. For given gas density and temperature structures, the collisional timescale is set by the mean grain size and the dust-to-gas mass ratio (${t}_{\mathrm{coll}}\propto \langle {a}^{3}\rangle /\langle {a}^{2}\rangle ,{({\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}})}^{-1}$). As the dust-to-gas mass ratio drops exponentially over height, gas molecules have significantly less frequent collisions with dust grains. This makes the collisional timescale orders of magnitude longer than the dynamical timescale in the surface layers.

The relaxation timescale t_relax, calculated as in Equation (21), is shown in Figure 2(f). The plot clearly shows that the assumption of thermal equilibrium between gas and dust is not necessarily valid in the outer and surface regions of the disk because of the long collisional timescale there. Note that this picture is in good agreement with previous studies of thermal relaxation in protoplanetary disks (e.g., Malygin et al. 2017; Barranco et al. 2018; Pfeil & Klahr 2019). In Section 5.1, we explore how different assumptions on the diffusion length scale, level of disk turbulence, grain size distribution, disk mass, and the existence of a gap in the disk can affect the cooling timescales. As we will show, the thermal relaxation of the gas in surface layers of protoplanetary disks is limited by infrequent gas–dust collision over a broad range parameter space.

As we mentioned in Section 1, buoyancy resonances require adiabatic responses to thermal perturbations to fully develop (i.e., ${t}_{\mathrm{relax}}\gtrsim {N}_{Z}^{-1}$). To see how t_relax compares with ${N}_{Z}^{-1}$, we present t_relax N_Z in Figure 3. As shown, when gas–dust collision is taken into account, ${t}_{\mathrm{relax}}\gtrsim {N}_{Z}^{-1}$ in the entire disk except at the midplane where N_Z is zero due to the symmetry across the midplane. This suggests that most parts of the disk, including the surface layers the CO lines probe, has favorable conditions for buoyancy resonances to develop. In contrast, if gas–dust collision is ignored, the surface layers have ${t}_{\mathrm{relax}}\ll {N}_{Z}^{-1}$, suggesting that buoyancy resonances are weak or unlikely to develop there.

Zoom In Zoom Out Reset image size

In addition to the cooling processes we considered above, atomic and molecular line cooling plays a role in the surface layers of protoplanetary disks (e.g., Gorti et al. 2011; Du & Bergin 2014; Kama et al. 2016; Facchini et al. 2018). The exact height beyond which line cooling dominates depends on the underlying thermal/chemical properties of the disk as well as the stellar/external irradiation. While including comprehensive thermo/photochemistry requires full thermochemical radiative transfer calculations, which is beyond the scope of the paper, here we test the potential effect line cooling would have on buoyancy resonances by assuming that the disk gas cools efficiently beyond a certain height of the disk Z = Z_line. Taking this into account, we adopt the following form for the relaxation timescale considering diffusion, radiation, gas–dust collision, and line cooling:

$$\begin{eqnarray}&&{t}_{\mathrm{relax}}=\left[{t}_{\mathrm{coll}}+\max ({t}_{\mathrm{rad}},{t}_{\mathrm{diff}})\right]\exp \left[-{\left(\displaystyle \frac{Z}{{Z}_{\mathrm{line}}}\right)}^{-12}\right].\end{eqnarray} \tag{ 27 }$$

The exponential term on the right-hand side of the equation is added to mimic the effect of line cooling. ⁵

In our fiducial models, we adopt Z_line = 4H_g, which corresponds to Z/R = 0.3 at the radial location of the planet in our simulations, 90 au. This choice is motivated by the thermochemical model of TW Hya presented in Kama et al. (2016), where the major atomic line emission, including [C i], [C ii], and [O i] lines, originates from Z/R ∼ 0.3–0.4 (see their Figure D2.). However, it is important to note that the exact line cooling rate is dependent upon various factors, including the gas-phase abundance of the coolants and electron number density, which is determined by UV flux as well as full chemical chains. To test the effect of line cooling, we ran additional simulations adopting Z_line = 2H_g (Z/R = 0.15 at 90 au; Section 5.2). In this model, the ¹²CO molecular line probes the layers where cooling is rapid, dominated by line cooling.

For the sake of reproducing the hydrodynamic simulations we will present in the following sections, we provide parameterized fits to the diffusion and collisional timescales:

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\mathrm{diff}}}{{t}_{\mathrm{orb}}}=2.4{\left(\displaystyle \frac{R}{90\mathrm{au}}\right)}^{-3.9}\exp \left[-{\left(\displaystyle \frac{Z}{{Z}_{\mathrm{diff}}}\right)}^{0.9}\right],\end{eqnarray} \tag{ 28 }$$

where Z_diff = 2.0 au (R/90 au)^1.4, and

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\mathrm{coll}}}{{t}_{\mathrm{orb}}}=0.4{\left(\displaystyle \frac{R}{90\mathrm{au}}\right)}^{-1.1}\exp \left[{\left(\displaystyle \frac{Z}{{Z}_{\mathrm{coll}}}\right)}^{2.2}\right],\end{eqnarray} \tag{ 29 }$$

where Z_coll = 11.7 au (R/90 au)^1.2.

We solve the hydrodynamic equations for mass, momentum, and energy conservation in the three-dimensional spherical coordinates (r, θ, ϕ) using FARGO3D (Masset 2000; Benítez-Llambay & Masset 2016):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{g}}}v)=0,\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}\left(\displaystyle \frac{\partial v}{\partial t}+v\cdot {\rm{\nabla }}v\right)=-{\rm{\nabla }}P-{\rho }_{{\rm{g}}}{\rm{\nabla }}({{\rm{\Phi }}}_{* }+{{\rm{\Phi }}}_{p})+{\rm{\nabla }}\cdot {\rm{\Pi }},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\rm{\nabla }}\cdot ({ev})=-P{\rm{\nabla }}\cdot v+{Q}_{\mathrm{relax}}.\end{eqnarray} \tag{ 32 }$$

In the above equations, ρ_g is the gas density, v is the velocity vector, P is the gas pressure, Φ_* = −GM_*/r is the gravitational potential of the central star having mass M_*, Φ_p is the gravitational potential of the planet, Π is the viscous stress tensor, e is the internal energy per unit volume, and Q_relax is the rate at which the gas thermally relaxes to the initial state (see Section 2.2). Note that Q_relax can be either a positive or a negative value depending on whether the gas is colder or hotter than the initial equilibrium temperature.

The gravitational potential of the planet is computed as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}(r,\theta ,\phi )=-\displaystyle \frac{{{GM}}_{p}}{{\left(| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{p}{| }^{2}+{s}^{2}\right)}^{1/2}},\end{eqnarray} \tag{ 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 cell in question and of the planet, and s is the smoothing length. Because the smoothing length in three-dimensional calculations is used only to avoid the singularity in the potential on the grid scale, we adopt the cell diagonal size at the position of the planet for the smoothing length. We insert the planet at R_p = 90 au with a fixed, circular orbit. We use three planet masses: M_p = 0.5, 1, and 2 M_Jup. We run simulations for 500 planetary orbits. The planet mass is linearly increased over the first five planetary orbits.

For the thermal evolution, we adopt an adiabatic equation of state with an adiabatic index γ = 1.4. The gas pressure and the internal energy are related as P = (γ − 1)e. In addition to the thermal energy evolution via PdV work, which is accounted for by the first term of the right-hand side of Equation (32), cooling/heating of the gas is realized through the relaxation of the temperature T_g toward the initial temperature T_g,init (described in Equations (4)–(6)) over the thermal relaxation timescale t_relax computed in Section 2. The thermal relaxation rate can be written as

$$\begin{eqnarray}&&{Q}_{\mathrm{relax}}=-{\rho }_{{\rm{g}}}{c}_{V}\displaystyle \frac{{T}_{{\rm{g}}}-{T}_{{\rm{g}},\mathrm{init}}}{{t}_{\mathrm{relax}}}.\end{eqnarray} \tag{ 34 }$$

In practice, we use $\max ({t}_{\mathrm{relax}},{\rm{\Delta }}{t}_{\mathrm{hydro}})$ in the denominator of Equation (34), where Δt_hydro is the hydrodynamic time step, in order to avoid the overrelaxation that can happen when t_relax < Δt_hydro.

For our standard model, we take into account radiation, diffusion, gas–dust collision, and line cooling (Section 3.3.1). In practice, this is done by adopting a prescribed relaxation timescale using Equation (27), along with the fits in Equations (28) and (29). As a comparison to the standard model, we additionally carry out simulations without gas–dust collision: i.e., ${t}_{\mathrm{relax}}\,=\max {({t}_{\mathrm{rad}},{t}_{\mathrm{diff}})\exp [-(Z/{Z}_{\mathrm{line}})}^{-12}]$ (Section 3.3.2; hereafter the gas–dust thermal equilibrium model). This model is similar to what is often used in three-dimensional protoplanetary disk simulations where it is implicitly assumed that the gas and dust have instantaneous energy balance.

The simulation domain extends from 30 au (=0.33 R_p ) to 210 au (=2.33 R_p ) in r, from π/2 − 0.4 to π/2 in θ, which covers 5.6 scale heights at the radial location of the planet, and from 0 to 2π in ϕ. We adopt 460 logarithmically spaced grid cells in the radial direction, 96 uniformly spaced grid cells in the meridional direction, and 1482 uniformly spaced grid cells in the azimuthal direction. With this choice, one gas scale height at the location of the planet is resolved with about 18 grid cells in all directions.

At the radial boundaries, we adopt a wave-damping zone to suppress wave reflection (de Val-Borro et al. 2006). At the lower meridional boundary, which is the disk midplane, we adopt the symmetric boundary condition for all variables but the meridional velocity for 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 turbulence observationally constrained for the TW Hya disk (Teague et al. 2016; Flaherty et al. 2018).

In order to ensure that no other hydrodynamic instabilities operate in the disk, we ran both the standard model and gas–dust thermal equilibrium model in the absence of a planet. We found that no hydrodynamic instabilities develop in these runs. The vertical shear instability (Urpin & Brandenburg 1998; Nelson et al. 2013) is suppressed due to the nonzero viscosity and a long thermal relaxation time, in agreement with previous studies (Nelson et al. 2013).

3.3.1. Standard Model

We start by discussing the results from our standard model where we consider radiation, diffusion, gas–dust collision, and line cooling for the thermal relaxation of the disk gas. Figure 4 shows the perturbed density, perturbed temperature, and vertical velocity in a ϕ−r plane at Z/R = 0.23 for the M_p = 0.5 M_Jup model. This height corresponds to about three scale heights above the midplane at the radial location of the planet and is close to the ¹²CO emission surface in the synthetic observation we will present in Section 4.

Zoom In Zoom Out Reset image size

As most clearly shown in the perturbed density distribution, the planet excites a pair of primary Lindblad spirals, one in the inner disk (the arc crossing the r = 40 au boundary at ϕ ∼ 3π/4) and one in the outer disk (the arc crossing the r = 140 au boundary at ϕ ∼ − π/4). The Lindblad spirals are nearly perpendicular to the azimuth axis near the planet, suggesting that they have a large pitch angle close to 90° there. In addition to the primary Lindblad spirals, the planet excites a secondary Lindblad spiral in the inner disk, which emerges at r ∼ 50 au and ϕ ∼ π. At Z/R = 0.23, density, temperature, vertical velocity perturbations created by the Lindblad spirals of the 0.5 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 a family of spirals via buoyancy resonances, which can be most clearly seen in the vertical velocity plot. Two main features that distinguish buoyancy spirals from Lindblad spirals are (1) the tightly wound morphology and (2) the large vertical motions, which we will explain in detail one by one.

In Figure 5, we show the measured pitch angle of Lindblad and buoyancy spirals with a 5 au interval in radius. For Lindblad spirals, we measure the pitch angle using the peak in the density perturbation. For buoyancy spirals, we measure the pitch angle using the peak in positive vertical velocities. We opt to use the vertical velocity instead of the density because (1) buoyancy spirals are more clearly identifiable with the vertical velocity and (2) we can make a more direct comparison to the TW Hya spiral detected in the velocity space (Section 4.1.1).

Zoom In Zoom Out Reset image size

We also present the pitch angle derived with linear theory. For Lindblad spirals, the phase angle ϕ_L (i.e., the azimuthal angle from the spiral to the planet) as a function of the cylindrical disk radius R is

$$\begin{eqnarray}&&\begin{array}{l}{\phi }_{{\rm{L}}}(R)=-\mathrm{sgn}(R-{R}_{p})\displaystyle \frac{\pi }{4m}\\ \quad \ -\ {\displaystyle \int }_{{R}_{m}^{\pm }}^{R}\displaystyle \frac{{\rm{\Omega }}(R^{\prime} )}{{c}_{s}(R^{\prime} )}{\left|{\left(1-\displaystyle \frac{R{{\prime} }^{3/2}}{{R}_{p}^{3/2}}\right)}^{2}-\displaystyle \frac{1}{{m}^{2}}\right|}^{1/2}{dR}^{\prime} \end{array}\end{eqnarray} \tag{ 35 }$$

following Bae & Zhu (2018a), where $m=(1/2)(H/R{)}_{p}^{-1}$ and ${R}_{m}^{\pm }={(1\pm 1/m)}^{2/3}{R}_{p}$.

For buoyancy spirals, the phase angle ϕ_B is

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{{\rm{B}}}(R) & = & \pm 2n\pi \displaystyle \frac{\left|{{\rm{\Omega }}}_{p}-{{\rm{\Omega }}}_{K}\right|}{{{\rm{\Omega }}}_{K}}\displaystyle \frac{{H}_{{\rm{g}}}}{Z}{\left(1+\displaystyle \frac{{Z}^{2}}{{R}^{2}}\right)}^{3/2}\\ & & \times \ {\left[\displaystyle \frac{\gamma -1}{\gamma }+\displaystyle \frac{2{c}_{s}}{g}\displaystyle \frac{\partial {c}_{s}}{\partial Z}\right]}^{-1/2},\end{array}\end{eqnarray} \tag{ 36 }$$

following Zhu et al. (2015), where n is a positive integer and g is the gravity from the star. Compared with the analytic form given in Zhu et al. (2015), note that Equation (36) has an additional term in the parentheses of the right-hand side, (2c_s /g) ∂c_s /∂Z, which takes the vertical temperature stratification into account. When there is no vertical temperature stratification, Equation (36) reduces to the phase angle derived in Zhu et al. (2015). The pitch angles of the Lindblad and buoyancy spirals are computed as $\beta \equiv -{\tan }^{-1}({dR}/{Rd}\phi )$, using the phase equations in Equations (35) and (36).

As shown in Figure 5, the buoyancy spirals' pitch angle is smaller than the Lindblad spirals' pitch angle over a broad range of distance from the planet. In particular, within a few scale heights from the planet, the Lindblad spirals' pitch angle increases to 90° toward the planet, whereas buoyancy spirals' pitch angle remains ≲ 10°. We note that Figure 5 shows pitch angles of the first-order buoyancy spirals only (i.e., n = 1 in Equation (36)) for both linear theory prediction and measurement from the simulations. For higher-order buoyancy spirals (i.e., n > 1), the pitch angle is smaller than that of first-order buoyancy spirals by a factor of ≃1/n, meaning that they are more tightly wound.

It is also worth pointing out that the buoyancy spirals' pitch angle varies continuously as a function of radius without a singularity at the radial location of the planet, unlike Lindblad spirals. From the observational point of view, this implies that the inner and outer buoyancy spirals can appear connected to each other as a single spiral, especially when the spatial resolution is insufficient.

Another characteristic that distinguishes buoyancy spirals from Lindblad spirals is the large vertical motions. In order to visualize the density and velocity perturbations at different heights, we present the perturbed density and vertical, radial, azimuthal velocities at the midplane, Z/R = 0.08, and Z/R = 0.23 in Figure 6. The azimuthal profiles of the perturbed density, perturbed temperature, and vertical, radial, and azimuthal velocities at R = R_p ± 2 H_g are presented in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Focusing on the Lindblad spirals first, both density and velocity perturbations driven by Lindblad spirals are the strongest at the midplane and decrease over height. The only exception is the near-zero vertical velocity at the midplane, which is due to the symmetry across the midplane. For the buoyancy spirals, on the other hand, perturbations are the smallest at the midplane because N_Z = 0 there. Between Z/R = 0.08 and 0.23, perturbations remain comparable to or become stronger than their Lindblad counterparts. In particular, we note that the vertical velocity perturbation associated with the buoyancy spirals becomes stronger and more extended in azimuth over height. Compared to Lindblad spirals, buoyancy spirals generally produce smaller perturbations, but the vertical motions associated with buoyancy spirals can be stronger especially as we move to the surface layers. These characteristics of buoyancy spirals suggest that the best strategy to observe buoyancy spirals is to look for vertical velocity perturbations in the surface layers of face-on disks.

Finally, we note that there is a π/2 phase shift between the vertical velocity and density/temperature perturbations driven by buoyancy resonances. At R = 103.5 au, for example, the vertical velocity has peaks at ϕ = −0.29, −0.91, −1.50, and −2.23 rad, and the density and temperature perturbation are close to zero at those azimuthal locations (see the red dashed lines in Figure 7).

3.3.2. Gas–Dust Thermal Equilibrium Model

We now turn our discussion to the gas–dust thermal equilibrium model in which the gas is assumed to be thermally coupled with the dust. With this assumption, the relaxation timescale in the surface layers is much shorter than the dynamical time (Figure 2(c)), and buoyancy resonances are expected to be weak or absent in the surface layers (Figure 3 right).

The perturbed density, perturbed temperature, and vertical velocity at Z/R = 0.23 are presented in Figure 8. In Figure 9, we present the azimuthal profiles of the perturbed density, perturbed temperature, and vertical, radial, and azimuthal velocities at R = R_p ± 2H_g from the planet in the radial direction. As apparent from the figures, buoyancy spirals do not generate clear density and temperature perturbations. Vertical velocity perturbations arising from buoyancy resonances are seen, but they are much weaker compared with the standard model. We note that the resonance is not completely suppressed in this model, as suggested by the vertical velocity perturbation, because the buoyancy frequency is not strictly zero with the imposed stratified disk temperature (see Equation (2)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Speaking of Lindblad spirals, we find that they produce smaller perturbations as we move toward the surface layers, in agreement with what is seen in the standard model. However, we note that the level of density perturbations at Z/R = 0.08 and 0.23 are larger than the standard model by about a factor of 2. Note also that the Lindblad spirals are more tightly wound—compare the azimuthal angles where the Lindblad spirals meet r = 40 and 140 au in Figures 4 and 8. This is because the gas behaves (nearly) isothermally, and thus, the sound speed is smaller than the standard model where the gas behaves adiabatically.

Before we present Keplerian-subtracted moment maps, we show the vertical velocity distribution from standard 0.5 M_Jup, 1 M_Jup, and 2 M_Jup models in Figure 10. Buoyancy resonances produce vertical motions of the order of 100 m s⁻¹, with a larger magnitude for more massive planets. As discussed earlier, Lindblad spirals produce much smaller vertical perturbations compared with buoyancy spirals, especially in the surface layers. We thus do not expect to detect them in ¹²CO in face-on disks.

Zoom In Zoom Out Reset image size

We generate Keplerian-subtracted moment maps using the quadratic method described in Teague & Foreman-Mackey (2018). Using the Python package eddy (Teague 2019), we fit a Keplerian rotation profile to the rotation maps, allowing the source center, the disk inclination and position angle, and the systemic velocity to vary. Given the low inclination of TW Hya (i ∼ 5°), we do not include any terms describing an elevated emission surface.

The Keplerian-subtracted moment maps are shown in Figure 11. We note that the residual maps clearly show a coherent, tightly wound spiral structure. It is also worth pointing out the excellent agreement between the input velocity field from hydrodynamic simulations and the retrieved velocities.

Zoom In Zoom Out Reset image size

4.1.1. Case Study: TW Hya

When the perturbations driven by buoyancy spirals are sufficiently large, the spirals can be seen in velocity channel maps. We present channel maps of the synthetic ¹²CO line observation from the standard 2 M_Jup model in Figure 12. Channel maps from standard 0.5 M_Jup and 1 M_Jup models are presented in Figures 22 and 23 in Appendix B.

Zoom In Zoom Out Reset image size

The main characteristic of buoyancy spirals in channel maps is wedge-like features standing out of the so-called butterfly pattern of the Keplerian disk. In Figure 12, these non-Keplerian features are most clearly seen close to the planet, on the northeast side of the disk at v = 0.2–0.3 km s⁻¹. Along the major axis of the disk, buoyancy spirals can appear as an arc bridging the Keplerian wings (e.g., southeast side at v = − 0.1 km s⁻¹, northwest side at v = 0.1 km s⁻¹) or as an arc beyond the inner disk (e.g., southeast side at v = 0.5–0.6 km s⁻¹).

In order to examine the observational appearance of buoyancy spirals in more inclined disks, we generate additional image cubes adopting 15°, 30°, 45°, and 60° inclinations. Keplerian-subtracted moment maps for standard 1 M_Jup model are shown in Figure 13.

Zoom In Zoom Out Reset image size

Although non-Keplerian velocity components arising from buoyancy resonances are still present, they appear more as an ellipse for more inclined geometry, making it challenging to identify buoyancy spirals. For the i = 15° and i = 30° cases, the redshifted buoyancy spirals appear stronger than the i = 5° case. This is because we are more sensitive to in-plane velocity perturbations for inclined disks. As the planet opens a gap around its orbit, it alters the gas pressure profile from the unperturbed one, such that the disk gas rotates at sub-Keplerian speed inside of the planet's orbit and at super-Keplerian speed outside of the planet's orbit (Teague et al. 2018). For an inclined disk, sub-Keplerian rotation would appear as a blue-/redshifted semiellipse on the red-/blueshifted side of the disk, while super-Keplerian rotation would appear as a red-/blueshifted semiellipse on the red-/blueshifted side of the disk (see Figure 5 of Teague et al. 2019). The pattern we see across the planet's orbit in the Keplerian-subtracted moment maps in Figure 13 exactly matches with the aforementioned expectation. If we first look at the i = 15° model, there is an elongated blueshifted arc on the redshifted side (southeast side) of the disk, right inward of the planet's orbit. On the blueshifted side of the disk, we see an elongated redshifted arc inward of the planet's orbit, along with the redshifted buoyancy spiral. Right outside of the planet's orbit on the redshifted side, the super-Keplerian rotation adds redshifted residuals to the buoyancy spiral, enhancing the overall magnitude of the residual. The same pattern is consistently observed in the i = 30° model, and the rotation modulation is stronger in this case.

It is thus reasonable to conclude that disentangling the vertical motions induced by buoyancy resonances and the modulation in rotational motions associated with the gap is generally more challenging for inclined disks. However, if a buoyancy spiral is sufficiently extended in azimuth such that the spiral crosses the minor axis of the disk, this offers a possibility to disentangle vertical motions from rotational motions. This is because vertical motions do not change their sign across the minor axis while rotational motions do change their sign.

We now turn our attention to velocity channel maps. Representative velocity channels for different inclination are shown in Figure 14. Channel maps for standard 0.5 M_Jup and 1 M_Jup models are presented in Appendix B. While the characteristic features of buoyancy spirals are still present, they are clearly weaker for inclined disks due to the $\sin i$ components of the projection, in agreement with what we see in Keplerian-subtracted moment maps. This suggests that disks with small inclinations of ≲ 30° offer the best opportunities to search for buoyancy spirals.

Zoom In Zoom Out Reset image size

As we have shown with hydrodynamic simulations in Section 3.3, the development of buoyancy resonances depends on the local thermal relaxation timescale. To examine the development and observability of buoyancy resonances under various disk conditions, we explore a broad range of parameter space and compute the relaxation timescale following the approach described in Section 2. Specifically, we vary (1) the diffusion length scale λ_diff, (2) the dust scale height, (3) the maximum dust grain size, (4) the power-law slope in the dust size distribution, and (5) the disk mass. The resulting thermal relaxation timescales are shown in Figure 15. We compare t_relax with ${N}_{Z}^{-1}$ in Appendix A and Figure 21. In what follows, we focus our discussion on the surface layers of the disk (Z/R ≳ 0.1) as we are mostly concerned about the observability of buoyancy resonances using optically thick CO lines which will trace elevated regions of the disk.

Zoom In Zoom Out Reset image size

Due to the steep dependency (${t}_{\mathrm{diff}}\propto {\lambda }_{\mathrm{diff}}^{2}$), having a short diffusion length scale of λ = 0.1 H_g results in a diffusion timescale that is much shorter than the dynamical timescale almost everywhere in the disk. Even in such a case, however, the relaxation of the gas is limited by long collisional timescales between gas molecules and dust grains. In particular, the relaxation timescale in the surface layers where the ¹²CO lines probe is determined by the collisional timescale and is insensitive to the choice of λ_diff.

The gas–dust collision rate is dependent upon the detailed spatial and size distribution of grains. We first vary the scale height of dust grains by changing α in Equation (22). This has two opposite effects on the collisional timescale. For increased scale height, the mean dust size increases, which would lengthen the collisional timescale (Equation (20)). At the same time, the dust-to-gas mass ratio also increases, which would shorten the collisional timescale. The two effects effectively cancel out, and the collisional timescale in the surface layers remains sufficiently long for buoyancy resonances to develop, as can be seen in Figures 15 and 21.

When ${a}_{\max }$ is decreased, the mean grain size decreases. The dust-to-gas mass ratio increases in the surface layers because, with a smaller ${a}_{\max }$, a larger fraction of the total dust mass is in the grains that can be lofted sufficiently high. Together, this shortens the collisional timescale, resulting in ${t}_{\mathrm{relax}}\lt {N}_{Z}^{-1}$ in a larger region of the disk as shown in Figure 21. Nevertheless, ${t}_{\mathrm{relax}}\gtrsim {N}_{Z}^{-1}$ near the ¹²CO emission surface, suggesting that buoyancy resonances likely develop there.

Changes in the power-law index of the dust size distribution work in a similar way. When the dust size distribution follows a steeper power-law distribution, a larger fraction of the total dust mass is in small grains that can be lofted up high. This will result in a smaller mean grain size and a larger dust-to-gas mass ratio, which would shorten the collisional timescale.

Next, we vary the disk mass by a factor of 3. Having a larger disk mass increases the number of colliders, increasing the gas–dust collision rate. The relaxation timescale would be therefore shortened, while ${t}_{\mathrm{relax}}\gtrsim {N}_{Z}^{-1}$ in the majority of the disk regions.

Lastly, we adopt the azimuthally averaged gas density from the standard 1 M_Jup model to examine the influence the gap has on the relaxation timescale. Because dust grains are depleted within the gap, the collisional timescale becomes longer, facilitating the development of buoyancy resonances there.

In summary, we argue that the relaxation timescale of the gas in the surface layers (Z/R ≳ 0.1) is limited by infrequent gas–dust collisions under a broad range of conditions applicable to protoplanetary disks (Figure 15). This is because the dust-to-gas mass ratio is small (ρ_d/ρ_g ≲ 10⁻³) in the surface layers due to the vertical settling of dust grains. The resulting thermal relaxation timescale is comparable to or longer than the timescale associated with buoyancy oscillations ${N}_{Z}^{-1}$ (Figure 21), suggesting that the surface layers have a favorable condition for buoyancy resonances to develop.

5.1.1. Implications for Planet-induced Gaps

While we focused our analysis mainly on the surface layers so far, it is worth pointing out the finite relaxation timescale near the disk midplane. Even at large radii for which it is typically thought that less absorbing materials would result in more efficient cooling, we find that the infrequent gas–dust collision could prevent the gas from cooling instantaneously.

The finite relaxation timescale in the main body of a disk can have important implications for the formation of gaps by planets. Adopting a constant dimensionless relaxation timescale β ≡ 2π t_relax/t_orb in vertically integrated two-dimensional disks, Miranda & Rafikov (2020b) and Zhang & Zhu (2020) independently showed that the isothermal assumption does not provide a good approximation when β ≳ 10⁻³—note again that we consistently find β ≳ 10⁻³ under various disk conditions (see Figure 15). In particular, when β ∼ 1, the radiative dissipation of Lindblad spirals becomes important, and the gap around the planet becomes narrower than it otherwise would be.

In Figure 16, we present the radial profiles of the gas surface density from the standard model and the gas–dust equilibrium model. We additionally ran a simulation adopting an isothermal equation of state and included the density profile from the isothermal simulation in the same figure. The gap widths measured at the half-maximum are Δ_gap = 27.4, 36.1, and 40.4 au for the three models. Note that the gap is wider in the isothermal simulation, consistent with previous two-dimensional simulations (Miranda & Rafikov 2020b; Zhang & Zhu 2020).

Zoom In Zoom Out Reset image size

There have been many attempts to infer the masses of planets responsible for the observed gaps, using planet–disk interaction simulations or empirical planet mass–gap width relations (see Bae et al. 2018; Disk Dynamics Collaboration et al. 2020 and references therein). We note that most of the simulations were carried out adopting an isothermal equation of state, implying that the inferred planet masses can underestimate the actual planet masses. Along the same lines, we argue that empirical gap depth/width–planet mass relations that are typically derived under the isothermal assumption need a revision.

5.1.2. Implications for Hydrodynamic Instabilities

As discussed in Section 2.2, exactly at which height atomic/molecular line cooling becomes the dominant cooling mechanism depends upon the balance between various heating and cooling mechanisms, which in turn is determined by various factors including the amount of dust grains, gas temperature, the abundance of coolant atoms/molecules, and electron number density. Here, we test the effect of efficient line cooling near the surface, by adopting Z_line = 2 H_g in Equation (27) (cf. Z_line = 4 H_g in the standard model).

Figure 17 presents the vertical velocity distribution. Not surprisingly, buoyancy resonances are weaker at Z = 3H_g (Z/R = 0.23) due to the rapid cooling in the surface layers. However, we note that buoyancy resonances are not completely suppressed. This is because even when cooling is rapid (effectively γ ≃ 1), the buoyancy frequency is not strictly zero due to the stratified temperature profile (see Equation (2)). At Z = 1H_g (Z/R = 0.08), buoyancy resonances develop at the level they develop in the fiducial model adopting Z_line = 4H_g (see Figure 10). This suggests that the development of buoyancy resonances depends on the local thermodynamic properties and that, as far as line cooling is not efficient all the way to the midplane, there will be regions where buoyancy resonances would develop.

Zoom In Zoom Out Reset image size

From the observational point of view, this suggests that we can choose optically thinner lines that probe the adequate heights. To support this argument, we generate Keplerian-subtracted moment maps from the models with efficient line cooling, which are shown in Figure 18. The residual velocities are smaller than the fiducial model in ¹²CO because the line probes the surface layers where cooling is rapid. On the other hand, the morphology and magnitude of the buoyancy spirals in the ¹³CO maps are nearly identical to the standard model (see Figure 11).

Zoom In Zoom Out Reset image size

Here we discuss potential ways to discriminate the non-Keplerian motions driven by buoyancy resonances from those driven by other origins, Lindblad resonance, corrugated vertical flows, and gas pressure changes.

Lindblad spirals: The vertical dependence of the perturbation driven by buoyancy and Lindblad spirals is the key to discriminate the two. Because the buoyancy frequency is strictly zero at the disk midplane, we expect no or weak buoyancy resonances there. On the other hand, perturbations driven by Lindblad spirals are the strongest 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 versus ¹³CO or C¹⁸O, will thus help discriminate between buoyancy spirals and Lindblad spirals.

Corrugated vertical flows: Various (magneto)hydrodynamic processes are known to create radially alternating corrugated vertical flow patterns, including vertical shear instability (Nelson et al. 2013), spiral wave instability (Bae et al. 2016a, 2016b), and magnetically driven zonal flows (Johansen et al. 2009; Flock et al. 2015; Riols et al. 2020). Vertical velocity perturbations driven by buoyancy resonances (and Lindblad resonance) are symmetric against the midplane, so at a given radius, the gas motion will be either toward the midplane or toward the surface. In contrast, vertical velocity perturbations associated with corrugated vertical flows are not symmetric against the midplane. Rather, instabilities develop into corrugation modes, and the entire column oscillates vertically (Nelson et al. 2013; Bae et al. 2016a, 2016b). This difference suggests that we can distinguish the two scenarios if we probe velocity perturbations at the upper and lower surfaces of the disk separately. Optically thick tracers (e.g., ¹²CO), a sufficiently high spatial resolution, and a moderately inclined disk geometry would provide the best chance to separate the upper and lower surface emission.

Gas pressure changes: Gas pressure changes across the radius can lead the sub-/super-Keplerian rotation of the gas to maintain the radial force balance (Teague et al. 2018). As we discussed in Section 4.3, distinguishing non-Keplerian motions associated with buoyancy spirals from rotation modulation can become a challenge in inclined disks as we become sensitive to both vertical and azimuthal velocities. As an example, we present channel maps of the ¹²CO line emission of HD 143006 in Figure 19. Morphologically, the arcs connecting the Keplerian wings (most clearly seen on the north side of the disk at v = 7.24 km s⁻¹ and on the south side of the disk at v = 8.20 km s⁻¹) show a good resemblance to those features expected from buoyancy resonance. However, these arcs can instead be interpreted as rotation modulation. The non-Keplerian features are most prominently seen as redshifted arcs in blueshifted channels (e.g., v = 7.24 km s⁻¹) and blueshifted arcs in redshifted channels (e.g., v = 8.20 km s⁻¹), which are consistent with the expected modulation associated with sub-Keplerian rotation (Teague et al. 2019). In fact, the inner arcs at ∼04 coincide with the outermost continuum ring, suggesting that the arcs could arise from the rotation modulation around the pressure peak. The outer arcs at ∼08 are close to where ¹²CO emission fades, suggesting they could be due to a rapid drop in the gas density at that radius. Because of the degeneracy between vertical and azimuthal velocities in channel maps, we recommend using Keplerian-subtracted moment maps rather than velocity channel maps when it comes to distinguishing buoyancy spirals and rotation modulation arising from gas pressure changes (see Section 4.3).

Zoom In Zoom Out Reset image size

How would buoyancy resonances appear in scattered light and (sub)millimeter continuum observations? Interestingly enough, van Boekel et al. (2017) reported a tightly wound spiral in the near-infrared polarized intensity map of the TW Hya disk. The spiral is located at the outer edge of an annular gap centered at about 93 au, ⁷ and extends about 90° in azimuth on the southwest side of the disk. While we opt out of making simulated scattered light observations from our models because our simulations do not include dust grains, it is interesting to point out that a buoyancy spiral produces positive density perturbations at the exact location where the scattered light spiral is revealed (see the lower-right quadrant of δ ρ/〈ρ〉_ϕ at Z/R = 0.08 and 0.23 in Figure 6). It is interesting to speculate that the velocity spiral in the CO observation and the density spiral in the scattered light observation are probing buoyancy resonances driven by a planet embedded within the gap at ∼90 au. Future simulations including dust particles will help further investigate this possibility.

On the other hand, given that buoyancy resonances are weak or absent near the midplane and that buoyancy spirals are confined in the corotating region where large (sub)millimeter-sized grains are expected to depleted due to radial drift, we believe observing buoyancy resonances in (sub)millimeter continuum observations is less likely.