© The Authors 2024

1 Introduction

The high spectral and spatial resolution of the Atacama Large Millimeter Array (ALMA) have made it possible to accurately probe temperatures and velocities at the τ = 1 surfaces of various molecular lines, such as those associated with ¹²CO, ¹³CO, and C¹⁸O (Pinte et al. 2023b). Observations in systems such as HD 163296 (Pinte et al. 2018), HD 97048 (Pinte et al. 2019), TW Hya (Teague et al. 2019, 2022), CQ Tau (Wölfer et al. 2021), and Elias 2-24 (Pinte et al. 2023a) – with the background temperature and (sub-)Keplerian velocity profiles subtracted off -have revealed localized velocity kinks, as well as large-scale spiral structures in temperature and velocity. Numerical (e.g., Pérez et al. 2018) and analytical (e.g., Bollati et al. 2021) studies indicate that such signatures are consistent with those caused by spiral wakes (Goodman & Rafikov 2001) launched at Lind-blad resonances in the disk, where the Doppler-shifted planetary forcing frequency equals the local epicyclic frequency (e.g., Goldreich & Tremaine 1978, 1979).

For computational efficiency and ease of interpretation, simulations of disk–planet interaction have historically used a 2D, vertically integrated approach, with the gas following a locally isothermal equation of state. But as the quality of observations improve, it has become increasingly necessary to account for more detailed disk structure and thermodynamics in order to make a meaningful comparison. Zhu et al. (2012) and Lubow & Zhu (2014) simulated 3D adiabatic disks, discovering additional spirals excited at distinct buoyancy resonances, where the Doppler-shifted planetary forcing frequency equals the Brunt-Väisälä frequency. Lobo Gomes et al. (2015) simulated disk-planet interaction with cooling in 2D, running their simulations to gap-opening timescales with an emphasis on vortex formation at the outer pressure bump formed by the planetary gap. Zhu et al. (2015) performed global 3D simulations of planet-driven spirals with so-called β cooling, in which the temperature in a grid cell is relaxed to a preset background value (with β being the ratio of the cooling time to the local dynamical time, ${\Omega }_{\text{K}}^{-1}$). Juhász & Rosotti (2018) performed 3D locally isothermal prescriptions, but with a vertically stratified temperature. Miranda & Rafikov (2020a,b) performed high-resolution, 2D simulations of spirals with β cooling, studying the details of angular momentum transport as cooling times varied from short (isothermal) to long (adiabatic) relative to the local dynamical time. The 3D simulations of Muley et al. (2021) incorporated β cooling along with a vertically stratified temperature structure obtained from radiative-transfer simulations, while Bae et al. (2021) went a step further and computed βcooling timescales at each point in the disk, based on radiative diffusion and gas–grain collision times. These studies conclude that Lindblad spirals propagate through the disk for β ≪ 1 (isothermal limit) or β ≪ 1 (adiabatic limit) but damp close to the planet location for β ≈ 1. Temperature stratification also causes the pitch angle and morphology of these spirals to deviate from the 2D, vertically averaged expectation, particularly at the high altitudes amenable to observations.

In addition to parametrized cooling, radiation-hydrodynamic techniques have also been used in the context of disk-planet interaction. A number of works have concentrated on planetary torques and migration (e.g., Kley et al. 2009; Lega et al. 2014; Fung et al. 2017; Chrenko & Nesvorný 2020; Yun et al. 2022), while others have focused on flows in the circumplanetary region (e.g., Klahr & Kley 2006; Szulágyi 2017; Chrenko & Lambrechts 2019). Such studies have typically employed a one-temperature (1T; Kley 1989) scheme – in which gas and radiation are assumed to have the same temperature – or a more involved two-temperature (2T; Bitsch et al. 2013b) scheme – in which matter and radiation have separate temperatures, coupled by opacity. For the transport of radiation, these simulations used flux-limited diffusion (FLD; Levermore & Pomraning 1981) approaches to solve the governing equations, including various parametrizations for the stellar irradiation component. The high-resolution 2D simulations of Ziampras et al. (2023), which used the 2T FLD scheme outlined in Ziampras et al. (2020), gave more emphasis to the morphologies of the spirals themselves. Their work finds that the transport of energy across the spiral shock (Ensman 1994; Commerçon et al. 2011) shifts and broadens the spiral temperature perturbation in ways that an inherently local, β cooling approximation cannot. Concurrently, Muley et al. (2023) ran 3D simulations with M1 radiation transport (Levermore 1984; Melon Fuksman & Mignone 2019) – which accounts for both the optically thick diffusion and optically thin free-streaming regimes – as well as ray-traced irradiation from the central star, as a test of their “three-temperature” (3T) approach, in which gas, grains, and radiation are coupled via collisions and opacity . In both 2T and 3T simulations, they find weaker pre-shock heating than in Ziampras et al. (2023), potentially attributable to a lower disk mass enabling efficient vertical cooling (Ziampras 2023, priv. comm.).

In this work, we build on these previous disk-planet simulations using our 3T method, with the aim of better connecting kinematic and thermal spiral signatures to the properties of the planets driving them. In Sect. 2, we describe our methods, including our thermodynamic prescriptions, treatment of planetary accretion luminosity, and disk initial conditions. In Sect. 3, we discuss the spiral structure, flow patterns, and background temperature created by a non-accreting, Saturn-mass planet orbiting at 40 au in the disk, testing disk–planet interaction under four different thermodynamic prescriptions (local isothermality, physically motivated β cooling, 2T radiation hydrodynamics, and 3T radiation hydrodynamics). In Sect. 4, we run 3T simulations only, measuring the effects of changing the mass and accretion luminosity. In Sect. 5, we present polar cuts of temperature and sky-projected velocity high in our simulated disk, and comment on the observational implications. Finally, in Sect. 6, we summarize and conclude our work.

2 Methods

2.1 Basic equations

For our study of spiral arms, we used a version of the PLUTO hydrodynamical code (Mignone et al. 2007), modified to solve the equations of radiation hydrodynamics (Melon Fuksman & Mignone 2019; Melon Fuksman et al. 2021) with an additional dust internal energy field (Muley et al. 2023). This field interacts thermally with the gas and radiation field, but passively traces the same velocity field as the gas without any back-reaction (implying a Stokes number St ≪ 1, as well as a globally constant dust-to-gas ratio f_d ≪ 1): $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathit{\upsilon })=0$$(1a) $$\frac{\partial (\rho \mathit{\upsilon })}{\partial t}+\nabla \cdot (\rho \mathit{\upsilon }\mathit{\upsilon })=-\nabla p-\rho \nabla \Phi +{\mathit{S}}_{\text{m}}+\mathit{G}$$(1b) $$\frac{\partial E_{\text{g}}}{\partial t}+\nabla \cdot \left(E_{\text{g}}\mathit{\upsilon }\right)=-\nabla \cdot ((p+\rho \Phi )\mathit{\upsilon })+S_{\text{m}}+X_{\text{gd}}+c{G}_{\text{g}}+S_{\text{g}}^{\text{irr}}$$(1c) $$\frac{\partial E_{\text{d}}}{\partial t}+\nabla \cdot \left(E_{\text{d}}\mathit{\upsilon }\right)=-X_{\text{gd}}+c{G}_{\text{d}}+S_{\text{d}}^{\text{irr}}$$(1d) $$\frac{\partial E_{\text{r}}}{\partial t}+\widehat{c}\nabla \cdot {\mathit{F}}_{\text{r}}=-\widehat{c}\left(G_{\text{g}}+G_{\text{d}}\right)$$(1e) $$\frac{\partial {\mathit{F}}_{\text{r}}}{\partial t}+\widehat{c}\nabla \cdot P_{\text{r}}=-\widehat{c}\mathit{G},$$(1f)

in which ρ, υ, and p represent the gas density, velocity, and pressure, respectively. ρ_d is the dust density, while f_d is the dust-to-gas ratio, Φ is the gravitational potential, and {E_g, E_d, E_r} are total energy densities for gas, dust, and radiation, respectively. S_m represents parabolic source terms such as α-viscosity, X_gd represents energy exchange between gas and dust. F_r is the radiative flux, while P_r is the radiation pressure. The G_g, G_d, and G terms represent opacity-mediated interaction between the gas, dust, and radiation, respectively. $S_{\text{g}}^{\text{irr}}$ and $S_{\text{d}}^{\text{irr}}$ represent, respectively, the gas and dust absorption of stellar irradiation. c represents the speed of light, whereas the c term is a “reduced speed of light” (Gnedin & Abel 2001), which enables longer timesteps but must nevertheless exceed all hydrodynamic velocities relevant to the problem (Skinner & Ostriker 2013).

