Observational Signatures of Planets in Protoplanetary Disks: Temperature Structures in Spiral Arms

We obtain an axisymmetric background temperature profile, stratified in both radial and vertical directions, using the radiative-transfer code HOCHUNK3D (Whitney et al. 2013). Because radiative transfer is not scale-free, we set fiducial parameters loosely inspired by TW Hya; the planetary location at R_p = 1 is scaled to 37 au (where one of the gaps in the system is located; see Tsukagoshi et al. 2016), and the disk mass set to 30M_J with a dust-to-gas ratio of 0.01 (all interstellar medium grains from Kim et al. 1994; well coupled to gas). We set the stellar radius to 2.09 R_⊕ and temperature to 4000 K, typical of T Tauri stars. To ensure that the vertical distribution of disk material reflects the temperature profile, we enable the HSEQ mode of HOCHUNK3D (Whitney et al. 2013), which iterates the vertical disk profile after each RT iteration until convergence is reached.

We fit the resulting HOCHUNK3D temperature as described in Appendix A, and plot the result in Figure 1. From here, we find the sound speed as ${c}_{s,0}^{2}({\boldsymbol{x}})=\gamma {k}_{B}T({\boldsymbol{x}})/\mu {m}_{{\rm{H}}}$, taking the mean molecular weight μ = 2.34 as in the minimum-mass solar nebula. At the location of the planet, the effective aspect ratio h_eff ≡ c_s,0( x )/(γ^1/2 RΩ) ≈ 0.07. For a vertically isothermal disk in hydrostatic equilibrium, this would equal the true aspect ratio h, but in a stratified disk we must define it explicitly as

$$\begin{eqnarray}&&h(R)\equiv {\left(\displaystyle \frac{1}{{R}^{2}{\rm{\Sigma }}(R)}{\int }_{-\infty }^{\infty }{z}^{2}\rho (z){dz}\right)}^{1/2}.\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

We initialize the disk with a surface density

$$\begin{eqnarray}&&{\rm{\Sigma }}={{\rm{\Sigma }}}_{0}{\left(\displaystyle \frac{R}{{R}_{p}}\right)}^{-1}.\end{eqnarray} \tag{ 7 }$$

Our science simulations use a resolution of 267(r) × 683(ϕ) × 42(θ), spanning a radial range of r = {0.4, 5}, a polar range of θ = {π/2, π/2 − 0.42} (covering ∼6 scale heights from the midplane at the planet location), ⁴ and the full 2π in azimuth. Cells are spaced logarithmically in the radial direction, but uniformly in the polar and azimuthal directions. This yields a resolution of roughly seven cells per effective scale height at the location of the planet. We use periodic boundary conditions in the azimuthal direction, and reflective boundary conditions in the polar direction; in the radial direction, we use outflow boundaries because appropriate fixed boundary conditions for density, velocity, etc for our stratified temperature structure lack an analytic expression.

We verify our implementation of beta cooling with tests at t_c = 10⁻⁴ and 10⁻⁶; these demonstrate typical fractional azimuthal temperature perturbations of order ∼ 10⁻⁵ and ∼ 10⁻¹⁴, respectively (not shown). In the vertical average, however, we see fluctuations in temperature deviating from the background at percent-level. This is because, even though temperature in any given grid cell deviates little from its background value, mass is redistributed upward and downward in the disk by the spiral density wave, changing the column-averaged temperature. This effect would not be captured in 2D simulations, so we cover it in more detail in the following section.

In the opposite limit of long cooling times, it is not immediately clear whether our simulations would reach any sort of steady state. To test whether they do, we run t_c = 10⁴ models for each planet mass, and find that they yield essentially identical temperature and density structures to our t_c = 10² runs. As an additional test, we run a M_p = 50M_⊕, t_c = 10⁴ simulation out to 400 orbits, and find that once spiral morphology and amplitude is established at t ≈ 10 orbits, it remains in steady state (up to the gap opening) throughout. This is because planet-induced spiral arms are patterns that gas enters (to be compressed and heated) and leaves (to expand and cool) over the course of a single orbit, rather than persistent accumulations of gas.

Throughout this work, we define the relative strength along the spine of a Lindblad spiral at a given cylindrical radius (ΔT/T)_spiral(R) as:

