Buoyancy torques prevent low-mass planets from stalling in low-turbulence radiative disks

Low-mass planets migrating inwards in laminar protoplanetary disks (PPDs) experience a dynamical corotation torque, which is expected to slow down migration to a stall. However, baroclinic effects can reduce or even reverse this effect, leading to rapid inward migration. In the radiatively inefficient inner disk, one such mechanism is the buoyancy response of the disk to an embedded planet. Recent work has suggested that radiative cooling can quench this response, but for parameters that are not necessarily representative of the inner regions of PPDs. We perform global three dimensional inviscid radiation hydrodynamics simulations of planet-disk interaction to investigate the effect of radiative cooling on the buoyancy-driven torque in a more realistic disk model. We find that the buoyancy response exerts a negative dynamical corotation torque -- albeit partially damped due to radiative cooling -- resulting in sustained, rapid inward migration. Models that adopt a local cooling prescription significantly overestimate the impact of the buoyancy response, highlighting the importance of a realistic treatment of radiation transport that includes radiative diffusion. Our results suggest that low-mass planets should migrate inwards faster than has been previously expected in radiative disks, with implications for the formation and orbital distribution of super-Earths and sub-Neptunes at intermediate distances from their host stars, unless additional physical processes that can slow down migration are considered.


INTRODUCTION
Planets are born in protoplanetary disks (PPDs), and their observed location after the disk has dispersed is highly sensitive to their interaction with it.In particular, for low-mass planets, addressing their origin is a long-standing problem given how rapidly they migrate through the disk (Ward 1997).A solution to this problem is crucial in order to explain the large population of super-Earths and sub-Neptunes in the observed exoplanet demographics.For a recent review, see Paardekooper et al. (2022).
Regardless of planet mass, different mechanisms have been proposed to halt or slow the inward migration of young planets or even reverse it.Massive enough planets that can open a deep gap around their orbit (type-II regime, Lin & Papaloizou 1986;Rafikov 2002) can migrate very slowly in disks that sustain very low levels of turbulent viscosity (Ward 1997;Lega et al. 2022) although interaction with the Rossby-wave unstable gap edge can lead to vortexassisted inward or outward migration (McNally et al. 2019a;Lega et al. 2022).For intermediate-mass planets, migration has previously been expected to eventually stall as the planet acts akin to a snowplough and shovels a wall of gas ahead of its corotating region, in a phenomenon dubbed the "inertial limit" (Hourigan & Ward 1984;Ward & Hourigan 1989;Ward 1997).However, more recent work ⋆ E-mail: a.ziampras@qmul.ac.uk has shown that vortex activity can help sustain inward migration (McNally et al. 2019a).
For low-mass planets (type-I regime), the disk remains relatively unperturbed by the planet's presence beyond the formation of spiral arms (Ogilvie & Lubow 2002) and the planet's corotating region (Goldreich & Tremaine 1980).Here, corotation torques have been shown to efficiently stop or even reverse migration when the turbulent viscosity is sufficiently high (Paardekooper et al. 2011), and thermal torques due to the planet's accretion luminosity can do the same for roughly Earth-mass planets (Benítez-Llambay et al. 2015;Masset 2017).For super-Earths and sub-Neptunes in laminar disks however, none of these effects operate efficiently enough to significantly slow down type-I migration.
In laminar disks, where accretion is driven through the disk surface in the form of a magnetothermal wind (Bai & Stone 2013;Gressel et al. 2015;Wang et al. 2019), the corotation torque described above vanishes and the planet is forced to drag its corotating region along as it migrates inwards (Paardekooper 2014).As long as the vortensity of the material enclosed in that region is conserved (i.e., in the absence of baroclinic effects), the associated drag takes the form of a dynamical corotation torque (DCT) and can slow the planet down to an effective halt (Paardekooper 2014) even though migration doesn't formally stop.
This last scenario could help reconcile theoretical modeling with exoplanet demographics by dramatically slowing down type-I migration, thereby preventing low-mass planets from migrating to the inner rim of the disk.However, as stated above, it relies on the conservation of vortensity in the corotating region.While this condition is satisfied for barotropic flows (e.g., isothermal disks), radiative effects have been shown to generate vortensity near the planet (Pierens 2015;Ziampras et al. 2024), resulting in a negative torque that can not only reduce the drag due to the DCT but even accelerate inward migration beyond the expected type-I speed.At the same time, the buoyancy response of the disk to a low-mass embedded planet has been shown to produce a similar vorticity-generating effect in radiatively inefficient, adiabatic disks (McNally et al. 2020), in addition to providing an apparently smaller and separate buoyancy torque that arises through direct interaction between the planet and the buoyant modes in the disk (Zhu et al. 2012;Lubow & Zhu 2014).Recent work has shown that, in the context of a specific disk model, radiative cooling can act to quench this buoyancy response, preventing vortensity generation and erasing its contribution to the total torque (Yun et al. 2022).
While hydrodynamical modeling with a treatment of radiation transport is far more realistic than a simplified isothermal or adiabatic equation of state, it is significantly more expensive and the results are highly sensitive to the underlying disk model, which determines the cooling timescale and therefore the efficiency of radiative processes.With that in mind, while Yun et al. (2022) have demonstrated that radiative cooling can damp the disk buoyancy response, and given that its operation can decide the fate of embedded, low-mass planets (Paardekooper 2014), it is crucial to investigate in which disk conditions their result holds true.As shown by Yun et al. (2022), this requires simulations that treat radiative cooling in a realistic way.
Therefore, in this study, we perform hydrodynamical simulations that feature the disk buoyancy response to an embedded planet, designed to represent the inner few au of the disk and with a realistic treatment of radiation transport.Our goal is to examine whether the buoyancy-associated torque is ever relevant in real disks, and, if so, to what extent it is capable of modifying the planet's migration rate.
We describe our physical framework in Sect. 2. We motivate our numerical parameters in Sect.3, present our results from numerical simulations in Sect.4, and discuss their implications in Sect. 5. We finally summarize our findings in Sect.6.

