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

Alexandros Ziampras1, Richard P. Nelson1, Sijme-Jan Paardekooper1,2
1Astronomy Unit, School of Physics and Astronomy, Queen Mary University of
London, London E1 4NS, UK
2Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg
1, 2600 AA Delft, The Netherlands
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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.

###### keywords:

planet–disc interactions — protoplanetary discs — hydrodynamics — radiation:
dynamics — methods: numerical

††pubyear: 2023††pagerange: Buoyancy torques prevent low-mass planets from
stalling in low-turbulence radiative disks–9

##  1 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
vortex-assisted 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 snow-plough 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 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.

##  2 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.

###  2.1 Equations of hydrodynamics

We consider a disk of ideal gas with adiabatic index
$\gamma=\nicefrac{{7}}{{5}}$ and mean molecular weight $\mu=2.353$, orbiting a
star with mass $M_{\star}=\mathrm{M}_{\odot}$ and luminosity
$L_{\star}=\mathrm{L}_{\odot}$. The volume density $\rho$, velocity field
$\bm{u}$, and pressure $P$ of the gas evolve according to the (inviscid) Euler
equations

$\displaystyle\frac{\partial{\rho}}{\partial{t}}+\bm{u}\cdot\nabla\rho=-\rho\nabla\cdot\bm{u},$
(1a)
$\displaystyle\frac{\partial{\bm{u}}}{\partial{t}}+(\bm{u}\cdot\nabla)\bm{u}=-\frac{1}{\rho}\nabla
P-\nabla\Phi,$ (1b)
$\displaystyle\frac{\partial{e}}{\partial{t}}+\bm{u}\cdot\nabla e=-\gamma
e\nabla\cdot\bm{u}+Q,$ (1c)

where $e=P/(\gamma-1)$ is the internal energy density given by the ideal gas
law, $\Phi=\Phi_{\star}=\text{G}M_{\star}/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_{\mathrm{s,iso}}=\sqrt{P/\rho}$, and the temperature is
$T=\mu c_{\mathrm{s,iso}}^{2}/\mathcal{R}$, with $\mathcal{R}$ denoting the
gas constant.

In a cylindrical coordinate system $(R,\varphi,z)$ with
$r=\sqrt{R^{2}+z^{2}}$, by assuming that the disk is non-accreting and
axisymmetric and further requiring that the midplane density
$\rho_{\mathrm{mid}}$ and vertically constant temperature follow radial power-
law profiles such that

$\rho_{\mathrm{mid}}(R)=\rho_{0}\left(\frac{R}{R_{0}}\right)^{p},\quad
T(R)=T_{0}\left(\frac{R}{R_{0}}\right)^{q},$ (2)

the disk structure in equilibrium can be described by (Nelson et al., 2013)

$\displaystyle\rho^{\text{eq}}(R,z)$
$\displaystyle=\rho_{\mathrm{mid}}(R)\,\exp\left[-\frac{1}{h^{2}}\left(1-\frac{R}{r}\right)\right],$
(3) $\displaystyle u_{\phi}(R,z)$
$\displaystyle=R\Omega_{\mathrm{K}}\left[1+(p+q)h^{2}+q\left(1-\frac{R}{r}\right)\right]^{1/2}.$
(4)

Here, $\Omega_{\mathrm{K}}=\sqrt{\text{G}M_{\star}/R^{3}}$ is the Keplerian
angular velocity, $h=H/R$ is the aspect ratio of the disk, and
$H=c_{\mathrm{s,iso}}/\Omega_{\mathrm{K}}$ is the pressure scale height.
Finally, we can define the surface density through
$\Sigma=\int_{-\infty}^{\infty}\rho\,\mathrm{d}z$.

###  2.2 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_{\mathrm{rad}}$ evolves as

$\frac{\partial{E_{\mathrm{rad}}}}{\partial{t}}+\nabla\cdot\bm{F}=\kappa_{\mathrm{P}}\rho
c\left(a_{\mathrm{R}}T^{4}-E_{\mathrm{rad}}\right),$ (5)

and is coupled to the gas energy density with a source term

$Q_{\mathrm{rad}}=-\kappa_{\mathrm{P}}\rho
c\left(a_{\mathrm{R}}T^{4}-E_{\mathrm{rad}}\right)$ (6)

in Eq. (1c). Here, $\kappa_{\mathrm{P}}$ is the Planck mean opacity, $c$ is
the speed of light, and $a_{\mathrm{R}}$ is the radiation constant. The
radiation flux $\bm{F}$ is given by

$\bm{F}=-\frac{\lambda c}{\kappa_{\mathrm{R}}\rho}\nabla E_{\mathrm{rad}},$
(7)

where $\kappa_{\mathrm{R}}$ is the Rosseland mean opacity and $\lambda$ is the
flux limiter following Kley (1989):

$\lambda=\begin{cases}\frac{2}{3+\sqrt{9+10x^{2}}},&x\leq 2,\\\
\frac{10}{10x+9+\sqrt{180x+81}},&x>2.\end{cases},\quad x:=\frac{|\nabla
E_{\mathrm{rad}}|}{\kappa_{\mathrm{R}}\rho E_{\mathrm{rad}}}.$ (8)

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\sim 4H$, see Chiang & Goldreich, 1997), this is not a
problem as we are mainly interested in the region $z\lesssim 2H$. We can
sidestep this issue by prescribing $E_{\mathrm{rad}}=a_{\mathrm{R}}T_{0}^{4}$
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
$\beta_{\mathrm{cool}}=t_{\mathrm{cool}}\Omega_{\mathrm{K}}$ can be
approximated following Flock et al. (2017b)