$$\begin{eqnarray}&&\,\begin{array}{l}{({\rm{\Delta }}T/T)}_{\mathrm{spiral}}(R)=\displaystyle \frac{1}{2{h}_{\mathrm{eff},{\rm{p}}}}\\ \quad \times {\displaystyle \int }_{{\phi }_{\mathrm{peak}}(R)-2{h}_{\mathrm{eff},{\rm{p}}}}^{{\phi }_{\mathrm{peak}}(R)+2{h}_{\mathrm{eff},{\rm{p}}}}d\phi ({\rm{\Delta }}T/{\left\langle T\right\rangle }_{\phi })(R,\phi )\end{array}\,\end{eqnarray} \tag{ 8 }$$

where ϕ_peak(R) is the azimuthal location of the spiral density peak and h_eff,planet = 0.07 is the effective aspect ratio at the location of the planet, and $({\rm{\Delta }}T/{\left\langle T\right\rangle }_{\phi })(R,\phi )$ is the mass-weighted, fractional azimuthal perturbation in vertically averaged temperature at a given 2D location:

$$\begin{eqnarray}&&{\left\langle T\right\rangle }_{\phi }(R)=\displaystyle \frac{{\int }_{0}^{2\pi }{\int }_{0}^{\infty }d\phi {dz}\rho T}{{\int }_{0}^{2\pi }d\phi {\rm{\Sigma }}(R,\phi )}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&({\rm{\Delta }}T/{\left\langle T\right\rangle }_{\phi })(R,\phi )=\displaystyle \frac{1}{{\rm{\Sigma }}(R,\phi )}\displaystyle \frac{{\int }_{0}^{\infty }{dz}\rho T}{{\left\langle T\right\rangle }_{\phi }(R,\phi )}-1.\end{eqnarray} \tag{ 10 }$$

We define (ΔΣ/Σ)_spiral(R) and ${\rm{\Delta }}{\rm{\Sigma }}/{\left\langle {\rm{\Sigma }}\right\rangle }_{\phi }$ analogously. This formulation is motivated by observational relevance—beam convolution would smear out a point estimate of spiral amplitude.

As an additional check, we run convergence tests in both grid resolution and smoothing length r_s at a fiducial M_p = 50M_⊕ and t_c = 10², and display our results in Figure 2. Our upper panel, plotting (ΔT/T)_spiral(R), shows a well-converged Lindblad spiral amplitude for radii within a factor of 2 of the planet location. Farther away, however, spiral strength is supported by acoustic propagation rather than resonant driving, and thus becomes degraded by numerical diffusion.

Zoom In Zoom Out Reset image size

Our lower panel—an azimuthal cut at r = 1.5—underscores that within a factor ∼ 2 of the planet position, most of the difference between resolutions arises from nonlinear wave steepening at the peak amplitude location ϕ_peak. While properly capturing spiral amplitude at this point would require resolutions far exceeding what is feasible in 3D global simulations (Dong et al. 2011a, 2011b), this only changes the integrated arm amplitude by ∼ 1.15 × . As for changing r_s , we find that it affects features in the immediate co-orbital region, but otherwise has negligible impact on our results.

Our analysis centers on a grid of 12 models, with cooling timescales t_c = {10⁻², 10⁰, 10²}—typical of real protoplanetary disks at tens of au (e.g., Miranda & Rafikov 2020; Zhang & Zhu 2020; Ziampras et al. 2020)—and planetary masses M_p = {50, 100, 200, 400}M_⊕ (corresponding to planet-star mass ratio q = {1.5 × 10⁻⁴, 3 × 10⁻⁴, 6 × 10⁻⁴ 1.2 × 10⁻³}, or q/q_th = {0.48, 0.96, 1.91, 3.83}, where ${q}_{\mathrm{th}}\equiv {h}_{\mathrm{eff}}^{3}$ is the thermal mass of the planet in the disk). At t = 0, our simulations are vertically isothermal, with a radially varying temperature given by Equation (A2) of Appendix A. We relax to the vertically stratified temperature profile (Equation (A1)) with a cooling timescale of t_cool,init = 0.1 × 2π over 10 planetary orbits. At that point, we set the cooling time to its notional value, and initialize the planet, growing it to its final mass over 1 orbit. We run each simulation for 15 more orbits, for a total of 25. We plot Cartesian density and temperature maps for one representative simulation—the M_p = 50M_⊕, t_c = 10² run used in our resolution test—in Figure 3.

Zoom In Zoom Out Reset image size