PHYSICAL FRAMEWORK
In this section we lay out the physical framework of our models.We describe in detail our approach to radiation transport, and motivate our numerical experiments by investigating the conditions under which buoyancy torques could be relevant.

Equations of hydrodynamics
We consider a disk of ideal gas with adiabatic index γ = 7 /5 and mean molecular weight µ = 2.353, orbiting a star with mass M ⋆ = M ⊙ and luminosity L ⋆ = L ⊙ .The volume density ρ, velocity field u, and pressure P of the gas evolve according to the (inviscid) Euler equations where e = P/(γ − 1) is the internal energy density given by the ideal gas law, Φ = Φ ⋆ = GM ⋆ /r is the gravitational potential of the star at distance r, G is the gravitational constant, and Q encapsulates any additional sources of heating or cooling.The isothermal sound speed is then c s,iso = P/ρ, and the temperature is T = µc 2 s,iso /R, with R denoting the gas constant.
In a cylindrical coordinate system (R, φ, z) with r = √ R 2 + z 2 , by assuming that the disk is non-accreting and axisymmetric and further requiring that the midplane density ρ mid and vertically constant temperature follow radial power-law profiles such that the disk structure in equilibrium can be described by (Nelson et al. 2013) Here, Ω K = GM ⋆ /R 3 is the Keplerian angular velocity, h = H/R is the aspect ratio of the disk, and H = c s,iso /Ω K is the pressure scale height.Finally, we can define the surface density through Σ = ∞ −∞ ρ dz.

Radiative cooling
We model the radiative cooling of the disk using the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981).In this approach, the radiation energy density E rad evolves as and is coupled to the gas energy density with a source term in Eq. (1c).Here, κ P is the Planck mean opacity, c is the speed of light, and a R is the radiation constant.The radiation flux F is given by where κ R is the Rosseland mean opacity and λ is the flux limiter following Kley (1989): We do not explicitly include the heating due to stellar irradiation in our models.Although this implies that our disk model is only accurate up to the disk optical surface (z ∼ 4H, see Chiang & Goldreich 1997), this is not a problem as we are mainly interested in the region z ≲ 2H.We can sidestep this issue by prescribing E rad = a R T 4 0 at the upper boundary of our disk which prevents the gas from rapidly cooling off (see also Sect. 4.1).
In the FLD approach, the cooling timescale β cool = t cool Ω K can be approximated following Flock et al. (2017b) Here, σ SB is the Stefan-Boltzmann constant, and is the specific heat at constant volume.
The above expression for β is valid under the assumption that the solid particles ("dust") in the disk, which facilitate cooling, are wellcoupled to the gas.In particular, we focus on the small grains with size a d = 0.1 µm, as they provide the most efficient cooling channel.We can then express the collision timescale β coll between the gas and small grains following Dullemond et al. (2022) Here, s d = 4πa 2 d and m d = 4π 3 ρ d a 3 d are the surface area and mass of a single dust particle with bulk density ρ d = 2.08 g/cm 3 , and ϵ is the dust-to-gas mass ratio.The combined cooling timescale β tot is then given by We note that the cooling timescale is a sensitive function of the disk's temperature and density structure, and depends strongly on the choice of opacity model.In our models, we use the density-and temperature-dependent opacity model of Lin & Papaloizou (1985) with κ R = κ P .

Buoyancy-driven torque
The gravitational interaction between disk and planet excites a buoyancy response in the disk, which can generate a torque on the planet.Zhu et al. (2012) showed with local shearing box simulations that the instantaneous torque due to the direct interaction between the planet and the buoyantly excited gas can become comparable to the Lindblad torque (Ward 1997) in the adiabatic limit (β tot → ∞).McNally et al. (2020) then carried out global adiabatic models with a migrating planet and showed that, in addition to the effect shown by Zhu et al. (2012), buoyancy oscillations result in a long-term build-up of vortensity (ϖ) in the planet's corotating region, which can erase or even reverse the stalling effect of the dynamical corotation torque (Paardekooper 2014;McNally et al. 2017) Here, a subscript "p" denotes quantities at the planet's location, u R is the radial velocity of the background disk, and x h is the horseshoe half-width (Paardekooper et al. 2010b) defined as The smoothing length ϵ is chosen to match 3D models (Müller & Kley 2012).In Eq. ( 12), the vortensity is given in the vertically integrated approximation by ϖ = (∇ × u)/Σ • ẑ, while Masset & Benítez-Llambay (2016) derived the equivalent quantity in 3D to be Then, ϖ(R p ) is the vortensity at R p for an unperturbed disk, and ϖ h is the characteristic vortensity enclosed in the planet's corotating region.As the planet migrates inwards it carries with it the corotating material that was present at its initial location, while preserving the value of ϖ h in the absence of turbulent viscosity or vortensitygenerating baroclinic effects.This in turns leads to evolution of the quantity ϖ(R p )/ϖ h , and a change in the corotation torque.
In the absence of vortensity-generating (baroclinic) effects the unperturbed vortensity in the background disk is given by ϖ 0 (R) = 1 2 Ω K /Σ eq , which would result in a continuously increasing positive torque as the planet migrates inwards for Σ(R) shallower than R −3/2 .However, the excitation of buoyancy oscillations in the planet's horseshoe region induces a baroclinic forcing (Ziampras et al. 2023), generates vortensity, and results eventually in a negative dynamical corotation torque (McNally et al. 2020).
Yun et al. ( 2022) have called into question the relevance of the buoyancy-driven torque in radiative disks, showing that radiative diffusion can quench the buoyancy response and erase the associated torque when the cooling timescale due to radiative diffusion (β diff = Ω K H 2 /η, see Eq. ( 9)) is comparable to or shorter than the buoyancy timescale β buoy .The latter is related to the Brunt-Väisälä frequency through which we can write In their simulations and for their disk model, Yun et al. (2022) found that β diff ≲ β buoy for z ≳ 2H, and concluded that the buoyancy-driven torque is quenched in radiative disks.However, as stated above, β diff is a sensitive function of the underlying disk model.It is therefore necessary to investigate the conditions under which the buoyancydriven torque is relevant in a more slowly cooling, passively irradiated disk model (e.g., Chiang & Goldreich 1997).

