© The Authors 2026

1 Introduction

In recent years, thermal torques on low-mass planets and planetary embryos embedded in regions of the gaseous protoplanetary disks where the flux is laminar (Lega et al. 2014; Masset 2017; Benítez-Llambay et al. 2015; Chametla & Masset 2021) have emerged as a possible solution to the problem of fast inward migration (see Paardekooper et al. 2023, for a review). Thermal torques on a planet are due to the asymmetric gas density structures that formed in the vicinity of the planet, which result from thermal diffusion in the disk. When a planet accretes pebbles, it releases heat into the environment, producing the formation of low-density, asymmetric hot lobes in the disk. These lobes, in turn, produce a torque on the planet (the heating torque; Masset 2017). If the planet does not accrete pebbles, thermal diffusion produces the formation of high-density cold lobes in the coorbital region of the planet (the cold thermal torque; Lega et al. 2014; Masset 2017). Both heating and cold thermal torques can exceed the Lindblad and corotation torques by up to an order of magnitude, such that the Lindblad and corotation torques become negligible. Heating thermal torques can slow down, or even reverse, the inward migration of growing low-mass planets (Benítez-Llambay et al. 2015).

Using global disk models, Guilera et al. (2019) showed that thermal torques can modify both migration maps and planet formation tracks, i.e., the evolution of the planetary mass as a function of semi-major axis (but see also Baumann & Bitsch 2020). In a more comprehensive framework of mass growth by pebble accretion, Guilera et al. (2021) found that thermal torques play a significant role in shaping planet formation tracks in low-viscosity disks, particularly when the disks are both massive and metal-rich. Under these conditions, planets forming beyond the ice line may undergo substantial outward migration.

Thermal forces can also excite eccentricity and inclination. In a model that additionally follows the evolution of eccentricity and inclination for a population of pebble-accreting embryos, Cornejo et al. (2023) report that the inclusion of thermal torques leads to a markedly different dependence of the final embryo mass spectrum on disk metallicity. The excitation of eccentricity by thermal torques may also have important implications for the trapping and formation of planets in pressure bumps (Chrenko & Chametla 2023; Pierens & Raymond 2024; Velasco-Romero et al. 2024).

All studies of thermal torques to date have been conducted under the assumption of an unmagnetized, laminar disk. However, there exist regions of protoplanetary disks in which turbulence is expected to dominate the transport of angular momentum. In such turbulent environments, it has been shown that the migration of planetary embryos can become stochastic (Nelson & Papaloizou 2004; Nelson 2005). This occurs, for instance, outside the dead zone, where strong magnetic coupling enables turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1991). In addition, hydrodynamical instabilities may also lead to turbulence even within the dead zone (see Lyra & Umurhan 2019). In this work, we assess whether thermal torques remain effective at driving the migration of low-mass planets or planetary embryos in turbulent regions of the disk. To address this, we conducted high-resolution magnetohydrodynamic (MHD) simulations of a planetary embryo embedded in a thermally diffusive, stratified protoplanetary disk threaded by a toroidal magnetic field.

The paper is laid out as follows. Section 2 describes the MHD equations, the numerical setup, and the code used in this work. In Section 3, we present our results. A brief discussion of possible improvements in more complex models is presented in Section 4. Finally, concluding remarks can be found in Section 5.

2 Physical and numerical setup

We consider a 3D magnetized, inviscid gas disk whose dynamical evolution is described by the MHD equations $${\partial }_t\rho +\nabla \cdot \left(\rho v\right)=0,$$(1) $${\partial }_t\left(\rho v\right)+\nabla \cdot \left(\rho v\otimes v+pI\right)=-\rho \nabla \Phi +J\times B,$$(2) $${\partial }_{t}e+\nabla \cdot \left(ev\right)=-p\nabla \cdot v-\nabla \cdot F_{\text{H}}+S_p,$$(3) $$\frac{\partial B}{\partial t}=\nabla \times \left(v\times B\right),$$(4)

where ρ and v denote the gas density and velocity, respectively, and J and B are the current density and magnetic field. Furthermore, I denotes the unit tensor and Φ is the gravitational potential. For the gravitational potential of the planetary embryo, we adopted a softening length comparable to our mesh resolution. For the pressure p, we consider the equation of state $$p=\left(\gamma -1\right)e,$$(5)

