On Wave Interference in Planet Migration: Dead Zone Torques Modified by Active Zone Forcing

MHD equations describing the gas flow in 3D-MHD disks, considering a frame rotating with a uniform angular velocity Ω, are given by the equation of continuity

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }})\rho -\rho {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 1 }$$

the gas momentum equation

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}=-({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}}-\displaystyle \frac{1}{\rho }{\boldsymbol{\nabla }}p-{\boldsymbol{\nabla }}{\rm{\Phi }}+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{\rho }\\ -\,{\boldsymbol{\Omega }}\times ({\boldsymbol{\Omega }}\times {\boldsymbol{r}})-2{\boldsymbol{\Omega }}\times {\boldsymbol{u}},\end{array}\end{eqnarray} \tag{ 2 }$$

and the induction equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times \left[{\boldsymbol{u}}\times {\boldsymbol{B}}-\eta {\boldsymbol{J}}\right],\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&\rho ={\rho }_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-{q}_{\rho }}\end{eqnarray} \tag{ 4 }$$

is the gas density, which is related to the surface density in our cylindrical disk approximation by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Sigma }}}{{L}_{z}}=\displaystyle \frac{1}{2\pi {L}_{z}}\int \rho d\phi {dz},\end{eqnarray} \tag{ 5 }$$

where L_z represents the vertical extent of our computational domain.

In Equations (1)–(2), u is the gas velocity, Φ is the gravitational potential, and p is the pressure given as

$$\begin{eqnarray}&&p=\rho {c}_{s}^{2},\end{eqnarray} \tag{ 6 }$$

with c_s the sound speed. Note that the phrase aspect ratio (see Table 1) in unstratified 3D-MHD disks merely implies that c_s = 0.1rΩ_Kep everywhere in the disk, so the disk scale height H ≡ 0.1r. In the induction equation ${\boldsymbol{J}}={\mu }_{0}^{-1}{\boldsymbol{\nabla }}\times {\boldsymbol{B}}$, is the current density, and η is the resistivity modeled as in Lyra et al. (2015):

$$\begin{eqnarray}&&\eta (r)={\eta }_{0}-\displaystyle \frac{{\eta }_{0}}{2}\left[\tanh \left(\displaystyle \frac{r-{r}_{1}}{{h}_{1}}\right)-\tanh \left(\displaystyle \frac{r-{r}_{2}}{{h}_{2}}\right)\right],\end{eqnarray} \tag{ 7 }$$

which leads to a smooth resistivity transition between the inner dead zone and the outer active zone.

Table 1. Initial Conditions and Main Parameters of Our Simulations (r₀ Denotes a Reference Radius ^a )

	Parameter	Value (code units)
	Central Star
Mass of the Star	M⋆	1.0
	Disk
Aspect ratio at r0	h	0.1
Density at r0	ρ0	6 × 10−4
Surface density slope	qρ	1.5
Temperature slope	qT	1.0
Initial resistivity	η0	{2.28, 5} × 10−4
Initial magnetic field b	B0	{4.47 × 10−3, 10−2}
Magnetic permittivity
of vacuum	μ0	1.0
Plasma parameter at
r0	β	{2 × 102, 103}
Kinematic viscosity (only 2D α-disks)	ν0	2.56 × 10−6
	Planet
Planet mass interval	Mp	[3 × 10−6, 10−2]
Planet location	(rp , ϕp , zp )	(1, 0, 0)
Softening length		{0.2, 0.6}H(rp )
	Global Mesh
Radial extension	r	[0.4, 9.6]
Azimuthal extension	ϕ	[0, 2π]
Vertical extension	z	[−0.1,0.1]
Radial resolution	Nr	518 cells
Azimuthal resolution	Nϕ	1024 cells
Vertical resolution	Nz	64 cells
Radial spacing	log	Logarithmic
	Additional
	Parameters
Reference frame	C	Corotating
Frame angular speed	Ω(rp )	Ωp
Gravitational constant	G	1.0
Orbital period at rp	T0	$2\pi {{\rm{\Omega }}}_{p}^{-1}$

Notes.

^a We consider r₀ = 10 au and M_⋆ = 1M_⊙ when scaling back to physical units. ^b We use a net vertical magnetic field.

Download table as:  ASCIITypeset image

We use a reference frame centered on the star and rotating with the angular frequency ${\boldsymbol{\Omega }}={{\rm{\Omega }}}_{p}\hat{{\boldsymbol{k}}}$, where ${{\rm{\Omega }}}_{p}\equiv \sqrt{{{GM}}_{\star }/{r}_{p}^{3}}$ is the angular velocity of the planet (with r_p the planet's orbital radius) and $\hat{{\boldsymbol{k}}}$ is a unit vector along the rotation axis. The unit of time adopted when discussing the results is ${T}_{0}=2\pi {{\rm{\Omega }}}_{p}^{-1}$. The gravitational potential Φ is given by

$$\begin{eqnarray}&&{\rm{\Phi }}={{\rm{\Phi }}}_{S}+{{\rm{\Phi }}}_{p}+{{\rm{\Phi }}}_{\mathrm{ind}},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{S}=-\displaystyle \frac{{{GM}}_{\star }}{r}\end{eqnarray} \tag{ 9 }$$

is the stellar potential, with M_⋆ the mass of the central star and G the gravitational constant;

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}=-\displaystyle \frac{{{GM}}_{p}}{\sqrt{r{{\prime} }^{2}+{\epsilon }^{2}}}\end{eqnarray} \tag{ 10 }$$