COOLING TIMESCALE MAPS
In this section we present maps of the cooling timescale β tot and the ratio β buoy /β tot for several different disk models.We use these maps to motivate our choice of numerical parameters for our hydrodynamical simulations.
In agreement with the findings of Yun et al. ( 2022), we find that β buoy /β tot ∼ 0.1 at z ∼ 2H in the inner few au, which is why buoyancy torques were quenched in their disk model.In general, with the exception of a small region around the water iceline at 5 au, we find that β buoy /β tot ≳ 0.01 between z ∼ 1-2H throughout the whole disk, and therefore expect that buoyancy torques are indeed not relevant for this disk model.

A more realistic disk model
The disk model used by McNally et al. (2020), while both simple and effective in demonstrating the effect of buoyancy torques in a scalefree environment, is not quite representative of a real PPD.Passively heated disks are flared (q ≈ −1/2, Chiang & Goldreich 1997) and significantly thinner (h ≈ 0.02-0.025at 1 au) than the constant aspect ratio disk model used above.Achieving such a comparatively high value for h would require considerable internal heating due to turbulence, with a corresponding α ∼ 10 −3 -10 −2 (Shakura & Sunyaev 1973).This is both difficult to reconcile with the low turbulence expected in the dead zone (Bai & Stone 2013)  With that in mind, we repeat the above exercise assuming a passively heated disk model.We assume that the aspect ratio h follows a power-law profile h(R) = 0.02 (R/au) 2/7 , or equivalently that q = −3/7 (Chiang & Goldreich 1997;Menou & Goodman 2004).We also assume that the surface density profile is given by Σ(R) = 1700 (R/au) −1 g/cm 2 , matching the Minimum Mass Solar Nebula model (Hayashi 1981) at R = 1 au.We then recompute the cooling timescale β tot and the ratio β buoy /β tot , and present the results in the middle panels of Fig. 1.
Given that β diff ∝ T −3 ∝ h −6 , such a significant drop in the aspect ratio results in dramatically longer cooling timescales in the inner disk, even though the surface density is lower.In fact, we find that β buoy /β tot ≲ 10 −3 up to z ∼ 2-3H for R ≲ 2 au, which could be sufficient for buoyancy torques to operate in the inner disk.

A computationally friendly hybrid model
The model computed in the previous section, while more realistic and certainly more optimistic for the operation of buoyancy torques, is also significantly more expensive to use in the global models we would like to execute.It would also result in much more rapid migra-tion as the Lindblad torque scales with h −2 (Goldreich & Tremaine 1979).Given that the buoyancy response of the disk is a subtle effect highly sensitive to the numerical resolution (Ziampras et al. 2023), achieving a high enough resolution (∼16 cells per scale height) to capture the effect of buoyancy torques in the inner disk (h ∼ 0.02) while also integrating the radiation hydrodynamics equations in Eqs. ( 1) & (5) for several hundred planetary orbits is computationally prohibitive.
We note here that our goal is not to measure the precise migration rate of a planet in a realistic disk model, but rather to investigate whether the buoyancy-driven torque is even active in the inner few au of a disk, where there still exists the possibility that radiative diffusion can damp the buoyancy response of the disk.With the optimistic model in Sect.3.2 in mind, we therefore choose to use a hybrid disk model that is both computationally friendly and still representative of the radiative environment in a real protoplanetary disk.
To that end, we assume a flared disk with q = −3/7 and Σ ∝ 1/R as before, but adjust Σ 0 and h 0 such that the values of Σ and h at the planet's initial radial location of 2 au are identical to those in the model of McNally et al. (2020) and Yun et al. (2022).In this way, our model represents similar conditions to the realistic model in Sect.3.2, is computationally feasible, and also maintains compatibility with previous work to an extent.We therefore choose Σ 0 = 4800 g/cm 2 (R/au) −1 and h = 0.041 (R/au) 2/7 .Of course, adjusting both Σ and h in this way will inevitably also change the cooling timescale compared to our previous model.To account for this difference and ensure that the disk buoyancy response is similar between the model in Sect.3.2 and this hybrid model, we artificially increase the dust opacity by a factor of ≈ 7. The result is a near-perfect match between the two models in terms of the ratio β buoy /β tot both radially and vertically, as shown in Fig. 2. We also show the cooling timescale maps for this model in the right panels of Fig. 1.We will be using this hybrid disk model in our hydrodynamical simulations in Sect. 4.