where e is the internal energy density and γ = 7/5 the gas adiabatic index.

In Eq. (3), Sp is a source termarising fromthe energy release of the planetary embryo defined as $$S_p=L\delta \left(r-r_p\right),$$(6)

with L and r_p being the luminosity and the location of the planet, respectively, and 6 being Dirac's delta function. While F_H is the heat flux, given by $$F_{\text{H}}=-\chi \rho \nabla \left(\frac{e}{\rho }\right),$$(7)

χ is the thermal diffusivity. This advection-diffusion approximation has been commonly adopted in previous studies (e.g., Masset 2017; Velasco Romero & Masset 2020; Chametla & Masset 2021; Chametla et al. 2025). Other authors instead employed the flux-limited diffusion scheme within a two-temperature framework (Benítez-Llambay et al. 2015). In the regime considered here, however, both approaches describe radiative transport in the diffusion limit, and therefore we do not expect significant differences in the resulting thermal torques.

2.1 Setup

We largely followed the setup of Chametla & Masset (2021). Only a brief summary is given here. We considered a 3D patch of a nonself-gravitating disk. The simulations were performed in a frame corotating with the planetary embryo, using a polar spherical grid centered on the star. The numerical domain extends for 0.9r_p to 1.1r_p in the radial direction, from −π/18 to π/18 in the azimuthal direction and between π/2 − h_p and π/2 + h_p in colatitude (here h_p = 0.05 is the aspect ratio at the position of the planetary embryo). The simulations employed (N_r, N_e, N,) = (880, 128, 1536) grid zones, corresponding to cell sizes of (2.3 × 10⁻⁴r_p, 7.8 × 10⁻⁴r_p, 2.3 × 10⁻⁴r_p).

Initially, the surface density and temperature followed constant radial profiles, i.e., Σ ∝ r^−σ and T ∝ r^−ζ, with σ = ζ = 0. The disk was initially threaded by a toroidal magnetic field given by $$B_{\varphi }=H\Omega \sqrt{\frac{8\pi \rho }{\beta }},$$(8)

whose strength was parameterized by the plasma-β parameter, with ß = 50 and 1000. Here, H denotes the pressure scale length.

We considered a planetary embryo of mass M_p, with a luminosity of L ∈ {0, L_c}, where L_c is the critical luminosity $$L_c=\frac{4\pi G{M}_{p\chi }{\rho }_0}{\gamma },$$(9)

with ${\rho }_0={\displaystyle {\sum }_0/\left(\sqrt{2\pi }H\right)}$ being the unperturbed midplane density at r = r_p (Masset 2017). For L > 0.94L_c, heat release dominates over thermal diffusion and the total thermal torque is positive, whereas for L < 0.94L_c, diffusion dominates and the total thermal torque is negative. At L = 0.94L_c, the thermal contribution vanishes (see Fig. 6 of Chametla & Masset 2021). We set thermal diffusivity to $\chi ={10}^{-5}{r}_p^2{\Omega }_p$.

Since planetary embryos are the seeds from which larger planets grow, a natural first step is to investigate whether the thermal torques that govern their migration in laminar disks (where they can slow down or even reverse migration and thereby favor their growth) remain effective in a turbulent medium. Therefore, unless otherwise stated, we adopted M_p = 0.33 M_♂, with M_♂ being the mass of Mars. In physical units, further adopting r_p = 5.2au and Σ₀ = 208 g cm⁻² yielded χ = 1.1 × 10²⁰ cm² s⁻¹ and L_c = 7.8 × 10²⁵ erg s⁻¹. The value of χ appears reasonable at this distance from the central star. In principle, the results of our simulations can be rescaled to a different orbital radius r_p. In that case, the thermal diffusivity scales as χ = 1.1 × 10²⁰(r/5.2au)^1/2cm² s⁻¹. At large radii, however, this simple scaling is expected to deviate from the true value because the radial dependencies of the opacity (Bell & Lin 1994), density, and temperature are unlikely to conspire to produce such a scaling of the diffusivity. Therefore, caution must be exercised when rescaling the simulation results to other orbital distances.

2.2 Code and boundary conditions