is the planetary potential; and

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{ind}}=\displaystyle \frac{{{GM}}_{p}}{{r}_{p}^{3}}{\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{p}+G{\int }_{V}\displaystyle \frac{{dm}({\boldsymbol{r}}^{\prime} )}{r{{\prime} }^{3}}{\boldsymbol{r}}\cdot {\boldsymbol{r}}^{\prime} \end{eqnarray} \tag{ 11 }$$

is the indirect potential arising from the planet and the gravitational force of the disk, respectively. Here, M_p is the planet mass, $r^{\prime} \equiv | {\boldsymbol{r}}-{{\boldsymbol{r}}}_{p}|$ is the distance to the planet, and is a softening length used to avoid a singularity of the potential in the vicinity of the planet. In 3D disks when the MHD is included, we use a softening length of = 0.2H(r_p ).

The equation of motion for the planet can be written as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{2}{{\boldsymbol{r}}}_{p}}{{{dt}}^{2}}=-\displaystyle \frac{G({M}_{\star }+{M}_{p})}{{r}_{p}^{3}}{{\boldsymbol{r}}}_{p}+{{\boldsymbol{f}}}_{\mathrm{disk}}\\ \,\,-\,2{\boldsymbol{\Omega }}\times {{\boldsymbol{v}}}_{p}-{\boldsymbol{\Omega }}\times ({\boldsymbol{\Omega }}\times {{\boldsymbol{r}}}_{p}),\end{array}\end{eqnarray} \tag{ 12 }$$

with the acceleration due to the gas disk

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{disk}}=-G{\int }_{V}\displaystyle \frac{\rho ({\boldsymbol{r}})({{\boldsymbol{r}}}_{p}-{\boldsymbol{r}})}{{\left(r{{\prime} }^{2}+{\epsilon }^{2}\right)}^{3/2}}{dV},\end{eqnarray} \tag{ 13 }$$

where v _p is the planet velocity. In Equations (11) and (13), the integrations are performed over the whole volume V of our computational domain.

In the case of our 2D α-disk models, the continuity and Navier–Stokes equations are solved on a polar grid centered on the star and in a frame corotating with the planet. A locally isothermal equation of state is used with the vertically integrated gas pressure P and the isothermal sound speed c_s satisfying $P={\rm{\Sigma }}{c}_{s}^{2}$, with c_s kept time independent. The source terms in the Navier–Stokes equation include the gravitational potential due to the star, the planet, and the indirect terms arising from the stellar acceleration imparted by the disk gas and the planet (now the integrals in Equations 11 and 13 are performed on a surface S of the 2D computational domain). In the planet potential, a softening length = 0.6H(r_p ) is included to mimic the effects of finite vertical thickness, with H ≡ c_s /Ω the disk's pressure scale height. A viscous acceleration is also included with a turbulent kinematic viscosity ν, which is related to the dimensionless alpha parameter by ν = α c_s H. It should be noted that the inclusion of this viscosity is motivated in principle by the results of the MHD simulations since, considering the resistivity transitions, it is found that the Reynolds stress is not null within the dead zone (see Lyra et al. 2015). The initial radial profiles of the surface density and temperature are power-law functions with exponents q_ρ and q_T , respectively.

The main purpose of the presence of the 2D HD simulations is only to compare the migration rate with our MHD models; therefore, details of planetary migration in viscous disks can be found elsewhere (see, for instance, Paardekooper et al. 2022 and references therein). We recall that the width considered for the transition zone (h₁ = h₂ = 0.8) does not produce any formation of vortices in the viscous simulation studied in Lyra et al. (2015). Therefore, we have not included any viscosity transition zone for these models.

To numerically solve Equations (1)–(6), we use the publicly available magnetohydrodynamic code FARGO3D ⁴ (Benítez-Llambay & Masset 2016) with MHD-orbital advection enabled (see Masset 2000; Benítez-Llambay & Masset 2016 for details). The FARGO3D code solves the HD equations with a time-explicit method, using operator splitting and upwind techniques on an Eulerian mesh. The update of the magnetic field governed by the induction equation (Equation (3)) is done by the method of characteristics (Stone & Norman 1992), and the constrained transport method (Evans & Hawley 1988) is used to preserve the divergence-free property of the magnetic field. Orbital evolution of the planets is integrated by a fifth-order Runge–Kutta method (Cash & Karp 1990) with a time step provided by the Courant–Friedrichs–Levy (CFL) condition of the MHD algorithm.

Table 1 summarizes the set of disk parameters and initial conditions used in our simulations. Although we consider the same extent of disk components as in Lyra et al. (2015) where one has r = [0.4, 9.6]r₀, ϕ = [ − π, π], and z = [ −0.1, 0.1] and the same spacing (logarithmic in radius and uniform in azimuthal and vertical ⁵ directions, respectively), we use a higher resolution (N_r , N_ϕ , Nz) = (518, 1024, 64). We note that the results presented here are valid also for a lower resolution.