The opacity source terms are given by $$\begin{array}{l}{G}_{\text{g}}\equiv -\rho {\kappa }_{\text{g}}\left(a_{\text{r}}{T}_{\text{g}}^4-E_{\text{r}}\right)\\ -\rho \left(2{\kappa }_{\text{g}}-{\chi }_{\text{g}}\right)\mathit{\beta }\times {\mathit{F}}_{\text{r}}-\rho {\chi }_{\text{g}}\mathit{\beta }\cdot \left(E_{\text{r}}\mathit{\beta }+\mathit{\beta }\times {\mathbf{P}}_{\text{r}}\right)\end{array}$$(2a) $$\begin{array}{l}{G}_{\text{d}}\equiv -\rho {\kappa }_{\text{d}}{f}_{\text{d}}\left(a_{\text{r}}{T}_{\text{d}}^4-E_{\text{r}}\right)\\ -\rho f_{\text{d}}\left(2{\kappa }_{\text{d}}-{\chi }_{\text{d}}\right)\mathit{\beta }\cdot {\mathit{F}}_{\text{r}}-\rho f_{\text{d}}{\chi }_{\text{d}}\mathit{\beta }\cdot \left(E_{\text{r}}\mathit{\beta }+\mathit{\beta }\cdot P_{\text{r}}\right)\end{array}$$(2b) $$\begin{array}{l}{\mathit{G}}_{\text{g}}\equiv \rho {\chi }_{\text{g}}{\mathit{F}}_{\text{r}}-\rho {\kappa }_{\text{g}}\left(a_{\text{r}}{T}_{\text{g}}^4-E_{\text{r}}\right)\mathit{\beta }\\ -2\rho {\kappa }_{\text{g}}\left(\mathit{\beta }\cdot {\mathit{F}}_{\text{r}}\right)\mathit{\beta }-\rho {\chi }_{\text{g}}\left(E_{\text{r}}\mathit{\beta }+\mathit{\beta }\cdot P_{\text{r}}\right)\end{array}$$(2c) $$\begin{array}{l}{\mathit{G}}_{\text{d}}\equiv \rho f_{\text{d}}{\chi }_{\text{d}}{\mathit{F}}_{\text{r}}-\rho f_{\text{d}}{\kappa }_{\text{d}}\left(a_{\text{r}}{T}_{\text{d}}^4-E_{\text{r}}\right)\mathit{\beta }\\ -2\rho f_{\text{d}}{\kappa }_{\text{d}}\left(\mathit{\beta }\cdot {\mathit{F}}_{\text{r}}\right)\mathit{\beta }-\rho f_{\text{d}}{\chi }_{\text{d}}\left(E_{\text{r}}\mathit{\beta }+\mathit{\beta }\cdot P_{\text{r}}\right),\end{array}$$(2d)

where G ≡ G_g + G_d and β ≡ υ/c (not to be confused with the β cooling timescale).

The solution strategy for these equations is described thoroughly in the aforementioned articles (Melon Fuksman & Mignone 2019; Melon Fuksman et al. 2021; Muley et al. 2023), and we recapitulate the most relevant details here. Equations (1) are divided into radiation and hydrodynamic subsystems, which are solved using a Strang split (half-step radiation, full-step hydro, half-step radiation). For the radiation subsystem, the Courant–Friedrichs–Levy (CFL) condition imposed by ĉ is far more stringent than that imposed on the hydrodynamic subsystem by the sound speed c_s, so the radiation half-step is in turn divided into a number of sub-steps. Each individual sub-step is handled using an implicit-explicit (IMEX) strategy (Pareschi & Russo 2005): radiation transport terms are computed explicitly using the M1 formalism to handle both the optically thick diffusion and optically thin beaming limits (Levermore 1984), whereas the stiff G_g, G_d, G, and X_gd source terms are solved implicitly with a Newton-Raphson method. The stellar irradiation source term $S_{\text{d}}^{\text{irr}}$ is computed explicitly by ray-tracing at the start of each hydrodynamic timestep (with $S_{\text{g}}^{\text{irr}}$ set to zero in the current study); because it depends only on the density distribution, it is sufficient to simply incorporate it into the initial guess of our Newton scheme, without updating it at each iteration.

We assumed that the gas follows an ideal equation of state, $$p={\rho }_{\text{g}}{k}_{\text{B}}{T}_{\text{g}}/\mu m_{\text{H}}\text{, }$$(3)

where k_B is the Boltzmann constant, T_g the gas temperature, µ the mean molecular weight, and m_H the mass of a hydrogen atom. We further specify that the adiabatic index γ ≡ (∂ ln p/∂ In ρ)_entropy is a constant, implying an internal energy density: $${\xi }_{\text{g}}=p/(\gamma -1).$$(4)

We note that realistic equations of state – in which γ varies with temperature as rotational, vibrational, and dissociation modes of para- and ortho-hydrogen are activated (Decampli et al. 1978; Boley et al. 2007, 2013) – would change Eqs. (3) and (4), as well as the gas–dust energy exchange term in Eq. (5), and the gas-dust coupling timescale in Eq. (6), in the following section.

2.1.1 Three-temperature simulations

For our 3T simulations, we solved the full set of Eq. (1) and defined the dust-gas collision term, $$X_{\text{gd}}\equiv t_{\text{c}}^{-1}\left(r_{\text{gd}}{\xi }_{\text{d}}-{\xi }_{\text{g}}\right),$$(5)

in which t_c is the dust-gas thermal coupling time, ξ_d = E_d represents the dust internal energy (equivalent in our scheme, but in general different if dust dynamics were to be accounted for), and r_gd = c_g/f_dc_d is the ratio of heat capacity per unit volume between gas and dust. c_d is the specific heat capacity of dust, while c_g ≡ k_B/µm_H(γ − 1) is that of the gas. We computed the gas cooling time, t_c, as a function of the dust-gas stopping time, t_s, calculated in the Epstein regime (Burke & Hollenbach 1983; Speedie et al. 2022): $$t_{\text{c}}=\frac{2/3}{\gamma -1}{f}_{\text{d}}^{-1}{t}_{\text{s}}{\eta }^{-1},$$(6)

where we set the “accommodation coefficient,” η, to unity.

2.1.2 Two-temperature simulations

In our studies of the 2T regime of traditional radiation hydrodynamics, in which dust and gas are perfectly coupled, we set the stopping time to an artificially low t_s = 10⁻¹⁰ yr.

2.1.3 β cooling simulations

For beta-cooling simulations, we replaced X_gd with a term of the form $$X_{\text{gd},\beta }=t_{\text{rel}}^{-1}\frac{\rho k_{\text{B}}}{\mu m_p(\gamma -1)}\left(T_{\text{g}}-T_{\text{g},0}(r,\theta )\right),$$(7)

where T_g,0 is the initial condition for gas temperature (see Sect. 2.3 for more details), x is a position in the disk and t_rel can in general be a function of any primitive variables. We ignored the G_g and $S_{\text{𝑔}}^{\text{irr}}$ terms and eliminated Eq. (1d), to (1f) entirely. Following Bae et al. (2021) and Melon Fuksman et al. (2024b), we set the cooling (more accurately, thermal relaxation) time as $$t_{\text{rel}}=t_{\text{c}}+t_{\text{rad}},$$(8)

where t_c is defined as in Eq. (6), and t_rad is a radiative cooling timescale incorporating both the optically thick diffusion and optically thin free-streaming limits: $$t_{\text{rad }}\equiv \max \left({\lambda }_{\text{thin }}^2,{\lambda }_{\text{diff }}^2\right)/D\text{, }$$(9)

where the radiative diffusion coefficient $D\equiv 4c{a}_{\text{r}}{T}_{\text{g}}^3/3c_{\text{g}}{\kappa }_{\text{R}}{\rho }^2$. The effective optically thin cooling length scale ${\lambda }_{\text{thin }}={\left(3{\kappa }_{\text{R}}{\kappa }_{\text{P}}{\rho }^2\right)}^{-1/2}$, while the thick diffusion length is assumed to equal the local scale height, λ_diff ≡ H = c_s,isoΩ⁻¹, where Ω is the local Keplerian orbital frequency and $c_{\text{s},\text{ iso }}\equiv \sqrt{p/\rho }$ the “isothermal sound speed”.

2.1.4 Locally isothermal simulations

To test the locally isothermal case, we simply ran a β cooling simulation with t_rel = 10⁻¹⁰ yr everywhere in the disk.

2.2 Implementation of accretion luminosity

We implemented accretion luminosity as part of the irradiation source term, $S_{\text{d}}^{\text{irr}}$ This allowed us to use the implicit-explicit strategy discussed in Sect. 2.1 to deposit, exchange, and transport large amounts of energy into each grid cell, avoiding the numerical instabilities that would arise from an explicit approach. In this work, we modeled only the luminosity, without removing any mass from the grid or adding any to the planet.

We interpolated the luminosity onto the grid using a triangle-shaped cloud approach (e.g., Eastwood 1986), in which the kernel takes the form $$\begin{array}{l}\psi \left(\mathit{x},{\mathit{x}}_{\text{p}}\right)=\max \left[\left(1-\left|\frac{r-r_{\text{p}}}{\Delta r_{i_{\text{p}}}}\right|\right)\left(1-\left|\frac{\theta -{\theta }_{\text{p}}}{\Delta {\theta }_{j_{\text{p}}}}\right|\right)\left(1-\left|\frac{\varphi -{\varphi }_{\text{p}}}{\Delta {\varphi }_{k_{\text{p}}}}\right|\right),0\right]\\ \times {\left(r^2\sin \theta \right)}^{-1},\end{array}$$(10)