To solve the continuity, momentum, and energy equations for the gas, we used the hydrodynamic public FARGO3D multifluids' code (Benítez-Llambay et al. 2019), which is an extended version of the original FARGO3D code (Benítez-Llambay & Masset 2016), with orbital advection enabled (Masset 2000). A thermal diffusion module was included in the code FARGO3D¹ multi-fluids. This module was already tested and used in Chametla & Masset (2021) for the numerical study of thermal torques in a gaseous disk without magnetic fields. Using an operator-splitting scheme, the thermal diffusion term in the energy equation was evolved separately. Specifically, during the diffusion update, the internal energy satisfies $${\partial }_{t}e=-\nabla \cdot F_{\text{H}}.$$(10)

The source energy term S_p, representing the heat release of the planetary embryo, was implemented by increasing the internal energy density over eight adjacent grid cells (Benítez-Llambay et al. 2015; Chametla & Masset 2021; Velasco-Romero et al. 2024) $$\Delta e_{i,j,k,}=L\left(1-\frac{{\varphi }_{i,j,k}^p}{\Delta \varphi }\right)\left(1-\frac{r_{i,j,k}^p}{\Delta r}\right)\left(1-\frac{{\theta }_{i,j,k}^p}{\Delta \theta }\right)\frac{\Delta t}{V_{i,j,k}},$$(11)

where L is the luminosity, V_i,j,k is the volume of the cell (i, j, k), ${\varphi }_{i,j,k}^p=\left|{\varphi }_p-{\varphi }_{i,j,k}\right|$, $r_{i,j,k}^p=\left|r_p-r_{i,j,k}\right|$, ${\theta }_{i,j,k}^p=\left|{\theta }_p-{\theta }_{i,j,k}\right|$, and Δt is the difference between each timestep.

The boundary conditions follow those implemented for the gas density and velocity components in Chametla & Masset (2021). We included magnetic resistivity buffer zones in the radial and vertical directions of the disk, covering the same computational domain as the hydrodynamical damping zones.

2.3 Parameters and relevant scales within thermal torques' theory

Linear theory predicts that the size of the thermal disturbance is given by $${\lambda }_c=\sqrt{\frac{\chi }{\left(3/2\right){\Omega }_p\gamma }}$$(12)

(Masset 2017). The level of asymmetry of the thermal perturbations was determined by the planetary embryo's corotation offset, which is given by $$x_p=\eta h_p^2{r}_p,$$(13)

where η is a function of the power-law exponents of the surface density and temperature, given by $$\eta =\frac{\sigma }{3}+\frac{\zeta +3}{6}.$$(14)

In our case, for σ = ζ = 0, the parameter η reduces to 1/2.

3 Results

3.1 The turbulent disk

For the adopted disk parameters, Eqs. (12) and (13) yield x_p = 1.25 × 10⁻³r_p, λ_c = 2.2 × 10⁻³r_p, and H = 0.05r_p. The ordering x_p ≤ λ_c ≪ H implies that the density perturbation induced by the planetary embryo's luminosity should be observable, at least in the absence of turbulence or prior to its full development, which occurs after approximately three orbital periods. To verify this expectation, Fig. 1 presents the density perturbation for L = L_c and β = 50 after t = 1 orbit, i.e., before turbulence has had time to grow. At this stage, the hot lobes are already fully formed and exhibit a pronounced asymmetry with respect to the planetary embryo's position. This behavior is characteristic of the evolution of hot lobes in a turbulence-free gas disk (Chametla & Masset 2021).

The onset of turbulence is illustrated in Fig. 2, which shows the temporal evolution of the stress-to-pressure ratio α, a measure of the strength of the turbulence in the disk, computed as $$\alpha \left(r,t\right)=\frac{{\displaystyle \int \left(\rho \delta v_r\delta v_{\varphi }-B_r{B}_{\varphi }/4\pi \right)r^2\sin \theta d\theta d\varphi }}{{\displaystyle \int p{r}^2\sin \theta d\theta d\varphi }},$$(15)

for models with L = L_c and for the two values β = 50 and β = 1000. The development of MRI-driven turbulence begins at t ~ 1.8 orbits, and by t = 3 orbits α already exceeds 10⁻⁴. Note that when turbulence is fully developed, α reaches the values of 10⁻³ and 10⁻², when β = 1000 and β = 50, respectively.