We use periodic boundary conditions in the azimuthal and vertical directions. To avoid reflections at the radial boundaries of our computational domain, we use damping boundary conditions as in de Val-Borro et al. (2006), the width of the inner damping ring being 3.9 × 10⁻² r₀ and that of the outer ring being 8.54 × 10⁻¹ r₀. The damping timescale at the edge of each damping ring equals one-twentieth of the local orbital period. Following Lyra et al. (2015), we consider a net vertical magnetic field and ramp it smoothly to zero inward of r = 2.5. The resistivity profile is given by Equation (7); the dead-to-active transition in the outer disk is located at r₀ = 10 au (Dzyurkevich et al. 2013), which in our scale-free units translates to r₀ = 1. We set r₁ = 2.5, r₂ = 10, and h₁ = h₂ = 0.8 in Equation (7) in all our 3D-MHD simulations. Figure 1 shows the resulting resistivity profile in which a smooth transition can clearly be seen in the range r = 1–4, centered at r = 2.5. The dashed vertical lines in Figure 1 mark the radial region around the transition center (±2H) where the RWI occurs.

Zoom In Zoom Out Reset image size

Figure 2 shows the radial density profiles as t = 250 T₀. Regardless of the planet mass (M_p = 1, 10, 30, and 100 M_⊕), a density bump develops at the center of the resistivity transition, which gives rise to the RWI. For MR2, the density bump is slightly flattened. For MR4 (intermediate-mass planet with M_p = 100 M_⊕), a second peak appears at r = 1.4r₀ due to the fact that the planet opens a partial gap.

Zoom In Zoom Out Reset image size

In Figure 3 we show the 2D map of the gas density in the midplane (z = 0) at t = 250 T₀ of our MR1–MR4 models. In all cases, we see spiral waves excited by the turbulence in the active zone and by the RWI vortex at the resistivity transition. Consequently, the planetary migration can be modified by the background spiral waves because they generate a gas flow in the horseshoe (HS) region of the planet, and additionally, they generate constructive (destructive) interference with the outer (inner) spiral waves emitted by the planet.

Zoom In Zoom Out Reset image size

Figure 4 shows the evolution of the semimajor axis a_p and the eccentricity e_p for the AR1–AR4 and MR1–MR4 models. We find that the migration of low-mass Earth-like and Neptune-like planets (M_p ≤ 30M_⊕), when MHD is included, is directed inwards and proceeds faster compared to the case of α-disks.

Zoom In Zoom Out Reset image size

3.1.1. Fast Inward Migration

For a locally isothermal 3D disk with a power-law surface density profile, the total torque (Lindblad plus corotation), Γ, acting on a planet was estimated by Tanaka et al. (2002). They found that

$$\begin{eqnarray}&&{\rm{\Gamma }}=(1.364+0.541{q}_{\rho }){{\rm{\Gamma }}}_{0},\end{eqnarray} \tag{ 14 }$$

where Σ_p and c_s are the surface density and the sound speed at the planet position, respectively. Here, the total torque scales with ${{\rm{\Gamma }}}_{0}\equiv {({M}_{p}/{M}_{\star })}^{2}{({r}_{p}{{\rm{\Omega }}}_{p}/{c}_{s})}^{2}{{\rm{\Sigma }}}_{p}{r}_{p}^{4}{{\rm{\Omega }}}_{p}^{2}$. From Equation (14), the type-I radial migration speed of the planet can be calculated from the conservation of angular momentum as

$$\begin{eqnarray}&&\displaystyle \frac{{{dr}}_{p}}{{dt}}=2{r}_{p}\displaystyle \frac{{\rm{\Gamma }}}{{L}_{p}},\end{eqnarray} \tag{ 15 }$$

where ${L}_{p}={M}_{p}\sqrt{{{GM}}_{\star }{r}_{p}}$ is the planet's angular momentum. Using Equations (14) and (15), we can estimate the migration timescale as

$$\begin{eqnarray}&&\tau ={(2.7+1.1{q}_{\rho })}^{-1}{M}_{p}{r}_{p}^{2}{{\rm{\Omega }}}_{p}{{\rm{\Gamma }}}_{0}^{-1}.\end{eqnarray} \tag{ 16 }$$

To investigate the migration rates further, we compared the obtained values with $\tau ={r}_{p}{\left(\tfrac{{{dr}}_{p}}{{dt}}\right)}^{-1}$ resulting from our simulations ⁶ with the values given by Equation (16).

As in 3D models of cylindrical turbulent disks (see Nelson & Papaloizou 2004), we have been careful to consider a smoothing length for the planet potential appropriate for the 2D simulation set. This allows us to make a proper comparison with Equation (16) for 3D isothermal disks and also gives us the possibility to make the comparison with our 3D-MHD models. The comparison of migration timescales is shown in Figure 5 as a function of M_p . As expected, the 2D α-disk simulations are in agreement with the migration time appropriate to 3D locally isothermal disks. In fact, the greatest difference (14%) occurs between the AR4 and MR4 models.

Zoom In Zoom Out Reset image size

