Migration of low-mass planets in inviscid disks: the effect of radiation transport on the dynamical corotation torque

Low-mass planets migrate in the type-I regime. In the inviscid limit, the contrast between the vortensity trapped inside the planet's corotating region and the background disk vortensity leads to a dynamical corotation torque, which is thought to slow down inward migration. We investigate the effect of radiative cooling on low-mass planet migration using inviscid 2D hydrodynamical simulations. We find that cooling induces a baroclinic forcing on material U-turning near the planet, resulting in vortensity growth in the corotating region, which in turn weakens the dynamical corotation torque and leads to 2-3x faster inward migration. This mechanism is most efficient when cooling acts on a timescale similar to the U-turn time of material inside the corotating region, but is nonetheless relevant for a substantial radial range in a typical disk (5-50 au). As the planet migrates inwards, the contrast between the vortensity inside and outside the corotating region increases and partially regulates the effect of baroclinic forcing. As a secondary effect, we show that radiative damping can further weaken the vortensity barrier created by the planet's spiral shocks, supporting inward migration. Finally, we highlight that a self-consistent treatment of radiative diffusion as opposed to local cooling is critical in order to avoid overestimating the vortensity growth and the resulting migration rate.


INTRODUCTION
More than 5500 exoplanets have been discovered so far 1 , and their population shows remarkable diversity in orbital parameters, masses, and multiplicity.The idea that these planets formed in circumstellar disks is supported by both the direct observation of planets embedded in the disk around PDS 70 with VLT (Keppler et al. 2018;Haffert et al. 2019), as well as kinematic signatures in observations of gas emission with ALMA (e.g., Teague et al. 2018).Understanding how planets formed and evolved to their current state is a key goal of modern astrophysics.
The final orbital configuration of a planet is directly affected by its interaction with the disk it is embedded in.The disk exerts a torque on the planet, which causes it to migrate through the disk (Goldreich & Tremaine 1979;Ward 1997a;Tanaka et al. 2002).The migration rate depends on the disk's properties as well as on the planet's mass, with low-mass planets migrating in the rapid type-I regime (Ward 1997a).For a recent review, see Paardekooper et al. (2022).
Modeling planet-disk interaction over long timescales is not an easy task.The population synthesis approach (e.g., Mordasini et al. 2009a,b) and N-body modeling (e.g., Coleman & Nelson 2014;Izidoro et al. 2017) are computationally inexpensive and cover a wide range of parameters in a relatively short time, but rely on analytical prescriptions for the torques acting on the embedded ⋆ E-mail: a.ziampras@qmul.ac.uk 1 https://exoplanetarchive.ipac.caltech.edu/planet(s).Hydrodynamical simulations, on the other hand, are significantly more computationally restrictive, but allow a full treatment of planet-disk interaction that yields more accurate results.
Hydrodynamical modeling is favorable for massive-enough planets, where planet-disk interaction becomes nonlinear (Rafikov 2002) and gap opening can occur (Crida et al. 2006).It is also crucial for low-viscosity disks, where the corotation torque can desaturate (Paardekooper et al. 2011) and dynamical corotation torques can come into play (McNally et al. 2017), as well as for non-isothermal disks, where thermal effects can even reverse the direction of migration (Kley & Crida 2008;Pierens 2015).In light of the recent paradigm shift towards low-turbulence disks, where accretion is thought to be powered by MHD-driven winds (Bai & Stone 2013), and the constraints on the radiative properties of protoplanetary disks through their dust distribution (Birnstiel et al. 2018), the need for understanding how radiative effects can affect planet migration in the nearly inviscid limit is more relevant than ever.
A simplified thermodynamics model that is often employed in simulations of planet-disk interactions is the locally isothermal equation of state, which assumes instant cooling of the gas towards an equilibrium temperature.However, recent models (Miranda & Rafikov 2019, 2020a;Ziampras et al. 2020), based on observations (e.g., Andrews et al. 2018;Öberg et al. 2021), highlight the importance of a more realistic treatment of thermodynamics and therefore radiative effects, even in regions where cooling occurs on dynamical time scales.Further work has showcased the impact of thermal diffusion (Ziampras et al. 2023;Miranda & Rafikov 2020b) in disk thermodynamics.It therefore becomes clear that an investigation of how radiation transport and in particular thermal diffusion influence planet migration in the marginally optically thin 20-40 au range is necessary.
In this study, we carry out inviscid radiation hydrodynamics simulations of planet-disk interaction, with the focus being the effect of radiation transport in the type-I regime of planet migration.Our goal is to investigate the source of torques related to radiative effects using a realistic prescription of radiation transport in 2D, extending works that have focused on viscous models (e.g., Kley & Crida 2008;Pierens 2015) but also providing an explanation on the origin of these "radiative" torques and constraining where in the disk they might be important.
We describe our physical and numerical setup in Sect. 2. We present our results in Sect.3, and discuss them in Sect. 4. We then summarize our findings in Sect. 5.