Lindblad spirals are formed by constructive interference of the Fourier modes excited in the disk by the planetary potential. They have a characteristic width of ∼ hR and a pattern speed of Ω(R_p ), with waves launched at roughly R_p ± hR_p , so one can estimate in the rotating frame that it takes gas parcels at most one dynamical time t_dyn = Ω⁻¹(R_p ) to cross through. For the shortest cooling time in our model grid (t_c = 10⁻²); therefore, PdV work by the Lindblad spiral on a gas parcel is dissipated much faster than it is performed, so compression and expansion are nearly isothermal processes. For t_c ≳ 1, PdV work is retained through an arm crossing, making the compression and expansion effectively adiabatic. We emphasize that, even when the simulation is in "steady state," only the spiral pattern is fixed; the underlying gas is not in vertical hydrostatic equilibrium (Figure 3 of Dong & Fung 2017), as gas in the inner spiral arms has significant vertical motion (Zhu et al. 2015). In Figure 4, we plot (ΔT/T)_spiral at r = 1.5 in our Lindblad spirals for all 12 simulations in our model grid.

Zoom In Zoom Out Reset image size

For t_c = 10⁻², therefore, the temperature in any given grid cell remains essentially unchanged from the prescribed background value; any nonaxisymmetry in temperature arises because kinematic effects redistribute material vertically. In the region of the Lindblad arm where a gas parcel is being compressed azimuthally—in anticipation of the spiral density peak—it becomes vertically inflated and thus heats up in vertical average. Following the density peak, gas expands azimuthally and shrinks to the midplane, cooling down in vertical average. This effect is stronger for inner spiral arms (0.3%–3%) than for the outer arms (0.1%–0.5%) whose amplitudes appear in Figure 5. Results from our test runs at t_c = 10⁻⁴, not shown here, show quantitatively similar results.

Zoom In Zoom Out Reset image size

By contrast, when t_c > 1, the PdV work done by a Lindblad spiral is no longer dissipated before a gas parcel fully crosses through the arm. As a result, the high-density central spine of the arm becomes hotter adiabatically, while the lower-density regions before and after expand and cool. With the rising cooling time, this effect grows to dominate the overall temperature perturbation, causing it to follow the density (rather than its rate of change in azimuth, as in t_c = 10⁻² runs). As shown in the lower panel of Figure 5, this picture holds at all altitudes in the disk, although increasing sound speed and distance from the planet widen and somewhat weaken the spiral perturbation.

Motivated by these considerations, we present the following fit for the vertically averaged spiral temperature amplitudes as a function of both M_p and t_c , plotted in the upper panel of Figure 4 for fiducial radius R/R_p = 1.5:

$$\begin{eqnarray}&&{({\rm{\Delta }}T/T)}_{\mathrm{spiral},\mathrm{fit}}={\left(\displaystyle \frac{{M}_{p}}{{M}_{\mathrm{th}}}\right)}^{c}\left[{f}_{{PdV}}{e}^{-{t}_{\mathrm{arm}}/{t}_{c}}+{f}_{\mathrm{rad}}\right]\end{eqnarray} \tag{ 11 }$$

where c = 0.6 is the power-law exponent of the temperature curves, f_PdV = 0.06 is the temperature perturbation arising from PdV work in the adiabatic limit, f_rad = 0.002 that from vertical material redistribution in the stratified temperature structure, and t_arm = 0.25 a characteristic timescale for compression and expansion of material in the spiral arms. Differences between fitted and simulated amplitudes are typically ≲ 0.05 × (ΔT/T)_spiral,fit. We note that a more extensive parameter survey in the future may help further refine this.

At any given planet mass, the temperature perturbation at any given planet mass is weakest at t_c = 10⁻² and strongest at t_c = 10². As visible in our Figure 6, however, density perturbation is nonmonotonic in regions at least several scale heights from the planet—strongest for t_c = 10⁻², but weakening at t_c = 1 before regaining some strength at t_c = 10² (Ziampras et al. 2020; Zhang & Zhu 2020). The analytical study of Miranda & Rafikov (2020) investigates this in detail, finding (in the linear limit) a radial separation between wavecrests set by

$$\begin{eqnarray}&&{\rm{Re}}({k}_{m})\approx \displaystyle \frac{{\left|{{\rm{\Omega }}}^{2}-{\tilde{\omega }}^{2}\right|}^{1/2}}{{c}_{{\rm{s}}}}\end{eqnarray} \tag{ 12 }$$

and an amplitude damping rate

