EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 705 (January 2026) A&A, 705 (2026) A194 Full HTML
Open Access
Issue
A&A
Volume 705, January 2026
Article Number A194
Number of page(s) 12
Section Planets, planetary systems, and small bodies
DOI https://doi.org/10.1051/0004-6361/202557399
Published online 23 January 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Several recent observational efforts have demonstrated that substructure in the form of rings and gaps in continuum emission is ubiquitous in disks around ~Myr-old young stellar objects (e.g., ALMA Partnership 2015; Andrews et al. 2018; Teague et al. 2025). While several different mechanisms have been proposed to explain such substructure, young exoplanets present a particularly attractive scenario, as the planet formation process likely begins quite early in these systems (Birnstiel 2024; Kleine et al. 2020), and many of the observed features are well reproduced with models of planet–disk interaction (e.g., Zhang et al. 2018; Bae et al. 2019; Toci et al. 2020).

With the number of confirmed observed planets during the disk phase standing at about three (Keppler et al. 2018; Haffert et al. 2019; van Capelleveen et al. 2025), theoretical modeling is by and large the main avenue to study planet–disk interaction. A young planet interacts gravitationally with its surrounding disk, driving spiral density waves that propagate through the disk (Goldreich & Tremaine 1979; Rafikov 2002a; Ogilvie & Lubow 2002). As they travel, these density waves evolve due to nonlinear effects, eventually turning into shocks, delivering angular momentum to the disk and carving gaps around the planet’s orbit (Goodman & Rafikov 2001; Rafikov 2002b).

Models of gap-opening planets have increased in sophistication over the years, highlighting the effects of radiative cooling (Miranda & Rafikov 2020a,b; Ziampras et al. 2020b, 2023b), dust dynamics (Weber et al. 2019; Bi et al. 2021, 2023), and even the dependence of cooling on the local dust content (Binkert et al. 2023; Krapp et al. 2024; Ziampras et al. 2025b). In particular, it has been shown that planet-driven gap opening is highly sensitive to the cooling timescale (Miranda & Rafikov 2020a; Zhang & Zhu 2020; Zhang et al. 2024) as well as radiative diffusion across the planetary spiral shocks (Miranda & Rafikov 2020b; Ziampras et al. 2023b). Combined with the complicated opacity models in the literature (Semenov et al. 2003; Birnstiel et al. 2018; Woitke et al. 2016), surveying a wide parameter space while maintaining accuracy when attempting to model an observed system with a substructure-forming planet can be a lost cause.

While the aforementioned models are of key importance to understanding planet–disk interaction, they are also both sensitive to several – often poorly constrained – physical parameters (e.g., gas density, dust composition, dust-to-gas ratio) and computationally expensive, such that sweeping the parameter space with 3D models is unfeasible. To that end, vertically integrated, 2D models are quite attractive, as their lower computational cost enables the execution of large parameter studies at high resolution.

Of course, the question remains whether such 2D models can accurately capture the physics of gap opening in the first place. The flow around a gap-opening planet is inherently three dimensional, featuring zonal flows and meridional circulation (Morbidelli et al. 2014; Fung & Chiang 2016; Bi et al. 2021; Wafflard-Fernandez & Lesur 2023). At the same time, while extensive work has been carried out to verify that planet–disk interaction can be mapped from 3D to 2D (e.g., Müller et al. 2012; Fung & Chiang 2016; Lega et al. 2021; Hammer & Lin 2023; Zier & Springel 2023; Brown & Ogilvie 2024; Cordwell et al. 2025), radiative processes complicate this picture, as they are very sensitive to the optical depth and cooling timescale, both of which are strong functions ofz (e.g., Bae et al. 2021; Ziampras et al. 2024a). As such, a comparison of 2D to 3D models while including radiation hydrodynamics is necessary to determine whether and under what conditions the 2D, vertically integrated cooling prescriptions can be reconciled with fully 3D models.

In this work, we perform 2D and 3D radiation hydrodynamical simulations of planet–disk interaction, comparing the gap-opening efficiency of the two frameworks for the same physical parameters. We re-used some of the 2D models in Ziampras et al. (2023b) but re-analyzed them with a more appropriate method following Cordwell & Rafikov (2024) and supplement them with fully 3D simulations, which we analyzed following Cordwell et al. (2025). Our goals are to both determine whether 2D models can accurately capture the gap-opening efficiency of 3D models, and to identify cooling prescriptions that can approximate the behavior of fully radiative models with good accuracy.

The paper structure is as follows: in Sect. 2 we detail our physical and numerical frameworks in both 2D and 3D, and we define effective cooling timescales that we found to be useful in approximating the behavior of fully radiative models. We also define our simulation metrics, which we used to compare the two frameworks. In Sect. 3 we present our results, where we compare the gap-opening efficiency of2D and 3D models. We then discuss the implications of our findings in Sect. 4, and conclude in Sect. 5.

2 Physics and numerics

In this section we describe our physical framework and numerical approach. While both are largely similar to Ziampras et al. (2023b) and Ziampras et al. (2024a) for 2D and 3D respectively, we reiterate their description here for clarity and to facilitate the derivation of equivalent cooling timescales between 2D and 3D in Sect. 2.3.

2.1 Three-dimensional framework

We consider an inviscid disk of gas with volume density ρ, internal energy density e and velocity field u around a star with mass M* = 1 M⊙ and luminosity L* = 1 L⊙. The gas has a mean molecular weight μ = 2.353 and an adiabatic index γ = 7/5. The inviscid Navier–Stokes equations read as (e.g., Tassoul 1978)

∂ρ∂t+u⋅∇ρ=−ρ∇⋅u,${{{\partial \rho } \over {\partial t}} + u \cdot \nabla \rho = - \rho \nabla \cdot u,}$(1a)

∂u∂t+(u⋅∇)u=−1ρ∇P−∇Φ,${{{\partial u} \over {\partial t}} + (u \cdot \nabla )u = - {1 \over \rho }\nabla P - \nabla {\rm{\Phi }},}$(1b)

∂e∂t+u⋅∇e=−γe∇⋅u+Q.${{\partial e} \over {\partial t}} + u \cdot \nabla e = - \gamma e\nabla \cdot u + Q.$(1c)