PHYSICS AND NUMERICS
In this section we lay out our physical and numerical framework.We list the vertically integrated hydrodynamical equations, introduce the sources of cooling used in our models, and describe our numerical setup.

Physics
We consider a disk of ideal gas with adiabatic index γ = 7/5 and mean molecular weight µ = 2.35 around a star with mass M ⋆ and luminosity L ⋆ .The inviscid Navier-Stokes equations then read where Σ, u and P denote the surface density, velocity vector, and pressure of the gas, and the internal energy density is given by the ideal gas law as e = P/(γ − 1).The stellar potential at distance R is Φ ⋆ = −GM ⋆ /R, where G is the gravitational constant.We can also define the gas isothermal sound speed c s = √ P/Σ = RT/µ, where T is the gas temperature and R is the ideal gas constant.The pressure scale height is then H = c s /Ω K , where Ω K = GM ⋆ /R 3 is the Keplerian angular velocity, and the disk aspect ratio is h = H/R.
The planetary potential follows a Plummer-like prescription with a smoothing length ϵ = 0.6H p similar to Müller & Kley (2012) to account for the vertical disk stratification: Finally, Q contains any additional radiative terms.In our models the disk is heated by stellar irradiation following the passive, irradiated disk model of Menou & Goodman (2004) where θ is the flaring angle with θ = 2h/7 (as T ∝ R −3/7 , see Chiang & Goldreich 1997), ε is the disk albedo (here ε = 1/2), and τ eff is an effective optical depth following Hubeny (1990) Here, τ R,P = 1 2 κ R,P Σ is the optical depth due to the Rosseland and Planck mean opacities κ R and κ P .We assume κ R = κ P = κ, following the opacity model of Lin & Papaloizou (1985).
We then treat thermal cooling through the disk surfaces as and in-plane radiation transport with a flux-limited diffusion approach (FLD, Levermore & Pomraning 1981) with the flux limiter λ following Kley (1989).The balance between Q irr , Q cool , and Q rad results in a temperature profile T ∝ R −3/7 .

Vortensity growth due to cooling
In this study we often use the vertical component of the gas vortensity, defined in 2D as ϖ = ∇×u/Σ • ẑ, to interpret our results.Taking the curl of Eq. ( 1b) and dividing by Σ yields the vortensity equation, which, for a 2D flow with u z = 0, reads where S ∝ ∇Σ × ∇P is a vortensity source term due to baroclinic forcing.For a barotropic flow (P = P(Σ)), ϖ is conserved along streamlines.
Let us now consider a gas parcel in the planet's corotating region.In the planet's corotating frame this parcel follows a streamline with a horseshoe pattern (see e.g., Kley & Nelson 2012), performing a U-turn once behind and once ahead of the planet before completing a closed loop.We can now analyze what happens to the vortensity of the gas parcel during the U-turn at R = R p , along the radial direction (see schematic in Fig. 1).

Locally isothermal limit
In the limit where cooling happens on timescales much shorter than the orbital timescale Ω −1 K , we can assume that T (R) ∝ R q is set by a balance between the radiative terms listed in Sect.2.1, and U-turning material is expected to drive vortensity growth for q 0 (Casoli & Masset 2009).The vortensity of gas U-turning behind the planet (ϖ b ) will change according to where a subscript 'p' denotes the value at the planet's radial location.Assuming a dense, circular envelope around the planet such that ∂ ϕ Σ| b > 0 and q < 0, this implies that S iso b > 0. In other words, ϖ will increase as material U-turns behind the planet.In a similar way, the vortensity of gas U-turning ahead of the planet (ϖ a ) will change as resulting in a decrease in ϖ since ∂ ϕ Σ| a < 0. Finally, assuming that the U-turns are symmetric such that ∂ ϕ Σ| b = −∂ ϕ Σ| a and x a = x b for x ∈ [P, Σ, T ], we have that S iso a + S iso b = 0.In other words, while ϖ will increase (decrease) as material Uturns behind (ahead of) the planet, the total change along a closed horseshoe loop is zero and ϖ is ultimately conserved.
1.0 0.5 0.0 0.5 1.0 Gas streamlines around an embedded planet (marked with a black dot) in the {R, ϕ, z = 0} plane.Material just outside the planet's corotating region (CR) feels a kick as it shears past the planet (green curves).Inside the CR, material follows closed orbits as it U-turns behind (orange) and ahead (blue) of the planet.Purple orbits mark the L4 and L5 Lagrange points.The bottom panel is a zoom-in of the top between the dashed lines, showing that material U-turning closer to the planet librates closer to the edge of the CR.