To ensure that the MRI is well resolved in our models, we computed the quality factor Q_ϕ (Sorathia et al. 2012; Flock et al. 2017), $$Q_{\varphi }=2\pi \sqrt{\frac{16}{15}}\frac{\left|B_{\varphi }\right|}{\sqrt{4\pi \rho }}\frac{1}{\Omega }\frac{1}{r\Delta \varphi },$$(16)

for the model with a weak magnetic field (ß = 1000). Fig. 3 shows the radial profile of Q_ϕ (temporally and spatially {θ, ϕ} averaged). The fact that Q_ϕ > 8 indicates that the MRI is well resolved.

Because turbulence causes fluid elements in the vicinity of the planetary embryo (and more generally in the disk) to follow different trajectories in different realizations, a direct subtraction between simulations with L = L_c and L = 0 does not provide useful information. We therefore considered the density maps without subtracting the cold case. As shown in Fig. 4, independently of the magnetic field strength, turbulent fluctuations generate cells of high and low density (i.e., strong density gradients) in the same area that would otherwise be occupied by hot lobes (see Fig. 1). Turbulence thus effectively eliminates the thermal lobes.

Fig. 1 Perturbation of density arising from heat release obtained by subtracting a hot run (L = Lc) and a cold run (L = 0) for the case ß = 50 at t = 1 orbit. The perturbation was integrated over colatitude and normalized to $\gamma \left(\gamma -1\right)L/\chi c_s^2$. The solid purple line indicates the corotation radius.

Fig. 2 Temporal evolution of α, defined by Eq. (15), for the runs with L = Lc and β = 50 (left panel) and β = 1000 (right panel).

Fig. 3 Quality factor averaged over space and time. The dotted orange line shows the eight-cell limit. The vertical gray lines delineate the buffer zones.

3.2 Turbulence suppresses heating torques

The activation time of thermal torques is about one orbital period of the planetary embryo (Chametla & Masset 2021), whereas, as discussed above, the onset of MRI-driven turbulence occurs after about two orbital periods. In line with this, Fig. 5 shows that, for times t ≤ 1.5 orbits, the total torque Γ follows the behavior expected for the thermal torque in a laminar gas disk. During this phase, the torque is negative for L = 0, as expected for a non-luminous planetary embryo, whereas it becomes positive when L = L_c². Importantly, the torque behavior is largely insensitive to the value of β over this time interval (see the zoomed region in the right column of Fig. 5), reflecting the fact that turbulence has not yet developed in the disk.

For t > 1.5 orbits, the situation changes drastically: the torque increases substantially in all cases and exhibits a strongly oscillatory behavior, reflecting the onset of the turbulence in the disk. A careful comparison of the torque in the cases where β = 50, L = 0, and L = L_c (see the first and second rows in Fig. 5) suggests that the torque falls into the stochastic regime (Nelson & Papaloizou 2004; Nelson 2005; Baruteau & Lin 2010; Comins et al. 2016). In the case L = L_c, the torque is shifted to higher (more positive) values, while the overall shape of Γ(t) remains largely unchanged, even up to approximately three orbits. Furthermore, when the value of L is held fixed and β is increased from 50 to 1000, the amplitude of the torque oscillations is reduced (see the bottom row of Fig. 5). This indicates that the strength of the magnetic field governs the stochastic nature of the torque acting on the planetary embryo. We expect the stochastic nature of the torque to be largely insensitive if the luminosity or mass of the planetary embryo is increased, since thermal torques become significantly cut off once the inequality (λ_c/H)M_th < M_p < M_th is satisfied, where $M_{\text{th}}=c_s^3/G{\Omega }_p$ is the thermal mass (Velasco Romero & Masset 2020; Chametla & Masset 2021).

3.3 No traces of active magnetic resonances

When the toroidal magnetic field is sufficiently strong (β ≲ 4), the flow remains in a laminar regime, and magnetic resonances can make a significant contribution to the total torque. These resonances occur where the Doppler-shifted perturbation frequency in the fluid frame matches that of a slow MHD wave propagating along the magnetic field lines (Terquem 2003). The relative distance of these magnetic resonances from the planetary embryo is determined by the strength of the toroidal magnetic field through the β parameter (Terquem 2003; Uribe et al. 2015), $$\left|r_M-r_p\right|=\frac{2H_p}{3\sqrt{1+{\beta }_p}},$$(17)