$\beta_{\mathrm{cool}}=\frac{\Omega_{\mathrm{K}}}{\eta}\left(H^{2}+\frac{l_{\mathrm{rad}}^{2}}{3}\right),\qquad\eta=\frac{16\sigma_{\mathrm{SB}}T^{3}}{3\kappa_{\mathrm{P}}\rho^{2}c_{\mathrm{v}}},\quad
l_{\mathrm{rad}}=\frac{1}{\kappa_{\mathrm{P}}\rho}.$ (9)

Here, $\sigma_{\mathrm{SB}}$ is the Stefan–Boltzmann constant, and
$c_{\mathrm{v}}=\mathcal{R}/(\mu(\gamma-1))$ is the specific heat at constant
volume.

The above expression for $\beta$ is valid under the assumption that the solid
particles (“dust”) in the disk, which facilitate cooling, are well-coupled to
the gas. In particular, we focus on the small grains with size
$a_{\mathrm{d}}=0.1\,\mathrm{\mu m}$, as they provide the most efficient
cooling channel. We can then express the collision timescale
$\beta_{\mathrm{coll}}$ between the gas and small grains following Dullemond
et al. (2022)

$\beta_{\mathrm{coll}}=\frac{c_{\mathrm{v}}e}{\epsilon\rho^{2}\bar{C}_{H}}\frac{m_{\mathrm{d}}}{s_{\mathrm{d}}}\Omega_{\mathrm{K}},\qquad\bar{C}_{H}=c_{\mathrm{v}}\frac{c_{\mathrm{s,iso}}}{\sqrt{2\pi}}$
(10)

Here, $s_{\mathrm{d}}=4\pi a_{\mathrm{d}}^{2}$ and
$m_{\mathrm{d}}=\frac{4\pi}{3}\rho_{\mathrm{d}}a_{\mathrm{d}}^{3}$ are the
surface area and mass of a single dust particle with bulk density
$\rho_{\mathrm{d}}=2.08\,\mathrm{g}/\mathrm{cm}^{3}$, and $\epsilon$ is the
dust-to-gas mass ratio. The combined cooling timescale $\beta_{\mathrm{tot}}$
is then given by

$\beta_{\mathrm{tot}}=\beta_{\mathrm{cool}}+\beta_{\mathrm{coll}}.$ (11)

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
$\kappa_{\mathrm{R}}=\kappa_{\mathrm{P}}$.

###  2.3 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 ($\beta_{\mathrm{tot}}\to\infty$). 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 ($\varpi$) 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)

$\Gamma_{h}=2\pi\left(1-\frac{\varpi(R_{\mathrm{p}})}{\varpi_{h}}\right)\Sigma_{\mathrm{p}}R_{\mathrm{p}}^{2}{x}_{\mathrm{h}}\Omega_{\mathrm{p}}\left(\frac{\text{d}{R_{\mathrm{p}}}}{\text{d}{t}}-u_{R}\right).$
(12)

Here, a subscript “p” denotes quantities at the planet’s location, $u_{R}$ is
the radial velocity of the background disk, and ${x}_{\mathrm{h}}$ is the
horseshoe half-width (Paardekooper et al., 2010b) defined as

${x}_{\mathrm{h}}=\frac{1.1}{\gamma^{1/4}}\left(\frac{0.4}{\epsilon/H_{\text{p}}}\right)^{1/4}\sqrt{\frac{M_{\mathrm{p}}}{h_{\mathrm{p}}M_{\star}}}R_{\mathrm{p}},\qquad\epsilon=0.6H_{\text{p}}.$
(13)

The smoothing length $\epsilon$ is chosen to match 3D models (Müller & Kley,
2012).

In Eq. (12), the vortensity is given in the vertically integrated
approximation by $\varpi=(\nabla\times\bm{u})/\Sigma\cdot\hat{z}$, while
Masset & Benítez-Llambay (2016) derived the equivalent quantity in 3D to be

$\varpi=\left[\int\limits_{-\infty}^{\infty}\frac{\rho}{(\nabla\times\bm{u})\cdot\hat{z}}\mathrm{d}z\right]^{-1}.$
(14)