$$\begin{eqnarray}&&{\rm{Im}}({k}_{m})=\displaystyle \frac{(\gamma -1){\rm{\Omega }}\tilde{\omega }{t}_{c}}{2({{\rm{\Omega }}}^{2}+\gamma {\tilde{\omega }}^{2}{t}_{c}^{2})}{\rm{Re}}({k}_{m})\end{eqnarray} \tag{ 13 }$$

where $\tilde{\omega }\equiv \pm m({\rm{\Omega }}-{{\rm{\Omega }}}_{p})$ is the Doppler-shifted forcing frequency of the planet, times the azimuthal wavenumber m.

Zoom In Zoom Out Reset image size

Physically, ${\rm{Im}}({k}_{m})$ represents the energy lost as beta cooling erodes the temperature component of the wave over the course of each oscillation. At all Lindblad resonances Ω_m = Ω_p (1 ∓ 1/m), these losses are maximized for a t_c,crit ≈ γ^−1/2 ≈ 1, meaning waves are suppressed at their launching points. Cooling times longer and shorter than this, however, allow Lindblad waves to propagate more freely through the disk.

In Appendix B, we present vertically averaged galleries of 2D ${\rm{\Delta }}T/{\left\langle T\right\rangle }_{\phi }$ (Figure 11) and ${\rm{\Delta }}{\rm{\Sigma }}/{\left\langle {\rm{\Sigma }}\right\rangle }_{\phi }$ maps (Figure 12), as well as (r, θ) cuts of both temperature and density normalized to the azimuthal average (Figure 13). These provide a visual and intuitive understanding of the spiral perturbation for all simulations in our model grid.

Buoyancy waves become apparent in temperature plots for t_c ≳ 1 and saturate in strength at t_c = 10². In each simulation, the vertically averaged temperature structure reveals at least three sets of buoyancy spirals, corresponding to different azimuthal tidal forcing wavenumber m (Zhu et al. 2012, 2015). Whereas for low planet masses and short cooling times these perturbations are weak (∼1%–2.5%), they become stronger and substantially distorted in the opposite limit, wrapping around the full 2π in azimuth and intersecting the Lindblad-resonance arm. In all cases with t_c ≥ 1, we find buoyancy perturbation strengths several orders of magnitude stronger than in Zhu et al. (2015), who did not use a stratified background temperature.

The Lindblad spiral is strongest close to the midplane, widening and weakening in the disk atmosphere where sound speed is higher. Buoyancy spirals, however, are weak near the midplane but stronger in the disk atmosphere, the region typically probed by gas tracers such as ¹²CO. This is made quantitatively clear in the lower panel of Figure 5. Buoyancy causes hot material to rise and expand while pushing cold material to fall and contract, leaving pressure constant (at fixed altitude); consequently, in the long-t_c limit the ratio between temperature and density perturbations is ≈ −1 for buoyancy spirals, whereas for Lindblad spirals (an adiabatic phenomenon) the ratio is ≈ (γ − 1) = 0.4.

Previous work (e.g., Zhu et al. 2015; McNally et al. 2020; Bae et al. 2021) has established that buoyancy resonances occur where the Doppler-shifted forcing frequency of the planet commensurates with the local Brunt–Väisälä frequency:

$$\begin{eqnarray}&&\pm m({\rm{\Omega }}-{{\rm{\Omega }}}_{p})=N\end{eqnarray} \tag{ 14 }$$

where in a protoplanetary disk with g = Ω² z,

$$\begin{eqnarray}&&{N}^{2}={{\rm{\Omega }}}^{2}\left[\displaystyle \frac{z}{T}\displaystyle \frac{\partial T}{\partial z}-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{z}{P}\displaystyle \frac{\partial P}{\partial z}\right].\end{eqnarray} \tag{ 15 }$$

Because our disks are vertically stratified, the first (temperature) term is always present, resulting in very weak buoyancy arms close to corotation even with a t_c = 10⁻². For longer cooling times (equivalently, the effective γ rising from 1 to its notional value of 1.4), the second (pressure) term makes resonances diagonal in the upper disk atmosphere, while the first creates an additional bend at the interface between disk midplane and atmosphere, where temperature is changing rapidly. For superthermal planets, buoyant temperature perturbations are sufficiently strong as to change N and alter the resonance location, creating coupling between different buoyancy modes that manifests as merging and splitting of buoyant arms.

We stress that while the Brunt–Väisälä frequency is well behaved as γ → 1 from above, it is undefined for a truly isothermal equation of state, for which γ = 1 and the Navier–Stokes energy equation is degenerate. Thus there are no buoyant perturbations in simulations like Juhász & Rosotti (2018), even though the background temperature is vertically stratified.