HYDRODYNAMICAL SIMULATIONS
In this section we present the results from our hydrodynamical simulations.We first describe our numerical setup, and then present the results from our models with a migrating planet.

Numerical setup
We use the numerical hydrodynamics Godunov code PLUTO v.4.4 (Mignone et al. 2007) with a radiation transport module following the implementation of the FLD closure by Kolb et al. (2013) and described in Appendix A. We also utilize the FARGO algorithm (Masset 2000), implemented in PLUTO by Mignone et al. (2012).This method alleviates the strict timestep limitation imposed by the rapidly rotating inner radial boundary by subtracting the Keplerian rotation profile of the background disk before advecting the gas.In doing so, it provides a substantial speedup while also reducing numerical diffusion.
We use the ENTROPY_SWITCH option to offset the susceptibility of PLUTO to the high-Mach problem (Trac & Pen 2004).We also use the HLLC Riemann solver (Harten 1983), a second-order (RK2) time marching scheme, a third-order weighted essentially non-oscillatory reconstruction (WENO3, Yamaleev & Carpenter 2009) and the flux limiter by Van Leer (1974).This combination of numerical options has been shown to be both accurate and robust in the context of resolving the disk buoyancy response to a low-mass planet (Ziampras et al. 2023).Ziampras et al. (2023) showed that, even at high resolution and with a sophisticated numerical setup, PLUTO struggles to capture the disk buoyancy response and the associated torque compared to finite-difference, upwind codes such as FARGO3D (Benítez-Llambay & Masset 2016).For this purpose, while we carry out our radiative simulations with PLUTO, we also perform a set of simulations in the adiabatic limit and with a simplified, local cooling prescription with both PLUTO and FARGO3D v.1.2.In this way, we can gauge the extent to which PLUTO might underestimate the disk's buoyancy response and correct for it, allowing us to estimate the planet's migration track in a scenario where the buoyancy torque is captured more appropriately.
Our FARGO3D setup is based on the hybrid MPI/OpenMP parallelism implemented by McNally et al. (2019b) and solves the specific entropy s instead of the energy equation in Eq. (1c): where T ref and ρ ref are arbitrary constants.The thermal relaxation source term Q ′ , when applicable, essentially represents the same physics as in Eq. ( 1c) but is rewritten in the specific entropy formulation.
In both codes our grid extends radially between 0.8-4 au with a logarithmic spacing, vertically from the midplane up to z = 5H at R 0 = 2 au, and spans the full azimuthal range 0-2π.We use a grid of  2006) and smooth the planet's potential with the cubic spline curve described in Klahr & Kley (2006).The radiation energy is initialized as E rad,0 = a R T 4 0 .At the radial and vertical edges of the domain, all parameters are set to their initial values.For E rad , this implies that the disk is shielded from direct stellar irradiation, which is satisfied in our models (see green lines in Fig. 1).
The planet is held fixed at R 0 for 200 orbits and is then allowed to migrate for an additional 250, for a total of 450 orbits at R 0 by the end of the simulations.The planet grows over 200 orbits to its final mass of M p = 2 × 10 −5 M ⋆ ≈ 6.7 M ⊕ using the formula by de Val-Borro et al. (2006).Since we do not consider the disk selfgravity, we subtract the azimuthal average of the gas surface density before computing the acceleration on the planet (Baruteau & Masset 2008).Finally, we include the indirect term due to the star-planet system orbiting about their center of mass, and further include the indirect term due to star-disk interaction during the migration phase (Crida et al. 2022).
In our adiabatic runs, the source terms Q and Q ′ in Eqs. ( 1c) & ( 17) are set to zero.In our runs with a simplified local β cooling prescription we implement a thermal relaxation term similar to Gammie (2001) and McNally et al. (2020) for PLUTO and FARGO3D, respectively: where β tot is given by Eqs. ( 9)-(11).Our fully radiative run in PLUTO solves Eqs. ( 5) & (6) instead.By comparing a fully radiative to a locally cooled model, we can assess the effect of radiation transport on the buoyancy torque as not only a cooling effect, but also a diffusion mechanism.
Finally, in order to isolate the effect of the buoyancy torque, we supplement our set of simulations with three vertically integrated (R, ϕ) models using PLUTO.Similar to the 3D runs, we employ the following three cooling prescriptions: • local cooling: β relaxation term following the prescription of surface and in-plane cooling in Ziampras et al. (2024), • radiative: stellar irradiation, surface cooling, and in-plane radiative diffusion following the implementation of FLD in Ziampras et al. (2020).
In these 2D models we use a Plummer potential for the planet's gravity, with a smoothing length ϵ = 0.6H p (Müller & Kley 2012).

Static planet phase
During the first 200 orbits, the planet amasses a vortensity excess around the edge of its horseshoe region in all 3D models, with an additional excess in the center of the horseshoe region for the adiabatic and β-cooled models.We show the azimuthally averaged perturbed vortensity in Fig. 3.This already suggests that models with local cooling (i.e., that do not account for the effect of radiative diffusion) do not correctly capture the dissipation of buoyancy oscillations in the planet's horseshoe region, and therefore overestimate the associated vortensity generation.
We then plot perturbed vortensity heatmaps for the adiabatic and radiative 3D models in Fig. 4. Models with β cooling are not shown as they are identical to their adiabatic counterparts.Similar to Fig. 3, we find that the vortensity excess is more centered to R = R p in the adiabatic models, suggesting that higher-order buoyancy modes, which are excited radially closer to the planet, are damped due to radiative diffusion-fully consistent with the findings of Yun et al. (2022).The steeper vortensity gradients induced by the stronger vortensity growth in the adiabatic models also result in the excitation of vortices, which orbit around the horseshoe region and are not present in the radiative models.
Overall, until the planet is allowed to migrate, our results are consistent with the findings of McNally et al. (2020) and Yun et al. (2022), at least to an extent: the buoyancy response is partially damped in radiative models, although a vortensity excess in the horseshoe region is still visible.This indicates that the longer cooling timescales in the inner disk are indeed sufficient for buoyancy torques to operate, albeit not at the same level as in the adiabatic models.