Adiabatic limit
When cooling is inefficient, the parcel contracts (expands) adiabatically as it approaches (recedes from) the planet.This process is reversible and entropy is conserved along the streamline.Each streamline therefore retains its own entropy.Defining the entropy function K = P/Σ γ , with Σ ∝ R s , we initially have that K(R) ∝ R ξ , with ξ = q + (1 − γ)s.Assuming, without loss of generality, that ξ < 0, a high-entropy streamline at R < R p U-turning behind the planet will advect into a low-entropy radial zone at R > R p , resulting in an entropy discontinuity near R ≈ R p and a nonzero S b .The lowentropy streamline will, of course, also U-turn ahead of the planet into the high-entropy zone at R < R p , resulting in S a = −S b and a net S adb = 0. Furthermore, after a few libration timescales (see Eq. ( 14)), entropy will be completely mixed in the horseshoe region and therefore In our models, ξ ≈ −0.03 and therefore vortensity growth in adiabatic models is negligible.We do, however, explore different values of ξ in Appendix B.

Finite cooling timescale
When cooling acts on a finite timescale, this calculation becomes complicated.We can however show that the two source terms ahead of and behind the planet (S a , S b ) are asymmetric such that vortensity always increases along a pair of streamlines U-turning in opposite directions.We provide a non-rigorous derivation in Appendix A as a proof of concept.

Dynamical corotation torque
In a sufficiently laminar disk, where mixing between the planet's horseshoe region and the background disk is negligible, a migrating planet will experience a dynamical corotation torque (McNally et al. 2017) where the subscript 'p' denotes quantities at the planet's radial position and x h is an estimate of the half-width of the corotating region in the 2D approximation following Paardekooper et al. ( 2010) In Eq. ( 11), ϖ(R p ) denotes the vortensity of the background disk at R p and ϖ h refers to the characteristic vortensity enclosed in the planet's horseshoe region, which would correspond to the vortensity of the initial radial location of the planet (R 0 ) if ϖ is conserved inside the corotation region.
Provided that the background disk radial velocity u R is small compared to the planet's migration speed, this dynamical corotation torque is positive for a planet migrating inwards and becomes stronger as ϖ(R p )/ϖ h increases (assuming Σ ∝ R −3/2 or shallower), eventually stalling the planet (Paardekooper 2014).Based on the discussion in Sects.2.2.2-2.2.3, a vortensity source term in the horseshoe region would cause ϖ h to evolve as ϖ h (t) ≈ ϖ(R 0 ) + t 0 S rad (t ′ ) dt ′ , weakening the dynamical corotation torque and speeding up inward migration, possibly preventing the planet from stalling.
Our initial conditions assume power-law profiles for Σ and T with with Σ ref = 56.7 g/cm 2 and T ref = 21 K, or equivalently h 0 = 0.05, obtained by a balance between Q irr and Q cool (see Eqs. ( 3) & ( 5)).With these parameters, our chosen resolution yields 25 cells per scale height at R = R 0 .The disk surface density corresponds to 1700 g/cm 2 at 1 au.The chosen set of parameters also results in a cooling timescale (Ziampras et al. 2023, see also Eqs. ( 15) & ( 16)).
The velocity field is initialized with a Keplerian profile corrected for pressure support, and the planet is placed at R p = R 0 with M p = 2 × 10 −5 M ⋆ ≈ 6.7 M ⊕ .The disk boundaries are periodic in azimuth and closed in the radial direction, and we apply the wave-damping prescription of de Val-Borro et al. (2006) in the radial zones R < 0.525 and R > 1.525 with a damping timescale of 0.1 orbits at the respective boundary.
The planet is integrated using an adaptation of the N-body module by Thun & Kley (2018), where the correction due to the neglect of disk self-gravity by Baruteau & Masset (2008b) is applied to models where the planet is allowed to migrate.The indirect term by the planet and star orbiting their common center of mass is included in all models, but the star feels the gravitational influence of the disk only in models where the planet can migrate (Crida et al. 2022).The planet grows to its final mass over 10 orbits using the formula in de Val-Borro et al. (2006).
We run two sets of models, where the planet is either fixed or allowed to migrate through the disk.In each set we carry out the following models: These models are sometimes tagged "rad", "adb", and "iso" in our plots.For consistency, we always use orange, blue, and green respectively to refer to each model in all plots.

RESULTS
In this section we present the results of our numerical models.We first focus on the vortensity evolution in models with fixed planets, before moving on to models where the planet is allowed to migrate.We then investigate how different cooling regimes affect the vortensity evolution.