where the subscript p indicates that the corresponding quantity should be evaluated at the position of the planetary embryo. As expected, when the magnetic field is very weak, β_p → ∞, the position of the magnetic resonances converges to the position of the planetary embryo. These resonances can dominate the Lindblad torques, as they are located closer to the embryo than the Lindblad resonances.

To assess whether magnetic resonances operate within the thermal lobes prior to the onset of turbulence, we compared the radial separation |r_M – r_p| with the characteristic cutoff scale λ_c. In our models, |r_M – r_p| = 4.7 × 10⁻³r_p and |r_M – r_p| = 1.05 × 10⁻³ r_p, for β = 50 and β = 1000, respectively. On the other hand, λ_c = 2.2 × 10⁻³ (from Eq. (12)). Thus, λ_c and |r_M – r_p| are comparable. In our models, however, we find no evidence that magnetic resonances are active in the vicinity of the planetary embryo. Note that the timescale for the development of large-scale density waves associated with magnetic resonances is on the order of a few orbital periods (Uribe et al. 2015). Consequently, additional perturbations beyond those produced by the thermal lobes may not appear within the single orbital period shown in Fig. 1.

In the case of a relatively weak magnetic field (ß ≫ 1), as in the models considered in this work, we have seen that the disk becomes turbulent. In principle, turbulence can suppress the heating torques, whereas the magnetic torques may persist, at least partially. Returning to Fig. 4, we find no evidence of MHD ring-like waves such as those reported by Fromang et al. (2005, see their Figure 3). We emphasize that the resolution used in our simulations is sufficient to resolve the characteristic density perturbation λ_c induced by planetary heating with at least ten grid cells, and therefore it is also adequate to resolve the radial separation |r_M – r_p|. This indicates that the absence of slow MHD waves is not a resolution effect, but rather a consequence of the turbulent state of the disk, in which the magnetic field fluctuates locally on timescales shorter than those required for slow MHD waves to develop (see also Terquem 2003; Comins et al. 2016). This result complements what was previously reported for an isothermal disk with a toroidal magnetic field, in which magnetic resonances were likewise found to be ineffective (Uribe et al. 2015).

Fig. 4 Gas density (logarithmic scale) at the midplane at t = 11 orbits. The zoomed regions correspond to an area similar to that shown in Fig. around the planetary embryo. In these regions, no low-density structures resembling hot lobes were observed. The solid purple lines indicate the corotation radii. The dashed yellow and dotted green lines mark the locations of the magnetic resonances for ß = 50 and ß = 1000, respectively. The white plus symbol denotes the position of the planetary embryo.

Fig. 5 Temporal evolution of the total torque Γ (in units of Γ0/γ) on the planetary embryo.

Fig. 6 Temporal evolution of the total torque Γ (in units of Γ0/γ) on a planet of mass Mp = 1 M⊕, for ß = 50 (blue lines) and ß = 1000 (purple lines).

3.4 Thermal torques for Earth-mass planetary cores

In the previous subsections, we focused on the survival of thermal torques for planetary embryos with masses close to that of Mars in a turbulent disk. Since thermal torques reach their maximum efficiency for planetary cores on the order of one Earth mass (Velasco Romero & Masset 2020), we ran two additional models with M_p = 1 M_⊕ to assess whether thermal torques can survive in this mass range. Fig. 6 shows the resulting torques for L = L_c and two values of ß (50 and 1000)³. As in the case of the planetary embryo, after 1.5 orbits, turbulence induces strongly oscillatory torques for both values of ß. We do not show maps of the gas density as they are very similar to those presented in Figure 4.

Large torque fluctuations, very similar to those shown in Figure 6, have also been observed in simulations of super-Earth planets (M_p > 3 M_⊕). In that case, they arise from vortices forming at the edge of the horseshoe region (Chametla & Masset 2021), as the gas dynamics become highly nonlinear with the Hill sphere for such planetary masses. Therefore, in regions where MRI turbulence develops, migration is expected to be stochastic for 0.03 M_⊕ ≲ M_p ≲ 1 M_⊕ as a result of MRI turbulence and for M_p ≳ 3 M⊕ owing to vortices. Finally, we stress that the turn-off effect on the torque reported here is distinct from the cut-off effect investigated in Velasco Romero & Masset (2020).