where {i_p, j_p, k_p} indicate the indices in the {r, θ, ϕ} directions of the cell in which the planet is located. The kernel intersects the cells {i_p ± 1, j_p ± 1, k_p ± 1}, for a total of 27 cells¹; the fraction of accretion energy deposited in each cell during a timestep is equal to the integral of ψ over the cell volume. The triangle-shaped cloud approach avoids creating spurious discontinuities in the deposited radiation field or its gradient, unlike the cloud-in-cell or nearest grid point methods, where moving across a cell boundary or even within a cell can cause sharp changes in these quantities.

Fig. 1 Initial conditions for density (above) and temperature (below), obtained using the iterative procedure described in Sect. 2.3. White lines indicate the corresponding cooling-time contours, computed using Eq. (8). Due to our constant dust-to-gas ratio and assumption of small dust grains throughout the disk, we obtain shorter radiative–diffusion and gas–grain coupling times than in Bae et al. (2021).

2.3 Setup and initial conditions

In all simulations presented here, we assumed a disk gas surface density profile of $${\Sigma }_{\text{g}}=200 \text{g} {\text{cm}}^{-2}{\left(\frac{R}{1\text{au}}\right)}^{-1}.$$(11)

which corresponds to an M_d = 0.012 M_⊙ within our domain. As in our previous works (e.g., Melon Fuksman & Klahr 2022; Melon Fuksman et al. 2024a,b), we approximated the behavior of the gas with an ideal equation of state, with adiabatic index γ = 1.41 and mean molecular weight µ = 2.3² The total dust-to-gas ratio was assumed to be 1%, of which we took 10% by mass (corresponding to f_d = 10⁻³) to be the small grains that we modeled. These consist of 62.5% silicate and 37.5% graphite, with a material density ρ_gr = 2.5 g cm⁻³ and a specific heat capacity c_d = 0.7 Jg⁻¹ K⁻¹ ; their sizes follow an MRN (Mathis et al. 1977) distribution, n(a) ∝ a^−3.5, with a_min = 5 nm and a_max = 250 nm. Using the frequency-dependent opacities given by Krieger & Wolf (2020, 2022), we created tables of Planck and Rosseland opacity for this grain distribution as a function of temperature.

To obtain initial conditions, we employed the same iterative technique used in Melon Fuksman & Klahr (2022), cycling hydrostatic-equilibrium and radiative transfer calculations on a timescale (t_iter = 0.1 y) much shorter than thermal diffusion times through the disk. We ran 1000 iterations in order to obtain a well-converged background profile and used this profile as the initial condition for all simulations, regardless of the physics prescription implemented.

For our planet, located at r_p = 40 au, θ_p = π/2, and ϕ_ρ = π/4, this profile yields a local temperature T_p = T_g.0(r = 40au,θ = π/2) = 26 K and a scale-height ratio (H/r)_p = 0.065. We plot our initial conditions in Fig. 1, along with contours for the gas-grain coupling time (on the upper density plot) and radiative cooling time (on the lower temperature plot), normalized by the local dynamical time to obtain effective /β-values; in our β cooling simulations we simply added β_dg + β_rad. Compared with the detailed calculations of Bae et al. (2021), who self-consistently computed grain settling, opacities, and radiative cooling rates, our assumption of globally constant dust-to-gas ratio and grain size distribution yields shorter cooling times throughout the disk. In the upper atmosphere, this is because our unsettled small grain distribution makes collisional coupling more efficient; in the midplane, the low opacities of these small grains at tens of Kelvin shortens the thermal diffusion timescale.

For all simulations presented here, we used WENO3 reconstruction, along with a van Leer limiter and shock flattening to ensure stability. Our grid resolution was 268(r) × 116(θ) × 919(ϕ), logarithmically spaced in the radial and evenly spaced in the polar and azimuthal. Our domain extended between r = {0.4,2.5} × r_p, θ = {−0.4,0.4} + π/2, and ϕ = {0,2π}. Given these dimensions, our resolution yields ~10 cells per scale height at the planet location, the same amount used in historical 2D and 3D simulations of planetary gap-opening (Fung et al. 2014; Fung & Chiang 2016), and higher than in the 3D spiral simulations of Muley et al. (2021), but somewhat lower than in the 2D simulations of Ziampras et al. (2023), which contains a more detailed analysis of the angular-momentum fluxes and multi-gap opening induced by spirals. We used a reduced speed of light ĉ = 10⁻⁴c, which Muley et al. (2023) found to work well for analogous simulations of spirals driven by planets tens of au from their host stars.

In the 2T case, our test simulations show that disk columns can oscillate around the midplane, with low-order azimuthal modes undergoing linear growth. Because these modes do not manifest in the isothermal and β cooling simulations, and appear even in the absence of a planet, we tentatively attribute them to some form of irradiation (Fung & Artymowicz 2014) or self-shadowing (e.g., Melon Fuksman & Klahr 2022, and references therein) instability. In the 3T case, this effect is naturally suppressed – we surmise due to the decoupling between gas and dust temperatures at the high altitudes where stellar irradiation is intercepted. As such, for the 2T simulation only, we imposed rapid wave-damping zones in the upper and lower regions of the domain where |θ − π/2| > 0.3, and defer further investigation of this phenomenon to future work.

3 Disk physics

3.1 Thermodynamic prescriptions

In Fig. 2, we plot perturbations in density (ρ), gas temperature (T_g), dust temperature (T_d), and perturbations in radial (in center-of-mass coordinates) velocity (υ_r) for each of our four thermodynamic prescriptions (to recapitulate, local isothermal-ity, β cooling, 2T radiation hydrodynamics, and 3T radiation hydrodynamics). We took cuts of these variables at r = 1 . 5r_p at θ = 0 and θ = 0.2 above the midplane, at t = 2500 yr³.

In the midplane, we find that perturbations in ρ and υ_r are closely aligned with one another in ϕ. In agreement with previous results (e.g., Zhu et al. 2015; Muley et al. 2021), we find these perturbations strongest in the locally isothermal simulations. Among our non-isothermal runs, perturbations in ρ, υ_r, and T_g are similar regardless of whether a β cooling, 2T, or 3T prescription is used. T_g perturbations capture the work done by compression and expansion (pdV) on a fluid parcel in the disk over one thermal relaxation time. In the midplane of our fiducial disk, this is typically comparable to or shorter than the spiral-crossing time (β ≲ 1), so the measured T_g is relatively weak and somewhat offset in ϕ with respect to the ρ perturbation (which traces the integrated compression and expansion over a fluid parcel’s orbit).

Owing to the very short dust-gas coupling time t_c in this region, T_d closely follows T_g, and likewise agrees between the 2T and 3T setups. The β cooling simulation, which assumes a constant background temperature, is inherently unable to capture changes in T_d. Unlike in the 2D FLD simulations of Ziampras et al. (2023), we do not observe substantial radiative heating of the pre-shock region by emission from hot, post-shock gas (e.g., Ensman 1994). This may be caused by efficient cooling through the disk surface (Ziampras 2023, priv. comm.) in our relatively low-mass, optically thin disk.

This picture becomes somewhat more complicated at θ = 0.2, where the cooling time lengthens, and becomes dominated by dust-gas decoupling rather than radiative diffusion. The ρ spiral is sharp in the locally isothermal case, but weaker in the other simulations. T_d agrees between the 2T and 3T simulations, but the gas temperature differs, indicating that observational signatures observable in the upper disk gas (in, e.g.,¹²CO) may not be reflected in the dust.

Fig. 2 Relative differences in disk density, dust temperature, gas temperature, and radial velocity at t = 2500 yr ≈ 10 orbits with respect to initial conditions for different physical prescriptions. We take azimuthal cuts over the whole 2π at fixed r = 1.5rp and θ = 0 (left) and 0.2 (right) above the midplane. Dust temperatures largely agree between the 2T and 3T cases, whereas gas temperatures largely agree between the 3T and β cooling cases.

3.2 Lindblad and buoyancy resonances

The observed disk signatures are, in large part, driven by resonant interactions between the disk and the planet’s gravitational potential. One such interaction, studied by decades of analytical and numerical work, occurs at the “Lindblad resonances,” where the Doppler-shifted planetary forcing frequency times an azimuthal integer wavenumber m equals the epicyclic frequency κ (e.g., Goldreich & Tremaine 1978, 1979, 1980; Kley & Nelson 2012; Bae & Zhu 2018a,b). Assuming that the disk is thin, and that the spirals are tightly wound (k_r ⪢ mr⁻¹ ) so the Wentzel-Kramers-Brillouin (WKB) approximation holds, we obtain the following dispersion relation: $$m^2{\left(\Omega -{\Omega }_{\text{p}}\right)}^2={\kappa }^2+c_{\text{s}}^2{k}_{\text{r}}^2,$$(12)

where Ω_p is the planet’s orbital frequency, and k_r is the radial wavenumber of the excitation. Given that in a Keplerian disk κ = Ω ∝ r^−3/2, we can rewrite the above in terms of k_r and m, and integrate to find the phase of the mth spiral mode: $${\varphi }_{L,m}(r)={\displaystyle {\int }_{r_{\text{m}}}^r{\left[\frac{H\left(r^{\prime }\right)}{r^{\prime }}\right]}^{-1}}{\left[{\left(1-\frac{{\Omega }_{\text{p}}}{\Omega \left(r^{\prime }\right)}\right)}^2-m^{-2}\right]}^{1/2} \text{d}{r}^{\prime }$$(13)