Migrating planet phase
We now allow the planet to migrate through the disk for an additional 250 orbits.The migration tracks for all models and the associated migration timescales τ mig = R p / Ṙp are shown in Fig. 5. Several takeaways can be drawn from this figure: • In 2D, the planet migrates inward slower than the type-I rate, indicating a positive dynamical corotation torque (DCT) for all radiative treatments.This is expected given the very long cooling timescale and the lack of a buoyancy response in these 2D models, and leads to nearly perfectly overlapping tracks in all 2D models.
• The buoyancy response in our adiabatic 3D models induces a negative DCT, resulting in rapid inward migration.This is consistent with the findings of McNally et al. (2020).
• The planet initially migrates at the type-I rate in our radiative 3D models, but accelerates after ≈ 80 orbits.This suggests that the buoyancy-driven torque is operating, and reaches full strength as the planet continues to migrate.We investigate this behavior further below, but it is clear that the buoyancy torque is indeed active in our radiative models.
• Models with local cooling show a migration rate practically identical to the adiabatic models for both codes, indicating that the buoyancy torque is higher than in the radiative case and highlighting the need for realistic radiative modeling in order to correctly capture the buoyancy response of the disk.
To address the "knee" in the migration track of the 3D radiative model, we plot vortensity heatmaps near the planet's corotating region for this model and for several snapshots in Fig. 6.While  shows that this section shrinks abruptly to approximately a third of its initial azimuthal extent after ≈ 80 orbits, at which point the planet accelerates inwards.
Observing the migration timescale in the bottom panel of Fig. 5, we find that this delayed transition to faster migration is present to an extent in all 3D models and is consistent with the libration timescale which is the typical time it takes a gas streamline at the edge of the corotating region to complete a closed orbit and therefore allow the torque on the migrating planet to stabilize.We note that the planet's migration rate is roughly constant for all 3D models after ≈ 80 orbits, in line with this argument.Consistent with the findings of Ziampras et al. (2023), we find that PLUTO underestimates the buoyancy torque compared to FARGO3D.In our models this translates to a factor of τ Pluto mig /τ Fargo mig ≈ 3.5 between the two codes for both adiabatic and β-cooled models.It is tempting to also scale the migration timescales of our radiative models by this factor, but we caution that the thermal diffusion associated with FLD results in some physical damping of the buoyancy response (see Ziampras et al. 2023).As a result, the final migration timescale of the planet in a radiative disk model could in principle be up to 3.5× shorter that what we find in our PLUTO models, but the difference between the two codes is likely much smaller in the radiative case.Regardless, we stress that the buoyancy response exerts a substantial negative DCT on the planet in our radiative model, in contrast to the findings of Yun et al. (2022) who concluded that the buoyancy torque is quenched in radiative disks.This does not invalidate their results for their disk model, but rather shows that whether the disk buoyancy response can exert a significant torque on a migrating planet is highly sensitive to the underlying disk conditions.In the next section we compare the excited buoyancy modes between the two models in an attempt to identify key differences in the buoyancy response of the disk.

Buoyancy mode excitation
In the previous section we showed that the disk buoyancy response can exert a significant torque on a migrating planet in favorable conditions, even when radiative cooling is considered.At the same time, Yun et al. (2022) showed that the buoyancy torque is already quenched when β buoy /β diff ≲ 0.1 at z ∼ 2H.Here, we attempt to quantify the difference in the buoyancy response between the two disk models.
To that end, we carry out an additional set of short-term, global simulations with a fixed planet using PLUTO, and compare the am-  plitude and phase of the buoyancy modes generated by the planet in the disk model presented in Yun et al. ( 2022) and our hybrid model (see Sect. 3).For these models the planet grows over one orbit and the simulation runtime is 10 orbits, similar to Ziampras et al. (2023).
Figure 7 shows an azimuthal slice of the gas vertical velocity u z normalized to the local sound speed at z = {H, 2H} and R = R p − x h , tracing the excited buoyancy oscillations at the inner edge of the horseshoe region.We find that: • Both models agree very well in the adiabatic case (blue curves), with our model showing very slightly stronger buoyancy modes on average.This difference is due to the disk flaring in our model, which results in a slightly smaller aspect ratio at that location.We find the opposite behavior at R = R p + x h .
• Cooling damps buoyancy modes, more so for fully radiative (green) than for β models (orange).This further highlights the need for realistic radiative modeling, and is consistent with the faster inward migration found for β models in Fig. 5.
• Buoyancy modes in our radiative models (green) are more spaced out in azimuth compared to the adiabatic and β-cooled cases.This is consistent with the findings of Yun et al. ( 2022), and is most likely due to the diffusive nature of FLD.
The latter in particular might be a secondary reason for the weaker buoyancy torque, as a combination of weaker and fewer buoyancy modes will quickly diminish the rate at which vortensity can be generated in the corotating region.We expect that a "cutoff" in the buoyancy response of the disk will occur for β buoy /β tot ≳ 0.01 at z ∼ 2H, but finding the exact value is beyond the scope of this work.