Here, the pressure is given by the perfect gas equation of state with P = (γ – 1)e. Through the latter we can define the isothermal sound speed cs=P/ρ${c_{\rm{s}}} = \sqrt {P/\rho } $, the temperature T=μcs2/ℛ$T = \mu c_{\rm{s}}^2/{\cal R}$, the pressure scale height H = cs/ΩK, and the aspect ratio h = H/R, with ΩK=GM⋆/R3${\Omega _{\rm{K}}} = \sqrt {{\rm{G}}{M_ \star }/{R^3}} $ being the Keplerian angular speed at cylindrical distance R from the star. The perfect gas and gravitational constants are ℛ and G, respectively. Finally, Φ = Φ* + Φp is the gravitational potential of the star and planet with mass Mp and radial location Rp, given by Φ* = -GM*/r and

 Φp={ −GMp∈[ d3∈3−2d2∈2+2 ],−GMp/d, d≥ ∈d< ∈,$\,{\Phi _{\rm{p}}} = \left\{ {_{ - {{G{M_p}} \over \in }\left[ {{{{d^3}} \over {{ \in ^{\rm{3}}}}} - 2{{{d^2}} \over {{ \in ^{\rm{2}}}}} + 2} \right],}^{ - {\rm{G}}{M_p}/d,}} \right.\matrix{ {} & {d \ge \, \in } \cr {} & {d < \, \in ,} \cr } $(2)

in order to avoid singularities, following Klahr & Kley (2006). Here, ∈ = 0.5 RH = 0.5 Rp ϵ=0.5RH=0.5RpMp/3M⋆3$ = 0.5{R_{\rm{H}}} = 0.5{R_{\rm{p}}}\root 3 \of {{M_{\rm{p}}}/3{M_ \star }} $, and d = r - Rp is the distance between a point at position r and the planet. This choice of smoothing length ∈, while larger than that used in Cordwell et al. (2025), is still small enough to not affect the flow in our region of interest. We accounted for the indirect acceleration arising from the star–planet system orbiting their mutual center of mass.

In Eq. (1c), Q represents additional source terms that can affect the thermal energy evolution. We defined four different thermodynamical treatments similar to Ziampras et al. (2023b):

  • Locally isothermal: Eq. (1c) is ignored altogether, and a temperature profile T0 (r) is prescribed such that P = ℛρT0 /μ. This is the most simplistic approach within our framework, effectively translating to immediate thermal relaxation.

  • Adiabatic: Eq. (1c) is evolved with Q = 0. This approach captures adiabatic compression while conserving entropy (except at shock fronts), but assumes no thermal relaxation.

  • Local cooling: thermal relaxation is parametrized with a cooling timescale tcool=β ΩK−1${t_{{\rm{cool}}}} = \beta \,\Omega _{\rm{K}}^{ - 1}$, such that in Eq. (1c)

    Q=Qrelax=−fβe−e0tcool, e0=ℛμ(γ−1)ρT0=cvρT0,$Q = {Q_{{\rm{relax}}}} = - {f_\beta }{{e - {e_0}} \over {{t_{{\rm{cool}}}}}},\quad {e_0} = {{\cal R} \over {\mu (\gamma - 1)}}\rho {T_0} = {c_{\rm{v}}}\rho {T_0},$(3)

    with fβ being a prefactor that differs from unity for temperature-dependent β (Dullemond et al. 2022; Ziampras et al. 2023b).

  • Fully radiative: the disk cools by interacting with the local radiation field Erad in the flux-limited diffusion approximation (FLD, Levermore & Pomraning 1981). Here, Erad is dynamically evolved and added to the set of equations in Eq. (1) with

    ∂Erad∂t+∇⋅F=−Qrad, F=−λcκRρ∇Erad,${{\partial {E_{{\rm{rad}}}}} \over {\partial t}} + \nabla \cdot F = - {Q_{{\rm{rad}}}},\quad F = - {{\lambda c} \over {{\kappa _{\rm{R}}}\rho }}\nabla {E_{{\rm{rad}}}},$(4)

    and the matter–radiation field coupling is given by defining in Eq. (1c)

    Q=Qrad=−κPρc(aRT4−Erad).$Q = {Q_{{\rm{rad}}}} = - {\kappa _{\rm{P}}}\rho c\left( {{a_{\rm{R}}}{T^4} - {E_{{\rm{rad}}}}} \right).$(5)

    In this context, κR and κP are the Rosseland and Planck mean opacities, c is the speed of light, aR is the radiation constant, and λ is a flux limiter following Kley (1989) that captures the transition between the optically thick diffusion limit (λ → 1/3) and the optically thin, free-streaming limit (F → cErad). This is the most realistic approach within our framework, capturing both the thermal relaxation of e through Qrad and the radiative diffusion of Erad through Eq. (4). To prevent the disk from cooling indefinitely, we prescribed Erad = aR T0(R)4 at the surfaces and radial ends of the disk as in Ziampras et al. (2024a), assuming that the disk is heated by reprocessed starlight from the disk surfaces.

2.2 Two-dimensional framework

Similar to the set of equations in Eq. (1) and by defining the surface density Σ=∫−∞∞ρdz$\Sigma = \int_{ - \infty }^\infty \rho {\rm{d}}z$ and vertically integrated pressure P2D=Σcs2${P_{2{\rm{D}}}} = \Sigma c_{\rm{s}}^2$, (assuming T is independent of z) the inviscid Navier- Stokes equations in the thin disk approximation read

∂Σ∂t+u⋅∇Σ=−Σ∇⋅u,${{\partial {\rm{\Sigma }}} \over {\partial t}} + u \cdot \nabla {\rm{\Sigma }} = - {\rm{\Sigma }}\nabla \cdot u,$(6a)

∂u∂t+(u⋅∇)u=−1Σ∇P2D−∇Φ,${{\partial u} \over {\partial t}} + (u \cdot \nabla )u = - {1 \over {\rm{\Sigma }}}\nabla {P_{2{\rm{D}}}} - \nabla {\rm{\Phi }},$(6b)

∂e∂t+u⋅∇e=−γe∇⋅u+Q.${{\partial e} \over {\partial t}} + u \cdot \nabla e = - \gamma e\nabla \cdot u + Q.$(6c)

In the above equations the thermal energy density has been redefined as e = P2D/(γ - 1). To capture the vertical structure of the disk, the gravitational potential of the planet is also replaced by a Plummer potential with

Φp=−GMpd2+∈2,∈=0.6Hp,$\matrix{ {{{\rm{\Phi }}_{\rm{p}}} = - {{{\rm{G}}{M_{\rm{p}}}} \over {\sqrt {{d^2} + { \in ^2}} }},} & { \in = 0.6{H_{\rm{p}}},} \cr } $(7)

following Müller et al. (2012). While Cordwell et al. (2025) showed that the newer potential form of Brown & Ogilvie (2024) functions more optimally than the traditional Plummer potential, the main purpose of this paper is to explore the role of disk thermodynamics rather than the planetary potential on gap opening. As we seek to identify how well the overall morphology, rather than quantitative details, of gap opening can be represented in 2D simulations, the Plummer approximation is an appropriate choice.

While Qrelax maintains its form as in Eq. (3) within the local cooling approach (with Σ instead of ρ), the radiative cooling term Qrad in Eq. (4) is now split among a vertical, surface cooling component Qsurf, an in-plane radiative diffusion term Qmid, and an irradiation heating term Qirr. The two cooling terms can be written as

Qsurf=−2σSBT4τeff, τeff=3τR8+34+14τp, τR,P=12κR,PΣ${Q_{{\rm{surf}}}} = - 2{{{\sigma _{SB}}{T^4}} \over {{\tau _{{\rm{eff}}}}}},\quad {\tau _{{\rm{eff}}}} = {{3{\tau _{\rm{R}}}} \over 8} + {{\sqrt 3 } \over 4} + {1 \over {4{\tau _{\rm{p}}}}},\quad {\tau _{{\rm{R}},{\rm{P}}}} = {1 \over 2}{\kappa _{{\rm{R}},{\rm{P}}}}{\rm{\Sigma }}$(8)

for surface cooling following the effective optical depth model of Hubeny (1990), and

Qmid=2πH∇⋅(16λσSBT3κRρmid∇T), ρmid≡12πΣH${Q_{{\rm{mid}}}} = \sqrt {2\pi } H\nabla \cdot \left( {{{16\lambda {\sigma _{{\rm{SB}}}}{T^3}} \over {{\kappa _{\rm{R}}}{\rho _{{\rm{mid}}}}}}\nabla T} \right),\quad {\rho _{{\rm{mid}}}} \equiv {1 \over {\sqrt {2\pi } }}{{\rm{\Sigma }} \over H}$(9)

for in-plane radiative diffusion, where we further assume that Erad = aRT4 at the disk midplane (one-temperature approximation). In the above, σSB is the Stefan–Boltzmann constant and ρmid is the volume density at the disk midplane, obtained assuming vertical hydrostatic equilibrium.

Finally, we included the prescription of irradiation heating by Menou & Goodman (2004) to balance surface losses via Eq. (8):

Qirr=2L⋆4πR2(1−ε)(dlogHdlogR−1)hτeff,${Q_{{\rm{irr}}}} = 2{{{L_ \star }} \over {4\pi {R^2}}}(1 - \varepsilon )\left( {{{{\rm{d}}\log H} \over {{\rm{d}}\log R}} - 1} \right){h \over {{\tau _{{\rm{eff}}}}}},$(10)

with dlogH dlogR=9/7${{{\rm{d}}\log H} \over {{\rm{d}}\log R}} = 9/7$ and a disk albedo ε = 1 /2.

2.3 Effective cooling timescales

One of the goals of this work is to identify cooling timescale prescriptions that can approximate the behavior of fully radiative models with good accuracy, capturing radiative effects through local thermal relaxation (Eq. (3)) rather than solving for the full radiative diffusion problem (Eqs. (4) and (5)). Following previous work (Flock et al. 2017; Ziampras et al. 2024a) we define the cooling timescale within the 3D framework as

β3D=ΩKη(H2+κPκRlrad23), η=16σSBT33κRρ2cv, lrad=1κPρ,${\beta ^{3{\rm{D}}}} = {{{{\rm{\Omega }}_{\rm{K}}}} \over \eta }\left( {{H^2} + {{{\kappa _{\rm{P}}}} \over {{\kappa _{\rm{R}}}}}{{l_{rad}^2} \over 3}} \right),\quad \eta = {{16{\sigma _{{\rm{SB}}}}{T^3}} \over {3{\kappa _{\rm{R}}}{\rho ^2}{c_{\rm{v}}}}},\quad {l_{{\rm{rad}}}} = {1 \over {{\kappa _{\rm{P}}}\rho }},$(11)

where η and lrad are the radiative diffusion coefficient and photon mean free path, respectively. This recipe captures cooling in both the optically thick, diffusion-limited regime (with the assumption that the diffusion length scale is ~H) and in the optically thin regime where the cooling rate is proportional to the emissivity κP. We note that the factor κP/κR is a correction to the formulas provided by the aforementioned studies that arises when relaxing the assumption that κR = κP = κ. A heuristic derivation of Eq. (11) is provided in Appendix B.

In the 2D framework, the above formula is still relevant as it represents cooling via diffusion through the disk plane (defining all parameters at the midplane), effectively capturing the effects of Qmid in Eq. (9). To further include the effects of surface cooling, assuming that heat is removed locally along the z direction, we can write the two components of surface and in-plane cooling as

βmid=β3D(z=0), βsurf=e| Qsurf |ΩK.${\beta _{{\rm{mid}}}} = {\beta ^{3{\rm{D}}}}(z = 0),\quad {\beta _{{\rm{surf}}}} = {e \over {\left| {{Q_{{\rm{surf}}}}} \right|}}{{\rm{\Omega }}_{\rm{K}}}.$(12)

The two can then be combined such that cooling is always limited by the fastest of the two channels (e.g., Miranda & Rafikov 2020b) to yield an effective cooling timescale:

β2D=[ βmid−1+βsurf−1 ]−1.${\beta ^{2{\rm{D}}}} = {\left[ {\beta _{{\rm{mid}}}^{ - 1} + \beta _{{\rm{surf}}}^{ - 1}} \right]^{ - 1}}.$(13)

Similar to Ziampras et al. (2023b), Eq. (13) can be algebraically rearranged to express β2D as

β2D=1f+1βsurf=1f+1e| Qsurf |ΩK, f=16πτpτeff6τRτp+π,${\beta ^{2{\rm{D}}}} = {1 \over {f + 1}}{\beta _{{\rm{surf}}}} = {1 \over {f + 1}}{e \over {\left| {{Q_{{\rm{surf}}}}} \right|}}{{\rm{\Omega }}_{\rm{K}}},\quad f = {{16\pi {\tau _{\rm{p}}}{\tau _{{\rm{eff}}}}} \over {6{\tau _{\rm{R}}}{\tau _{\rm{p}}} + \pi }},$(14)

with f → π and f → 4 in the optically thick and thin limits (τ → ∞ and τ → 0, respectively) if τR = τP = τ.

We note that Eqs. (11) and (13) can be used to both estimate the cooling timescale in radiative models and compute a relaxation rate Qrelax when using a local cooling approach. In the latter case, Dullemond et al. (2022) showed that the prefactor in Eq. (3) takes the form fβ = 1/(4 ± b) in the optically thick (+) and thin (-) limits, where τeff ∝ T±b for temperature-dependent opacities (see Eq. (8)). Nevertheless, in the range τ ~ 0.1–10, where τeff is rather flat due to the constant term 3/4$\sqrt 3 /4$ in Eq. (8), this simplifies to fβ ≈ 1/4 (Ziampras et al. 2023b). We use this value of fβ for our models with local cooling. In Sect. 3, we will show that models using this local cooling approach with our analytic formulas in Eqs. (11) and (13) perform very well when compared to their fully radiative counterparts in the context of planet-driven gap opening.

2.4 Numerical setups

We use the numerical (magneto)hydrodynamics finite-volume code PLUTO (Mignone et al. 2007) in 3D spherical {r, θ, ϕ} or 2D cylindrical {R, ϕ} polar geometry, depending on framework. We use the flux-limited diffusion radiation transport modules detailed in Ziampras et al. (2024a) and Ziampras et al. (2020a) for our 3D and 2D runs, respectively1. The models are run in a frame corotating with the planet’s orbital speed Ωp=G(M⋆+Mp)/Rp3${\Omega _{\rm{p}}} = \sqrt {{\rm{G}}\left( {{M_ \star } + {M_{\rm{p}}}} \right)/R_{\rm{p}}^3} $ and make use of the FARGO transport algorithm (Masset 2000) implemented in PLUTO by Mignone et al. (2012) to both minimize numerical diffusion and alleviate the strict timestep limitation of the rapidly rotating inner radial boundary. Finally, we employ a third-order weighted essentially non-oscillatory reconstruction scheme (WENO, Yamaleev & Carpenter 2009) and third-order Runge–Kutta timestepping, the slope limiter by Van Leer (1974), and the HLLC Riemann solver (Toro et al. 1994).

Our numerical grid spans R or r ∈ [0.4, 2.5] Rp with logarithmic radial spacing, covers the full ϕ ∈ [0, 2π] in azimuth, and, in our 3D models, extends between θ ∈ [π/2 - 0.15, π/2] in the polar direction. For our choice of temperature profile with h = 0.05 (R/Rp)1/4 (see text below), this translates to zmax = 3H at R = Rp. Our boundaries are periodic in ϕ, symmetric about the midplane, and we impose wave-damping regions for R < 0.5 Rp, R > 2.1 Rp following de Val-Borro et al. (2006) with a damping timescale of 0.1 local orbits.

After carrying out a resolution study (see Appendix A), we choose a grid resolution of NR × Nϕ = 1200 × 4096 and Nr × Nθ × Nϕ = 600 × 48 × 2048 cells for our 2D and 3D models, which translates to 32 and 16 cells per scale height, respectively. This choice does not sacrifice quality for performance, however. Rather, it shows that our 3D models converged within 95% for half the effective resolution compared to our 2D runs with respect to our metrics of choice.

Similar to Ziampras et al. (2023b), we aim to build a set of models where τ, Σ, and β are constant within our simulation domain (for 3D models, we seek to satisfy this at the midplane). This is achieved by setting κR = κP = κ = const., Σ0 = Σ(t = 0) = const., and T0 = T(t = 0) = Tp (R/Rp)-1/2. For a solar-type star, thermal equilibrium between Qsurf and Qirr in Eqs. (8) and (10) translates to Tp = 21 K at Rp = 30 au, or hp = 0.05. We then use an iterative approach to compute κ and Σ0 to achieve different values of β2D, or β3D at the disk midplane while maintaining an optical depth τ=12κΣ=10$\tau = {1 \over 2}\kappa \Sigma = 10$. In sections where we compare our findings as a function of the cooling timescale, we use β2D as a reference value for the disk cooling time, regardless of the equation of state. This is done because it represents the total cooling timescale in our 2D runs, whereas in 3D the effective cooling timescale might be much shorter than β3D(z = 0) due to the entire column from the midplane to the vertical boundary contributing to cooling.

We note that the actual values of κ and Σ0, however unrealistic, are completely arbitrary. The goal is to establish numerical experiments where τ and β are (initially) constant throughout the disk midplane in order to measure the dissipation of planet-driven spirals as a function of the cooling timescale while maintaining a marginally optically thick disk (τmid = 10). While this is trivial for models with local cooling, where the value of β can be simply prescribed, running equivalent radiative models requires defining our initial conditions in this way such that β evaluates to its target value through Eq. (11) or (13) in 3D and 2D, respectively.

Finally, we chose a planet mass of Mp = 3.75 × 10-5 M* = 12.5 M⊕. For our choice of parameters, we defined the thermal mass following (Goodman & Rafikov 2001) as Mth=23h3M⋆≈28M⊕${M_{{\rm{th}}}} = {2 \over 3}{h^3}{M_ \star } \approx 28{{\rm{M}}_ \oplus }$, such that Mp = 0.42 Mth2. This places the planet in a quasi-nonlinear regime, where it is expected to carve a gap due to the lack of viscous diffusion, but its spiral wakes will not steepen into shocks before traveling a radial distance of (Goodman & Rafikov 2001)

xsh≈0.93Hp(γ+112/5MpMth)−2/5${x_{{\rm{sh}}}} \approx 0.93{H_{\rm{p}}}{\left( {{{\gamma + 1} \over {12/5}}{{{M_{\rm{p}}}} \over {{M_{{\rm{th}}}}}}} \right)^{ - 2/5}}$(15)

away from the planet. This allows us to measure differences in gap opening between isothermal, adiabatic, and radiative models, as no angular momentum should be deposited in the region |R - Rp| < xsh in the absence of radiative damping (Rafikov 2002a). The planet is ramped to its full mass during its first orbit in the simulation following the formula in de Val-Borro et al. (2006). This ensures that the spiral patterns have been established throughout the computational domain and have reached a steady state within a few (~5) planetary orbits.

2.5 Simulation metrics

As previous studies have already shown, the gap-opening process is a sensitive function of the cooling timescale. To quantify this effect, we focused on the planet’s primary gap, even though multiple gaps can in principle be opened in inviscid disks (e.g., Zhang et al. 2018). We therefore need a metric to measure the efficiency with which the planet removes material from its local orbit. Rather than computing an “angular momentum budget” as in Ziampras et al. (2023b), who integrated the total spiral angular momentum flux FJ throughout the disk, we instead computed the angular momentum deposition function ∂Fdep /∂R following Cordwell & Rafikov (2024). This can be defined as

∂Fdep∂R=∂T∂R−∂FJ∂R,${{\partial {F_{{\rm{dep}}}}} \over {\partial R}} = {{\partial T} \over {\partial R}} - {{\partial {F_J}} \over {\partial R}},$(16)

where the angular momentum flux FJ and gravitational torque density ∂T /∂R are given in the 2D framework by

FJ(R)=R2∮ΣδuRδuϕdϕ, δux=ux−〈 ux 〉,${F_J}(R) = {R^2}\oint {\mkern 1mu} {\mkern 1mu} {\rm{\Sigma }}\delta {u_R}\delta {u_\phi }{\rm{d}}\phi ,\quad \delta {u_x} = {u_x} - \left\langle {{u_x}} \right\rangle ,$(17a)

∂T∂R=R∮Σ∂Φp∂ϕdϕ,${{\partial T} \over {\partial R}} = R\oint {\mkern 1mu} {\mkern 1mu} {\rm{\Sigma }}{{\partial {{\rm{\Phi }}_{\rm{p}}}} \over {\partial \phi }}{\rm{d}}\phi ,$(17b)

and in the 3D framework computed after first integrating along z (Cordwell et al. 2025):

FJ(R)=R2∮[ ∫−∞∞ρδuRδuϕdz ]dϕ,${F_J}(R) = {R^2}\oint {\left[ {\int_{ - \infty }^\infty {\rho \delta {u_R}\delta {u_\phi }{\rm{d}}z} } \right]} {\mkern 1mu} {\mkern 1mu} {\rm{d}}\phi ,$(18a)

∂T∂R=R∮[ ∫−∞∞ρ∂Φp∂ϕdz ]dϕ.${{\partial T} \over {\partial R}} = R\oint {\mkern 1mu} {\mkern 1mu} \left[ {\int_{ - \infty }^\infty {\rho {{\partial {{\rm{\Phi }}_{\rm{p}}}} \over {\partial \phi }}{\rm{d}}z} } \right]{\rm{d}}\phi .$(18b)

In the above equations, 〈ux〉 denotes density-weighted azimuthal averages such that

〈 ux 〉=12πρ¯∮uxρdϕ,$\left\langle {{u_x}} \right\rangle = {1 \over {2\pi \bar \rho }}\oint {\mkern 1mu} {\mkern 1mu} {u_x}\rho {\rm{d}}\phi ,$(19)

with ρ¯$\bar \rho $ being the azimuthally averaged density (ρ and ρ¯$\bar \rho $ are replaced with Σ and Σ¯$\bar \Sigma $ in 2D). After computing ∂Fdep/∂R, we can then reconstruct the azimuthally averaged perturbed density ΔΣ¯/Σ0$\Delta \bar \Sigma /{\Sigma _0}$ as a function of time by computing the surface density evolution rate ∂ς/∂t using ∂Fdep/∂R, following Eqs. (C1)–(C12) in Appendix C of Cordwell & Rafikov (2024)3:

ΔΣ¯Σ0(t)≈t1Σ0∂Σ∂t.${{{\rm{\Delta \bar \Sigma }}} \over {{{\rm{\Sigma }}_0}}}(t) \approx t{1 \over {{{\rm{\Sigma }}_0}}}{{\partial {\rm{\Sigma }}} \over {\partial t}}.$(20)

Here, Σ0 is the initial surface density profile. We note that this approach assumes that the gap-opening process operates in its linear regime, where ∂ς/∂t is constant in time. This is a very accurate approximation for gap depths less than 20% (Cordwell & Rafikov 2024), a condition met by the circumstances considered in this paper. As we are interested in the inviscid regime of gap opening we are unable to make direct comparisons to standard viscous steady state models such as in Kanagawa et al. (2015), Zhang et al. (2018), or Duffell (2020).

In all figures shown below, ∂Fdep /∂R is normalized to

FJ,0Rp=(MpM⋆)2hp−3ΣpRp3Ωp2,${{{F_{J,0}}} \over {{R_{\rm{p}}}}} = {\left( {{{{M_{\rm{p}}}} \over {{M_ \star }}}} \right)^2}h_{\rm{p}}^{ - 3}{{\rm{\Sigma }}_{\rm{p}}}R_{\rm{p}}^3{\rm{\Omega }}_{\rm{p}}^2,$(21)

where FJ0 is the characteristic one-sided torque due to the planet (Goldreich & Tremaine 1980).

For our analysis in the following sections, we will focus on the angular momentum deposition due to spiral shocks rather than dynamics within the planet’s corotating region. The latter should not contribute to spiral-driven gap opening in the absence of radiative damping. However, dynamical processes such as the vertical shear instability (VSI, Nelson et al. 2013; Stoll et al. 2017; Ziampras et al. 2023a), buoyancy modes (Zhu et al. 2012; McNally et al. 2020), or transient features while the horseshoe region is undergoing phase mixing (e.g., the vortensity striping seen in Cordwell et al. (2025)) can still drive correlated velocity fluctuations in that region. In particular, the VSI will operate for very short cooling timescales (β ≲ 0.1, Lin & Youdin 2015), while the disk’s buoyancy response can drive substantial motion for nearly adiabatic conditions (β ≫ 100, Yun et al. 2022; Ziampras et al. 2024a). To avoid the effects of either of these processes, we set ∂Fdep /∂R = 0 within the planet’s corotating region |R - Rp| < xh, with the horseshoe half-width xh defined following Paardekooper et al. (2010):

xh=1.1γ1/4(0.4ϵ/H)1/4MphM⋆Rp, ϵ=0.6Hp.${x_{\rm{h}}} = {{1.1} \over {{\gamma ^{1/4}}}}{\left( {{{0.4} \over {/H}}} \right)^{1/4}}\sqrt {{{{M_{\rm{p}}}} \over {h{M_ \star }}}} {R_{\rm{p}}},\quad = 0.6{H_{\rm{p}}}.$(22)

The radial profile of ∂Fdep /∂R is then smoothed with a rolling average over a width of H to eliminate any numerical noise, and ∂ς/∂t is computed using this smoothed data. We found that this smoothing has a negligible effect on the results, but improves the readability of the figures. Nevertheless, we show the raw data in faded colors in Figs. 1-3 for transparency.

3 Results

In this section we present the results of our numerical experiments. We begin with a demonstration of our approach, continue to an analysis of the angular momentum deposition function and the resulting gap structures as a function of the cooling timescale, and conclude with a comparison between our 2D and 3D models.

3.1 A fiducial model

To illustrate our approach, we show a snapshot of the surface density structure after 10 planetary orbits for our locally isothermal 3D model in Fig. 1. On the left panel we show the entire simulation domain, with the planet located at R = R0 and spirals permeating the disk. The top right panel shows a zoomed-in view of the region around the planet (indicated with an orange box on the left), highlighting the density perturbations around the spiral shocks. The bottom right panel shows the angular momentum deposition function ∂Fdep /∂R in normalized units, and the surface density evolution rate in arbitrary units computed based on this ∂Fdep /∂R according to Cordwell & Rafikov (2024). Dashed vertical lines mark the shock distance xsh following Eq. (15), and dotted lines mark the edge of the planet’s corotating region xh following Eq. (22).

From the bottom right panel of Fig. 1 we can identify several key features of planet-driven gap opening by noticing that the angular momentum deposition function is – within noise levels – zero between the edge of the planet’s corotating region and the shock distance. This is expected, as the spirals have not yet steepened into shocks and travel as linear waves within this region4. Beyond xsh, however, we see a rise in ∂Fdep /∂R and a corresponding increase in the surface density evolution rate ∂Δς/∂t. The resulting “double-trough” gap shape is entirely consistent with expectations from literature for low-viscosity, locally isothermal models (e.g., Rafikov 2002b; Miranda & Rafikov 2020a; Zhang et al. 2024; Cordwell & Rafikov 2024).

This behavior no longer necessarily holds when a finite cooling timescale is introduced, as radiative damping can now drive angular momentum deposition well before xsh is reached. This is demonstrated in Fig. 2, where we show the angular momentum deposition function ∂Fdep /∂R for our fiducial 2D and 3D models (dashed and solid lines, respectively) and for four different equations of state: adiabatic (black), locally isothermal (orange), β cooling (purple), and fully radiative (green). Our models with cooling (β, radiative) have β2D ≈ 2. Here, it becomes clear that a finite cooling timescale leads to a nonzero angular momentum deposition within xsh. In particular, by comparing the purple and orange curves in the bottom panel of Fig. 2, we see that ∂Fdep /∂R not only peaks much closer to the planet, but also remains high up to ~2xsh away, suggesting the formation of a single, deep gap around the planet. This is fully in line with the findings of Miranda & Rafikov (2020a).

In Fig. 3 we then show the expected surface density profile after 100 planetary orbits by proxy of the surface density evolution rate through Eq. (20) for the same models as in Fig. 2. In this figure, and especially in the bottom panel, we can see that the gap profile is significantly different between the different models: the locally isothermal and adiabatic models (orange and black lines, respectively) feature a shallow, double-trough gap profile, while the models with cooling (purple and green lines) show a single, deep gap, which is slightly deeper for β cooling (purple). This is once again fully consistent with findings from previous studies (Miranda & Rafikov 2020b; Zhang & Zhu 2020; Zhang et al. 2024).

We note that the profiles shown in Fig. 3 are not the result of dynamical evolution over 100 planetary orbits, but rather extrapolations using the surface density evolution rates computed using ∂Fdep/∂R and the framework highlighted in Cordwell & Rafikov (2024) and Cordwell et al. (2025). Nevertheless, the aforementioned studies have shown that this approximation remains valid during the early phases of gap opening (or ΔΣ/Σ0 of a few percent). We refer the reader to these studies for further details.

thumbnail Fig. 1

Snapshot of the perturbed surface density after ten planetary orbits for our locally isothermal 3D model. Left: entire simulation domain. Spiral shocks are weaker near either radial end of the disk due to wave damping. Top right: a zoom-in into the orange box on the left, highlighting the region where gap opening is driven. Bottom right: angular momentum deposition function ∂Fdep/∂R and the surface density evolution rate ∂ς/∂t in arbitrary units computed according to Cordwell & Rafikov (2024). Saturated and faded colors indicate smoothed or raw data, respectively. Dashed and dotted vertical lines mark the isothermal shock distance xsh (Eq. (15) with γ = 1) and the planet’s corotating region (Eq. (22) with γ = 1), respectively.
thumbnail Fig. 2

Angular momentum deposition ∂Fdep/∂R for our fiducial 2D and 3D models (dashed and solid lines, respectively) and for four different equations of state. Models with cooling (β, radiative) have β2D ≈ 2. As in Fig. 1, faded and saturated colors indicate raw and smoothed data, respectively, showing that our smoothing procedure has negligible effects beyond the corotating region. The bottom panel compares the four equations of state for our 3D models, highlighting the role of cooling.
thumbnail Fig. 3

Estimates of the perturbed surface density expected after 100 planetary orbits for the models featured in Fig. 2 (using ∂Fdep/∂R from Fig. 2, Eq. (20), and the method in Cordwell & Rafikov 2024). The adiabatic and isothermal models show a double-trough gap structure, while radiative models show a single, deep gap. An apparent third dip in the adiabatic model at R = Rp is related to the perturbations driven by buoyancy modes near the edge of the planet’s horseshoe region.
thumbnail Fig. 4

Angular momentum deposition function ∂Fdep/∂R for our 3D (top) and 2D models (bottom) as a function of the cooling timescale defined through Eqs. (11) and (13) for 3D and 2D, respectively. The quantity approaches zero within xsh (dashed vertical lines) for very short or very long cooling timescales (blue and red lines, respectively), but peaks around β ~ 1 (gray dots).

3.2 Gap opening as a function of cooling timescale

We now carry out a suite of 2D and 3D models where we vary the cooling timescale β2D while keeping the optical depth constant at τ = 10, similar to Ziampras et al. (2023b). We then measure both the angular momentum deposition function ∂Fdep/∂R and the expected gap structure ΔΣ¯/Σ0$\Delta \bar \Sigma /{\Sigma _0}$ (via the surface density evolution rate ∂ς/∂t) as a function of β.

Our results for our fully radiative models in both 2D and 3D are shown in Figs. 4 and 5. Specifically, in Fig. 4 we see how very little angular momentum is deposited within |R - Rp| < xsh for very short or very long cooling timescales (blue or red limits, respectively), but the angular momentum deposition increases as we approach β ~ 1, highlighted with gray dots. Similar trends can be seen in Fig. 5, where the gap structure transitions from double- to single-trough as β increases from β ≪ 1 to β ~ 1, and then back to a double-trough gap for β ≫ 1. A secondary gap is also visible for R ~ 0.6 Rp for β ≪ 1 and β ≫ 1, induced by the planet’s secondary spiral arm in the inner disk. This spiral is strongly suppressed and cannot deposit a substantial amount of angular momentum for β ~ 1, resulting in a lack of additional gaps by a single planet, in line with the findings of Ziampras et al. (2020b) and Miranda & Rafikov (2020a).

Finally, we would like to quantify the contribution of radiative damping to gap opening as a function of β. To do this, we work under the assumption that the only mechanism of angular momentum deposition that contributes to gap opening is the one due to planetary spirals, and that it can be decomposed into two components: one due to shocks beyond xsh, and one due to radiative damping within xsh. Since ∂Fdep/∂R switches signs about Rp (and would therefore sum to zero over the gap-opening region), we choose ∂ς/∂t as a better indicator of gap opening near the planet. With these assumptions, we define the gap-opening efficiency factor G as

G=1ΩpRp∫Rp−xshRp+xsh1Σ∂Σ∂tdR.$G = {1 \over {{{\rm{\Omega }}_{\rm{p}}}{R_{\rm{p}}}}}\mathop \smallint \limits_{{R_{\rm{p}}} - {x_{sh}}}^{{R_{\rm{p}}} + {x_{sh}}} {1 \over {\rm{\Sigma }}}{{\partial {\rm{\Sigma }}} \over {\partial t}}{\rm{d}}R.$(23)

In other words, G reflects the mass flux exiting the gap region within the shock distance xsh from the planet. Since we are interested in relative changes in G due to cooling, we normalize it to the value for our adiabatic 2D model, Gadb2D$G_{{\rm{adb}}}^{2{\rm{D}}}$ when plotting.

We then present our results in Fig. 6, where we show that G changes substantially as a function of β, peaking for β ~ 1 and recovering the isothermal and adiabatic values for β ≪ 1 and β ≫ 1, respectively. These findings are consistent with those of Ziampras et al. (2023b), even though a different metric was used to quantify the gap-opening efficiency. Compared to that study, here we can obtain a more appropriate measurement of the gapopening efficiency in the vicinity of the planet rather than the entire disk.

From Fig. 6 we also observed that the 2D and 3D models agree very well in terms of the shape and amplitude of G across different β, albeit with a slight offset in the value of β at which the peak occurs. More importantly, however, we find that models with local, β cooling (purple lines) consistently overestimate the contribution of radiative damping to gap opening compared to fully radiative models (green lines). We interpret this as a result of the local cooling approach not properly capturing the effects of radiative diffusion, but rather assuming that the cooling is entirely localized around the spiral arms.

thumbnail Fig. 5

Estimated gap structures after 100 planetary orbits similar to Fig. 3 as a function of β. The gap structure transitions from double- to single-trough as β increases towards β ~ 1, and then back to a doubletrough gap for β ≫ 1. A secondary gap is also visible for R ~ 0.6 Rp for β ≪ 1 and β ≫ 1, induced by the planet’s secondary inner spiral arm.
thumbnail Fig. 6

Gap-opening efficiency factor G, normalized to the adiabatic 2D model value Gadb2D$G_{{\rm{adb}}}^{2{\rm{D}}}$, as a function of the cooling timescale β for our fully radiative models (green) and those with local, β cooling (purple), finding a good agreement between 2D and 3D models. Local β cooling overestimates G compared to fully radiative models.

3.3 Comparison between 2D and 3D models

When discussing our findings in previous sections, we noted a very good agreement between our 2D and 3D models in terms of the metrics we used to quantify gap opening. Both the angular momentum deposition function ∂Fdep/∂R and the expected surface density evolution rate ∂ς/∂t show both qualitatively and quantitatively similar radial behavior in Figs. 2 and 3, and the gap-opening efficiency factor G(β) in Fig. 6 also behaves very similarly between the two frameworks. This result paints an optimistic picture for the use of 2D radiative models in the context of planet- driven gap opening (see also Fung & Chiang 2016; Cordwell et al. 2025, for a complementary focus on isothermal models). They can be run with a much lower computational cost than their 3D counterparts while recovering the same gap structure, its dependence on the cooling timescale, and even the overestimation of G in models with β cooling compared to fully radiative ones.

Another noteworthy similarity can be found in the amplitude and shape of the surface density and midplane temperature deviations across the shock fronts themselves, as shown in Fig. 7 for our fiducial set with β2D ≈ 2 and at R = 1.15Rp, where ∂Fdep/∂R due to shocks is maximized (see Fig. 2). In this figure, we find that the 2D and 3D models agree remarkably well across all equations of state. This is particularly interesting since 3D shocks are expected to be weaker than their 2D counterparts (Lyra et al. 2016), but this is not straightforward to infer from this figure.

Nevertheless, some differences between the two frameworks are expected. The clearest deviation between 2D and 3D is found in the angular momentum deposition profiles for our adiabatic models in Fig. 2, where the 3D model shows a noticeably lower ∂Fdep/∂R beyond xsh but also a nonzero contribution within the planet’s corotating region. The latter is most likely due to the disk’s buoyancy response (Zhu et al. 2012), which can drive substantial radial and vertical motion in the region between the horseshoe edge and the shock distance (McNally et al. 2020; Ziampras et al. 2023c) and is therefore not filtered out by our smoothing procedure. In Fig. 8 we show heatmaps of the vertical velocity and the density-weighted correlations between the radial and azimuthal velocity components for our adiabatic 3D model at z = 2H, highlighting both the presence of buoyancy modes and a nonzero angular momentum flux in the vicinity of the planet. In contrast, no such features can be captured in a 2D model.

We note that this is different from the vortensity striping reported in Cordwell et al. (2025), which referred to a transient process during the horseshoe mixing phase in their isothermal model rather than the active vortensity evolution due to buoyancy modes in our adiabatic run. In fact, it is unclear whether the contribution to ∂Fdep/∂R due to buoyancy modes will persist as the gap-opening process develops. For this reason, the third, gaplike feature at R = Rp in Figs. 3 and 4 should be interpreted with caution. Clarifying this could be the subject of future work.

Finally, while the agreement between 2D and 3D models is excellent in several aspects in terms of the gap-opening efficiency factor G in Fig. 6, we note that the 3D models show a peak at β ≈ 2 rather than β ≈ 1 as in the 2D models. A difference between the two frameworks is expected (at least in the case of β cooling), as in 3D models the disk cools over the entire vertical column at each point along the disk plane, whereas 2D models attempt to capture this process in a vertically integrated fashion. Nevertheless, an agreement within a factor of two in terms of the value of β at which gap opening is most efficient is still very good, considering the significantly higher computational cost of 3D models and the excellent agreement between 2D and 3D regarding the behavior of G as a function of β in both amplitude and shape.

thumbnail Fig. 7

Azimuthal structure of the perturbed surface density (left) and midplane temperature (right) across the shock fronts at R = 1.15 Rp in our fiducial models for all cooling treatments. The 2D and 3D models (dashed and solid lines, respectively) agree very well with each other.
thumbnail Fig. 8

Heatmaps of the vertical velocity (top) and the density-weighted correlation between the radial and azimuthal velocity components (bottom) for our adiabatic 3D model at z = 2H. Buoyancy modes are clearly visible as diagonal stripes in the top panel, while the bottom panel shows a nonzero angular momentum flux around the planet (black dot).

4 Discussion

In this section, we discuss the implications of our findings in the context of planet-driven gap opening. We also address the limitations of 2D models in capturing the full complexity of the thermodynamic aspects of planet–disk interaction.

4.1 Cooling in parameter studies of planet–disk interaction

Modeling radiative processes in a satisfactory manner has been historically notoriously difficult in the context of protoplanetary disks. For one, the cooling timescale is highly sensitive to both gas and dust parameters such as the opacity, the dust-to-gas ratio, the gas surface density, and the temperature. Given the massive uncertainties on these parameters, and the sensitivity of the gap-opening process to the cooling timescale, it may even seem counterproductive to include radiative effects in models of planet– disk interaction. At the same time, the high computational cost of radiative simulations makes them impractical for parameter studies, especially in 3D.

However, several groundbreaking advancements have been made in recent years to both constrain disk parameters and develop more efficient numerical methods for radiative transfer. On the observational side, high-resolution ALMA observations have provided stringent constraints on the disk temperature structure in 3D (e.g., Law et al. 2021; Galloway-Sprietsma et al. 2025), as well as on the overall disk masses from which surface density profiles can be inferred (e.g., Zhang et al. 2025; Trapman et al. 2025; Longarini et al. 2025). Multi-wavelength observations of the same systems have also been utilized to infer dust properties such as the maximum grain size and dust size distribution (e.g., Doi & Kataoka 2023; Doi et al. 2024). These constraints, in combination with improved models of dust evolution (e.g., Birnstiel et al. 2012; Drążkowska et al. 2019; Stammler & Birnstiel 2022), dust opacities (e.g., Woitke et al. 2016; Birnstiel et al. 2018; Dominik et al. 2021), and radiative transfer (e.g., Ercolano et al. 2003; Dullemond et al. 2012), can be used to infer the properties of dust grains and constrain disk properties to the extent that a treatment of radiative cooling can be both motivated and yield realistic results. Finally, on the computational side, the rising usage of GPUs as tools for accelerating hydrodynamical simulations (e.g., Benítez-Llambay & Masset 2016; Lesur et al. 2023b; Stone et al. 2024) in combination with highly optimized linear algebra libraries that can be used for the express purpose of solving stiff operators such as the source terms in FLD (see Eq. (5); e.g., Balay et al. 2018; Trilinos Project Team 2020), have made it feasible to carry out large-scale radiative hydrodynamical simulations in 3D.

From the above it becomes clear that while modeling planet– disk interaction with a treatment for radiative transfer is quite involved and challenging, it is becoming increasingly more feasible and necessary to capture the full complexity of the problem. In fact, simple thermodynamic assumptions (e.g., locally isothermal) fail to capture the intricacies of planet–disk interaction in several ways including gap opening (e.g., Miranda & Rafikov 2020a; Ziampras et al. 2023b), planet-induced vortex dynamics (e.g., Fung & Ono 2021; Rometsch et al. 2021), circumplanetary disk dynamics (e.g., Krapp et al. 2024), gap edge dynamics (e.g., Muley et al. 2024), and planet migration (e.g., Lega et al. 2014; Ziampras et al. 2024b,a, 2025a).

Nevertheless, our results show that in the context of planet- induced gap opening it is both possible and accurate to use 2D models with a simple yet physically motivated treatment of radiative cooling via Eqs. (12) and (13) or (14). This opens up the possibility of carrying out large parameter studies of planet–disk interaction in 2D, while still capturing the key features of fully radiative models.

4.2 The role of viscous diffusion in gap opening

Viscosity acts as both a damping mechanism for spiral shocks and as a diffusive process across the gap region. As such, the angular momentum deposition function ∂Fdep/∂R is expected to be different in disks with significant levels of viscous diffusion. This could be important when comparing 2D to 3D results, as several sources of turbulence that require the inclusion of the vertical direction have been investigated in the context of protoplanetary disks but cannot be inherently captured in 2D models. Such mechanisms include the magnetorotational instability (MRI, Balbus & Hawley 1991; Hawley & Balbus 1991), the vertical shear instability (VSI, Nelson et al. 2013), the streaming instability (SI, Youdin & Johansen 2007; Johansen & Youdin 2007), the spiral wave instability by the planet’s spiral shocks (SWI, Bae et al. 2016b,a), the convective overstability (COS, Klahr & Hubbard 2014; Lyra 2014), and the zombie vortex instability (ZVI, Marcus et al. 2015, Marcus et al. 2016). Nevertheless, as Miranda & Rafikov (2020b) have shown, the effects of viscosity on the spiral angular momentum flux are negligible for α ≲ 10-2 (see Fig. 12 therein) which is generally the case for the aforementioned mechanisms and in realistic disk conditions (see Pfeil & Klahr 2019; Lyra & Umurhan 2019; Lesur et al. 2023a, for further details).

It should be noted, however, that while the general shape of the gap is expected to be similar between 2D and 3D models, the detailed structure of the gap edges may differ (see also Fung & Chiang 2016). This is particularly important in the context of vortices driven by the Rossby-wave instability (RWI, Lovelace et al. 1999; Chang & Youdin 2024), as their formation is sensitive to the shape of the gap edges and they can act as a source of turbulent diffusion during their lifetime (Lega et al. 2021). This effect is especially relevant for low-viscosity models such as the inviscid ones presented here, and could be the subject of future work.

5 Summary

We performed high-resolution two- and three-dimensional hydrodynamical simulations featuring embedded gap-opening planets. Our intent was to compare the 2D and 3D frameworks in the context of planet-driven gap opening across different thermodynamical treatments: locally isothermal, adiabatic, local β cooling, and full radiation transport including the effects of radiative diffusion. In doing so, we could both examine the accuracy of 2D models with respect to 3D, and identify simplified cooling prescriptions that could reproduce the key features of fully radiative models reasonably well.

We found that the overall angular momentum deposition profiles ∂Fdep (R)/∂R – which in turn determine the gap-opening efficiency and resulting gap shape – were very similar between 2D and 3D models for the same thermodynamical treatment. This is in good agreement with the findings of Cordwell et al. (2025) and Fung & Chiang (2016), who focused on locally isothermal models. Given the profiles of ∂Fdep (R)/∂R we could then synthesize the gap structure at its early phase following Cordwell & Rafikov (2024), and demonstrated that the expectation of a “single trough” gap shape for β ~ 1 and a “double trough” gap shape for locally isothermal and adiabatic models holds, in excellent agreement with previous studies (e.g., Miranda & Rafikov 2020a).

We then investigated the dependence of the gap-opening efficiency as a function of the cooling timescale by examining the total angular momentum deposition ∂Fdep (R)/∂R in the vicinity of the planet for different cooling prescriptions and disk parameters. Here, we found a continuous transition from double- to single- to double-trough gap shapes as β varied from ≪1, to ~1, and finally to ≫1, in line with the above expectations. We also found our metric for the gap-opening efficiency G(β) to peak for β ~ 1 as a direct consequence of that.

We found a very good match between 2D and 3D models in terms of the dependence of G(β) on the cooling timescale, albeit with a few quantitative differences. In particular, models with local β cooling slightly overestimated the overall gap-opening efficiency compared to their fully radiative counterparts, and the peak in G(β) for 2D models was found to be slightly offset towards shorter β by a factor of ≈2.

Nevertheless, the very good overall agreement between 2D and 3D models across all thermodynamical prescriptions is encouraging for two reasons. First, it suggests that 2D simulations can be used to accurately capture the gap structure around embedded planets, even when radiative effects are important. This is particularly relevant in the context of modeling the substructures found in observed systems with embedded planets, as 2D models are significantly less computationally expensive than their 3D counterparts. Second, it indicates that the use of a local, “β cooling” prescription can be a viable alternative to the more complex fully radiative approach when studying gap opening by planets. This makes it much easier to explore and/or constrain the parameter space before executing dedicated – possibly 3D – simulations with full radiation transport.

We believe that our results provide a solid foundation for the use of radiative processes in future studies of planet–disk interaction, and encourage the inclusion of radiation transport, even with simplified prescriptions, in such studies.

Data availability

Data from our numerical models are available upon reasonable request to the corresponding author.

Acknowledgements

This research utilized Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT (http://doi.org/10.5281/zenodo.438045). This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. A.Z. acknowledges funding from the European Union under the European Union’s Horizon Europe Research and Innovation Programme 101124282 (EARLYBIRD). A.J.C. is funded by the Royal Society of New Zealand Te Apārangi and the Cambridge Trust through the CambridgeRutherford Memorial Scholarship. R.R.R. acknowledges financial support through the STFC grant ST/T00049X/1 and the IAS. R.P.N. is supported by STFC grant ST/X000931/1. Views and opinions expressed are those of the authors only. All plots in this paper were made with the Python library matplotlib (Hunter 2007). Typesetting was expedited with the use of GitHub Copilot, but without the use of AI-generated text.

References

  1. ALMA Partnership (Fomalont, E. B., et al.) 2015, ApJ, 808, L1 [Google Scholar]
  2. Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bae, J., Nelson, R. P., & Hartmann, L. 2016a, ApJ, 833, 126 [Google Scholar]
  4. Bae, J., Nelson, R. P., Hartmann, L., & Richard, S. 2016b, ApJ, 829, 13 [Google Scholar]
  5. Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJ, 884, L41 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56 [NASA ADS] [CrossRef] [Google Scholar]
  7. Balay, S., Abhyankar, S., Adams, M. F., et al. 2018, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
  8. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
  9. Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11 [Google Scholar]
  10. Bi, J., Lin, M.-K., & Dong, R. 2021, ApJ, 912, 107 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bi, J., Lin, M.-K., & Dong, R. 2023, ApJ, 942, 80 [NASA ADS] [CrossRef] [Google Scholar]
  12. Binkert, F., Szulágyi, J., & Birnstiel, T. 2023, MNRAS, 523, 55 [NASA ADS] [CrossRef] [Google Scholar]
  13. Birnstiel, T. 2024, ARA&A, 62, 157 [NASA ADS] [CrossRef] [Google Scholar]
  14. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [CrossRef] [Google Scholar]
  16. Brown, J. J., & Ogilvie, G. I. 2024, MNRAS, 534, 39 [Google Scholar]
  17. Chang, E., & Youdin, A. N. 2024, ApJ, 976, 100 [Google Scholar]
  18. Cordwell, A. J., & Rafikov, R. R. 2024, MNRAS, 534, 1394 [Google Scholar]
  19. Cordwell, A. J., Ziampras, A., Brown, J. J., & Rafikov, R. R. 2025, MNRAS, 543, 4198 [Google Scholar]
  20. de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 Doi, K., & Kataoka, A. 2023, ApJ, 957, 11 [Google Scholar]
  21. Doi, K., Kataoka, A., Liu, H. B., et al. 2024, ApJ, 974, L25 [Google Scholar]
  22. Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library [record ascl:2104.010] [Google Scholar]
  23. Drążkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91 [Google Scholar]
  24. Duffell, P. C. 2020, ApJ, 889, 16 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
  26. Dullemond, C. P., Ziampras, A., Ostertag, D., & Dominik, C. 2022, A&A, 668, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X. W. 2003, MNRAS, 340, 1136 [Google Scholar]
  28. Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2017, ApJ, 835, 230 [CrossRef] [Google Scholar]
  29. Fung, J., & Chiang, E. 2016, ApJ, 832, 105 [NASA ADS] [CrossRef] [Google Scholar]
  30. Fung, J., & Ono, T. 2021, ApJ, 922, 13 [NASA ADS] [CrossRef] [Google Scholar]
  31. Galloway-Sprietsma, M., Bae, J., Izquierdo, A. F., et al. 2025, ApJ, 984, L10 [Google Scholar]
  32. Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [Google Scholar]
  33. Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [Google Scholar]
  34. Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793 [NASA ADS] [CrossRef] [Google Scholar]
  35. Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
  36. Hammer, M., & Lin, M.-K. 2023, MNRAS, 525, 123 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223 [NASA ADS] [CrossRef] [Google Scholar]
  38. Hubeny, I. 1990, ApJ, 351, 632 [NASA ADS] [CrossRef] [Google Scholar]
  39. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  40. Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [Google Scholar]
  41. Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJ, 806, L15 [NASA ADS] [CrossRef] [Google Scholar]
  42. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21 [Google Scholar]
  44. Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Kleine, T., Budde, G., Burkhardt, C., et al. 2020, Space Sci. Rev., 216, 55 [CrossRef] [Google Scholar]
  46. Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
  47. Krapp, L., Kratter, K. M., Youdin, A. N., et al. 2024, ApJ, 973, 153 [NASA ADS] [CrossRef] [Google Scholar]
  48. Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4 [NASA ADS] [CrossRef] [Google Scholar]
  49. Lega, E., Crida, A., Bitsch, B., & Morbidelli, A. 2014, MNRAS, 440, 683 [Google Scholar]
  50. Lega, E., Nelson, R. P., Morbidelli, A., et al. 2021, A&A, 646, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Lesur, G., Flock, M., Ercolano, B., et al. 2023a, ASP Conf. Ser., 534, 465 [Google Scholar]
  52. Lesur, G. R. J., Baghdadi, S., Wafflard-Fernandez, G., et al. 2023b, A&A, 677, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
  54. Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17 [NASA ADS] [CrossRef] [Google Scholar]
  55. Longarini, C., Lodato, G., Rosotti, G., et al. 2025, ApJ, 984, L17 [Google Scholar]
  56. Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
  57. Lyra, W. 2014, ApJ, 789, 77 [Google Scholar]
  58. Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001 [Google Scholar]
  59. Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102 [NASA ADS] [CrossRef] [Google Scholar]
  60. Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87 [NASA ADS] [CrossRef] [Google Scholar]
  61. Marcus, P. S., Pei, S., Jiang, C.-H., & Barranco, J. A. 2016, ApJ, 833, 148 [NASA ADS] [CrossRef] [Google Scholar]
  62. Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. McNally, C. P., Nelson, R. P., Paardekooper, S.-J., Benítez-Llambay, P., & Gressel, O. 2020, MNRAS, 493, 4382 [NASA ADS] [CrossRef] [Google Scholar]
  64. Menou, K., & Goodman, J. 2004, ApJ, 606, 520 [NASA ADS] [CrossRef] [Google Scholar]
  65. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
  66. Mignone, A., Flock, M., Stute, M., Kolb, S. M., & Muscianisi, G. 2012, A&A, 545, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Miranda, R., & Rafikov, R. R. 2020a, ApJ, 892, 65 [Google Scholar]
  68. Miranda, R., & Rafikov, R. R. 2020b, ApJ, 904, 121 [Google Scholar]
  69. Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [Google Scholar]
  70. Muley, D., Melon Fuksman, J. D., & Klahr, H. 2024, A&A, 690, A355 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
  73. Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950 [Google Scholar]
  74. Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. B. 2010, ApJ, 725, 146 [Google Scholar]
  75. Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150 [NASA ADS] [CrossRef] [Google Scholar]
  76. Pfeil, T., & Klahr, H. 2021, ApJ, 915, 130 [NASA ADS] [CrossRef] [Google Scholar]
  77. Rafikov, R. R. 2002a, ApJ, 569, 997 [Google Scholar]
  78. Rafikov, R. R. 2002b, ApJ, 572, 566 [Google Scholar]
  79. Rometsch, T., Ziampras, A., Kley, W., & Béthune, W. 2021, A&A, 656, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Stammler, S. M., & Birnstiel, T. 2022, ApJ, 935, 35 [NASA ADS] [CrossRef] [Google Scholar]
  82. Stoll, M. H. R., Picogna, G., & Kley, W. 2017, A&A, 604, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Stone, J. M., Mullen, P. D., Fielding, D., et al. 2024, arXiv e-prints [arXiv:2409.16053] [Google Scholar]
  84. Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton University Press) [Google Scholar]
  85. Teague, R., Benisty, M., Facchini, S., et al. 2025, ApJ, 984, L6 [Google Scholar]
  86. Toci, C., Lodato, G., Christiaens, V., et al. 2020, MNRAS, 499, 2015 [NASA ADS] [CrossRef] [Google Scholar]
  87. Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [Google Scholar]
  88. Trapman, L., Zhang, K., Rosotti, G. P., et al. 2025, ApJ, 989, 5 [Google Scholar]
  89. Trilinos Project Team, T. 2020, The Trilinos Project Website (accessed May 22, 2020) [Google Scholar]
  90. van Capelleveen, R. F., Ginski, C., Kenworthy, M. A., et al. 2025, ApJ, 990, L8 [Google Scholar]
  91. Van Leer, B. 1974, J. Comput. Phys., 14, 361 [NASA ADS] [CrossRef] [Google Scholar]
  92. Wafflard-Fernandez, G., & Lesur, G. 2023, A&A, 677, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178 [NASA ADS] [CrossRef] [Google Scholar]
  94. Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Yamaleev, N. K., & Carpenter, M. H. 2009, J. Comput. Phys., 228, 3025 [NASA ADS] [CrossRef] [Google Scholar]
  96. Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [Google Scholar]
  97. Yun, H.-G., Kim, W.-T., Bae, J., & Han, C. 2022, ApJ, 938, 102 [NASA ADS] [CrossRef] [Google Scholar]
  98. Zhang, S., & Zhu, Z. 2020, MNRAS, 493, 2287 [NASA ADS] [CrossRef] [Google Scholar]
  99. Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47 [Google Scholar]
  100. Zhang, M., Huang, P., & Dong, R. 2024, ApJ, 961, 86 [NASA ADS] [CrossRef] [Google Scholar]
  101. Zhang, K., Pérez, L. M., Pascucci, I., et al. 2025, ApJ, 989, 1 [Google Scholar]
  102. Zhu, Z., Stone, J. M., & Rafikov, R. R. 2012, ApJ, 758, L42 [NASA ADS] [CrossRef] [Google Scholar]
  103. Ziampras, A., Ataiee, S., Kley, W., Dullemond, C. P., & Baruteau, C. 2020a, A&A, 633, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Ziampras, A., Kley, W., & Dullemond, C. P. 2020b, A&A, 637, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Ziampras, A., Kley, W., & Nelson, R. P. 2023a, A&A, 670, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Ziampras, A., Nelson, R. P., & Rafikov, R. R. 2023b, MNRAS, 524, 3930 [Google Scholar]
  107. Ziampras, A., Paardekooper, S.-J., & Nelson, R. P. 2023c, MNRAS, 525, 5893 [Google Scholar]
  108. Ziampras, A., Nelson, R. P., & Paardekooper, S.-J. 2024a, MNRAS, 532, 351 [Google Scholar]
  109. Ziampras, A., Nelson, R. P., & Paardekooper, S.-J. 2024b, MNRAS, 528, 6130 [Google Scholar]
  110. Ziampras, A., Nelson, R. P., & Paardekooper, S.-J. 2025a, MNRAS, 542, 1685 [Google Scholar]
  111. Ziampras, A., Sudarshan, P., Dullemond, C. P., et al. 2025b, MNRAS, 536, 3322 [Google Scholar]
  112. Zier, O., & Springel, V. 2023, MNRAS, 520, 3097 [NASA ADS] [CrossRef] [Google Scholar]
1

For FLD, we use the PetSc library (Balay et al. 2018) in 3D, and a simple successive over-relaxation technique in 2D.

2

We note that this definition of the thermal mass is slightly different from that of Mth = h3 M* in, for example, Cordwell et al. (2025).

3

Tools for the calculation of ΔΣ/Σ0 from ∂Fdep/∂R are available at https://github.com/cordwella/disc_planet_analysis

4

Horseshoe dynamics leads to a nonzero ∂Fdep/∂R within xh (see Eq. (22)), but as mentioned above we set this quantity to zero in this region.

Appendix A Resolution study

To ensure that our results are appropriately converged with respect to numerical resolution, we carried out a series of runs for our fiducial 3D model with β2D ≈ 2 where we varied the resolution from 4 to 32 cells per scale height (cps) at the planet’s location. The same resolution check for 2D models was carried out in Ziampras et al. (2023b), where it was found that 32 cps were sufficient for convergence, so we did not repeat it here.

thumbnail Fig. A.1

Gap-opening efficiency metric G as a function of the number of grid cells per scale height at the planet’s location. Top: absolute value of G. Bottom: relative difference with respect to the highest-resolution run with 32 cps. We find that G is converged to within ≈ 5% for 16 cps.

Figure A.1 shows our gap-opening efficiency metric G (see Eq. (23)) as a function of the number of grid cells per scale height at the planet’s location. We find that G is converged to within ≈ 5% for 16 cps with respect to the 32 cps run. Given the high computational cost of our fully radiative 3D models, we chose to run our main suite of 3D models with 16 cps, which we consider to be sufficiently converged for our purposes.

Appendix B Derivation of the cooling timescale

Here we aim to provide a heuristic derivation of the dimensionless cooling timescale β used in Eq. (11). Our approach relies on the following key assumptions:

  • If cooling happens via dust thermal emission, the collisional coupling timescale between gas and dust is far shorter than the cooling timescale (see Dullemond et al. 2022, for details).

  • The radiation field Erad is isotropic and constant in time, and approximately satisfies Erad ≈ aRT4 (one-temperature approximation within the FLD framework).

  • Cooling due to advection or adiabatic decompression (terms u∇e and -γe∇ · u in Eq. (1c)) is negligible.

Equation (1c) can then be written as

∂e∂t≈−∇⋅F⇒ρcvdTdt≈∇⋅(λcκρ∇aRT4),${{\partial e} \over {\partial t}} \approx - \nabla \cdot F \Rightarrow \rho {c_{\rm{v}}}{{{\rm{d}}T} \over {{\rm{d}}t}} \approx \nabla \cdot \left( {{{\lambda c} \over {\kappa \rho }}\nabla {a_{\rm{R}}}{T^4}} \right),$(B.1)

We now consider the two limits where cooling is limited by radiative diffusion due to the medium being opaque to radiation (optically thick, λ → 1/3) or by the emissivity of the medium (optically thin, F → caRT 4).

For simplicity, we further assumed that quantities other than temperature are constant such that they can be pulled out of the gradient operators. In practice, our results can be corrected with a prefactor of order unity to account for a typical radial and/or temperature dependence of ρ or κ.

In the optically thick limit, Eq. (B.1) transforms into a diffusion equation with

ρcv∂T∂t≈∇⋅(caR4T33κρ∇T)⇒∂T∂t∼η∇2T,$\rho {c_{\rm{v}}}{{\partial T} \over {\partial t}} \approx \nabla \cdot \left( {{{c{a_{\rm{R}}}4{T^3}} \over {3\kappa \rho }}\nabla T} \right) \Rightarrow {{\partial T} \over {\partial t}} \sim \eta {\nabla ^2}T,$(B.2)

where we have used that aRc = 4σSB and the definition of the radiative diffusion coefficient η from Eq. (11). We can now estimate the diffusion timescale from the above equation as tthick =βthick ΩK−1=ldiff 2/η${t_{{\rm{thick }}}} = {\beta _{{\rm{thick }}}}\Omega _{\rm{K}}^{ - 1} = l_{{\rm{diff }}}^2/\eta $, where ldiff is a characteristic diffusion lengthscale. Assuming that ldiff ~ H (Flock et al. 2017; Ziampras et al. 2023b), we arrive at

βthick=ΩKηH2.${\beta _{{\rm{thick}}}} = {{{{\rm{\Omega }}_{\rm{K}}}} \over \eta }{H^2}.$(B.3)

In the optically thin limit, Eq. (B.1) instead becomes

ρcv∂T∂t≈∇⋅(caRT4)⇒∂T∂t≈V⋅∇T, V=16σSBT3ρcv,$\rho {c_{\rm{v}}}{{\partial T} \over {\partial t}} \approx \nabla \cdot \left( {c{a_{\rm{R}}}{T^4}} \right) \Rightarrow {{\partial T} \over {\partial t}} \approx V \cdot \nabla T,\quad V = {{16{\sigma _{SB}}{T^3}} \over {\rho {c_{\rm{v}}}}},$(B.4)

which corresponds to an advection of freely streaming radiation at a characteristic velocity V. Here, the characteristic lengthscale is the photon mean free path lrad (see Eq. (11)), and we can estimate the advection timescale as tthin =βthin ΩK−1=lrad/V${t_{{\rm{thin }}}} = {\beta _{{\rm{thin }}}}\Omega _{\rm{K}}^{ - 1} = {l_{{\rm{rad}}}}/V$. We then arrive at

βthin=ΩKρcvlrad16σSBT3=ΩKηκPκRlrad23.${\beta _{{\rm{thin}}}} = {{{{\rm{\Omega }}_{\rm{K}}}\rho {c_{\rm{v}}}{l_{{\rm{rad}}}}} \over {16{\sigma _{{\rm{SB}}}}{T^3}}} = {{{{\rm{\Omega }}_{\rm{K}}}} \over \eta }{{{\kappa _{\rm{P}}}} \over {{\kappa _{\rm{R}}}}}{{l_{rad}^2} \over 3}.$(B.5)

Realistically, cooling will be primarily limited by the slowest of the two channels. We can therefore estimate the total cooling timescale as

β3D=βthick+βthin=ΩKη(H2+κPκRlrad23),${\beta _{3{\rm{D}}}} = {\beta _{{\rm{thick}}}} + {\beta _{{\rm{thin}}}} = {{{{\rm{\Omega }}_{\rm{K}}}} \over \eta }\left( {{H^2} + {{{\kappa _{\rm{P}}}} \over {{\kappa _{\rm{R}}}}}{{l_{{\rm{rad}}}^2} \over 3}} \right),$(B.6)

which matches Eq. (11) in the main text (a similar, although slightly different in its final form, result was obtained in Miranda & Rafikov 2020b). Alternatively, one can write that β3D = max (βthick, βthin) (e.g., Pfeil & Klahr 2021) which behaves similarly but is not differentiable at the transition between the two regimes.

All Figures

thumbnail Fig. 1

Snapshot of the perturbed surface density after ten planetary orbits for our locally isothermal 3D model. Left: entire simulation domain. Spiral shocks are weaker near either radial end of the disk due to wave damping. Top right: a zoom-in into the orange box on the left, highlighting the region where gap opening is driven. Bottom right: angular momentum deposition function ∂Fdep/∂R and the surface density evolution rate ∂ς/∂t in arbitrary units computed according to Cordwell & Rafikov (2024). Saturated and faded colors indicate smoothed or raw data, respectively. Dashed and dotted vertical lines mark the isothermal shock distance xsh (Eq. (15) with γ = 1) and the planet’s corotating region (Eq. (22) with γ = 1), respectively.
In the text
thumbnail Fig. 2

Angular momentum deposition ∂Fdep/∂R for our fiducial 2D and 3D models (dashed and solid lines, respectively) and for four different equations of state. Models with cooling (β, radiative) have β2D ≈ 2. As in Fig. 1, faded and saturated colors indicate raw and smoothed data, respectively, showing that our smoothing procedure has negligible effects beyond the corotating region. The bottom panel compares the four equations of state for our 3D models, highlighting the role of cooling.
In the text
thumbnail Fig. 3

Estimates of the perturbed surface density expected after 100 planetary orbits for the models featured in Fig. 2 (using ∂Fdep/∂R from Fig. 2, Eq. (20), and the method in Cordwell & Rafikov 2024). The adiabatic and isothermal models show a double-trough gap structure, while radiative models show a single, deep gap. An apparent third dip in the adiabatic model at R = Rp is related to the perturbations driven by buoyancy modes near the edge of the planet’s horseshoe region.
In the text
thumbnail Fig. 4

Angular momentum deposition function ∂Fdep/∂R for our 3D (top) and 2D models (bottom) as a function of the cooling timescale defined through Eqs. (11) and (13) for 3D and 2D, respectively. The quantity approaches zero within xsh (dashed vertical lines) for very short or very long cooling timescales (blue and red lines, respectively), but peaks around β ~ 1 (gray dots).
In the text
thumbnail Fig. 5

Estimated gap structures after 100 planetary orbits similar to Fig. 3 as a function of β. The gap structure transitions from double- to single-trough as β increases towards β ~ 1, and then back to a doubletrough gap for β ≫ 1. A secondary gap is also visible for R ~ 0.6 Rp for β ≪ 1 and β ≫ 1, induced by the planet’s secondary inner spiral arm.
In the text
thumbnail Fig. 6

Gap-opening efficiency factor G, normalized to the adiabatic 2D model value Gadb2D$G_{{\rm{adb}}}^{2{\rm{D}}}$, as a function of the cooling timescale β for our fully radiative models (green) and those with local, β cooling (purple), finding a good agreement between 2D and 3D models. Local β cooling overestimates G compared to fully radiative models.
In the text
thumbnail Fig. 7

Azimuthal structure of the perturbed surface density (left) and midplane temperature (right) across the shock fronts at R = 1.15 Rp in our fiducial models for all cooling treatments. The 2D and 3D models (dashed and solid lines, respectively) agree very well with each other.
In the text
thumbnail Fig. 8

Heatmaps of the vertical velocity (top) and the density-weighted correlation between the radial and azimuthal velocity components (bottom) for our adiabatic 3D model at z = 2H. Buoyancy modes are clearly visible as diagonal stripes in the top panel, while the bottom panel shows a nonzero angular momentum flux around the planet (black dot).
In the text
thumbnail Fig. A.1

Gap-opening efficiency metric G as a function of the number of grid cells per scale height at the planet’s location. Top: absolute value of G. Bottom: relative difference with respect to the highest-resolution run with 32 cps. We find that G is converged to within ≈ 5% for 16 cps.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.