Fixed planet
We first compare the vortensity evolution in our models with a fixed planet.Even though there is initially a radial vortensity gradient through the disk (ϖ 0 (R) ≈ Ω K /2), the vortensity is mixed in the planet's corotating region.This process happens over the libration timescale After a few τ lib , ϖ is flat inside the corotating region and the corotation torque vanishes.In Fig. 2 we compare the vortensity evolution in our fiducial models with a fixed planet after 1 τ lib .In all models ϖ undergoes phase mixing and is roughly constant in the corotating region after a few libration timescales.In the radiative model, however, a vortensity excess is visible (in red) at R ≈ R p ± x h , indicating that vortensity is not conserved.We plot the azimuthally averaged vortensity in Fig. 3, where we show that the excess grows with time and affects the vortensity profile in the entire corotating region.
This excess, driven by cooling and unrelated to vortensity growth due to the Rossby Wave Instability (Lovelace et al. 1999), is due to the source term S in Eq. ( 7), which acts to increase ϖ as material U-turns near the planet (see Sect. 2.2).To verify this, we plot S for all models at t = τ lib in Fig. 4. Our findings agree with the expectations in Sect.2.2: S is zero in the adiabatic model, nonzero but antisymmetric about the planet in the locally isothermal model such that it averages to zero, and asymmetric in the radiative model such that S > 0 in the corotating region.As a result, ϖ in the corotating region only increases in models with cooling, and in particular on streamlines that U-turn very close to the planet.Since these streamlines librate at the edges of the corotating region, the excess piles up at R ≈ R p ± x h (see Figs. 2, 3).
Given that the dynamical corotation torque for an inwardly migrating planet is proportional to ϖ p /ϖ h − 1, where ϖ p is the background disk vortensity at the planet's location and ϖ h is the characteristic vortensity trapped in the corotating region (McNally et al. 2017), the vortensity growth we found in our radiative models should result in faster inward migration.We investigate this in Sect.3.2.

Migrating planet
We now repeat the models in Sect.3.1, allowing the planet to migrate after 10 orbits.The migration tracks for these models are shown in Fig. 5.
Initially, the planet migrates inwards faster in the radiative model than in the adiabatic.This can be attributed to a combination of a weaker dynamical corotation torque due to the vortensity growth mechanism we described in Sects.2.2 & 3.1, as well as an effectively smaller γ due to cooling (Paardekooper et al. 2010).
After 200 orbits, migration in the radiative model slows down to match the adiabatic model.This becomes clear when comparing the migration timescales τ mig = a/ȧ, which become equal for the radiative and adiabatic models for t ∼ 200-600 orbits.This can be interpreted as the result of a combination of effects.For one, the dynamical corotation torque scales with the planet's migration speed, see Eq. ( 11), therefore partially slowing the planet down.In addition, the vortensity excess via baroclinic forcing becomes less important as the planet migrates inwards, because the contrast to the local vortensity (ϖ p /ϖ h ) increases faster than the vortensity input via baroclinic forcing in the corotating region when the planet migrates quickly.
To support this argument we show that vortensity growth is not as efficient in the case of a migrating planet.In part, this happens because the mechanism relies on material U-turning near the planet repeatedly.As the planet migrates inwards, however, material Uturning behind the planet is no longer trapped in the corotating region.Instead, the latter is now a tadpole-shaped zone in front of the planet (Masset & Papaloizou 2003;Papaloizou et al. 2007), and the vortensity-rich streamlines behind the planet are simply left behind the planet's wake, unable to sustain the feedback loop that led to substantial vortensity growth in the fixed case in Sect.3.1.We illustrate this picture in Fig. 6.
However, after t ∼ 600 orbits, the planet speeds up in the radiative model, migrating roughly 3× and 2× faster compared to the adiabatic and locally isothermal models, respectively.We interpret this as a combination of two effects: a rapidly migrating planet will not capture as many U-turning streamlines, weakening vortensity growth due to cooling.This results in the planet slowing down over time as the vortensity contrast to the background disk increases and the dynamical corotation torque becomes more positive (Γ h ∝ dR p /dt), which is why the radiative model briefly matches the adiabatic migration timescale at t ∼ 200 orbits.By slowing down, however, the planet can capture more vortensity-rich U-turning streamlines, enhancing ϖ h and weakening the dynamical corotation torque, causing the planet to speed up once again.The two effects seemingly reach a balance after ∼ 1000 orbits, with the planet migrating inwards without slowing down.We nevertheless stress that this is only a qualitative explanation, and as such we cannot rule out the possibility of runaway migration in the radiative model (e.g., Pierens 2015).14)), showing that ϖ is being mixed throughout in the CR but an excess develops in the radiative run.The rightmost panel shows that the excess grows with time in the radiative model.Vertical dashed lines mark the edges of the CR at ±x h (see Eq. ( 12)).