where r_m = (1 + m⁻¹)^2/3 r_p is the resonance location. A more detailed analysis, going beyond the WKB approximation (Artymowicz 1993; Papaloizou et al. 2007) shows that the peak mode strength occurs at approximately m ≳ 2/(H/r), above which wave excitation becomes inefficient.

One can also perform a similar analysis on the buoyancy resonances, where the m-multiplied, Doppler-shifted forcing frequency equals the (vertical) Brunt-Väisälä frequency

N_ɀ = 𝑔γ⁻¹ [ln P/ρ^γ]^1/2 (Lubow & Zhu 2014). These are an inherently 3D phenomenon, depending sensitively on the disk’s thermal physics prescription. Applying the thin-disk approximation that vertical gravity 𝑔 ≈ Ω²ɀ and assuming a background in hydrostatic equilibrium, one can simplify the above to $$m^2{\left(\Omega -{\Omega }_{\text{p}}\right)}^2=N_z^2={\Omega }^{2}z\left[\frac{1}{T_{\text{g},0}}\frac{\text{d}{T}_{\text{g},0}}{ \text{d}z}+\frac{\gamma -1}{\gamma }\frac{{\Omega }^{2}z}{c_{\text{s},\text{iso},0}^2}\right],$$(14)

where $c_{\text{s},\text{iso},0}\equiv \sqrt{p_0/{\rho }_0}$ is the isothermal sound speed in the initial condition. Following (Zhu et al. 2015; Bae et al. 2021), this can be used to obtain a phase angle for the spiral: $${\varphi }_{\text{B},\text{m}}=\pm 2\pi m\left(1-{\Omega }_{\text{p}}/\Omega \right){\left[\frac{z}{T_{\text{g},0}}\frac{\text{d}{T}_{\text{g},0}}{ \text{d}z}+\frac{\gamma -1}{\gamma }\frac{{\Omega }^2{z}^2}{c_{\text{s},\text{iso},0}^2}\right]}^{-1/2}.$$(15)

In Fig. 3, we plot the phase angles for our 3T simulations, density for Lindblad spirals and meridional velocities for buoyancy spirals. Lindblad spirals are clearly visible at all altitudes, but are somewhat less tightly bound than the linear WKB theory of Eq. (13) would predict, given the fiducial planet’s “thermal mass” $q_{\text{th }}\equiv \left(M_{\text{p}}/M_{\odot }\right)h_{\text{p}}^{-3}$ of 1.12. The deviation becomes stronger, and the spiral arm becomes more open, with increasing altitude (an effect also seen in the locally isothermal, temperature-stratified simulations of Juhász & Rosotti 2018). The vertical velocity component of the Lindblad spirals also becomes more pronounced in the disk upper atmosphere, especially in the inner disk where scale height is lower than at the wave-launching location near the planet.

As expected, buoyancy spirals are only visible in the upper layers of the disk, but likewise deviate somewhat from the linear phase prediction from Eq. (15). In comparison to previous works, the fact that the thermal relaxation time $t_{\text{rel}}\lesssim N_{\text{Z}}^{-1}$, even at high altitudes, acts to damp the buoyant oscillations. We surmise that they could be strengthened through substantial dust growth, settling, or depletion, which would increase the gas–grain coupling time, t_c, or alternatively (especially at θ ≈ 0.1 from the midplane, and in the inner disk) by a high dust-to-gas ratio or disk mass, which would increase the thermal diffusion time t_diff.

Fig. 3 Density, with the planetary spirals indicated (left), and meridional velocity (right) for various altitudes in the disk. Phase predictions for the primary Lindblad arm (Eq. (13), thick line) and m = 3 buoyancy arm (Eq. (15), thin line) are overplotted. In all panels, the cylindrical radius, R = r sin θ, is plotted on the x-axis and the azimuthal angle, θ, on the y-axis. Due to the relatively rapid thermal relaxation in this simulation, buoyancy spirals are markedly weaker than Lindblad spirals in all cases.

3.3 Flow analysis

We next turned our attention to the flow patterns in the disk. In Figs. 4–6, we plot streamlines of the fluid flow in various disk cuts, along with quivers indicating the instantaneous flow direction (see the captions for more details); in the background, we include the density perturbation for reference. Close to the mid-plane, velocities are essentially restricted to the r –ϕ direction as in 2D simulations; however, at the temperature transition, the spiral pattern weakens and bends, with the velocity acquiring a θ-component (see also, e.g., Boley & Durisen 2006, for a more general discussion of vertical spiral velocities).

What gives rise to these flow patterns, and in turn, how do they influence the observed spiral perturbations in density, temperature, and kinematics? Our system is initialized in a background state ρ₀, T_g,0, υ₀, where ${\mathit{\upsilon }}_0={\upsilon }_{\varphi ,0}\widehat{\mathit{\varphi }}$ is a divergence-less, axisymmetric, quasi-Keplerian flow set by the stellar potential Φ_* and the initial pressure profile P₀. All quantities (ρ, T_g, T_d, υ_r, υ_θ, and υ_ϕ) can be expressed in terms of the initial condition and a perturbation, for example ρ = ρ₀ + ρ′; both the initial condition and perturbation are, in general, dependent on space Working in the rotating frame of the planet, and given that we do not include gas opacities, we can write the evolution of the perturbations as follows: $$\frac{\partial {\rho }^{\prime }}{\partial t}=-\left[\mathit{\upsilon }\cdot \nabla {\rho }^{\prime }-{\mathit{\upsilon }}^{\prime }\cdot \nabla {\rho }_0\right]-\left[\rho \nabla \cdot {\mathit{\upsilon }}^{\prime }\right]$$(16a) $$\frac{\partial T_{\text{g}}^{\prime }}{\partial t}=-\left[\mathit{\upsilon }\cdot \nabla T_{\text{g}}^{\prime }+{\mathit{\upsilon }}^{\prime }\cdot \nabla T_{\text{g},0}\right]-(\gamma -1)\left[T_{\text{g}}\nabla \cdot {\mathit{\upsilon }}^{\prime }\right]-t_{\text{c}}^{-1}\left[T_{\text{g}}^{\prime }-T_{\text{d}}^{\prime }\right]$$(16b) $$\frac{\partial {\mathit{\upsilon }}^{\prime }}{\partial t}=-\left[\mathit{\upsilon }\cdot \nabla {\mathit{\upsilon }}^{\prime }-{\mathit{\upsilon }}^{\prime }\cdot \nabla {\mathit{\upsilon }}_0\right]-\left[\nabla P/\rho -\nabla P_0/{\rho }_0+\nabla {\Phi }_{\text{P}}\right]-2{\Omega }_{\text{p}}\times {\mathit{\upsilon }}^{\prime }.$$(16c)

In these equations, the partial time-derivative terms represent overall evolution of the quantity at a fixed location; given that the spiral is well developed at t = 2500 yr, and that its pattern speed equals the frame rotation speed, the system is approximately in steady state and these terms net to a small number. The terms in the first set of square brackets, including velocities projected along gradients of various quantities, represent advec-tion of the flow. We considered both the vertical (formally, along the $\widehat{\theta }$ direction) and in-plane (formally, along the $\widehat{r}$ and $\widehat{\theta }$ directions) transport of perturbed density, temperature, and velocity within our stratified disk; we write explicit expressions for these terms in Appendix B.

In the second set of square brackets are source terms representing gas compressibility, pdV work, and pressure-gravity balance in Eqs. (16a), (16b), and (16c), respectively. For simplicity, we aggregated the last of these quantities into $${\left(\frac{\partial {\mathit{\upsilon }}^{\prime }}{\partial t}\right)}_{\text{source }}=-\left[\nabla p/\rho -\nabla p_0/{\rho }_0+\nabla {\Phi }_{\text{p}}\right]\text{. }$$(17)

The final bracketed term in Eq. (16b) represents the exchange of energy between gas and dust. The term 2Ω_P × υ′ in Eq. (16c) is the Coriolis acceleration arising from frame rotation. In the second set of square brackets are the source terms, involving divergences of velocities for the scalar quantities, and gradients of pressure/gravitational potential for the velocity.

We plot all of these terms and vertical-advective terms in a cut at r = 1.5r_p and θ = 0.2 in Fig. 7. As fluid in this upper-disk region enters the spiral density wave, its quasi-Keplerian orbit is perturbed downward (due to a change in the pressure-gravity balance, υ_θ) toward the midplane – against ∇ρ₀, but along ∇T_g,0 – decreasing the local T_g while increasing ρ (θ advection). Because this flow pattern (which also includes a perturbation in υ_r) has a nonzero divergence, it also increases ρ and T_g via compression and pdV work, respectively. These terms are balanced by quasi-Keplerian transport of gas through the spiral pattern (“in-plane” advection), and additionally for T_g, colli-sional relaxation to the background dust temperature (cooling), which hold the system in an approximately steady state.