Therefore, α-disks as well as Equation (16) can provide reference values when quantifying how τ differs in the dead zone of disks when the MHD is included. Here, we find that the migration time for M_p < 100M_⊕ in α-disks is ≈1.6 times larger than the migration time of planets in disks where there are background spiral waves propagating through the dead zone.

It is worth mentioning that the migration time of the planets does not depend on the strength of the initial magnetic field. We have analyzed the orbital evolution and the migration time for the MR2b model (those results are not shown here), and we have not found significant changes when the value of the initial magnetic field is increased.

In the case of the intermediate-mass planet with M_p = 100M_⊕, we find that the migration is only slightly faster in the simulation when MHD is included for t ≤ 800 T₀ (τ in the α-disk simulation is ≈1.15 times larger compared to the 3D-MHD disk simulation). At t > 800 T₀, the migration becomes stagnant at a distance from the central star that is larger than in the α-disk simulation (see lower right panel in Figure 4).

3.1.2. Time Evolution of the Eccentricity

Furthermore, we find that when the MHD is included in our disk simulations, the eccentricity exhibits oscillatory behavior, and the amplitude of oscillations varies with time. For MR1 and MR2, there are episodes when the eccentricity suddenly becomes pumped up to e_p = 0.02–0.03 at t = 750T₀ and t = 1050T₀, respectively.

This increase in eccentricity is due to transient vortices migrating from the outer edge of the dead zone inwards as they interact both directly and indirectly with the planet by modifying the background spiral waves. Once these vortices vanish, the eccentricity begins to decrease as can be seen in the upper right panel of Figure 4. For MR3 and MR4 models, the amplitude of oscillations can reach values around 10⁻³.

Although for planets embedded in α-disks with masses greater than or equal to 20M_⊕, the damping time of the eccentricity differs a little from the linear theory, it is known that the eccentricity should drop exponentially to zero in a timescale τ_ecc ∝ (H/r)² τ (Cresswell et al. 2007). However, we find that in our MHD-disk models, the eccentricity of low-mass planets is not damped in the same way. In fact, it shows a slightly increasing behavior on average (at least for M_p ≤ 30M_⊕; see inset boxes in Figure 4).

In addition, Nelson (2005) found that turbulence is the main source of eccentricity driving with a strong dependence on the planet mass (see his Figures 5, 6, 8 and 10). For instance, he found that a mass planet of M_p = 1M_⊕ can reach an eccentricity between 0.02 and 0.08, and for a planet with M_p = 10M_⊕, the value of the eccentricity oscillates in the range 0.02 ≤ e ≤ 0.03, over the whole simulation time (∼500 orbital periods). On the other hand, in this study we find that for M_p ≤ 30M_⊕, the values of the eccentricities driven by the background spiral waves, although these also show an oscillatory behavior, are an order of magnitude smaller within the same time period than in Nelson (2005).

Zoom In Zoom Out Reset image size

In this section, we study gap opening by massive planets in a fixed circular orbit, embedded in the dead zone. Since we are interested in analyzing the effect of background spiral waves on the shape and depth of the gap formed by a massive planet, we make the following considerations (in both MHD and α-disk models):

• 1.  
  The planet is not allowed to migrate during the entire simulation time. One reason for this (apart from that related to viscous relaxation) is to see the evolution of the vortices that can form at the edges of the gap.
• 2.  
  To avoid spurious shocks that would arise if the planet were initialized with its final mass, we introduce the mass of the planet smoothly as a function of time given by

  $$\begin{eqnarray}{M}_{p}(t)=\left\{\begin{array}{lll}\displaystyle \frac{1}{2}\left[1-\cos \left(\displaystyle \frac{\pi t}{{t}_{\mathrm{ramp}}}\right)\right] & \mathrm{if} & t\lt {t}_{\mathrm{ramp}}\\ 1 & & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 17 }$$

  Here, we consider t_ramp = 200T₀, which is a suitable timescale for α-disks (Lega et al. 2021).

With the above in mind, we consider the AR5–AR7 and MR5–MR7 models (see Table 2). In Figure 7 we show the vertically and azimuthally averaged density profiles of the dead zone for disks when the MHD is included, as well as azimuthally averaged density profiles for the case of 2D α-disks at t = 250 T₀. Focusing on AR5 and MR5 models first, a higher-density bump is formed at r = 0.6 in MR5. The reason is that at this orbital time, the gap depth is already considerable since the gap bottom has been emptied by a little more than 70%. On the other hand, at radii r > r_p , there is no well-defined density bump in the disk when the MHD is included compared to the α-disk. Additionally, we find that the width of the gap is larger in the former case, and we also notice that an increase in the density profile occurs from r ≈ 1.8 outwards. In summary, the planet in the dead zone when the MHD is included is more efficient at removing gas from the gap region, that is, the density peaks at the edges of the gap are spread out, and the bottom of the gap is deeper.

Zoom In Zoom Out Reset image size