Different cooling timescales
In the previous sections we showed how cooling can speed up inward migration.We also showed that in the limit of both inefficient and rapid cooling (adiabatic and isothermal, respectively), this effect is absent.We therefore expect that the effect depends on the cooling timescale β cool = t cool Ω K , and is maximized for a critical value.
To measure this dependency of vortensity growth on β cool , we run a set of radiative models where we vary the reference surface density Σ ref to control β cool .For these models, the planet is kept fixed at R 0 .We compute the cooling timescale following the prescription in Miranda & Rafikov (2020b) where β surf and β mid are the cooling timescales for surface and in-   7), which leads to vortensity growth.This term is zero throughout the corotating region in the adiabatic model, averages to zero in the locally isothermal model, but is asymmetric in the radiative model such that vortensity grows over time.Black curves mark the position of the planetary spirals following Rafikov (2002).
plane cooling and are given by (Ziampras et al. 2023) Figure 5. Migration tracks for the models with migrating planets.The radiative model (orange) initially migrates at the locally isothermal rate (green) for t ≲ 200 P 0 , then matches the adiabatic rate (blue) until t ∼ 600 P 0 , and afterwards diverges significantly as the planet speeds up.This is also evident in the migration timescales τ mig (bottom), which match between the radiative and adiabatic models for t ∼ 200-600 P 0 , but then show that the planet migrates roughly 2-3× as fast in the radiative model compared to the adiabatic.The dashed curve marks the migration rate without corotation torques.The purple curve (discussed in detail in Sect.3.4) shows a model where thermal diffusion was omitted, resulting in unphysically rapid migration at the rate given by Lindblad torques only.Given that we use a wave damping zone for R < 0.53 R 0 , we exclude data for a p < 0.55 R 0 .The damping zone interface is also the reason why migration slows down in all models for a p ≲ 0.6 R 0 .
We then compute the quantity S similar to Fig. 4 after 20 planetary orbits and integrate it over the corotating region for each model to obtain a proxy for the vortensity growth rate.We plot the results in Fig. 7 (orange curve), where we show that this integral peaks at β cool ≈ 4.3 and approaches zero in the adiabatic and locally isothermal limits.
The peak agrees very well with the estimate of the U-turn timescale β U-turn ≈ hτ lib Ω K ≈ 3.65 in Baruteau & Masset (2008a), highlighted with a vertical line in Fig. 7.This suggests that the vortensity growth rate is maximized when the cooling timescale is comparable to the U-turn timescale of material in the corotating region.In other words, the effect of cooling is strongest when it happens over a timescale similar to the time over which U-turning material interacts dynamically with the planet.

The effect of thermal diffusion
Our radiative models discussed above include a treatment of thermal diffusion along the disk plane through the flux-limited diffusion (FLD) approximation.Thermal diffusion serves to smooth out temperature gradients in the disk, and therefore could affect the vortensity growth rate.
A more simplistic (but less accurate) approach would be to implement radiation transport as a local cooling term in the energy equation, similar to the method in Miranda & Rafikov (2020b) (see    7), integrated over the corotating region, as a function of the cooling timescale β cool or the corresponding distance R for our models.The function peaks where β cool ≈ β U-turn (vertical line), or R ∼ 10-35 au for our setup.The orange curve shows the results for our radiative models, while the purple curve shows models without thermal diffusion (discussed in Sect.3.4).Ziampras et al. 2023, for a comparison).This approach, however, misses the diffusive component of in-plane radiation transport.To investigate the difference, we run a set of models where we replace Q irr + Q cool + Q rad in Eq. ( 6) with a local cooling term Q relax given by (e.g., Gammie 2001) where c v is the heat capacity at constant volume and β = β cool is the cooling timescale for both surface and in-plane cooling (see Eqs. ( 15) & ( 16)).The factor of 4 is a correction following Ziampras et al. ( 2023), obtained by linearizing Q cool in Eq. ( 5).We note that Dullemond et al. (2022) cite a correction of 4 + b with b = dlog κ dlog T ≈ 2 instead, but this applies only in the optically thin limit where τ eff ∝ τ −1 ∝ κ −1 (see Eq. ( 4)).In the optically thick limit one would instead use 4 − b since τ eff ∝ τ, but since we have τ ∼ 2 ⇒ τ eff ∼ 1 at R = R 0 in our models we ignore b.
In Fig. 5 we also plot the migration tracks for a run with local, β cooling.We find that the planet migrates much faster when thermal diffusion is omitted, leading to drastically faster migration.In addition, we compute the integrated baroclinic term as a function of β cool similar to Sect.3.3 and include it in Fig. 7 (purple curve).While the function also peaks at β cool ≈ β U-turn , it consistently overestimates the vortensity growth rate compared to our radiative models where thermal diffusion is included.This suggests that, while cooling results in vortensity growth due to baroclinic forcing and therefore faster inward migration, thermal diffusion acts to partially suppress this effect.