DISCUSSION
In this section we discuss the implications and limitations of our findings.We outline the current challenges in planet migration modeling, and possible pathways to reconciling with the observed exoplanet population.

Vortensity growth due to a radial entropy gradient
In this work we focused on the vortensity generation due to the disk buoyancy response.We note that, in principle, a radial entropy gradient can also lead to vortensity spikes around the edges of the horseshoe region (Masset & Casoli 2009;Paardekooper et al. 2010a), although this has only been investigated in 2D.In our models, the entropy function is P/ρ γ ∝ R 0.48 and P 2D /Σ γ ∝ R −0.03 in 3D and 2D, respectively.Thus, we can expect that a fraction of the vortensity growth shown for 3D models in Figs. 3 & 4 could be due to the contribution of this mechanism, rather than the disk buoyancy response (see also Fig. 14 in McNally et al. 2020, whose model satisfied P/ρ γ ∝ R −0.4 ).
Nevertheless, as found by McNally et al. (2020), the contribution by buoyancy modes-when they are active-dominates the vortensity growth in the horseshoe region, evident by the vortensity excess being centered on the planet in the adiabatic models while the planet is held fixed (see Fig. 4).Once the planet is allowed to migrate, the contribution by the buoyancy response dominates in the radiative model as well (see Fig. 6).Overall, we expect that the vortensity growth due to the disk buoyancy response is the dominant mechanism in accelerating migration in our models, and that the contribution by the radial entropy gradient is secondary.

Observational constraints on the cooling timescale
In Sects. 3 & 4 we showed that, for the disk buoyancy response to drive meaningful vortensity generation, a very long cooling timescale (technically, a large ratio β tot /β buoy ) is required.This condition is expected to be met in the inner, optically thick region of protoplanetary disks, but it is worth pointing out that β tot will depend sensitively on the disk density and temperature as well as the properties of dust particles, which facilitate cooling.
Assuming a passive, irradiated disk model allows us to constrain its temperature structure (Chiang & Goldreich 1997), but the gas surface density and the fraction of small, efficiently-cooling dust grains are both largely unknown.With recent observational evidence1 pointing towards relatively low disk masses of ⟨M disk ⟩ ≲ 0.01 M ⋆ , and assuming a median disk cutoff radius of ⟨R out ⟩ ∼ 30 au, a median stellar mass of ⟨M ⋆ ⟩ ≈ 0.5 M ⊙ , and Σ ∝ R −1 , we find that the surface density at 1 au for a typical star should be of the order of ≲ 400 g/cm 2 .At the same time, dust growth models predict that a significant fraction of the submicron-micron-sized grains should be depleted within ∼ 100 local orbits as small grains coagulate to larger grains that no longer contribute to cooling, significantly lowering the opacity (for a review, see Birnstiel 2023).
As a result, we can expect that we could easily be overestimating the cooling timescale by one, two, or more orders of magnitude, and that buoyancy-driven torques should not be active for the typical planet-forming disk, or be limited to a very narrow radial range near the magnetically active inner rim of the disk (see e.g., Flock et al. 2017a).

Additional physical mechanisms
While the previous paragraph suggests that buoyancy torques might not be ubiquitous, this does not imply that the stalling effect of the dynamical corotation torque (DCT) can efficiently slow down inward migration.Ziampras et al. (2024) showed that for intermediate cooling timescales (β tot ∼ 0.3-100) the DCT is significantly weaker due to cooling-induced vortensity generation in the corotating region of a low-mass planet.As a result, a low-mass planet can be efficiently transported to the optically thick, slowly-cooling inner disk where β tot ∼ 1000 and buoyancy torques are active.The combined effect of these two mechanisms highlights the importance of radiation transport in hydrodynamical modeling and suggests that additional physical processes are necessary to prevent rapid inward migration and thus explain the population of super-Earths and sub-Neptunes at separations of ≲ 2 au from their host stars.
Directly relevant to our discussion of dust coagulation in Sect.5.2, Guilera et al. (2023) showed that a low-mass planet can efficiently trap a large amount of dust in its corotating region.The torque exerted by this concentration of solids was shown to be positive and comparable to or even larger in magnitude than the Lindblad torque for Stokes numbers of St ∼ 0.1, able to stop or even reverse migration.Given that St ∼ Σ −1 , achieving St ∼ 0.1 at 10 au would require the presence of ≳cm-sized grains, even with the lower surface densities discussed in Sect.5.2.This requirement is in tension with both observational evidence that point to a typical maximum grain size a max ≲ 3 mm (Sierra et al. 2021), as well as dust evolution models (a max ≲ 1 mm, see e.g., Birnstiel 2023;Dominik & Dullemond 2024).Nevertheless, dust dynamics could help slow down planet migration to an extent even for St ∼ 0.01.
Moreover, thermal torques due to the planet's accretion luminosity could further slow down or reverse the direction of migration (Benítez-Llambay et al. 2015;Masset 2017).Combining dust dynamics with an accretion luminosity due to a pebble flux, Cummins et al. (2022) showed that the picture complicates significantly as the edge of the horseshoe region becomes unstable to the Rossby Wave Instability (RWI, Lovelace et al. 1999), leading to vortices which can induce a stochastic element to the planet's migration (Lega et al. 2021).We note, however, that none of these models included both a dust component and a self-consistent treatment of radiative cooling.Combining these mechanisms in one model while maintaining control over all processes involved is a challenge in itself, and will be the focus of follow-up work.