In fact, it can be seen in Figure 7 that the shape of the density profile for MR5 is somewhat similar to the density profile formed in AR6. In other words, a Jupiter-mass planet in the dead zone of a disk threaded by background density waves creates a gap with features that resemble a gap formed by a planet of 3 Jupiter masses in a disk governed by the classical viscous evolution. In the case of the MR6 model, the density profile in the inner part (r < r_p ), does not show any density bumps in comparison with the AR6 model. Similar to the previous case, for the outer part of the disk (r > r_p ) in model MR6, the density profile shows an increase from r ≈ 1.5; however, the width and depth of the gap are greater than in the AR6 model.

On the other hand, for AR7 and MR7 models, the density at the bottom of the gap is very low. In this last case, we find that the shape of the gap in both models only differs in the outer part of the disk since the gap in the α-disk has a smaller width and shows a density bump at r = 2. Note that in the AR7 model there is again an increase in the density profile of the disk starting at r ≈ 2.7.

In Figure 8 we show a comparison of the 2D-vortensity maps at the midplane (z = 0), between the AR5–AR7 and MR5–MR7 models at t = 250T₀. Let us focus on the comparison of the MR5 and AR5 models. In the case of the AR5 model, a large elongated vortex can be seen to form at the outer edge of the gap, while in the MR5 model in the same region, a very weak vortex of much shorter length is observed (which may be the result of the formation of a transient vortex). We argue that at this orbital time (t = 250T₀), the background spiral waves have removed enough density from the outer edge of the gap (and within the gap itself), which is accumulated on the inner edge of the gap, producing a density bump where the inner vortex forms. This vortex approaches its steady-state strength at t = 450T₀ and survives until the end of the simulation time. Note that in both models there is formation of a vortex in the coorbital region behind the planet. Lastly, for more massive planets like those considered in the MR6–MR7 models, we did not find vortex formation at either of the two edges of the gap (see Figure 8), which may be due to flattening of the gap edges by the BSWs.

Zoom In Zoom Out Reset image size

Although the spiral waves generated by a vortex shows a spiral wave pattern, which is reminiscent of an embedded planet, they do not have the same origin. The former is the result of Rossby waves, which are vorticity waves propagating on the gradient of vortensity (Lovelace et al. 1999; Li et al. 2000; Meheut et al. 2013). From linear theory, their dispersion relation is given as (Meheut et al. 2013)

$$\begin{eqnarray}&&\tilde{\omega }=\displaystyle \frac{2{\kappa }^{2}{c}_{s}^{2}}{{c}_{s}^{2}({k}_{r}^{2}+{m}^{2}/{r}^{2})+{\kappa }^{2}}\displaystyle \frac{m{ \mathcal L }^{\prime} }{r{\rm{\Sigma }}},\end{eqnarray} \tag{ 18 }$$

with a wave frequency ω. In Equation (18), $\tilde{\omega }\equiv \omega -m{\rm{\Omega }}$ is the wave frequency in a rotating frame; κ is the epicyclic frequency; m and k_r are the azimuthal and radial numbers, respectively; and ${ \mathcal L }^{\prime}$ represents the radial derivative of the quantity ${ \mathcal L }$ (which is related to the inverse of vortensity). Due to differential rotation, the Rossby waves are coupled to the density waves (Tagger 2001; Bodo et al. 2005), and the latter can only propagate outside of Lindblad resonance positions given by

$$\begin{eqnarray}&&\omega -m{\rm{\Omega }}({r}_{\mathrm{LR}})=\pm \kappa ({r}_{\mathrm{LR}}).\end{eqnarray} \tag{ 19 }$$

In fact, when the linear phase of the Rossby waves saturates, the amplitude of the density spiral waves is higher at the Lindblad resonances closer to the corotation radius (defined for the vortex; Meheut et al. 2013).

Note that we are interested in the amplitude of these density waves when they interfere (constructively or destructively) with the spiral waves emitted by the planet. Therefore, a study to parameterize the amplitude and nonlinear saturation criteria of background spiral waves (created by RWI and from turbulence in the active region of the disk) are beyond the scope of this work. With all of the above in mind, we now focus on discussing how the background spiral waves interfere with the (inner/outer) spiral arms of the planet.

In Figure 9 we show the shape of the background spiral waves in the dead zone of the disk (top row) and their physical amplitudes, δΣ/Σ₀ ≡ (Σ − Σ₀)/Σ₀, at the planet position r = r_p (bottom row) for the MR1–MR4 models. In the 2D density maps of this figure, it can be seen that when the mass of the planet increases, the spiral waves of the planet emerge in the density contrast and suffer interference from the background spiral waves.

Zoom In Zoom Out Reset image size

Since the peak in the δΣ/Σ₀ curve located at r = r_p and ϕ = 0 can be identified as the primary arm of the planet (Miranda & Rafikov 2019), we can compare the amplitude of the background spiral waves with the spiral waves emitted by the planet for models MR1–MR4. From Figure 9 (bottom row) it can be seen that for the MR1 model, the amplitude of the background spiral waves is much greater than the density waves emitted by the planet, which explains why they cannot be seen in the 2D density maps. For the MR2 model, both amplitudes are approximately equal. In the case of the MR3 model, the amplitude of the density waves emitted by the planet is greater than the amplitude of the background spiral waves, and the interference of the background spiral waves can be easily observed in the 2D map in Figure 9. Lastly, for the MR4 model, the amplitude of the density waves generated by the planet is much greater than that of the background spiral waves, and because the planet begins to open a gap, the shape of the background spiral waves is modified although its effect of interference remains. Note that, as the amplitude of the waves emitted by the planet increases, the amplitude of the background spiral waves decreases.

