Three-temperature radiation hydrodynamics with PLUTO: Thermal and kinematic signatures of accreting protoplanets

In circumstellar disks around young stars, the gravitational influence of nascent planets produces telltale patterns in density, temperature, and kinematics. To better understand these signatures, we first performed 3D hydrodynamical simulations of a 0.012 $M_{\odot}$ disk, with a Saturn-mass planet orbiting circularly in-plane at 40 au. We tested four different disk thermodynamic prescriptions (in increasing order of complexity, local isothermality, $\beta$-cooling, two-temperature radiation hydrodynamics, and three-temperature radiation hydrodynamics), finding that $\beta$-cooling offers a reasonable approximation for the three-temperature approach when the planet is not massive or luminous enough to substantially alter the background temperature and density structure. Thereafter, using the three-temperature scheme, we relaxed this assumption, simulating a range of different planet masses (Neptune-mass, Saturn-mass, Jupiter-mass) and accretion luminosities (0, $10^{-3} L_{\odot}$) in the same disk. Our investigation revealed that signatures of disk-planet interaction strengthen with increasing planet mass, with circumplanetary flows becoming prominent in the high-planet-mass regime. Accretion luminosity, which adds pressure support around the planet, was found to weaken the midplane Doppler-flip, potentially visible in optically thin tracers like C$^{18}$O, while strengthening the spiral signature, particularly in upper disk layers sensitive to thicker lines, like those of $^{12}$CO.


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 12 CO, 13 CO, and C 18 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(Teague et al. , 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 offhave 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 Lindblad 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, Ω −1 K ).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; A213, page 1 of 22 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.

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): in which ρ, u, and p represent the gas density, velocity, and pressure, respectively.ρ d is the dust density, while f d is the dustto-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 irr g and S irr d represent, respectively, the gas and dust absorption of stellar irradiation.c represents the speed of light, whereas the ĉ 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 where G ≡ G g + G d and β ≡ u/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, fullstep 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 A213, page 2 of 22 Muley, D., et al.: A&A, 687, A213 (2024) implicitly with a Newton-Raphson method.The stellar irradiation source term S irr d is computed explicitly by ray-tracing at the start of each hydrodynamic timestep (with S irr g 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, 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/∂ ln ρ) entropy is a constant, implying an internal energy density: 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. 2007Boley et al. , 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.

Three-temperature simulations
For our 3T simulations, we solved the full set of Eq. ( 1) and defined the dust-gas collision term, 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 d c 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): where we set the "accommodation coefficient," η, to unity.

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 −10 yr.

β cooling simulations
For beta-cooling simulations, we replaced X gd with a term of the form 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 irr g terms and eliminated Eq. (1d), to (1f) entirely.Following Bae et al. (2021) andMelon Fuksman et al. (2024b), we set the cooling (more accurately, thermal relaxation) time as 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: where the radiative diffusion coefficient D ≡ 4ca r T 3 g /3c g κ R ρ 2 .The effective optically thin cooling length scale λ thin = (3κ R κ P ρ 2 ) −1/2 , while the thick diffusion length is assumed to equal the local scale height, λ diff ≡ H = c s,iso Ω −1 , where Ω is the local Keplerian orbital frequency and c s,,iso ≡ p/ρ the "isothermal sound speed".

Locally isothermal simulations
To test the locally isothermal case, we simply ran a β cooling simulation with t rel = 10 −10 yr everywhere in the disk.

Implementation of accretion luminosity
We implemented accretion luminosity as part of the irradiation source term, S irr d .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 triangleshaped cloud approach (e.g., Eastwood 1986), in which the kernel takes the form 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 1 ; the fraction of accretion energy deposited in each cell during a timestep is equal to the integral of ψ over the cell volume.The triangleshaped cloud approach avoids creating spurious discontinuities in the deposited radiation field or its gradient, unlike the cloudin-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).