For velocity, vertical gradients are much weaker, so the θ velocity does not play a significant role. Instead, perturbations in the υ_r and υ_θ components are governed primarily by the planet-driven density wave (pressure-gravity) and counterbalanced by in-plane transport across the spiral. Along the ϕ-component, the pressure gradient is weak, and υ_ϕ is governed instead by a balance between in-plane advection and Coriolis terms. The fact that υ_r and υ_phi are somewhat out of balance reflects the fact that over the long term, planet-driven spirals open a gap in the disk.

Fig. 4 Midplane flow pattern in our disk, in the corotating frame of the planet. Streamlines proceed from bottom to top inside the planet radius rp = 40 au, and from top to bottom outside the planet radius.

Fig. 5 Flow pattern at r = 1.5rp = 60au in an azimuthal cut of the disk. Streamlines flow from right to left; unlike in Fig. 4, vector arrows represent velocity differences from the local initial (quasi-)Keplerian value, rather than from the planet’s Keplerian speed.

Fig. 6 Flow pattern at r = 0.66rp = 26.6au in an azimuthal cut of the disk. Streamlines flow from left to right; as in Fig. 5, vector arrows represent velocity differences from the local initial (quasi-)Keplerian value.

Fig. 7 Source (red) and advective (blue) terms at a fixed r = 1.5rp and θ = 0.2, plotted over azimuthal angle ϕ. From top to bottom, we plot these terms for ρ (in units of Ωpρ0,p), Tg (ΩpTg,0,p), υr, υθ, and υϕ (Ωpcs,iso,p), where the p subscript indicates quantities taken at the planet location in the initial condition. A thin dashed gray line passes through the azimuthal peak of the density spiral, showing the significant offset between terms driving spiral perturbations in each quantity.

3.4 Equilibrium temperature

Disk-planet interaction not only heats the disk locally through pdV work, but non-locally by changing the background radiation field. Spiral arms, for instance, push disk material to higher altitudes, intercepting direct stellar irradiation while gently shadowing the regions behind them. Closer to the midplane, the accumulation of material in circumplanetary regions leads to clearer, radially directed shadowing . The gas heating in both regions, whether through spiral compression or gas accretion, also makes them weak⁴ sources of radiation. All this impacts the equilibrium temperature T_eq, defined by the equation $$0=S_{\text{d}}^{\text{irr}}+{\rho }_{\text{d}}{\kappa }_{\text{d}}\left(T_{\text{eq}}\right)\left(E_{\text{r}}-a_{\text{r}}{T}_{\text{eq}}^4\right),$$(18)

which incorporates both stellar irradiation and the thermalized radiation field, and whose solution we found using Newton-Raphson iterations. Because these effects depend on the (inherently non-local) transport of radiation, they cannot be accounted for using a β cooling approach; conversely, the deviation of T_d,eq from the initial condition T_d,0 provides a good measure of the suitability of a β cooling prescription for a particular problem.

In Fig. 8, we plot deviations of Teq, T_d, and T_g in our 3T, Saturn-mass simulation; we fixed r = 1.5r_p and θ = {0,0.1,0.2}, and plotted the deviations as a function of azimuthal angle. We find that deviations in equilibrium temperature are strongest at ϕ ≈ π/4, the angular position of the planet, resulting from mid-plane shadowing and upper-atmosphere exposure to radiation. This corresponds to the bright “pseudo-arm” observed in Monte Carlo radiative-transfer modeling of near-infrared scattered light from simulated planetary spirals (Muley et al. 2021). The deviations in T_eq located approximately π/4 rad from the spiral arms are much weaker; given that they are not centered on the gas temperature bump, we attribute them to rearrangement of disk density affecting transport of stellar and reprocessed radiation, rather than emission from the spiral itself found by Ziampras et al. (2023).

As in Muley et al. (2023), we find that dust and gas temperatures largely agree with each other at lower disk layers, whereas in the upper atmosphere, longer dust-gas coupling mean that the gas temperature reflects pdV work from the spirals, while the dust temperature closely tracks the equilibrium temperature. The generally small deviation of T_eq explains why the β cooling approach discussed in Sect. 2.1.3 provides generally accurate results for non-accreting, Saturn-mass planets. However, for other setups – such as those we discuss in the following sections – this need not be the case, and radiation hydrodynamics are essential to obtaining physically consistent results.

4 Parameter study