Therefore, we conjecture that the fast migration found in 3D-MHD disks is mainly due to the fact that the background spiral waves interact strongly with the spiral arms of the planet (and with the planet itself), substantially modifying the Lindblad's torque (and also corotation torque when present, as we discuss in the next section). Specifically, these propagating toward the central star may eventually overlap with the planet's outer spiral arm (i.e., the pitch angle of the latter coincides with the pitch angle of the background spiral waves) and produce a more negative differential Lindblad torque. In other words, the outer spiral arm of the planet suffers from constructive wave interference.

As a demonstration of this hypothesis, let us consider the MR3 model. Figure 10 shows the background spiral waves and spiral waves of the planet. In this case, it can be clearly seen that the outer planet's spiral arm suffers from constructive interference (since their pitch angles match), while the planet's inner spiral arm suffers from negative interference. In addition, it can be seen from this figure that the interference in the outer arm of the planet occurs exactly at some external Lindblad resonances, contrary to the inner arm, where the interference does not occur at the internal Lindblad resonances.

Zoom In Zoom Out Reset image size

On the other hand, when the resonances are closer to the planet (i.e., in the radial region r_p ± H, which contributes more to the Lindblad torque), we find that while constructive interference from the outer arm of the planet with the background spiral waves do not make a total match, the latter produce a redistribution of density in this region that drives the strengthening of the outer wake of the planet continuously.

In other words, the effect of the background spiral waves builds the density toward the planet when r > r_p and pushes it away from the planet for r < r_p in addition to producing a greater flux of angular momentum in the outer part of the dead zone (r > r_p ). To verify the latter, we have calculated the angular momentum flux in the dead zone of the disk for the MR3 model at t = 210T₀ by means of the expression

$$\begin{eqnarray}&&{F}_{\mathrm{BSW}}={\int }_{0}^{2\pi }({\rm{\Sigma }}{{rv}}_{\phi }^{\prime} ){v}_{r}^{{\prime} }{rd}\phi ,\end{eqnarray} \tag{ 20 }$$

where ${v}_{r}^{{\prime} }={v}_{r}-{\bar{v}}_{r}$ and ${v}_{\phi }^{{\prime} }={v}_{\phi }-{\bar{v}}_{\phi }$, with ${\bar{v}}_{r}$ and ${\bar{v}}_{\phi }$ the radial and azimuthal components of the velocity averaged azimuthally. The result of this calculation is shown in Figure 11 where we also included the case of the flux of angular momentum in the disk without a planet. A direct comparison of the angular momentum flux in the disk for the MR model with the case of M_p = 0 (no planet in the disk) shows a larger angular momentum flux on the outer part of the dead zone when the interference effect of BSWs takes place.

Zoom In Zoom Out Reset image size

Here it is important to emphasize that the angular momentum flux is carried by a set of spiral waves launched "continuously" by the vortex at the transition zone and the turbulence of the active zone. Therefore, the planet's spiral arms can interact many times with the same or different background spiral waves within a planet's orbital period, unlike an interaction of a planet with the spiral arms of a second planet for example, where the interaction occurs in the HS region (see Cui et al. 2021).

As we mentioned earlier in Section 1, several studies on turbulent disks have shown that the total torque on a migrating planet can present a strongly oscillatory behavior, leading to a random walk migration (Laughlin et al. 2004; Nelson & Papaloizou 2004; Papaloizou et al. 2004; Nelson 2005).

The main cause of random walk migration is the density fluctuations in the disk, which in turn causes in most cases that the wakes generated by the planet are also turbulent. Furthermore, density perturbations generated in the turbulent zone of the layers in a disk can propagate into the dead zone and drastically change the torque on the planets (Oishi et al. 2007), which again can be stochastic as in a fully turbulent disk.

In our case, we show that although the torques on the planet present both negative and positive oscillations (see Figure 6), on average they result in a negative torque in most cases, thus generating a faster convergent migration toward the central star (see Figure 4).

Recently, McNally et al. (2020) showed that the migration of low-mass planets in near-inviscid 3D disks accelerates the inward migration due to the buoyancy torques arising as a consequence of the disk response through vertical motions to the planet potential at buoyancy resonances, which leads to evolution of the vortensity of the libration region.

For a slow radial accretion flow, planet migration would still be accelerated by the evolution of the corotation torque since the buoyancy response of the disk causes the vortensity of the coorbital disk material to evolve too quickly.

Although, McNally et al. (2020) found that the gas flow rate in the disk surface layers (which corresponds to the fast flow case studied in McNally et al. 2018b) has minimal impact on inward migration, a laminar radial gas flow in the midplane induced by Hall effect (see McNally et al. 2017, 2018b) may be a way to counteract the influence of the buoyancy response of the disk on the torque of the corotation.