Radiative damping of spiral shocks
While the focus of this study has been how radiative cooling can induce vortensity growth via baroclinic effects, it is worth highlighting that an additional source of vortensity exists in the planet's vicinity.The planet's spiral arms will steepen into shocks as they propagate through the disk (Rafikov 2002), resulting in a vortensity excess (Lin & Papaloizou 2010).This creates a vortensity pileup about ±x sh away from the planet, where x sh is the shock distance following Zhu et al. (2015) (see also Goodman & Rafikov 2001) x sh ≈ 0.93 γ + 1 12/5 Here, M th is the thermal mass at the location of the planet (Goodman & Rafikov 2001).However, radiative cooling is expected to weaken the spiral shocks (Miranda & Rafikov 2020a;Zhang & Zhu 2020;Ziampras et al. 2023), possibly resulting in weaker vortensity growth.This would imply that, while a planet in an adiabatic or isothermal model would also need to overcome the dynamical corotation torque due to the vortensity "barrier" created by its own spiral shocks, this effect is mitigated or even absent in a radiative model.The result would be faster inward migration when β ∼ 1, where the damping effect of radiative cooling on spiral shocks is maximized (Miranda & Rafikov 2020a;Ziampras et al. 2023).This condition is met in our models, as we have β ∼ 1.5 at R = 30 au (see Fig. 7).
Figure 8 shows a comparison of the azimuthally averaged perturbed vortensity for our fiducial models similar to Fig. 3, but covering a wider radial range and normalized to the Keplerian vortensity to highlight any excess.Here, it becomes clear that radiative cooling interferes with the damping of spiral shocks, resulting in a vortensity excess that is weaker and closer to the planet.Finally, in Fig. 9 0.3 0.2 0.1 0.0 0.1 0.  we show a heatmap of the perturbed vortensity for models with migrating planets, where we find that the excess due to spiral shocks is practically absent in models with radiative cooling.
We expect that this effect is of secondary importance compared to the baroclinic forcing we discuss in this paper.Nevertheless, we stress that this can further enhance the effect of radiative cooling on inward migration such that baroclinic forcing in the horseshoe region (β ∼ β U-turn ∼ 10) and radiative damping of shocks (β ∼ 1) can operate efficiently together for β ∼ 1-10.Such an investigation is beyond the scope of this paper, but could be the subject of future work.

DISCUSSION
In this section we discuss our findings in the context of previous studies and more realistic physics.We also highlight the relevance of our results in planet migration.Lega et al. (2014) showed that radiative cooling results in the "cold finger" effect in the corotating region, as material is colder and denser after a U-turn compared to an adiabatic case.This would formally result in vortensity generation due to the effect discussed in this paper, but only if the cooling timescale is right (see Fig. 7).Pierens (2015) also investigated the effect of radiation transport on planet migration, but did not report on the mechanism described here.Yun et al. (2022) similarly carried out 3D radiative hydrodynamic simulations of planet migration, showcasing the "cold finger" effect but not finding evidence of cooling-induced vortensity growth.This is likely due to the fact that all of these studies probed the inner 1-5 au, where the cooling timescale is β cool ∼ 10 3 -10 4 ≫ β U-turn , such that the effect we describe was too inefficient (see Fig. 7).For a solar-type star, we expect cooling-induced vortensity growth to be most efficient at R ∼ 10-35 au (see Fig. 7).

More realistic physics
In this paper we aimed to isolate the effect of radiative cooling on the migration rate of low-mass planets.In reality, however, many more processes that we did not include and which can directly or indirectly affect migration could come into play.
For example, accretion heating can give rise to thermal torques (Benítez-Llambay et al. 2015;Masset 2017, who also included thermal diffusion) or even result in vortex formation (Cummins et al. 2022) that could affect migration (Lega et al. 2021).The Hall effect in non-ideal MHD can lead to a torque in the disk midplane, which strongly affects migration (McNally et al. 2017(McNally et al. , 2018)).In addition, the torque by an MHD-driven wind could reverse the direction of migration (Kimmig et al. 2020).Finally, including the vertical direction can give rise to buoyancy-related torques (Zhu et al. 2012;McNally et al. 2020) as well as introducing a vertical dependence of the MHD wind-driven torques (Wafflard-Fernandez & Lesur 2023).
We expect that a more realistic picture should include all such effects, but also note that a treatment of cooling-and in particular radiative diffusion-would likely regulate their contribution to the total torque.This could be an interesting topic for future work.