Transition to type-II regime and the "inertial limit"
In the above, we used the term "low-mass" in a rather liberal way, referring in fact to planets well below the gap-opening or thermal mass M th ≈ h 3 M ⋆ and below the disk feedback mass (Goodman & Rafikov 2001;Rafikov 2002).Of course, in a cold, passively heated protoplanetary disk, where we typically expect h ∼ 0.02 at 1 au, even an Earth-mass-sized planet can imprint on its disk and open a partial gap, slowing down migration.This concept of an "inertial limit" (Hourigan & Ward 1984;Ward & Hourigan 1989;Ward 1997) could in theory prevent the rapid inward migration discussed above, although McNally et al. (2019a) showed that intermittent vortex formation and dissipation can sustain an inward track.Based on the above paragraphs, investigating the effect of dust dynamics and vortex activity in disk models with realistic hydrodynamics is necessary in approaching a definitive answer to how low-mass planets actually interact with their surrounding disk.

SUMMARY
We have investigated the presence of a negative dynamical corotation torque driven by the buoyancy response of a protoplanetary disk in the presence of a migrating planet (Zhu et al. 2012;Lubow & Zhu 2014;McNally et al. 2020).Our approach involved constructing a quasi-realistic yet computationally feasible disk model and then carrying out global, inviscid, radiation hydrodynamics simulations in the flux-limited diffusion (FLD) approximation with an embedded migrating planet.
We first computed the cooling timescale β tot for a disk model used by McNally et al. (2020) andYun et al. (2022), and confirmed that the ratio β buoy /β tot is ≲ 0.01 at z ∼ 2H for R ∼ 2 au, consistent with the quenching of buoyancy torques in their work.We then constructed a more realistic, passive, colder and flared disk model, and found that the conditions for the buoyancy torque to operate are more favorable at the same radius.However, due to the prohibitive computational cost of this model, we reach a compromise by constructing a hybrid setup that combines the reference parameters of the model used by Yun et al. (2022) with the radial structure and longer cooling timescales of our realistic setup.We used this model in our hydrodynamical simulations.
Our numerical calculations showed that the planet migrates inwards significantly faster than in equivalent 2D models, indicating that the buoyancy torque is indeed active in the inner disk even when radiative diffusion is considered.In our radiative models, the planet's migration track exhibits a "knee" during its orbital evolution.Further analysis showed that this happens over a libration timescale, as the corotating region shrinks to a tadpole-shaped libration island and the torque on the planet stabilizes to a more negative value sustained by the vortensity generation due to the dissipation of buoyancy modes inside the planet's corotating region.
At the same time, models with a local, β cooling prescription matched the migration rate of our adiabatic models for both PLUTO and FARGO3D.Comparing individual buoyancy modes between our models revealed that FLD not only damps the disk buoyancy response, but also stretches the modes in azimuth, resulting in fewer and weaker oscillations and therefore much less efficient vortensity growth in the planet's corotating region (consistent with Yun et al. 2022).These findings highlight that the diffusive component of realistic radiation transport is necessary to correctly capture the buoyancy response of the disk, and further that the buoyancy response of the disk is likely to be rapidly quenched for larger β buoy /β tot , as weaker and fewer buoyancy modes will be excited.
Comparing PLUTO with FARGO3D, we found that the former underestimates the buoyancy torque by a factor of ≈ 3.5, consistent with the findings of Ziampras et al. (2023).While in the adiabatic and locally cooled models this discrepancy is of numerical origin, the thermal diffusion due to FLD should result in a physical damping of the buoyancy response in our radiative models.It is therefore plausible that the true migration timescale of the planet in a radiative model with FARGO3D could be up to 3.5× shorter than what we find in our PLUTO models, although we expect the two codes to largely agree in the radiative case due to the damping associated with radiative diffusion.
Overall, our results show that the buoyancy torque can indeed reverse the dynamical corotation torque acting on a migrating planet in a realistic protoplanetary disk, and highlight the sensitivity of the buoyancy torque to the disk model and treatment of radiation transport.This suggests that planet migration in the low-turbulence dead zone might not ever slow down due to dynamical corotation torques (Paardekooper 2014), and that low-mass planets (i.e., below the feedback regime, see Hourigan & Ward 1984;Ward & Hourigan 1989;Ward 1997) will likely migrate rapidly to the inner edge of the dead zone unless additional physical processes are considered.With that in mind, investigating the interplay of dust dynamics, radiation transport, and nonlinear processes such as gap opening and vortex formation will be critical to understanding the origin of super-Earths and sub-Neptunes in the disk dead zone.

APPENDIX A: IMPLEMENTATION OF FLD SOLVER
Our FLD module closely follows the numerical implementation of Kolb et al. (2013).We nevertheless provide a brief description of our implementation, and carry out a supplementary validation test to ensure that it is correct.

A1 Numerics
Similarly to the implementation of Kolb et al. (2013) (see Sect. 3.2 therein), we write Eq. ( 5) in implicit form for E rad , taking into account the finite-volume nature of PLUTO.Following Commerçon et al. (2011), the nonlinear term T n+1 4 in Q rad is approximated as We then solve the resulting linear system of equations using the BiCGStab method, implemented in the sparse matrix solver module of the PetSc library (Balay et al. 1997(Balay et al. , 2024a,b),b).After obtaining the updated E rad , we update the gas temperature as