Setup and initial conditions
In all simulations presented here, we assumed a disk gas surface density profile of 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.32 The total dustto-gas ratio was assumed to be 1%, of which we took 10% by mass (corresponding to f d = 10 −3 ) 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 −3 and a specific heat capacity c d = 0.7 J g −1 K −1 ; 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 ϕ p = π/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 gasgrain 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 −4 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 selfshadowing (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.

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 (v r ) for each of our four thermodynamic prescriptions (to recapitulate, local isothermality, β 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 3 .
In the midplane, we find that perturbations in ρ and v 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 ρ, v 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 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.5r p 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.
disk, this is typically comparable to or shorter than the spiralcrossing 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., 12 CO) may not be reflected in the dust.

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 −1 ) so the Wentzel-Kramers-Brillouin (WKB) approximation holds, we obtain the following dispersion relation: 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: where r m = (1 ∓ m −1 ) 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 A213, page 5 of 22 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.
N z ≡ gγ −1 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 g ≈ Ω 2 z and assuming a background in hydrostatic equilibrium, one can simplify the above to where c s,iso,0 ≡ p 0 /ρ 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: 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 th ≡ (M p /M ⊙ )h −3 p 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 rel ≲ N −1 z , 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 .

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 midplane, 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 ρ 0 , T g,0 , u 0 , where u 0 = v ϕ,0 φ is a divergenceless, axisymmetric, quasi-Keplerian flow set by the stellar potential Φ * and the initial pressure profile P 0 .All quantities A213, page 6 of 22  (ρ, T g , T d , v r , v θ , and v ϕ ) can be expressed in terms of the initial condition and a perturbation, for example ρ = ρ 0 + ρ ′ ; 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: ∂u 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 advection of the flow.We considered both the vertical (formally, along the θ direction) and in-plane (formally, along the r and θ 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 The final bracketed term in Eq. ( 16b) represents the exchange of energy between gas and dust.The term 2Ω p × u ′ 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, v θ ) toward the midplane -against ∇ρ 0 , but along ∇T g,0 -decreasing the local T g while increasing ρ (θ advection).Because this flow pattern (which also includes a perturbation in v r ) has a nonzero divergence, it also increases ρ and T g via compression and pdV work, respectively.These terms are A213, page 7 of 22 Fig. 7. Source (red) and advective (blue) terms at a fixed r = 1.5r p and θ = 0.2, plotted over azimuthal angle ϕ.From top to bottom, we plot these terms for ρ (in units of Ω p ρ 0,p ), T g (Ω p T g,0,p ), v r , v θ , and v ϕ (Ω p c s,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.
balanced by quasi-Keplerian transport of gas through the spiral pattern ("in-plane" advection), and additionally for T g , collisional 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 v r and v θ components are governed primarily by the planetdriven density wave (pressure-gravity) and counterbalanced by in-plane transport across the spiral.Along the ϕ-component, the pressure gradient is weak, and v ϕ is governed instead by a balance between in-plane advection and Coriolis terms.The fact that v r and v phi are somewhat out of balance reflects the fact that over the long term, planet-driven spirals open a gap in the disk.

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 weak4 sources of radiation.All this impacts the equilibrium temperature T eq , defined by the equation 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 T eq , 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 midplane 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.3provides 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.

Parameter study
We tested three different planet masses (M p = {5 × 10 −5 , 3 × 10 −4 , 1 × 10 −3 M ⊙ , corresponding to Neptune-, Saturn-, and Jupiter-mass, respectively) 5 with no accretion luminosity, and two planetary accretion luminosities (L acc,p = {0, 1 × 10 −3 }) 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.

Planet mass
In Fig. 9, we plot the spiral-averaged perturbation amplitude ⟨∆x⟩ ϕ,spiral -where x can be any one of various normalized 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.5r p .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.T d agrees well with T g at θ = 0.0 and 0.1 (and is covered by the line for T g ), and with T eq at θ = 0.2 (covering the line for T eq ).
quantities (ln ρ/ρ 0 , ln T g /T g,0 , ln T d /T d,0 , v r /c s,iso,p , v θ /c s,iso,p , v ϕ / 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 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 v r /c s,iso,p perturbations follow each other closely; both weaken in the upper atmosphere, but the density perturbation much more so.v θ /c s,iso,p is nearly zero in the disk midplane, but rises with increasing altitude (see Sect. 3.3).v ϕ /c s,iso,p is relatively weak in most cases, but is stronger for highmass 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 v 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 |v r | ≲ 0.1c s,iso,p -corresponding to a physical velocity of ≲30m s −1 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 v 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. 2015Fung et al. , 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 circumplanetary 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 (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 Saturnmass, but less so at Jupiter-mass, where the T g Lindblad spiral is distorted by the effects of circumplanetary and turbulent flows.

Accretion luminosity
In recent years, the effect of planetary accretion luminosity has been studied in parametrized 2D global disk simulations (Montesinos et al. 2015(Montesinos et al. , 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 −3 L ⊙ we used corresponds to a mass accretion rate of Ṁ = 7 M J Myr −1 , given the fiducial planet mass M p = 4 × 10 −4 M ⊙ and assuming a typical radius of A213, page 9 of 22 Fig. 9. Measurement of the azimuthal perturbation in each normalized quantity x, averaged over an azimuthal range of ϕ peak ± 2h p 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 T dust amplitude reflects the radial "pseudo-arm" (see the discussion in Sect.3.4) rather than the Lindblad spiral.
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 v 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 singlebanded 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.

Observational implications
Of the cases we tested, we find that (non-accreting) Saturnmass 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  For an 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, v ′ ⊥ increasingly reflects v ′ r (and to a lesser extent v ′ ϕ ).Because the typical magnitude of v ′ r in spirals is typically much larger than v ′ θ , 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 v ′ ⊥ does not change with projection.Furthermore, as discussed in Sects.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 hydrodynamical 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 v ′ ϕ 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 themalong with Monte Carlo radiative-transfer post-processing and synthetic observations -in a forthcoming work.

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, T eq , 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 12 CO -thermal and kinematic signatures are driven largely by local source terms, rather than transport across vertical gradients within the disk.
A213, page 12 of 22 Fig. 12. Difference from the initial condition in v r (above) and T g (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.
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. 2015Fung et al. , 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 idealgas assumption and compute γ as a function of temperature, as well as hydrogen (para-, ortho-, and atomic), helium, and metal fractions (Boley et al. 2007(Boley et al. , 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(Bi et al. , 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.A213, page 21 of 22 Muley, D., et al.: A&A, 687, A213 (2024)

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 axisymmetric in all variables, and that v r (t = 0) = v θ (t = 0) ≡ 0, we can write the evolution of the density perturbation ρ 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: 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 (∂u ′ /∂t) adv,in−plane = (∂u ′ /∂t) advr + (∂u ′ /∂t) advϕ .

Fig. 4 .
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 r p = 40 au, and from top to bottom outside the planet radius.

Fig. 5 .
Fig. 5. Flow pattern at r = 1.5r p = 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 .
Fig. 6.Flow pattern at r = 0.66r p = 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. 10 .
Fig. 10.Radial velocities, v r , from 3T, zero-accretion-luminosity simulations with Neptune-mass (top row), Saturn-mass (middle row), and Jupitermass (bottom row) planets, at various altitudes in the disk (left, middle, and right columns).Gray circles indicate the planetary Hill radius.All v r values are expressed as a function of initial isothermal sound speed at the location of the planet, c iso,p = p p,0 /ρ 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 T g,0 .
Fig. 13.Density perturbation in a cut through ϕ p = π/4 as a function of cylindrical radius, R, and vertical position, z, 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.

Fig
Fig. A.2. Sky-projected velocity perturbation, v ′ ⊥ , for a Saturn-mass planet with an accretion luminosity L acc,p = 10 −3 L ⊙ .Compared to the fiducial case, the outer spiral becomes significantly stronger and extends through a larger radial range of the disk.