Sustaining, retriggering, or exaggerating the effect
While the dynamical corotation torque will not fully halt inward migration, it can slow the planet down enough for it to open a partial gap and slowly transition into the much slower type-II migration (Ward 1997b).Alternatively, migration can halt as the planet approaches the inner edge of its gap and the (negative) outer Lindblad torque is suppressed (though see McNally et al. 2019, for a discussion on the effect of vortices).
In both cases the corotating region will span the full disk azimuth, trapping trailing vortensity-rich streamlines and retriggering rapid inward migration if the planet can break the stall (e.g., McNally et al. 2019).This could lead to a sustained faster inward migration, or a series of migration stalls and restarts.Investigating this scenario could be the subject of future work.
In the case of a more traditionally viscous disk, one or more zerotorque regions ("migration traps") could exist in the disk, such as near icelines where the background disk density and temperature profiles change (e.g., Coleman & Nelson 2016).In such a scenario, the (static, outward) corotation torque remains unsaturated and balances against the (inward) Lindblad torque, allowing the planet to remain stationary.Here, the cooling-induced vortensity growth we discuss could result in the planet migrating either inwards or outwards (depending on the accretion velocity of the disk, see Eq. ( 11)) at a rate that will strongly depend on the disk properties (Pierens 2015).
Our analysis in Sect.3.4 further showed that thermal diffusion plays a critical role in the vortensity generation and subsequent torque by smoothing temperature gradients.It is therefore important to treat radiation transport appropriately (using e.g., FLD instead of a local cooling prescription for the in-plane cooling component) to avoid exaggerating the related torque.
In the more realistic scenario of a planet migrating through the disk as it grows, we expect that the effect discussed in this work will be less important during the earlier stages of migration, when the planet is still small and the corotating region is not yet fully developed.However, as the planet grows, the corotating region widens, and the cooling-induced vortensity growth becomes stronger, its impact on migration will quickly become relevant.Finding a critical mass beyond which this happens is beyond the scope of our work, however.
Finally, as we have discussed, it becomes apparent that holding the planet in a fixed radial location would both hide the contribution of this effect to the total torque, as it mainly affects the dynamical corotation torque, but also artificially speed up inward migration once the planet is released due to the pileup of vortensity in the corotating region.For this reason, we recommend that numerical models of planet migration allow the planet to migrate as soon as possible after the start of the simulation.

Connecting to population synthesis models
As discussed in Sect. 1, we expect that dynamical corotation torques will be important for low-to moderate-mass planets in lowturbulence disks.In order to investigate its effects on planet populations, however, a recipe for the dynamical corotation torque is required.McNally et al. (2018) provided such a recipe for a planet migrating in a disk where the Hall effect drives a Lorentz acceleration on the disk, resulting in a background radial velocity u R and a vortensity source term on the planet given by dϖ The vortensity source term S and the associated torque that we discuss in this paper is not related to magnetic fields as in the work of McNally et al. (2018), but it would be useful to parametrize S as a function of planet and disk properties such that its effects can be included in population synthesis models.This will be the focus of followup work.At the same time, it is important to note that our models focus on disks with a single planet.In a multiple-planet configuration, on top of planet-planet interactions, it is possible that the vortensityrich streamlines left behind by the innermost, inwardly-migrating planet can interfere with the migration of planets further out by increasing the local disk vortensity above its nominal Keplerian value.This could enhance the dynamical corotation torque on the remaining planets, slowing them down.

SUMMARY
We carried out numerical simulations of planet-disk interaction in inviscid disks, and investigated the effect of cooling on the vortensity evolution in the planet's corotating region.We summarize our findings below.
We found that radiative cooling interferes with the otherwise adiabatic compression and expansion of material in the planet's corotating region as it U-turns near the planet.The result is a net increase in vortensity in the corotating region, which weakens the stalling effect of the dynamical corotation torque (McNally et al. 2017) and substantially speeds up inward migration.We identified the reason for this vortensity growth as baroclinic forcing along streamlines Uturning very close to the planet, causing the vortensity excess to pile up at the edges of the corotating region.
We then compared models of migrating planets with and without radiative cooling, and found that this baroclinic forcing results in slightly faster inward migration for the first 200 orbits.As the planet continues to migrate inwards, a competition between the dynamical corotation torque becoming stronger for a rapidly-migrating planet and baroclinic forcing enhancing vortensity growth for a slowlymigrating planet prevents the planet from stalling and instead allows it to sustain a migration rate 2-3× faster than the isothermal or adiabatic rates.We also found that the vortensity growth rate (and therefore inward migration rate) is maximized when the cooling timescale is comparable to the U-turn timescale of material in the corotating region, making this mechanism relevant in a quite broad radial range in the disk (R ∼ 5-50 au).
Furthermore, we showed that thermal diffusion acts to suppress vortensity growth due to baroclinic forcing.We therefore conclude that the effect of radiative cooling on inward migration is twofold: it acts to speed up inward migration via vortensity growth, but this growth is partly suppressed by thermal diffusion.
Finally, we discussed how radiative damping of spiral shocks can further enhance the effect of cooling on inward migration by weakening the vortensity barrier created by the planet's own spiral shocks.While this effect is of secondary importance compared to the baroclinic forcing we discuss in this paper, we expect that it will be relevant in the radial range where β ∼ 0.1-10 (Miranda & Rafikov 2020a; Ziampras et al. 2023), and note that the two effects could operate efficiently together for β ∼ 1-10.
Our results highlight the importance of radiative effects in planetdisk interaction, and in particular the role of thermal diffusion.An important take-home message is that type-I migration in lowturbulence, radiative disks might be more challenging to understand than previously thought, and mechanisms that can naturally lead to outward migration (e.g., magnetic fields, planetary accretion luminosity) will become necessary from a modeling perspective in order to explain the population of low-mass planets at R ∼ 1-10 au.We also expect that the effect of thermal diffusion will be even more central when additional thermal effects are considered.In the case with cooling, we assume that thermal emission and in-plane cooling happen alongside adiabatic expansion.The heated gas parcel will then cool faster (and expand less) than it would have adiabatically.We then write that (A6) As we have defined τ ϕ and σ ϕ to be positive, we find that S cool > 0 during a U-turn both ahead of and behind the planet.In other words, vortensity will always increase along streamlines during a U-turn.We note, however, that this is a crude approach as, in principle, the assumption that σ adb ϕ = σ cool ϕ and τ adb ϕ = τ cool ϕ will not hold.We also did not take into account the presence of a background temperature and density gradient, which would further complicate this calculation.We can nevertheless expect that, in any case, the two terms in Eqs.(A5) & (A6) will sum to a positive value such that there is a net vortensity growth.