It should be mentioned that although in this study we consider the density perturbations generated from the active zone of the disk, it can be said that they are density perturbations that show a well-defined pattern that induces a different effect on the torques on the planet compared to the constant radial gas flow studied in McNally et al. (2017, 2018b) and the bouyancy torques studied in McNally et al. (2020). Note that latter effect is only present in 3D stratified disks. Therefore, a direct comparison with the background spiral wave effect presented here is not possible.

Since we find that the effect of the spiral waves propagating in the dead zone of the disk mainly modifies the Lindblad torque when the corotation torque is not active, a natural question that arises is what happens when the corotation torque is also present on the disk. To Section 4.2.1 we redirect this question.

4.2.1. Null Vortensity Gradient in the Dead Zone

In the case of the corotation torque, it is well known that it can depend on the radial gradients of vortensity, entropy, and temperature, as well as on the viscous production of vortensity (Casoli & Masset 2009; Masset & Casoli 2010; Paardekooper et al. 2010; Baruteau & Masset 2013; Jiménez & Masset 2017). The latter arises from the contact discontinuities on the downstream separatrices if there is an entropy gradient (Jiménez & Masset 2017).

On the other hand, in an isothermal disk, it holds that, if the surface density and angular frequency are power laws of the distance to the star, the vortensity gradient can become null for a particular value of the density slope (Jiménez & Masset 2017).

So, the choice of the value of the surface density slope q_ρ = 3/2 (and that we also consider a locally isothermal disk) gives us the possibility of being able to analyze in more detail the effect of the background spiral waves on the total torque on the planet since the vortensity gradient (equal to q_ρ − 3/2) is null.

Note that if the vortensity gradient is null, the contribution of the corotation torque to the total torque may be negligible, and therefore, the total torque is governed only by the Lindblad torque.

Finally, to support our hypothesis that the fast inward migration is mainly due to the effects produced by spiral waves propagating in the dead zone and thus affecting the density in the HS region and the spiral arms created by the planet, we performed a set of numerical simulations (MHD and α disks) considering a surface density of the disk that decays radially with q_ρ = 1/2.

In this way, the vortensity gradient is no longer zero, and the corotation torque can make an additional contribution to the total torque, which could modify the migration of the planet. Figure 12 shows the temporal evolution of the semimajor axis of the planet for the AR2 and MR2 models, in both cases with two different slopes for the radial density of the disk (q_ρ = 1/2 and q_ρ = 3/2).

Zoom In Zoom Out Reset image size

It is clear that when q_ρ = 1/2, the corotation torque has a positive contribution to the total torque, so the migration of the planet is slower, but in both cases the relative difference in the migration most of the time is very similar for both disk models (α-disk and MHD disks). However, particular attention should be paid to time t ≲ 50T₀ in Figure 12 for the MR2a model. Since within this time interval, the migration is the same as for the AR2a model, this behavior can be explained by means of two effects on the corotation torque due to background spiral waves:

• 1.  
  The corotation torque still does not saturate, and the perturbations of the background spiral waves are very weak.
• 2.  
  This time interval is how long it takes for the background spiral waves to deactivate the corotation torque by destroying the dynamics of the HS region. In other words, from t > 50T₀ the dynamics of the HS region no longer exists.

We think that this behavior in the migration of the planets in both disk models is due to point (ii) (which should be studied in detail in a subsequent study; R. O. Chametla et al. 2023, in preparation). Let us show a hint that the dynamics of the HS region breaks down after a certain orbital period. In Figure 13 we show the streamlines overlaying the gas density at the midplane, calculated at t = 5, t = 60 and t = 80 orbits. In a short time (t = 5T₀), it can be seen that the streamlines define a clear pattern of the HS region in the corotation region of the planet. At a time of t = 60T₀, there is no well-defined behavior of the streamlines that form the HS region since much greater disturbance is seen from the BSW. Lastly, at t = 80T₀, the streamlines do not define an HS pattern.

Zoom In Zoom Out Reset image size

To further support our statement that the fast migration found in the MR2a model with respect to the AR2a model is due to the destruction of the dynamics of the HS region, we calculated the radial torque density. In Figure 14 we show the results of this calculation for different orbital times. At t = 5T₀ the torque is fully unsaturated since it shows a well-defined peak close to r = 1, which is characteristic of the contribution of the corotation torque (see Kley et al. 2012). We find that behavior of the torque density at this short time in both models is basically the same. On the contrary, for orbital times of t = 60T₀ and t = 80T₀ we find that the torque density in the MR2a model shows a strongly oscillating behavior near of the planet. This oscillatory behavior of the torque density shows a lower positive contribution near the planet than in the case of the AR2a model. Note that for a time t = 80T₀ the torque density for the AR2a model still remains in the unsaturated regime.

Zoom In Zoom Out Reset image size

On the other hand, if it were due to the fact that the perturbations of the spiral waves are very weak, in the case when the corotation torque is not active (q_ρ = 3/2), both planets should have a similar behavior at the beginning of their migration, which is ruled out because the semimajor axes are separate initially.

The migration of massive planets embedded in gas disks is mainly governed by the opening of a gap. Recently, Lega et al. (2021) studied the migration of massive planets in adiabatic gas disks considering a low viscosity. They found that the formation of a vortex at the outer end of the gap is feasible and conclude that two migration modes are possible that differ from classical type-II migration in the sense that they are not proportional to the disk's viscosity.