A2 Test problem: radiatively damped linear wave
To test our implementation, we design a 1D test problem similar to Benítez-Llambay et al. (2019) where we excite small perturbations that correspond to eigenvectors of the coupled system of equations ( 1), ( 5) & ( 6) and measure their oscillation frequency and damping rate.This both ensures that our implementation of the FLD solver for E rad in Eq. ( 5) is correct, and that the thermal energy and radiation energy field are coupled correctly in Eq. (6).

A2.1 Dispersion relation and eigenvectors
We assume that λ, κ P , and κ R are constant, and that all quantities q ∈ {ρ, u, P, E rad } can be expressed in the form q = q 0 + q 1 , with q 0 being a constant background state and q 1 = Aq ′ e ikx−ωt a wavelike perturbation with A ≪ 1.By plugging them into the set of Eqs.

A2.2 Numerical solution
In our numerical setup we consider a periodic Cartesian domain x ∈ [0, L] with L = 1 au and 1000 uniformly spaced grid cells, and define k = 2π/L.We set µ = 2.353, γ = 7/5, λ = 1/3, and T 0 = 30 K. The opacity κ R = κ P = κ is set to a constant value within each setup and varies between 0.003-300 cm 2 /g gas throughout the set of models we run.The background state is given by ρ 0 = {10 −13 , 10 −12 } g/cm 3 , u 0 = 0, and the perturbed quantities are initialized as q 1 = A Re (q ′ ) cos(kx) − Im (q ′ ) sin(kx) , where A = 10 −4 is the amplitude of the perturbation.We then integrate for a time corresponding to 150 years.Finally, we extract the oscillation frequency ω osc and damping rate ω damp from the perturbed quantities by fitting the gas pressure at x = 0 as a function of time with a function of the form P(t) = P 0 + c 0 P 0 sin(ω osc t + c 1 )e ω damp t , (A9) where c 0 and c 1 are additional, free fitting parameters to account for the amplitude and phase of the excited eigenmode.An example of the fitted oscillation is shown in Fig. A1 for all quantities.
We then present our results in Fig. A2, showing that we recover the expected damping and oscillation frequencies to very good degree for all values of κ and both reference values of ρ 0 .We therefore conclude that our FLD module is suitable for use in our models.
This paper has been typeset from a T E X/L A T E X file prepared by the author.We normalize ω to ω iso = kc s,iso 0 .The oscillation frequency fluctuates between the isothermal (ω iso ) and adiabatic expectation ( √ γω iso , dashed black line), depending on the coupling between the thermal and radiation energy fields.

Figure 1 .
Figure 1.Cooling timescale β tot (top) and the ratio to the buoyancy timescale β buoy /β tot (bottom) for different disk models.Left: the model of McNally et al. (2020) and Yun et al. (2022).Middle: a more realistic model of a passive, irradiated disk.Right: a hybrid model that combines the computational cost and reference parameters of model A with the radiative environment of model B. Green contours in the upper panels denote the "optical surface" of the disk, where τ irr (r) = r R⊙ κ P (T ⊙ )ρ dr = 1.Solid, dashed, and dotted black contours in the lower panels denote β buoy /β tot = 1, 0.1, and 0.01, respectively.

Figure 2 .
Figure2.Comparison of the ratio β buoy /β tot for the realistic and hybrid disk models in our region of interest, where buoyancy torques could operate (R ≲ 4 au).Top: radial profiles of β buoy /β tot , showing a near-perfect match for R ≈ 1-2 au.Bottom: vertical profile of β buoy /β tot at R = 2 au.
which achieves a resolution of 16 cells per scale height at R 0 in all directions.Similar to McNally et al. (2020) and Ziampras et al. (2023), we employ wave-damping zones near the radial boundaries of the domain following de Val-Borro et al. (

Figure 3 .
Figure 3. Azimuthally averaged perturbed vortensity for our models at t = 200 orbits.A vortensity excess inside the horseshoe region (marked by the dashed lines) is indicative of the dissipation of buoyancy oscillations in 3D runs.

Figure 4 .
Figure 4. Perturbed vortensity for a selection of our 3D models at t = 200 orbits.The vortensity excess in the adiabatic models is more centered to the planet's location, whereas the excess traces the edge of the horseshoe region in the radiative model.The planet is located at ϕ = 0 and moves towards the left.

Figure 5 .
Figure 5. Migration tracks for all models.Top: semimajor axis of the planet as a function of time.Bottom: the migration timescale τ mig = R p / Ṙp .Different colors represent different radiative prescriptions.Faint curves denote raw data, while solid curves denote a 20-orbit rolling average.Black dotted lines show the type-I migration track and cooling timescale.All 2D models overlap nearly perfectly.

Figure 6 .
Figure6.Vortensity heatmaps for the radiative 3D model at several snapshots during the planet's migration, highlighting the section of the corotating region where gas maintains its original vortensity.That section initially spans the full width of the horseshoe region, but shrinks to approximately a third of its initial size as the flow changes while the planet migrates.The planet then accelerates once this transition is complete.

Figure 7 .
Figure 7. Azimuthal slice of the gas vertical velocity u z /c s at R = R p − x h for the disk models of Yun et al. (2022) and our hybrid model.A vertical line marks the planet's azimuthal location.Buoyancy modes are more strongly damped in radiative models, and have a longer azimuthal wavelength for the model in Yun et al. (2022).

Figure A1 .Figure A2 .
Figure A1.An example of fitted oscillations for our model with ρ 0 = 10 −12 g/cm 3 and κ = 1 cm 2 /g at x = 0.All quantities have been normalized and shifted for clarity.Gray dots denote the analytical solutions.