We tested three different planet masses (M_p = {5 × 10⁻⁵, 3 × 10⁻⁴, 1 × 10⁻³ M_⊙, corresponding to Neptune-, Saturn-, and Jupiter-mass, respectively)⁵ with no accretion luminosity, and two planetary accretion luminosities (L_acc,p = {0,1 × 10⁻³}) with a Saturn-mass planet. In order to capture changes to the background radiation field – which become especially pronounced for high-mass or accreting protoplanets – we used the full 3T scheme for all of these simulations.

4.1 Planet mass

In Fig. 9, we plot the spiral-averaged perturbation amplitude ${\langle \Delta x\rangle }_{\varphi \text{,spiral }}$ – where x can be any one of various normalized quantities (lnρ/ρ0, ln T_g/T_g,0, ln T_d/T_d,0, υ_r/c_s,iso,p, υ_θ/c_s,iSo,p, υ_ϕ/ c_S,iso,p), and ∆x its deviation from the initial condition – as a function of planet mass. As in previous figures, we fixed r = 1.5r_p and tested θ = {0.0,0.1,0.2} above the midplane. We defined $${\langle \Delta x\rangle }_{\varphi ,\text{ spiral }}(r,\theta )=\frac{1}{4h_{\text{p}}}{\displaystyle {\int }_{{\varphi }_{\text{peak }}-2h_{\text{p}}}^{{\varphi }_{\text{peak }}+2h_{\text{p}}}\Delta }x(r,\theta ,\varphi )\text{d}\varphi ,$$(19)

which is analogous to the definition in Muley et al. (2021), with the important distinction that in the present work the amplitude is not vertically averaged. Open circles indicate simulations without accretion luminosity, while filled circles correspond to our simulation with it (Sect. 4.2).

As found in previous works (Fung & Dong 2015; Dong & Fung 2017; Muley et al. 2021), spiral amplitude increases substantially with planet mass, irrespective of the measure used. At θ = 0.0 and θ = 0.1, ρ and υ_r/c_S,iSo,p perturbations follow each other closely; both weaken in the upper atmosphere, but the density perturbation much more so. υ_θ/c_s,iso,p is nearly zero in the disk midplane, but rises with increasing altitude (see Sect. 3.3). υ_ϕ/c_S,iso,p is relatively weak in most cases, but is stronger for high-mass cases at θ = 0.1,0.2 thanks to the distortions introduced by the circumplanetary flow. At every altitude and planet mass, temperature deviations tend to be weaker than those in density and velocity. At θ = 0.0,0.1 where dust and gas are well coupled, they reflect both the Lindblad spiral and the radial pseudo-arm instead (Sect. 3.4), but at θ = 0.2, gas temperature primarily reflects the Lindblad spiral while dust temperature reflects the pseudo-arm.

For a more qualitative understanding, we show 2D polar cuts of υ_r in Fig. 10. Sub-thermal, Neptune-mass planets excite well-behaved kinematic spirals that are clearly distinct from the laminar background; however, these spirals are rather weak, with typical |υ_r| ≲ 0.1c_S,iso,p – corresponding to a physical velocity of ≲30m s⁻¹ in our disk model. In the Saturn-mass case, velocity spirals have a similar shape, but are significantly stronger – by a factor of ~4 in the outer disk (Fig. 9), and even more so in the inner disk. Here, flows at the corotation radius and in the vicinity of the planet, which can be identified with the velocity kinks observed in channel maps (e.g., Pinte et al. 2019), become significantly more prominent.

In the Jupiter-mass case, Lindblad spirals become somewhat stronger, with higher-order spirals (Bae & Zhu 2018a,b) particularly visible in the inner and upper disks. However, the observed υ_r signature is dominated by flows near the planet, which become larger and faster thanks to the greater planet mass and gravitational sphere of influence. In the upper atmosphere, this flow takes on a wing-like shape, and acts to bring material directly above/below the planet, where is funneled toward the planet’s Hill radius through the poles (e.g., Fung et al. 2015, 2019). At all altitudes, the background flow also becomes unsteady and more turbulent.

To complement our understanding from the kinematic signatures, we present perturbations in the gas temperature T_g in Sect. 3.4, taken at the same 2D polar cuts as in Fig. 10. In all simulations, we observe a radial shadow from the circumplane-tary region (see Sect. 3.4) whose strength increases with the size of the circumplanetary region and thus, with planet mass. Especially for higher-mass cases, the fact that the integrated one-sided Lindblad torque scales as ${\mathcal{𝒯}}_{L,+}=C_{L,+}{\left(M_{\text{p}}/M_{*}\right)}^2{h}^{-3}\propto T_{\text{g}}^{-3/2}$ (Tsang 2011) may contribute to strengthening the outer spiral. In the upper atmosphere, the funneling of gas into a flow toward the midplane along the planetary pole leads to compressional heating, manifesting as a hot spot in T_g above the planet location.

Thermal relaxation, first by gas–grain collision and then by thermal emission from heated dust grains, dissipates pdV work to the radiation field. For the disk setup we chose, the effective relaxation timescale (β = 0.1−1) is approximately equal to or less than the spiral crossing time. As a result, the T_g spiral has a multiband structure, reflecting the initial compression when the gas strikes the spiral (high-temperature band), expands behind the spiral (low-temperature band), and compresses again to return to the background density (high-temperature band). The first high-temperature band is dominant in the midplane, but both roughly even in the disk upper atmosphere. This stands in contrast to the effectively adiabatic situation where thermal relaxation is much slower than spiral crossing (Miranda & Rafikov 2020b; Muley et al. 2021), where the T_g spiral reflects the accumulated total of compression and expansion rather than individual short phases of it. This picture is clear up to Saturn-mass, but less so at Jupiter-mass, where the T_g Lindblad spiral is distorted by the effects of circumplanetary and turbulent flows.

Fig. 8 Background equilibrium, gas, and dust temperatures in our 3T, Saturn-mass simulations at selected altitudes above the midplane for a fixed r = 1.5rp. Background temperatures deviate by only a few percent from the initial condition, most strongly in the region of the planetary shadow (ϕ ≈ π/4) and somewhat more weakly near the Lindblad spiral. Td agrees well with Tg at θ = 0.0 and 0.1 (and is covered by the line for Tg), and with Teq at θ = 0.2 (covering the line for Teq).

4.2 Accretion luminosity

In recent years, the effect of planetary accretion luminosity has been studied in parametrized 2D global disk simulations (Montesinos et al. 2015, 2021; Gárate et al. 2021), and in local, 3D simulations with full radiative transfer (Szulágyi 2017) that emphasize the circumplanetary disk. We add to this body of work with our 3D global simulations, including 3T radiation transport and a luminous planet. The L_acc,p = 10⁻³L_⊙ we used corresponds to a mass accretion rate of $\dot{M}=7M_{\text{J}}{\text{Myr}}^{-1}$, given the fiducial planet mass M_p = 4 × 10⁻⁴ M_⊙ and assuming a typical radius of 2 R_J for the forming planet. For this mass, such an accretion rate can only be sustained over brief periods. We thus believe that this scenario, along with the no-accretion base-case, bracket the range of accretion luminosities that the planet might experience during its growth.

We revisit Fig. 9, in which spiral amplitudes from our simulations with accretion luminosity are plotted with filled circles. At r = 1.5r_p and θ = 0, accretion luminosity somewhat weakens the ρ and υ_r perturbations, while leaving the others unchanged. At higher altitudes, by contrast, most perturbations are significantly strengthened by accretion luminosity – sometimes by factors of 2 or more – with respect to the nonluminous planet case.

Along the lines of Figs. 10 and 11, we also made qualitative plots of kinematic and temperature perturbations at various disk cuts, and present these in Fig. 12. With planetary accretion luminosity, the circumplanetary envelope becomes more pressure-supported, with the hot core overflowing the planetary Hill radius. Close to the midplane, this disrupts the orderly Doppler-flip kinematic signature and enhances the Lindblad spirals in both temperature and velocity, without changing their overall morphology. It also significantly strengthens the outer radial shadow (Montesinos et al. 2021) through its impact on the envelope’s vertical and azimuthal structure.

In the upper atmosphere, however, accretion luminosity intensifies the spiral and changes its shape. The radial velocity shows a double-banded structure – as opposed to the single-banded one expected without accretion luminosity – while the temperature shows a cold-hot-cold band structure, rather than hot-cold-hot. With accretion, the spiral also remains strong much farther away from the planet than without.

To understand these results better, we plot the density perturbation in Fig. 13, at a vertical cut through the planet’s azimuth at ϕ_p, = π/4. In this view, it is clear that accretion luminosity puffs the circumplanetary region vertically. This changes the flow pattern significantly, with accretion primarily happening by material diagonally striking the edge of the envelope, rather than being funneled downward into a narrow polar flow perpendicular to the disk midplane, as is classically the picture (e.g., Fung et al. 2015). This is responsible for many of the observed changes to kinematic and thermal signatures, particularly in the upper atmosphere. Moreover, this means that material flowing into the Hill sphere has higher angular momentum in the accreting than non-accreting case and, as such, preferably gets expelled outward.

Fig. 9 Measurement of the azimuthal perturbation in each normalized quantity x, averaged over an azimuthal range of ϕpeak ± 2hp according to Eq. (19). Densities and temperatures are plotted in the left panel and velocities in the right. Open circles indicate non-accreting planets, and closed circles indicate those with accretion luminosity (see Sect. 4.2 for a discussion). At z/r ≈ θ = 0.2, where dust and gas are not well coupled, the plotted Tdust amplitude reflects the radial “pseudo-arm” (see the discussion in Sect. 3.4) rather than the Lindblad spiral.

5 Observational implications

Of the cases we tested, we find that (non-accreting) Saturn-mass planets, accreting Saturn-mass planets, and Jupiter-mass (non-accreting) planets are best at driving thermal and kinematic signatures that are amenable to observation. In each of these cases, we plot the planet-induced total velocity perturbation – projected along the line of sight ($v_{\perp }^{\prime }$) – in Figs. A.1, A.2, and A.3, respectively. In Figs. A.4, A.5, and A.6, we likewise plot the perturbations in gas temperature T_g. In these figures, we take a cut at θ = 0.3 above the midplane; this corresponds to ~4.6 scale heights at r_p, and roughly aligns with the expected ¹²CO J = (2−1) emission layer (Barraza-Alfaro et al. 2024). (The θ = 0.2 surface shown prominently in previous figures, at ~3.1 scale heights above the midplane, is closer to the ¹2CO J = (3−2) layer.) We tested disk inclinations i_d = {0°, 30°, 60°} and planetary position angles ϕ_p = {0°, 90°, 180°, 270°}.

For an $i_{\text{d}}=0^{°},{\upsilon }_{\perp }^{\prime }$, is dominated by ${\upsilon }_{\theta }^{\prime }$. In the non-accreting, Saturn-mass fiducial case, the observed velocity spiral is driven by the local pressure-gradient term, as identified in Sect. 3.3; the redshifted spot above the planet location reflects the polar planetary flow shown in Fig. 13. For the Saturn-mass planet with accretion luminosity, the outer Lindblad spiral is noticeably stronger and extends over a greater distance, while the planetary flow – no longer oriented in a polar direction – aligns with it, forming a double-banded velocity profile as discussed in Sect. 2.2. In the Jupiter-mass case, the polar accretion flow is more prominent than in the fiducial case, with the spiral becoming subdominant. In observations, these features could help distinguish between prominent spirals created by high-mass planets and those created by highly luminous ones.

As disk inclination increases from zero, ${\upsilon }_{\perp }^{\prime }$ increasingly reflects ${\upsilon }_{\text{r}}^{\prime }$ (and to a lesser extent ${\upsilon }_{\varphi }^{\prime }$). Because the typical magnitude of ${\upsilon }_{\text{r}}^{\prime }$ in spirals is typically much larger than ${\upsilon }_{\theta }^{\prime }$, this strengthens the observed spiral signature, both in the inner and outer disks. Inclination also interacts with the elevation angle θ = 0.3 of the plotted disk surface to distort the sky-projected areas of surface elements in an azimuthally dependent way, with disk patches on the near side of the star appearing smaller than those on the far side. Depending on the planetary position angle,ϕ_p these properties can emphasize or deemphasize the circumplanetary region, inner, or outer spirals.

In each case, temperature spirals generally follow the same structure as the kinematic spirals; however, the magnitude of the observed temperature perturbation at a given point in the disk is a scalar quantity, and unlike ${\upsilon }_{\perp }^{\prime }$ does not change with projection. Furthermore, as discussed in Sects. 3.1 and 4.1, the nonzero, finite cooling time introduces a slight offset between the temperature and velocity perturbations. Conversely, the offset between kinematic and thermal spirals – as well as the prominence of buoyancy spirals – in different tracers probing various disk layers (e.g., as studied in the MAPS program; Calahan et al. 2021) could help put limits on cooling time at various vertical positions in the disk. Planet-induced radial midplane shadows could also be used to measure cooling rates – by measuring their deviation from a straight line going through planet and star – provided that the planet’s position angle is well constrained. However, fully exploring these mechanisms would require a large hydrodynam-ical parameter that covers a range of disk masses and dust size distributions.

Our simulations were run for a relatively short period of time, emphasizing the development of spiral density waves. This means that the effects of gap opening - in particular, changes to disk illumination and equilibrium temperature (Jang-Condell & Turner 2012) and changes to ${\upsilon }_{\varphi }^{\prime }$ due to the pressure gradient at gap edges (Armitage 2020, Sect. 2.4) – have not had the time to fully develop. Especially in the Jupiter-mass case, these effects can be substantial. Simulations to gap-opening timescales are currently in progress, and we intend to present them – along with Monte Carlo radiative-transfer post-processing and synthetic observations – in a forthcoming work.

Fig. 10 Radial velocities, υr, from 3T, zero-accretion-luminosity simulations with Neptune-mass (top row), Saturn-mass (middle row), and Jupiter-mass (bottom row) planets, at various altitudes in the disk (left, middle, and right columns). Gray circles indicate the planetary Hill radius. All υr values are expressed as a function of initial isothermal sound speed at the location of the planet, $c_{\text{iso},\text{p}}=\sqrt{p_{\text{p},0}/{\rho }_{\text{p},0}}$

Fig. 11 Gas temperature perturbation for Neptune-mass (top row), Saturn-mass (middle row), and Jupiter-mass (bottom row) planets, at various altitudes in the disk (left, middle, and right columns). Gray circles indicate the planetary Hill radius. All temperature values are expressed with respect to the initial condition Tg,0.

6 Conclusion

We have run 3D, 3T radiation hydrodynamical simulations with the aim of better understanding the kinematic signatures that would be generated by forming protoplanets. Our fiducial setup, with a non-accreting Saturn-mass planet located at r_p = 40au in a 15 M_J disk, draws inspiration from the TW Hya system, where spiral arms – potentially excited by a forming planet (e.g., Muley et al. 2021; Bae et al. 2021) – have already been observed via temperature and velocity measurements (e.g., Teague et al. 2019). First, we studied the physical properties of the simulated planet-driven spirals and compared our 3T approach to several thermodynamic prescriptions commonly used in the literature. Thereafter, we investigated the effects that changing planetary masses and accretion luminosities have on the strength and morphology of planet-driven disk features.

For our fiducial disk, we find that the results of our 3T simulations agree well with those from physically motivated β cooling. This naturally follows from the fact that for this setup the background equilibrium temperature, Teq, does not change much from the initial condition. As expected from previous works (e.g., Juhász & Rosotti 2018; Muley et al. 2021; Bae et al. 2021), Lindblad spirals become more open with higher altitudes in the stratified temperature background; buoyancy spirals are weak in strength due to the relatively short cooling times we used. In the upper disk – most amenable to observation in ¹²CO – thermal and kinematic signatures are driven largely by local source terms, rather than transport across vertical gradients within the disk.

Our parameter survey shows that thermal and kinematic Lindblad spirals become stronger at high planet masses. Especially in the super-thermal mass regime, deeper planetary potential wells and larger Hill radii also enhance the signatures of circumplanetary flows. Planetary accretion luminosity adds pressure support to the circumplanetary region and reorients the classic polar accretion flows (e.g., Fung et al. 2015, 2019; Fung & Chiang 2016) slightly outward. The separate effects of increasing planet masses and accretion luminosities, relative to our fiducial case, are clearly visible in our sky-projected, inclined, and rotated plots of the upper disk layers and, thus, potentially in ALMA line observations.

Future areas of research could include testing different disk masses and dust-to-gas ratios, increasing resolution to allow for the development of hydrodynamical instabilities, which may impact the visibility of spiral signatures (e.g., Barraza-Alfaro et al. 2024), or running simulations to gap-opening timescales (e.g., Fung et al. 2014; Fung & Chiang 2016) and post-processing the results to compare to observations. The last work is currently in progress.

More sophisticated numerical approaches would expand the scope of these comparisons. To more closely reproduce spiral pitch angles and morphologies in real disks (Miranda & Rafikov 2019), particularly for the 20–200 K temperature range spanned by our simulations, it would significantly help to relax the ideal-gas assumption and compute γ as a function of temperature, as well as hydrogen (para-, ortho-, and atomic), helium, and metal fractions (Boley et al. 2007, 2013; Bitsch et al. 2013a). Incorporating multiple dust species (including momentum exchange and turbulent diffusion) would enable modeling of widened, thickened rings of millimeter grains at planetary gap edges (Bi et al. 2021, 2023), which are visible in the dust continuum. On the other hand, a short-characteristics approach to radiation transport (e.g., Davis et al. 2012), which would allow for beam-crossing between accretion luminosity and reprocessed stellar radiation, could more accurately simulate strengths for radial shadows from luminous planets. Coupling such methods with our 3T scheme, however, is a computationally expensive and technically ambitious undertaking that we defer to future work.

Fig. 12 Difference from the initial condition in υr (above) and Tg (below) for 3T simulations of Saturn-mass planets without and with accretion luminosity. The strongly luminous planet alters the vertical and azimuthal structure of the circumplanetary region, causing strong shadows behind the planet and greatly enhancing the kinematic signatures of accretion in the outer disk.

Fig. 13 Density perturbation in a cut through ϕρ = π/4 as a function of cylindrical radius, R, and vertical position, ɀ, for 3T simulations of Saturn-mass planets without (above) and with (below) accretion luminosity. Gray arrows indicate the velocity field, and the gray circle encloses the planetary Hill radius. When accretion luminosity is included, the circumplanetary envelope becomes larger and the flow around it is altered.

Acknowledgements

We thank Myriam Benisty and Kees Dullemond for the suggestion to include sky-projected maps of velocity and temperature, as well as Marcelo Barraza Alfaro and Richard Teague for detailed comments on the manuscript. We also acknowledge useful discussions with Xue-Ning Bai, Jiaqing Bi, Remo Burn, Can Cui, Mario Flock, Gabriel-Dominique Marleau, Dominik Ostertag, Jiahan Shi, Jess Speedie, and Alexandros Ziampras. We thank the anonymous referee for a helpful report, and particularly for encouraging us to comment on the equation of state. Numerical simulations were run on the Cobra and Raven clusters of the Max-Planck-Gesellschaft and the Vera Cluster of the Max-Planck-Institut für Astronomie, both hosted by the Max Planck Computing and Data Facility (MPCDF) in Garching bei München. The research of J.D.M.F. and H.K. is supported by the German Science Foundation (DFG) under the priority program SPP 1992: “Exoplanet Diversity” under contract KL 1469/16-1/2.

Appendix A Sky-projected figures

Fig. A.1 Sky-projected velocity perturbation, ${\upsilon }_{\perp }^{\prime }$, for our fiducial Saturn-mass planet, in a cut at θ = 0.3 radians above the midplane.

Fig. A.2 Sky-projected velocity perturbation, ${\upsilon }_{\perp }^{\prime }$, for a Saturn-mass planet with an accretion luminosity Lacc,p = 10−3L⊙. Compared to the fiducial case, the outer spiral becomes significantly stronger and extends through a larger radial range of the disk.

Fig. A.3 Sky-projected velocity perturbation, ${\upsilon }_{\perp }^{\prime }$, for a Jupiter-mass planet with no accretion luminosity. Relative to the fiducial case, circumplan-etary and in-gap flow patterns are significantly stronger.

Fig. A.4 Perturbation in gas temperature, Tg − Tg,0, at θ = 0.3 above the midplane for our fiducial setup with a Saturn-mass, non-accreting planet. The double-armed structure of the Lindblad spiral in the upper atmosphere is clearly visible here.

Fig. A.5 Gas temperature perturbation for an accreting Saturn-mass planet. As with the kinematic spiral, the thermal spiral extends over a larger radial range and is more prominent than in the non-accreting case.

Fig. A.6 Gas temperature perturbation for a Jupiter-mass planet. As with kinematics, the circumplanetary region and gap make a much larger contribution to the thermal signature than in the fiducial Saturn-mass, non-accreting case.

Appendix B Advective terms in spherical coordinates

For a scalar quantity such as density or temperature, the rate of change at a specific location due to advection in any given direction is given by the (negative) gradient of the scalar, projected along the velocity in that direction. The advection of perturbations is given by the advective term in the current state, minus advection in the initial condition.

Making use of the fact that our initial condition is axisym-metric in all variables, and that υ_r(t = 0) = υ_θ(t = 0) = 0, we can write the evolution of the density perturbation ρ′ as $${\frac{\partial {\rho }^{\prime }}{\partial t}|}_{\text{advr}}=-{\upsilon }_{\text{r}}^{\prime }{\partial }_r\rho$$(B.1a) $${\frac{\partial {\rho }^{\prime }}{\partial t}|}_{\text{adv}\theta }=-\frac{{\upsilon }_{\theta }^{\prime }}{r}{\partial }_{\theta }\rho$$(B.1b) $${\frac{\partial {\rho }^{\prime }}{\partial t}|}_{\text{adv}\varphi }=-\frac{{\upsilon }_{\varphi }}{r\sin \theta }{\partial }_{\varphi }\rho ,$$(B.1c)

and that of the gas temperature perturbation, $T_{\text{g}}^{\prime }$, as $${\frac{\partial T_{\text{g}}^{\prime }}{\partial t}|}_{\text{advr}}=-{\upsilon }_{\text{r}}^{\prime }{\partial }_r{T}_{\text{g}}$$(B.2a) $${\frac{\partial T_{\text{g}}^{\prime }}{\partial t}|}_{\text{adv}\theta }=-\frac{{\upsilon }_{\theta }^{\prime }}{r}{\partial }_{\theta }{T}_{\text{g}}$$(B.2b) $${\frac{\partial T_{\text{g}}^{\prime }}{\partial t}|}_{\text{adv}\theta }=-\frac{{\upsilon }_{\varphi }^{\prime }}{r\sin \theta }{\partial }_{\varphi }{T}_{\text{g}}.$$(B.2c)

For vectors, such as velocity, the advection term includes not only partial derivatives of each component in each coordinate, but also geometric connection terms arising from the change in the basis vectors themselves as a function of coordinate. This yields the following advective terms in each direction: $${\frac{\partial {\mathit{\upsilon }}^{\prime }}{\partial t}|}_{\text{advr}}=-{\upsilon }_{\text{r}}^{\prime }\left[{\partial }_{\mathit{r}}{\upsilon }_{\text{r}}^{\prime }\widehat{\mathit{r}}+{\partial }_r{\upsilon }_{\mathit{\theta }}^{\prime }\widehat{\mathit{\theta }}+{\partial }_r{\upsilon }_{\varphi }\widehat{\mathit{\varphi }}\right]$$(B.3a) $${\frac{\partial {\mathit{\upsilon }}^{\prime }}{\partial t}|}_{\text{adv}\theta }=-\frac{{\upsilon }_{\theta }^{\prime }}{r}\left[\left({\partial }_{\theta }{\upsilon }_{\text{r}}^{\prime }-{\upsilon }_{\theta }^{\prime }\right)\widehat{\mathit{r}}+\left({\partial }_{\theta }{\upsilon }_{\theta }^{\prime }+{\upsilon }_r^{\prime }\right)\widehat{\mathit{\theta }}+\left({\partial }_{\theta }{\upsilon }_{\varphi }^{\prime }+{\partial }_{\theta }{\upsilon }_{0,\varphi }\right)\widehat{\mathit{\varphi }}\right]$$(B.3b) $$\begin{array}{l}{\frac{\partial {\mathit{\upsilon }}^{\prime }}{\partial t}|}_{\text{adv}\varphi }=-\frac{{\upsilon }_{\varphi }}{r\sin \theta }\left[{\partial }_{\varphi }{\upsilon }_{\text{r}}^{\prime }\widehat{\mathit{r}}+{\partial }_{\varphi }{\upsilon }_{\theta }^{\prime }\widehat{\mathit{\theta }}+\left({\partial }_{\varphi }{\upsilon }_{\varphi }+{\upsilon }_{\text{r}}\sin \theta +{\upsilon }_{\theta }\cos \theta \right)\widehat{\mathit{\varphi }}\right]\\ +\frac{{\upsilon }_{\varphi }^2-{\upsilon }_{\varphi ,0}^2}{r\sin \theta }(\sin \theta \widehat{\mathit{r}}+\cos \theta \widehat{\mathit{\theta }}).\end{array}$$(B.3c)

We emphasize that despite considering only advection in one given coordinate direction, each of the expressions in Equation B.3 has three components, arising from the advection of velocity components orthogonal to the advection direction, as well as the aforementioned geometric connection terms. For ease of interpretation, in Figure 7, we define the in-plane advection as the sum (∂υ′/∂t)_{adv,in−plane} = (∂υ′/∂t)_advr + (∂υ′/∂t)_advϕ.

References

All Figures

Fig. 1 Initial conditions for density (above) and temperature (below), obtained using the iterative procedure described in Sect. 2.3. White lines indicate the corresponding cooling-time contours, computed using Eq. (8). Due to our constant dust-to-gas ratio and assumption of small dust grains throughout the disk, we obtain shorter radiative–diffusion and gas–grain coupling times than in Bae et al. (2021).

In the text

Fig. 2 Relative differences in disk density, dust temperature, gas temperature, and radial velocity at t = 2500 yr ≈ 10 orbits with respect to initial conditions for different physical prescriptions. We take azimuthal cuts over the whole 2π at fixed r = 1.5rp and θ = 0 (left) and 0.2 (right) above the midplane. Dust temperatures largely agree between the 2T and 3T cases, whereas gas temperatures largely agree between the 3T and β cooling cases.

In the text

Fig. 3 Density, with the planetary spirals indicated (left), and meridional velocity (right) for various altitudes in the disk. Phase predictions for the primary Lindblad arm (Eq. (13), thick line) and m = 3 buoyancy arm (Eq. (15), thin line) are overplotted. In all panels, the cylindrical radius, R = r sin θ, is plotted on the x-axis and the azimuthal angle, θ, on the y-axis. Due to the relatively rapid thermal relaxation in this simulation, buoyancy spirals are markedly weaker than Lindblad spirals in all cases.

In the text

	Fig. 4 Midplane flow pattern in our disk, in the corotating frame of the planet. Streamlines proceed from bottom to top inside the planet radius rp = 40 au, and from top to bottom outside the planet radius.
In the text

	Fig. 5 Flow pattern at r = 1.5rp = 60au in an azimuthal cut of the disk. Streamlines flow from right to left; unlike in Fig. 4, vector arrows represent velocity differences from the local initial (quasi-)Keplerian value, rather than from the planet’s Keplerian speed.
In the text

	Fig. 6 Flow pattern at r = 0.66rp = 26.6au in an azimuthal cut of the disk. Streamlines flow from left to right; as in Fig. 5, vector arrows represent velocity differences from the local initial (quasi-)Keplerian value.
In the text

Fig. 7 Source (red) and advective (blue) terms at a fixed r = 1.5rp and θ = 0.2, plotted over azimuthal angle ϕ. From top to bottom, we plot these terms for ρ (in units of Ωpρ0,p), Tg (ΩpTg,0,p), υr, υθ, and υϕ (Ωpcs,iso,p), where the p subscript indicates quantities taken at the planet location in the initial condition. A thin dashed gray line passes through the azimuthal peak of the density spiral, showing the significant offset between terms driving spiral perturbations in each quantity.

In the text

Fig. 8 Background equilibrium, gas, and dust temperatures in our 3T, Saturn-mass simulations at selected altitudes above the midplane for a fixed r = 1.5rp. Background temperatures deviate by only a few percent from the initial condition, most strongly in the region of the planetary shadow (ϕ ≈ π/4) and somewhat more weakly near the Lindblad spiral. Td agrees well with Tg at θ = 0.0 and 0.1 (and is covered by the line for Tg), and with Teq at θ = 0.2 (covering the line for Teq).

In the text

Fig. 9 Measurement of the azimuthal perturbation in each normalized quantity x, averaged over an azimuthal range of ϕpeak ± 2hp according to Eq. (19). Densities and temperatures are plotted in the left panel and velocities in the right. Open circles indicate non-accreting planets, and closed circles indicate those with accretion luminosity (see Sect. 4.2 for a discussion). At z/r ≈ θ = 0.2, where dust and gas are not well coupled, the plotted Tdust amplitude reflects the radial “pseudo-arm” (see the discussion in Sect. 3.4) rather than the Lindblad spiral.

In the text

Fig. 10 Radial velocities, υr, from 3T, zero-accretion-luminosity simulations with Neptune-mass (top row), Saturn-mass (middle row), and Jupiter-mass (bottom row) planets, at various altitudes in the disk (left, middle, and right columns). Gray circles indicate the planetary Hill radius. All υr values are expressed as a function of initial isothermal sound speed at the location of the planet, $c_{\text{iso},\text{p}}=\sqrt{p_{\text{p},0}/{\rho }_{\text{p},0}}$

In the text

	Fig. 11 Gas temperature perturbation for Neptune-mass (top row), Saturn-mass (middle row), and Jupiter-mass (bottom row) planets, at various altitudes in the disk (left, middle, and right columns). Gray circles indicate the planetary Hill radius. All temperature values are expressed with respect to the initial condition Tg,0.
In the text

	Fig. 12 Difference from the initial condition in υr (above) and Tg (below) for 3T simulations of Saturn-mass planets without and with accretion luminosity. The strongly luminous planet alters the vertical and azimuthal structure of the circumplanetary region, causing strong shadows behind the planet and greatly enhancing the kinematic signatures of accretion in the outer disk.
In the text

Fig. 13 Density perturbation in a cut through ϕρ = π/4 as a function of cylindrical radius, R, and vertical position, ɀ, for 3T simulations of Saturn-mass planets without (above) and with (below) accretion luminosity. Gray arrows indicate the velocity field, and the gray circle encloses the planetary Hill radius. When accretion luminosity is included, the circumplanetary envelope becomes larger and the flow around it is altered.

In the text

	Fig. A.1 Sky-projected velocity perturbation, ${\upsilon }_{\perp }^{\prime }$, for our fiducial Saturn-mass planet, in a cut at θ = 0.3 radians above the midplane.
In the text

	Fig. A.2 Sky-projected velocity perturbation, ${\upsilon }_{\perp }^{\prime }$, for a Saturn-mass planet with an accretion luminosity Lacc,p = 10−3L⊙. Compared to the fiducial case, the outer spiral becomes significantly stronger and extends through a larger radial range of the disk.
In the text

	Fig. A.3 Sky-projected velocity perturbation, ${\upsilon }_{\perp }^{\prime }$, for a Jupiter-mass planet with no accretion luminosity. Relative to the fiducial case, circumplan-etary and in-gap flow patterns are significantly stronger.
In the text

	Fig. A.4 Perturbation in gas temperature, Tg − Tg,0, at θ = 0.3 above the midplane for our fiducial setup with a Saturn-mass, non-accreting planet. The double-armed structure of the Lindblad spiral in the upper atmosphere is clearly visible here.
In the text

	Fig. A.5 Gas temperature perturbation for an accreting Saturn-mass planet. As with the kinematic spiral, the thermal spiral extends over a larger radial range and is more prominent than in the non-accreting case.
In the text

	Fig. A.6 Gas temperature perturbation for a Jupiter-mass planet. As with kinematics, the circumplanetary region and gap make a much larger contribution to the thermal signature than in the fiducial Saturn-mass, non-accreting case.
In the text