The first migration mode occurs when the gap opened by the planet is not very deep. This occurs when a big vortex forms at the outer edge of the planetary gap, diffusing material into the gap. The second migration mode occurs when the gap is deep so that the planet's eccentricity grows to a value 0.2, which can produce a slower migration and even reverse it. With all this in mind, we have investigated the background spiral waves effect on the opening of the gap and the migration of a massive planet in the dead zone of an unstratified 3D-MHD disk with weak resistivity gradient.

Our results show that the background spiral waves propagating inward can generate a greater flow of gas within the gap as well as induce the formation of a vortex at the inner edge of the gap for the MR5 model (see Figure 8), and for all the models studied (MR5–MR7), the background spiral waves prevent the formation of vortices at the outer edges of the gaps. Therefore, this could modify the migration of a massive planet with respect to the standard type-II migration and also regarding the migration found by Lega et al. (2021).

Remarkably, we find that the gap formed by a massive planet in the dead zone (where there are spiral waves propagating, and therefore a radial flow of gas occurs toward the inner part of the disk) is asymmetric.

Considering such asymmetry in the gap and the spiral waves passing through the gap could definitely unlock the planet from the viscous evolution of the disk to give rise to a different behavior in their migration.

Although the fact that massive planets embedded in turbulent disks can open a wider and deeper gap than in the case of a laminar disk was already reported by Nelson & Papaloizou (2003) and also by Zhu et al. (2013), the novelty of our results is that a massive planet can also open a wider and deeper gap in the dead zone of the disk, in addition to the fact that said gap can have an asymmetric shape.

As we mentioned in Section 1, in McNally et al. (2019) a "city map" diagram of the different mechanisms that can modify planetary migration was introduced. We think that the effect of the background spiral waves could be located in the last three blocks of the upper left corner of this figure. Since, as we have discussed before, our results suggest that the background spiral waves effect can modify the redistribution of the gas density in the dead zone, it can therefore change the total torque on low- to intermediate-mass planets (see Figure 6). For massive planets, it can modify the gas flow through the gap formed by the planet and also induce/prevent the formation of vortices at the edges of the gap depending on the planet mass. It is important to note that the effect of background spiral waves is in part the result of RWI in the inner/outer transition regions of the disk dead zone as several previous studies have shown (Gammie 1996; Varnière & Tagger 2006; Turner & Drake 2009; Lyra & Mac Low 2012; Lyra et al. 2015; Faure & Nelson 2016). Therefore, it can be considered a natural and frequent process in the temporal evolution of the protoplanetary disk.

Nelson (2005) has studied the orbital evolution of low-mass planets (M_p ≤ 30M_⊕) embedded in unstratified, full turbulent magnetized disks although we cannot make a direct comparison between our results and those presented in said study since there are minor differences. This work allows us to show the main differences in the orbital evolution of a planet migrating in the dead zone of the disk, perturbed by background spiral waves, with respect to its orbital evolution in a full turbulent disk. Nelson (2005) found in all his models that the planets undergo a stochastic migration, whereas we find (at least for the same range of planet masses) a well-defined pattern of rapid migration toward the central star (see Figure 4).

Another marked difference between results presented in Nelson (2005) and our results is that the long-term torque fluctuations dominate over the type-I migration torques. In other words, the torque felt by a migrating planet in a turbulent disk is less negative (even positive) than in a viscous laminar disk (see his Figure 12). While we find that although the torque on the planet when the MHD is included in the disk has an oscillatory behavior, the said torque on average does not exceed the torque in an α-disk (the latter is less negative), as can be seen in Figure 6.

Note that we have not included in this comparison the case when the planet has a mass of M_p = 30M_⊕ since Nelson (2005) includes three planets in the same simulation, and although they do not get too close to each other, the turbulent density fluctuations are of similar amplitude to the spiral wakes generated by the planets in this case, and therefore, the spiral waves emitted by each of the three planets can also modify the total torque and cause a different migration (Baruteau et al. 2014; Cui et al. 2021; Chametla & Chrenko 2022).

It is important to note that we have neglected the dependence of the gravitational potential on z along with the vertical stratification of the disk for reasons of computational cost. Thus, our simulations are of cylindrical disks (e.g., Armitage 1998; Hawley 2000, 2001; Steinacker & Papaloizou 2002). These simulations have been used in previous studies to show the effects of the magnetic field on gas dynamics and for planet–disk interaction in global simulations (Nelson & Papaloizou 2003, 2004; Nelson 2005; Lyra et al. 2015).

Although these types of simulations have the disadvantage that they can overestimate the results of the total torque on a planet (which applies to our models), they also give us the possibility to run our numerical models over a larger orbital time.

The question of how much the fast migration toward the star can change if the vertical stratification is included must be redirected to a subsequent work (R. O. Chametla et al. 2023, in preparation). However, it is reasonable to infer that the perturbations in density reported here can modify the total torque on the planet even if disk stratification and the potential dependence on z are included since the interference produced by spiral waves propagating in the disk on the spiral arms of the planet can produce an oscillating behavior of the torque even in 2D simulations (Chametla & Chrenko 2022).