In Figure 7, we plot temperature cuts in the disk at θ = 1, 2, and 3H; this is intended to qualitatively mimic the different layers of the disk probed by different gas tracers. As before, we normalize and subtract away the average in each azimuthal ring; however, we extend the radial range of the plots to {0.5, 5} to emphasize that ${({\rm{\Delta }}T/{\left\langle T\right\rangle }_{\phi })}_{\mathrm{spiral}}$ is only ≳ 0.025 within a factor of 2–3 orbital radii from the perturber. We present results from the t_c = 1, 50 M_⊕ and 400 M_⊕ models as representative cases. In Figure 8 we compare our simulated results (scaled to physical units) with those of Teague et al. (2019).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For both planet masses, the primary Lindblad spiral widens, and a secondary spiral emerges (Bae & Zhu 2018a, 2018b), with increasing altitude in the disk. Buoyancy spirals likewise grow in strength and number with increasing altitude. Both Lindblad spirals wind more tightly with increasing distance from the planet, as expected from Equation (12)—which, in the limit r ≫ r_p , reduces to ${\rm{R}}{\rm{e}}({k}_{m})\approx (m/\gamma {h}_{{\rm{e}}{\rm{f}}{\rm{f}}}R)({{\rm{\Omega }}}_{p}/{\rm{\Omega }})$. Convolved with a beam, these two arms emanating from a single point may resemble observations of the CQ Tau (Wölfer et al. 2020). On the other hand, buoyancy spirals are tightly wound even close to the planet. We find, in line with previous isothermal simulations (e.g., Fung et al. 2014; Fung & Chiang 2016), that planet-carved gaps are circular in the midplane, and widen with increasing altitude; but, in a departure from previous work, at higher altitudes our nonisothermal gaps exhibit a substantial "break" at the planet location.

For a better observational understanding, we postprocess our hydrodynamical simulations in radiative transfer with HOCHUNK3D, using 10⁹ photons to obtain background temperatures and scattered-light images at face on for each of our t_c = 1 runs. In these simulations, PdV contributions from the hydrodynamics are obviated, and only explicit irradiation effects on ΔT are included. We plot both in Figure 9.

Zoom In Zoom Out Reset image size

For our temperature panels, we plot the θ = 0.21 ≈ 3H layer studied by gas observations, as in Figure 7. Taking a cut in this fashion allows us to isolate the effects of spiral shadowing from those of vertical redistribution of material (see Section 3.1, paragraph 2, for a more detailed discussion). In the inner disk, we find azimuthal asymmetries that trace the primary and secondary Lindblad arms, ranging from 1% for our 50 M_⊕ case to 8% for the 400 M_⊕ case whose arms strongly perturb the disk surface. The outer spirals cause no shadowing effect on temperature; but, for superthermal planets, the reduction of disk scale height near the planet location exposes the outer disk at ϕ = π to increased stellar irradiation, giving the impression of a radial "arm" in temperature. We note that this effect is expected to be transient and disappear once a gap is opened.

While spiral shadowing can be noticeable at high altitudes, it has a markedly lower impact on the disk as a whole. In mass-weighted vertical average, the azimuthal temperature perturbation obtained with HOCHUNK3D is typically <1%–3% interior to the planet, but <1% in the outer disk. As these differences are consistent with those plotted in Figures 4 and 6 for our t_c = 10⁻² runs, we surmise that they arise primarily from vertical redistribution of disk material in the preexisting stratified temperature structure, already accounted for in our hydrodynamics. In any case, PdV work in Lindblad spirals (at high altitudes and large planet masses, buoyancy spirals as well) overwhelms any radiative effects.

In our near-infrared scattered-light images, the inner Lindblad spiral becomes prominent for the superthermal M_p = 200, 400M_⊕ runs, and the outer arms are less clearly observable. This aligns with expectations that Lindblad spirals ought to be visible only when planets substantially alter the disk scattering surface, and with the simulations of Dong et al. (2016), who test substantially higher thermal masses (albeit in isothermal disks with temperature constant along cylinders). In scattered light, the inner primary and secondary Lindblad arms are nearly as strong as each other, but the outer secondary Lindblad and buoyancy arms are barely visible. As in the temperature, there is a radial pseudo"arm" in scattered light for superthermal planets that reshape disk structure around them.