APPENDIX B: VORTENSITY GROWTH IN ADIABATIC MODELS WITH A RADIAL ENTROPY GRADIENT
In Sect.2.2.2, we argued that a radial entropy gradient K ∝ R ξ should result in a vortensity source term about the planet's azimuthal location.This term should however average to zero for each pair of streamlines U-turning ahead of and behind the planet, and only exist until the horseshoe is completely phase-mixed to a constant entropy state.
In Fig. B1 we show the vortensity source term similar to Fig. 4 at t = 20 orbits for three adiabatic models with different surface density radial profiles Σ ∝ R s , where s ∈ {−2, −1, 0}.We find that the source term is zero throughout the corotating region for our fiducial model with s = −1 and ξ ≈ 0, but nonzero for the models with ξ 0, where the source term is antisymmetric about the planet such that it averages to zero.This is consistent with our expectations in Sect.2.2.2.
We further run a set of models with different surface density exponents s to investigate the long-term evolution of S adb for different

Figure 2 .
Figure 2. Two-dimensional heatmaps showing the vortensity ϖ in the planet's corotating region (CR).The first three panels from left to right compare different models at t = τ lib (see Eq. (14)), showing that ϖ is being mixed throughout in the CR but an excess develops in the radiative run.The rightmost panel shows that the excess grows with time in the radiative model.Vertical dashed lines mark the edges of the CR at ±x h (see Eq. (12)).

Figure 3 .
Figure 3. Azimuthally averaged perturbed vortensity in the planet's corotating region (CR) for the models in Fig. 2. The vortensity excess in the radiative model grows with time and affects the entire CR.

Figure 4 .
Figure 4.The baroclinic forcing term S in Eq. (7), which leads to vortensity growth.This term is zero throughout the corotating region in the adiabatic model, averages to zero in the locally isothermal model, but is asymmetric in the radiative model such that vortensity grows over time.Black curves mark the position of the planetary spirals followingRafikov (2002).

Figure 6 .
Figure 6.Similar to Fig. 2, but for the models with migrating planets.The corotating region is now tadpole-shaped, with a vortensity excess around it (resembling a flame) in the radiative model.The vortensity-rich streamlines U-turning behind the planet are no longer confined to the corotating region, as the planet has drifted inwards by the time they would have reached the planet from the front.

Figure 7 .
Figure7.The baroclinic source term S in Eq. (7), integrated over the corotating region, as a function of the cooling timescale β cool or the corresponding distance R for our models.The function peaks where β cool ≈ β U-turn (vertical line), or R ∼ 10-35 au for our setup.The orange curve shows the results for our radiative models, while the purple curve shows models without thermal diffusion (discussed in Sect.3.4).

Figure 8 .
Figure 8. Azimuthally averaged perturbed vortensity for static planet models at t = 2 τ lib , similar to Fig. 3 but normalized to the Keplerian vortensity.Colored bands mark the location where spiral shocks deposit vortensity in the disk.The radiative model shows a weaker vortensity excess due to spiral shocks, closer to the planet.

Figure 9 .
Figure 9. Similar to Fig. 6, with a focus on the vortensity profile ahead of (inside) the planet's orbit.Models with radiative cooling (right) show practically no vortensity excess due to spiral shocks.
. 2.2, we now assume a hot and dense circular envelope around the planet such that ∂ϕ Σ| b = −∂ ϕ Σ| a = σ ϕ and ∂ ϕ T | b = −∂ ϕ T | a = τ ϕ .We also assume that the gas parcel in a disk with cooling contracts and heats up in the same way that it would in an adiabatic disk, except that the expansion and cooling phases are different.In other words, σ adb ϕ + |δT | σ ϕ .