Then, $\varpi(R_{\mathrm{p}})$ is the vortensity at $R_{\mathrm{p}}$ for an
unperturbed disk, and $\varpi_{\mathrm{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 $\varpi_{h}$ in the absence of
turbulent viscosity or vortensity-generating baroclinic effects. This in turns
leads to evolution of the quantity ${\varpi(R_{\mathrm{p}})}/{\varpi_{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
$\varpi_{0}(R)=\frac{1}{2}\Omega_{\mathrm{K}}/\Sigma^{\text{eq}}$, which would
result in a continuously increasing positive torque as the planet migrates
inwards for $\Sigma(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
($\beta_{\text{diff}}=\Omega_{\mathrm{K}}H^{2}/\eta$, see Eq. (9)) is
comparable to or shorter than the buoyancy timescale $\beta_{\mathrm{buoy}}$.
The latter is related to the Brunt–Väisälä frequency

$N=\sqrt{-\frac{1}{\gamma}\frac{\partial{\Phi_{\star}}}{\partial{z}}\frac{\partial{}}{\partial{z}}\left[\ln\left(\frac{P}{\rho^{\gamma}}\right)\right]}=\Omega_{\mathrm{K}}\sqrt{\frac{\gamma-1}{\gamma}}\frac{z}{H}\left(\frac{R}{r}\right)^{3},$
(15)

through which we can write

$\beta_{\mathrm{buoy}}=\frac{\Omega_{\mathrm{K}}}{N}=\sqrt{\frac{\gamma}{\gamma-1}}\frac{H}{z}\left(1+\frac{z}{R}\right)^{3/2}.$
(16)

In their simulations and for their disk model, Yun et al. (2022) found that
$\beta_{\text{diff}}\lesssim\beta_{\mathrm{buoy}}$ for $z\gtrsim 2H$, and
concluded that the buoyancy-driven torque is quenched in radiative disks.
However, as stated above, $\beta_{\text{diff}}$ is a sensitive function of the
underlying disk model. It is therefore necessary to investigate the conditions
under which the buoyancy-driven torque is relevant in a more slowly cooling,
passively irradiated disk model (e.g., Chiang & Goldreich, 1997).

##  3 Cooling timescale maps

In this section we present maps of the cooling timescale
$\beta_{\mathrm{tot}}$ and the ratio
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$ for several different disk
models. We use these maps to motivate our choice of numerical parameters for
our hydrodynamical simulations.

###  3.1 Connecting to previous work

We first compute the cooling timescale $\beta_{\mathrm{tot}}$ for the disk
model of McNally et al. (2020) and Yun et al. (2022). Their model assumes a
disk with
$\Sigma(R)=3400\,\left(R/\text{au}\right)^{-1/2}\,\text{g}/\text{cm}^{2}$, a
constant aspect ratio $h=0.05$ (i.e., $T\propto 1/R$) and a constant dust-to-
gas mass ratio $\epsilon=0.01$ for submicron grains. Combining Eqs. (2), (3),
(9), and (10) we compute the cooling timescale
$\beta_{\mathrm{tot}}=\beta_{\mathrm{cool}}+\beta_{\mathrm{coll}}$ and the
ratio $\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$ for this disk model, and
present the results in the left panels of Fig. 1.

Figure 1: Cooling timescale $\beta_{\mathrm{tot}}$ (top) and the ratio to the
buoyancy timescale $\beta_{\mathrm{buoy}}/\beta_{\mathrm{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
$\tau_{\mathrm{irr}}(r)=\int_{R_{\odot}}^{r}\kappa_{\mathrm{P}}(T_{\odot})\rho\,\mathrm{d}r=1$.
Solid, dashed, and dotted black contours in the lower panels denote
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}=1$, 0.1, and 0.01, respectively.

In agreement with the findings of Yun et al. (2022), we find that
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}\sim 0.1$ at $z\sim 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 $\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}\gtrsim 0.01$
between $z\sim 1\text{--}2H$ throughout the whole disk, and therefore expect
that buoyancy torques are indeed not relevant for this disk model.

###  3.2 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 scale-free environment,
is not quite representative of a real PPD. Passively heated disks are flared
($q\approx-1/2$, Chiang & Goldreich, 1997) and significantly thinner
($h\approx 0.02\text{--}0.025$ at 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
$\alpha\sim 10^{-3}\text{--}10^{-2}$ (Shakura & Sunyaev, 1973). This is both
difficult to reconcile with the low turbulence expected in the dead zone (Bai
& Stone, 2013) and inconsistent with the assumption that the models by both
McNally et al. (2020) and Yun et al. (2022) reflect laminar disks, where
buoyancy torques can operate in the first place.

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\,\left(R/\text{au}\right)^{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
$\Sigma(R)=1700\,\left(R/\text{au}\right)^{-1}\,\text{g}/\text{cm}^{2}$,
matching the Minimum Mass Solar Nebula model (Hayashi, 1981) at $R=1$ au. We
then recompute the cooling timescale $\beta_{\mathrm{tot}}$ and the ratio
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$, and present the results in the
middle panels of Fig. 1.

Given that $\beta_{\text{diff}}\propto T^{-3}\propto 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 $\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}\lesssim
10^{-3}$ up to $z\sim 2\text{--}3H$ for $R\lesssim 2\,\text{au}$, which could
be sufficient for buoyancy torques to operate in the inner disk.

###  3.3 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 migration 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 ($\sim$16 cells per scale
height) to capture the effect of buoyancy torques in the inner disk ($h\sim
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 $\Sigma\propto 1/R$ as
before, but adjust $\Sigma_{0}$ and $h_{0}$ such that the values of $\Sigma$
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
$\Sigma_{0}=4800~{}\text{g}/\text{cm}^{2}~{}(R/\text{au})^{-1}$ and
$h=0.041~{}(R/\text{au})^{2/7}$.

Of course, adjusting both $\Sigma$ 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 $\approx 7$. The result is a near-perfect match
between the two models in terms of the ratio
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{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.

Figure 2: Comparison of the ratio $\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$
for the realistic and hybrid disk models in our region of interest, where
buoyancy torques could operate ($R\lesssim 4$ au). Top: radial profiles of
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$, showing a near-perfect match for
$R\approx 1$–2 au. Bottom: vertical profile of
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$ at $R=2$ au.

##  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.

###  4.1 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):

$\frac{\partial{s}}{\partial{t}}+\bm{u}\cdot\nabla s=Q^{\prime},\qquad s\equiv
c_{\mathrm{v}}\log\left[\frac{T}{T_{\text{ref}}}\left(\frac{\rho}{\rho_{\text{ref}}}\right)^{1-\gamma}\right],$
(17)

where $T_{\text{ref}}$ and $\rho_{\text{ref}}$ are arbitrary constants. The
thermal relaxation source term $Q^{\prime}$, 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\pi$. We use a grid of $N_{r}\times
N_{\theta}\times N_{\phi}=528\times 80\times 2048$ cells, 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.
(2006) and smooth the planet’s potential with the cubic spline curve described
in Klahr & Kley (2006). The radiation energy is initialized as
$E_{\mathrm{rad,0}}=a_{\mathrm{R}}T_{0}^{4}$. At the radial and vertical edges
of the domain, all parameters are set to their initial values. For
$E_{\mathrm{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_{\mathrm{p}}=2\times 10^{-5}\,M_{\star}\approx 6.7\,\mathrm{M}_{\oplus}$
using the formula by de Val-Borro et al. (2006). Since we do not consider the
disk self-gravity, 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^{\prime}$ in Eqs. (1c) &
(17) are set to zero. In our runs with a simplified local $\beta$ cooling
prescription we implement a thermal relaxation term similar to Gammie (2001)
and McNally et al. (2020) for PLUTO and FARGO3D, respectively:

$Q_{\mathrm{relax}}=-\rho
c_{\mathrm{v}}\frac{T-T_{0}}{\beta_{\mathrm{tot}}}\Omega_{\mathrm{K}},\qquad
Q_{\mathrm{relax}}^{\prime}=-\frac{s(\rho,T)-s(\rho,T_{0})}{\beta_{\mathrm{tot}}}\Omega_{\mathrm{K}},$
(18)

where $\beta_{\mathrm{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,\phi)$ models
using PLUTO. Similar to the 3D runs, we employ the following three cooling
prescriptions:

  * •

_adiabatic_ : no cooling ($Q=0$),

  * •

_local cooling_ : $\beta$ 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 $\epsilon=0.6H_{\mathrm{p}}$ (Müller & Kley, 2012).

###  4.2 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 $\beta$-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.

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.

We then plot perturbed vortensity heatmaps for the adiabatic and radiative 3D
models in Fig. 4. Models with $\beta$ 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_{\mathrm{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.

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 $\phi=0$ and moves towards
the left.

###  4.3 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 $\tau_{\mathrm{mig}}=R_{\mathrm{p}}/\dot{R}_{\mathrm{p}}$ are shown
in Fig. 5.

Figure 5: Migration tracks for all models. Top: semimajor axis of the planet
as a function of time. Bottom: the migration timescale
$\tau_{\mathrm{mig}}=R_{\mathrm{p}}/\dot{R}_{\mathrm{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.

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 $\approx 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.

Figure 6: 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.

While the corotating region of the planet initially spans the full azimuthal
range, it shrinks to a tadpole-shaped libration island (Masset & Papaloizou,
2003; Papaloizou et al., 2007) ahead of the planet as the latter migrates
inwards (see also Fig. 18 in McNally et al., 2020). At the same time,
vortensity is continuously generated due to the dissipation of buoyancy modes
close to the planet. Tracing the section of the corotating region where the
gas maintains its original vortensity $\varpi_{\mathrm{ref}}=\varpi(R_{0})$
shows that this section shrinks abruptly to approximately a third of its
initial azimuthal extent after $\approx 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

$\tau_{\mathrm{lib}}=\frac{8\pi
R_{\mathrm{p}}}{3\Omega_{\mathrm{p}}{x}_{\mathrm{h}}}\approx
72~{}\mathrm{orbits},$ (19)

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 $\approx 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
$\tau^{\mathrm{Pluto}}_{\mathrm{mig}}/\tau^{\mathrm{Fargo}}_{\mathrm{mig}}\approx
3.5$ between the two codes for both adiabatic and $\beta$-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$\times$ 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.

###  4.4 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
$\beta_{\mathrm{buoy}}/\beta_{\text{diff}}\lesssim 0.1$ at $z\sim 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 amplitude 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_{\mathrm{p}}-{x}_{\mathrm{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_{\mathrm{p}}+{x}_{\mathrm{h}}$.

  * •

Cooling damps buoyancy modes, more so for fully radiative (green) than for
$\beta$ models (orange). This further highlights the need for realistic
radiative modeling, and is consistent with the faster inward migration found
for $\beta$ models in Fig. 5.

  * •

Buoyancy modes in our radiative models (green) are more spaced out in azimuth
compared to the adiabatic and $\beta$-cooled cases. This is consistent with
the findings of Yun et al. (2022), and is most likely due to the diffusive
nature of FLD.

Figure 7: Azimuthal slice of the gas vertical velocity $u_{z}/c_{\mathrm{s}}$
at $R=R_{\mathrm{p}}-{x}_{\mathrm{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).

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 $\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}\gtrsim 0.01$ at $z\sim
2H$, but finding the exact value is beyond the scope of this work.

##  5 Discussion

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

###  5.1 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/\rho^{\gamma}\propto R^{0.48}$ and $P_{\text{2D}}/\Sigma^{\gamma}\propto
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/\rho^{\gamma}\propto 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.

###  5.2 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 $\beta_{\mathrm{tot}}/\beta_{\mathrm{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 $\beta_{\mathrm{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 evidence111Based on unpublished data from
the AGE-PRO ALMA large program, private communication. pointing towards
relatively low disk masses of $\langle M_{\mathrm{disk}}\rangle\lesssim
0.01\,M_{\star}$, and assuming a median disk cutoff radius of $\langle
R_{\mathrm{out}}\rangle\sim 30$ au, a median stellar mass of $\langle
M_{\star}\rangle\approx 0.5\,\mathrm{M}_{\odot}$, and $\Sigma\propto R^{-1}$,
we find that the surface density at 1 au for a typical star should be of the
order of $\lesssim 400$ g/cm2. At the same time, dust growth models predict
that a significant fraction of the submicron–micron-sized grains should be
depleted within $\sim 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).

###  5.3 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
($\beta_{\mathrm{tot}}\sim 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 $\beta_{\mathrm{tot}}\sim
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 $\lesssim 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
$\mathrm{St}\sim 0.1$, able to stop or even reverse migration. Given that
$\mathrm{St}\sim\Sigma^{-1}$, achieving $\mathrm{St}\sim 0.1$ at 10 au would
require the presence of $\gtrsim$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_{\mathrm{max}}\lesssim 3$ mm (Sierra et al., 2021), as well as dust
evolution models ($a_{\mathrm{max}}\lesssim 1$ mm, see e.g., Birnstiel, 2023;
Dominik & Dullemond, 2024). Nevertheless, dust dynamics could help slow down
planet migration to an extent even for $\mathrm{St}\sim 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.

###  5.4 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_{\mathrm{th}}\approx h^{3}M_{\star}$ 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\sim 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.

##  6 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 $\beta_{\mathrm{tot}}$ for a disk
model used by McNally et al. (2020) and Yun et al. (2022), and confirmed that
the ratio $\beta_{\mathrm{buoy}}/\beta_{\mathrm{tot}}$ is $\lesssim 0.01$ at
$z\sim 2H$ for $R\sim 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, $\beta$ 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
$\beta_{\mathrm{buoy}}/\beta_{\mathrm{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 $\approx 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$\times$ 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.

## Acknowledgements

AZ would like to thank Will Béthune for their help in implementing and testing
the FLD module, and Mario Flock for their suggestions and helpful discussions.
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.
AZ and RPN are supported by STFC grants ST/T000341/1 and ST/X000931/1. This
project has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No 101054502). All plots in this paper were made with the Python
library matplotlib (Hunter, 2007).

## Data Availability

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

## References

  * Bai & Stone (2013) Bai X.-N., Stone J. M., 2013, ApJ, 769, 76
  * Balay et al. (1997) Balay S., Gropp W. D., McInnes L. C., Smith B. F., 1997, in Arge E., Bruaset A. M., Langtangen H. P., eds, Modern Software Tools in Scientific Computing. Birkhäuser Press, pp 163–202
  * Balay et al. (2024a) Balay S., et al., 2024a, PETSc Web page, https://petsc.org/, https://petsc.org/
  * Balay et al. (2024b) Balay S., et al., 2024b, Technical Report ANL-21/39 - Revision 3.21, PETSc/TAO Users Manual. Argonne National Laboratory, doi:10.2172/2205494
  * Baruteau & Masset (2008) Baruteau C., Masset F., 2008, ApJ, 678, 483
  * Benítez-Llambay & Masset (2016) Benítez-Llambay P., Masset F. S., 2016, ApJS, 223, 11
  * Benítez-Llambay et al. (2015) Benítez-Llambay P., Masset F., Koenigsberger G., Szulágyi J., 2015, Nature, 520, 63
  * Benítez-Llambay et al. (2019) Benítez-Llambay P., Krapp L., Pessah M. E., 2019, ApJS, 241, 25
  * Birnstiel (2023) Birnstiel T., 2023, arXiv e-prints, p. arXiv:2312.13287
  * Chiang & Goldreich (1997) Chiang E. I., Goldreich P., 1997, ApJ, 490, 368
  * Commerçon et al. (2011) Commerçon B., Teyssier R., Audit E., Hennebelle P., Chabrier G., 2011, A&A, 529, A35
  * Crida et al. (2022) Crida A., Griveaud P., Lega E., Masset F., Morbidelli A., Kloster D., Marques L., Minker K., 2022, in Richard J., et al., eds, SF2A-2022: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. Eds.: J. Richard. pp 315–317
  * Cummins et al. (2022) Cummins D. P., Owen J. E., Booth R. A., 2022, MNRAS, 515, 1276
  * Dominik & Dullemond (2024) Dominik C., Dullemond C. P., 2024, A&A, 682, A144
  * Dullemond et al. (2022) Dullemond C. P., Ziampras A., Ostertag D., Dominik C., 2022, A&A, 668, A105
  * Flock et al. (2017a) Flock M., Fromang S., Turner N. J., Benisty M., 2017a, ApJ, 835, 230
  * Flock et al. (2017b) Flock M., Nelson R. P., Turner N. J., Bertrang G. H. M., Carrasco-González C., Henning T., Lyra W., Teague R., 2017b, ApJ, 850, 131
  * Gammie (2001) Gammie C. F., 2001, ApJ, 553, 174
  * Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, ApJ, 233, 857
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Goodman & Rafikov (2001) Goodman J., Rafikov R. R., 2001, ApJ, 552, 793
  * Gressel et al. (2015) Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015, ApJ, 801, 84
  * Guilera et al. (2023) Guilera O. M., Benitez-Llambay P., Miller Bertolami M. M., Pessah M. E., 2023, ApJ, 953, 97
  * Harten (1983) Harten A., 1983, Journal of Computational Physics, 49, 357
  * Hayashi (1981) Hayashi C., 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Hourigan & Ward (1984) Hourigan K., Ward W. R., 1984, Icarus, 60, 29
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Klahr & Kley (2006) Klahr H., Kley W., 2006, A&A, 445, 747
  * Kley (1989) Kley W., 1989, A&A, 208, 98
  * Kolb et al. (2013) Kolb S. M., Stute M., Kley W., Mignone A., 2013, A&A, 559, A80
  * Lega et al. (2021) Lega E., et al., 2021, A&A, 646, A166
  * Lega et al. (2022) Lega E., Morbidelli A., Nelson R. P., Ramos X. S., Crida A., Béthune W., Batygin K., 2022, A&A, 658, A32
  * Levermore & Pomraning (1981) Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321
  * Lin & Papaloizou (1985) Lin D. N. C., Papaloizou J., 1985, in Black D. C., Matthews M. S., eds, Protostars and Planets II. pp 981–1072
  * Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
  * Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
  * Lubow & Zhu (2014) Lubow S. H., Zhu Z., 2014, ApJ, 785, 32
  * Masset (2000) Masset F., 2000, A&AS, 141, 165
  * Masset (2017) Masset F. S., 2017, MNRAS, 472, 4204
  * Masset & Benítez-Llambay (2016) Masset F. S., Benítez-Llambay P., 2016, ApJ, 817, 19
  * Masset & Casoli (2009) Masset F. S., Casoli J., 2009, ApJ, 703, 857
  * Masset & Papaloizou (2003) Masset F. S., Papaloizou J. C. B., 2003, ApJ, 588, 494
  * McNally et al. (2017) McNally C. P., Nelson R. P., Paardekooper S.-J., Gressel O., Lyra W., 2017, MNRAS, 472, 1565
  * McNally et al. (2019a) McNally C. P., Nelson R. P., Paardekooper S.-J., Benítez-Llambay P., 2019a, MNRAS, 484, 728
  * McNally et al. (2019b) McNally C. P., Nelson R. P., Paardekooper S.-J., 2019b, MNRAS, 489, L17
  * McNally et al. (2020) McNally C. P., Nelson R. P., Paardekooper S.-J., Benítez-Llambay P., Gressel O., 2020, MNRAS, 493, 4382
  * Menou & Goodman (2004) Menou K., Goodman J., 2004, ApJ, 606, 520
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, The Astrophysical Journal Supplement Series, 170, 228
  * Mignone et al. (2012) Mignone A., Flock M., Stute M., Kolb S. M., Muscianisi G., 2012, A&A, 545, A152
  * Müller & Kley (2012) Müller T. W. A., Kley W., 2012, A&A, 539, A18
  * Nelson et al. (2013) Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610
  * Ogilvie & Lubow (2002) Ogilvie G. I., Lubow S. H., 2002, MNRAS, 330, 950
  * Paardekooper (2014) Paardekooper S. J., 2014, MNRAS, 444, 2031
  * Paardekooper et al. (2010a) Paardekooper S. J., Baruteau C., Crida A., Kley W., 2010a, MNRAS, 401, 1950
  * Paardekooper et al. (2010b) Paardekooper S.-J., Lesur G., Papaloizou J. C. B., 2010b, ApJ, 725, 146
  * Paardekooper et al. (2011) Paardekooper S. J., Baruteau C., Kley W., 2011, MNRAS, 410, 293
  * Paardekooper et al. (2022) Paardekooper S.-J., Dong R., Duffell P., Fung J., Masset F. S., Ogilvie G., Tanaka H., 2022, arXiv e-prints, p. arXiv:2203.09595
  * Papaloizou et al. (2007) Papaloizou J. C. B., Nelson R. P., Kley W., Masset F. S., Artymowicz P., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. p. 655 (arXiv:astro-ph/0603196), doi:10.48550/arXiv.astro-ph/0603196
  * Pierens (2015) Pierens A., 2015, MNRAS, 454, 2003
  * Rafikov (2002) Rafikov R. R., 2002, ApJ, 569, 997
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Sierra et al. (2021) Sierra A., et al., 2021, ApJS, 257, 14
  * Trac & Pen (2004) Trac H., Pen U.-L., 2004, New Astron., 9, 443
  * Van Leer (1974) Van Leer B., 1974, Journal of computational physics, 14, 361
  * Wang et al. (2019) Wang L., Bai X.-N., Goodman J., 2019, ApJ, 874, 90
  * Ward (1997) Ward W. R., 1997, Icarus, 126, 261
  * Ward & Hourigan (1989) Ward W. R., Hourigan K., 1989, ApJ, 347, 490
  * Yamaleev & Carpenter (2009) Yamaleev N. K., Carpenter M. H., 2009, Journal of Computational Physics, 228, 3025
  * Yun et al. (2022) Yun H.-G., Kim W.-T., Bae J., Han C., 2022, ApJ, 938, 102
  * Zhu et al. (2012) Zhu Z., Stone J. M., Rafikov R. R., 2012, ApJ, 758, L42
  * Ziampras et al. (2020) Ziampras A., Ataiee S., Kley W., Dullemond C. P., Baruteau C., 2020, A&A, 633, A29
  * Ziampras et al. (2023) Ziampras A., Paardekooper S.-J., Nelson R. P., 2023, MNRAS, 525, 5893
  * Ziampras et al. (2024) Ziampras A., Nelson R. P., Paardekooper S.-J., 2024, arXiv e-prints, p. arXiv:2402.00125
  * de Val-Borro et al. (2006) de Val-Borro M., et al., 2006, MNRAS, 370, 529

##  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.

###  A.1 Numerics

Similarly to the implementation of Kolb et al. (2013) (see Sect. 3.2 therein),
we write Eq. (5) in implicit form for $E_{\mathrm{rad}}$, taking into account
the finite-volume nature of PLUTO. Following Commerçon et al. (2011), the
nonlinear term $\left(T^{n+1}\right)^{4}$ in $Q_{\mathrm{rad}}$ is
approximated as

$\left(T^{n+1}\right)^{4}\approx
4\left(T^{n}\right)^{3}T^{n+1}-3\left(T^{n}\right)^{4}.$ (20)

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, 2024a, 2024b). After obtaining the updated
$E_{\mathrm{rad}}$, we update the gas temperature as

$T^{n+1}=\frac{\kappa_{\mathrm{P}}c\left(3a_{\mathrm{R}}T^{4}+E_{\mathrm{rad}}^{n+1}\right)\Delta
t+c_{\mathrm{v}}T}{c_{\mathrm{v}}+4\kappa_{\mathrm{P}}ca_{\mathrm{R}}T^{3}\Delta
t}.$ (21)

###  A.2 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_{\mathrm{rad}}$ in Eq. (5) is correct,
and that the thermal energy and radiation energy field are coupled correctly
in Eq. (6).

####  A.2.1 Dispersion relation and eigenvectors

We assume that $\lambda$, $\kappa_{\mathrm{P}}$, and $\kappa_{\mathrm{R}}$ are
constant, and that all quantities $q\in\\{\rho,u,P,E_{\mathrm{rad}}\\}$ can be
expressed in the form $q=q_{0}+q_{1}$, with $q_{0}$ being a constant
background state and $q_{1}=Aq^{\prime}e^{ikx-\omega t}$ a wave-like
perturbation with $A\ll 1$. By plugging them into the set of Eqs. (1), (5) &
(6) and discarding terms of $\mathcal{O}(A^{2})$ or higher we obtain

$\displaystyle\frac{\partial{\rho_{1}}}{\partial{t}}+\rho_{0}\frac{\partial{u_{1}}}{\partial{x}}=0,$
(22a)
$\displaystyle\frac{\partial{u_{1}}}{\partial{t}}=-\frac{1}{\rho_{0}}\frac{\partial{P_{1}}}{\partial{x}},$
(22b) $\displaystyle\frac{\partial{P_{1}}}{\partial{t}}=-\gamma
P_{0}\frac{\partial{u_{1}}}{\partial{x}}+(\gamma-1)Q_{\mathrm{rad}}^{\prime}$
(22c)
$\displaystyle\frac{\partial{{E_{\mathrm{rad}}}_{1}}}{\partial{t}}-\frac{\lambda
c}{\kappa_{\mathrm{R}}\rho_{0}}\frac{\partial^{2}{{E_{\mathrm{rad}}}_{1}}}{\partial{x}^{2}}=-Q_{\mathrm{rad}}^{\prime},$
(22d)

where $Q_{\mathrm{rad}}^{\prime}$ is the linearized radiation source term
after expanding $T^{4}\propto P^{4}/\rho^{4}$ to first order and using that
${E_{\mathrm{rad}}}_{0}=a_{\mathrm{R}}T_{0}^{4}$:

$Q_{\mathrm{rad}}^{\prime}=-\kappa_{\mathrm{P}}\rho_{0}c\left(4\frac{P_{1}}{P_{0}}-4\frac{\rho_{1}}{\rho_{0}}-\frac{{E_{\mathrm{rad}}}_{1}}{{E_{\mathrm{rad}}}_{0}}\right).$
(23)

By combining Eqs. (22a)–(23) and seeking a solution for $\omega$ we obtain the
dispersion relation

$\begin{split}\omega^{4}&\left[\chi\omega_{\mathrm{R}}\right]+\omega^{3}\left[-\chi\omega_{\mathrm{c}}^{2}-(4+\chi)\omega_{\mathrm{P}}\omega_{\mathrm{R}}\right]+\omega^{2}\left[\chi\omega_{\mathrm{P}}\omega_{\mathrm{c}}^{2}+4\omega_{\mathrm{P}}^{2}\right]\\\
+\omega&\left[-\gamma\chi\omega_{\mathrm{c}}^{2}\omega_{\mathrm{s}}^{2}-(4+\gamma\chi)\omega_{\mathrm{P}}\omega_{\mathrm{R}}\omega_{\mathrm{s}}^{2}\right]+4\omega_{\mathrm{P}}\omega_{\mathrm{c}}^{2}\omega_{\mathrm{s}}^{2}=0,\end{split}$
(24)

where we have defined for convenience
$\omega_{\mathrm{s}}=k{c_{\mathrm{s}}}_{0}$,
$\omega_{\mathrm{R}}=\kappa_{\mathrm{R}}\rho_{0}c$,
$\omega_{\mathrm{P}}=\kappa_{\mathrm{P}}\rho_{0}c$,
$\omega_{\mathrm{c}}=\sqrt{\lambda}kc$, and
$\chi=e_{0}/{E_{\mathrm{rad}}}_{0}$. This equation can be solved numerically
to obtain the oscillation frequency and damping rate of radiatively damped
sound waves, provided that $\text{Im}(\omega)\neq 0$. The components of the
related eigenvectors are then

$\begin{split}\frac{u^{\prime}}{{c_{\mathrm{s}}}_{0}}&=-\frac{\omega}{\omega_{\mathrm{s}}}\frac{\rho^{\prime}}{\rho_{0}},\\\
\frac{P^{\prime}}{P_{0}}&=-\frac{\omega^{2}}{\omega_{\mathrm{s}}^{2}}\frac{\rho^{\prime}}{\rho_{0}},\\\
\frac{E_{\mathrm{rad}}^{\prime}}{{E_{\mathrm{rad}}}_{0}}&=-\left(\gamma\chi\frac{\omega}{\omega_{\mathrm{P}}}-4\frac{\omega^{2}}{\omega_{\mathrm{s}}^{2}}-4+\chi\frac{\omega^{3}}{\omega_{\mathrm{P}}\omega_{\mathrm{s}}^{2}}\right)\frac{\rho^{\prime}}{\rho_{0}}.\\\
\end{split}$ (25)

####  A.2.2 Numerical solution

In our numerical setup we consider a periodic Cartesian domain $x\in[0,L]$
with $L=1$ au and 1000 uniformly spaced grid cells, and define $k=2\pi/L$. We
set $\mu=2.353$, $\gamma=7/5$, $\lambda=1/3$, and $T_{0}=30$ K. The opacity
$\kappa_{\mathrm{R}}=\kappa_{\mathrm{P}}=\kappa$ is set to a constant value
within each setup and varies between $0.003$–300
$\text{cm}^{2}/\text{g}_{\text{gas}}$ throughout the set of models we run. The
background state is given by

$\begin{split}\rho_{0}&=\\{10^{-13},10^{-12}\\}~{}\text{g}/\text{cm}^{3},&\qquad&u_{0}=0,\\\
P_{0}&=\frac{\mathcal{R}}{\mu}\rho_{0}T_{0},&\qquad&{E_{\mathrm{rad}}}_{0}=a_{\mathrm{R}}T_{0}^{4},\end{split}$
(26)

and the perturbed quantities are initialized as

$\begin{split}q_{1}=A\left[\text{Re}\left(q^{\prime}\right)\cos(kx)-\text{Im}\left(q^{\prime}\right)\sin(kx)\right],\end{split}$
(27)

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 $\omega_{\text{osc}}$ and damping rate $\omega_{\text{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(\omega_{\text{osc}}t+c_{1})e^{\omega_{\text{damp}}t},$
(28)

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. 8 for all quantities.

We then present our results in Fig. 9, showing that we recover the expected
damping and oscillation frequencies to very good degree for all values of
$\kappa$ and both reference values of $\rho_{0}$. We therefore conclude that
our FLD module is suitable for use in our models.

Figure 8: An example of fitted oscillations for our model with
$\rho_{0}=10^{-12}$ g/cm3 and $\kappa=1$ cm2/g at $x=0$. All quantities have
been normalized and shifted for clarity. Gray dots denote the analytical
solutions. Figure 9: Results from the linear perturbation test, showing the
damping $\omega_{\mathrm{damp}}$ (top) and oscillation frequency
$\omega_{\mathrm{osc}}$ (bottom) as a function of the opacity $\kappa$ for two
reference values of the background density $\rho_{0}$. We normalize $\omega$
to $\omega_{\text{iso}}=k{c_{\mathrm{s,iso}}}_{0}$. The oscillation frequency
fluctuates between the isothermal ($\omega_{\text{iso}}$) and adiabatic
expectation ($\sqrt{\gamma}\omega_{\text{iso}}$, dashed black line), depending
on the coupling between the thermal and radiation energy fields.