4 Toward more complex models

In our simulations, we adopted a constant thermal diffusivity $\chi ={10}^{-5}{r}_p^2{\Omega }_p$. This value was chosen to allow the formation of hot thermal lobes around the luminous planet, enabling us to quantify the heating torques and construct a thermally controlled environment prior to the onset of MRI. In real protoplanetary disks, however, χ is expected to vary spatially as it depends on local physical conditions, notably the gas density, temperature, and the opacity of the medium (Jiménez & Masset 2017, and references therein). These quantities are, in turn, influenced by the local dust concentration and by dust properties such as the grain size distribution and composition (Chametla et al. 2025). Therefore, more sophisticated models that include non-constant thermal diffusivity and an explicit treatment of dust are warranted and will be explored in future work.

5 Conclusions

In this study, we report the results of high-resolution 3D MHD simulations of a planetary embryo with a mass one-third of that of the planet Mars and for an Earth-mass core, embedded in a stratified gaseous disk with thermal diffusivity and a toroidal magnetic field. Our aim was to assess whether thermal torques can drive planetary migration in turbulent regions of a protoplanetary disk, such as the inner and outer regions of the disk located outside the so-called dead zone. To this end, we have considered two different luminosities for both the planetary embryo and core: L = 0 corresponding to a cold thermal torque, and L = L_c, corresponding to a heating torque. In all cases, we find that as turbulence develops in the disk, the thermal lobes are disrupted for both ß = 50 and ß = 1000. In this turbulent regime, the total torque acting on the Mars planetary embryo and Earth-sized core displays the characteristic pattern of a stochastic torque, implying that its migration should proceed in the stochastic migration regime, similar to that observed in Nelson & Papaloizou (2004) and Comins et al. (2016). Given the results of Chametla & Masset (2021), which show that larger planetary masses also experience stochastic torques, we conclude that planets with masses 0.03 M_⊕ ≲ M_p ≲ 1 M_⊕ and M_p ≳ 3 M_⊕ are insensitive to thermal torques in a disk where turbulence is driven by MRI or vortices. Lastly, in agreement with the results of Comins et al. (2016), we find that the level of turbulence attained for ß ≥ 50 renders magnetic resonances ineffective.

Acknowledgements

We thank the referee for the careful reading of the manuscript and for the constructive comments that helped improve the paper. The work of R.O.C, was supported by the Czech Science Foundation (grant 25-16507S). Computational resources were available thanks to the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254). Also, R.O.C. thanks Neal J. Turner for the fruitful discussions on this manuscript.

References

All Figures

	Fig. 1 Perturbation of density arising from heat release obtained by subtracting a hot run (L = Lc) and a cold run (L = 0) for the case ß = 50 at t = 1 orbit. The perturbation was integrated over colatitude and normalized to $\gamma \left(\gamma -1\right)L/\chi c_s^2$. The solid purple line indicates the corotation radius.
In the text

	Fig. 2 Temporal evolution of α, defined by Eq. (15), for the runs with L = Lc and β = 50 (left panel) and β = 1000 (right panel).
In the text

	Fig. 3 Quality factor averaged over space and time. The dotted orange line shows the eight-cell limit. The vertical gray lines delineate the buffer zones.
In the text

Fig. 4 Gas density (logarithmic scale) at the midplane at t = 11 orbits. The zoomed regions correspond to an area similar to that shown in Fig. around the planetary embryo. In these regions, no low-density structures resembling hot lobes were observed. The solid purple lines indicate the corotation radii. The dashed yellow and dotted green lines mark the locations of the magnetic resonances for ß = 50 and ß = 1000, respectively. The white plus symbol denotes the position of the planetary embryo.

In the text

	Fig. 5 Temporal evolution of the total torque Γ (in units of Γ0/γ) on the planetary embryo.
In the text

	Fig. 6 Temporal evolution of the total torque Γ (in units of Γ0/γ) on a planet of mass Mp = 1 M⊕, for ß = 50 (blue lines) and ß = 1000 (purple lines).
In the text