© The Authors 2023

1 Introduction

The question of accretion and its origin in the context of protoplanetary disks is an utmost step if we want to develop reliable and realistic models, whether it be of the evolution of the disk itself (Pascucci et al. 2023), of its constituents such as the gas, the dust, and the planets, but also of their interactions (see, e.g., Rabago & Zhu 2021). Observational evidence of an ultraviolet excess in the spectral energy distribution gives, via continuum fitting methods, a quantitative estimate of the stellar accretion rate, with a typical value of 10⁻⁸ M_⊙ yr⁻¹ (Venuti et al. 2014). In order to explain such high values of the stellar accretion rate, one needs to rely on mass and angular momentum transport models. The conventional viscous model relies on a turbulent transport of angular momentum. It aims to parameterize the disk’s turbulent viscosity v with an α_v parameter (Shakura & Sunyaev 1973). Turbulence in disks is in general supposed to be the consequence of the nonlinear saturation of the magneto-rotational instability (MRI, Balbus & Hawley 1991; Lesur et al. 2022). Although the MRI is the most promising instability to explain turbulence-driven accretion, recent global simulations of protoplanetary disks with nonideal magneto-hydrodynamic (MHD) effects have shown that the MRI may not occur in most regions of protoplanetary disks (see, e.g., Perez-Becker & Chiang 2011b,a; Bai & Stone 2013; Lesur et al. 2014). In addition to the MRI, a wide range of hydrodynamic instabilities can play a role in the production of turbulence and therefore accretion, such as the gravitational instability (GI, Kratter & Lodato 2016; Béthune & Latter 2022) or the vertical shear instability (VSI, Nelson et al. 2013; Stoll & Kley 2014), but they rely on either very massive disks (GI) or on very fast cooling timescales (VSI) which are probably not always applicable to Class II objects (Lesur et al. 2022).

The diversity of increasingly resolved observations provides valuable constraints on disk accretion models. In particular, the turbulence-driven accretion model agrees with the inferred disk lifetimes and measured stellar accretion rates if α_v lies in the range 10⁻⁴−10⁻². Moreover, such hydrodynamic models are able to reproduce flows falling from the disk surface toward the disk midplane observed via ¹²CO emission at the radial locations of dark rings, as in Teague et al. (2019). Such fast meridional flows (≃ 0.1 c_s, with c_s being the local sound speed) are expected to occur in planet-induced gaps, and would result from the refill of material at the gap edges (see, e.g., Szulágyi et al. 2014; Morbidelli et al. 2014; Fung & Chiang 2016). However, other observational constraints tend to challenge this turbulence-driven accretion picture, when trying to estimate the level of turbulence in protoplanetary disks. Indeed, both the turbulent broadening of CO ro-vibrationnal line emissions (Flaherty et al. 2015, 2017) and dust settling measures from rings (Pinte et al. 2016) or edge-on disks (Villenave et al. 2020) all tend to suggest very low values of turbulence with α_v < 10⁻⁴−10⁻³ close to the midplane, a value much smaller than expected from full-blown MRI-driven turbulence.

Knowing that the gas mass should be efficiently transported from the outer parts of the disk toward the star and yet that the level of turbulence is (much) lower than expected has triggered the development of another paradigm of disk evolution, different from the viscous accretion model. In the MHD wind-driven accretion model, angular momentum is evacuated from the surface of the disk by a magnetized wind (Blandford & Payne 1982), which induces a radial laminar transport of the gas mass following the idea of Wardle & Koenigl (1993). This model reconciles at once the high stellar accretion rate, the a priori low turbulence, and also part of the disk dispersal in the latest stages of protoplanetary disks’ evolution (Bai & Stone 2013; Lesur et al. 2014). Both accretion models (viscous and MHD wind) are not mutually exclusive, and they both come from angular momentum conservation. Wind-driven accretion is an even more promising mechanism as it is supported by the detection of outflows coming from the inner parts of a protoplanetary disk that are probably too massive to be explained by thermal winds only (e.g., Louvet et al. 2018; de Valon et al. 2020). Readers can refer to Pascucci et al. (2023) for a comprehensive review of recent observations probing gaseous outflows.

An important challenge is to have a thorough understanding of the MHD wind-driven accretion mechanisms in order to be able to prescribe the impact of these winds on the gas in phenomenological models in 1D or 2D. It corresponds to the same process that made it possible to carry out 2D hydrodynamic simulations with an α_v prescription in viscous disks. For instance, Tabone et al. (2022) derived 1D equations modeling mass and momentum transport by a wind, and parameterized by an α_v-like dimensionless parameter, noted here as α_dw. Such a prescription is necessary in order to study long-term disk evolution or to perform systematic parameter space exploration, as highly resolved 3D MHD global simulations are time and energy consuming.

Recently, more and more studies have investigated disk and planet interactions as well as planetary migration in wind-driven accretion disks, using approaches similar to that of Tabone et al. (2022). For example, Kimmig et al. (2020) have explored how the wind-driven accretion process affects the migration of massive planets (Saturn-mass and Jupiter-mass planets) via 2D hydrodynamic simulations, with prescribed wind torques and ejection rates. They found, in particular, that planets can undergo episodes of runaway type-III-like outward migration (see, e.g., Masset & Papaloizou 2003; Pepliński et al. 2008a) if the wind ejection rate is sufficiently large. Concerning type-I planet migration with a wind-driven mass loss (Ogihara et al. 2015) and with inward gas flows induced by a wind torque at the midplane (Ogihara et al. 2017), studies have shown via simplified 1D models that type-I migration can also be slowed down and even reversed depending on the wind efficiency, with a positive unsaturated corotation torque that prevents the infall of super-Earths onto the central star. McNally et al. (2020) have focused on the migration of low-mass (6.7 Earth-mass) planets in a wind-driven accretion disk, via inviscid 3D hydrodynamic simulations, with a prescribed wind torque localized on the disk surface. They found a decoupling between the accretion layers and the passive planet-bearing midplane, leading to a negligible impact of the wind torque on the migration of the embedded planet. Lega et al. (2022) have conducted a similar study, but focusing on massive planets (M_j and M_s). They were able to identify two phases in the migration of giant planets, related to the formation and evolution of a vortex at the gap’s outer edge. Elbakyan et al. (2022) have examined via 1D and 2D hydrodynamic simulations the ability of fixed massive planets to open gaps in wind-driven accretion disks parameterized with a prescribed wind torque. They have reported that opening deep gaps is easier when angular momentum transport is dominated by magnetized winds, and therefore conclude that gap-opening planets may be less massive than currently believed. Finally, Aoyama & Bai (2023) have recently shown via 3D global non-ideal MHD simulations that magnetic flux tends to accumulate in the corotation region of massive planets, associated with nearly sonic accretion streamers and an enhanced angular momentum extraction from the gap region. They have shown that the gaps obtained in such MHD simulations are in essence deeper and wider than those obtained in viscous and inviscid simulations. Last, they found that the Lindblad torques are reduced, and the total gravitational torque exerted by the gas onto the planet is negative and fluctuates stochastically. In this paper, we revisit these processes in a more systematic manner and on longer timescales, varying the field strength and including low-mass planets which do not open gaps.

Most of these approaches rely on prescribed MHD wind torque (and in some cases also ejection rates) with very simple functional dependencies. However, it is well known that the wind torque and ejection rate strongly depend on disk magnetization, which is set by the disk surface density and the poloidal magnetic field strength (Lesur 2021). In turn, the surface density and the distribution of a poloidal field is expected to evolve as a result of the disk dynamics (Guilet & Ogilvie 2012, 2013, Guilet & Ogilvie 2014; Leung & Ogilvie 2019). Hence, the planet wake and gap necessarily affect the disk magnetization which in turn is expected to have an impact on the wind properties in a way that is not captured by phenomenological models. Our aim is to circumvent these difficulties by carrying out self-consistent 3D nonideal MHD simulations of a disk-planet-wind system where the magnetic field is free to evolve. More specifically, we want to focus on the opening of a gap by a fixed planet in a wind-driven disk, and reciprocally consider the accretion behavior in the vicinity of such planets. We also want to provide some leads to improve the 1D and 2D hydrodynamic models that use a wind torque prescription, as well as verify if planet gaps as formed with a wind can reproduce meridional flows.

The plan of this paper is the following. In Sect. 2, we describe the physical model and numerical methods of the nonideal MHD simulations. The results are then presented in Sect. 3, where we focus on the opening (Sect. 3.1) and evolution (Sects. 3.2.2 and 3.2.3) of the planet gap, as well as meridional flows (Sect. 3.2.1) and gravitational torques (Sect. 3.3), important for planet migration. A summary and a discussion follow in Sect. 4.

2 Numerical methods and setups

We performed MHD simulations using the IDEFIX¹ code (Lesur et al. 2023) that integrates the compressible MHD equations via a finite-volume method with a Godunov scheme. For this particular problem, IDEFIX was run on the Jean Zay cluster at IDRIS (France) on Nvidia V100 GPUs (graphics processing units) using a second order Runge–Kutta scheme, second order spatial reconstruction and the HLLD Riemann solver. The solenoidal condition was enforced using the constrained transport method and using an electromotive force reconstruction based on the contact wave upwinding strategy of Gardiner & Stone (2005). Diffusion terms (ambipolar diffusion and Ohmic resistivity) were integrated using the second order Runge-Kutta Legendre (RKL) scheme (Meyer et al. 2014). We note that resistivity is used here only as an inner boundary damping, as indicated in Sects. 2.4.1 and 2.4.3. The RKL scheme helps reduce the time to solution by at least a factor 2.5. We do not use the FARGO orbital advection scheme because the timestep is limited by the Alfvén velocity, not the azimuthal velocity of the gas at the inner radius.

2.1 Grid and resolution

We adopt a cylindrical coordinate system (R, ϕ, z), with R the radial cylindrical coordinate measured from the central star, ϕ the azimuthal angle, and z the vertical direction. We also make use of the spherical coordinate system (r, θ, ϕ), with r being the radial spherical coordinate and θ the colatitude.

In the azimuthal direction, the grid extends uniformly from 0 to 2π with 2048 cells. In the radial direction, we use two blocks with 512 evenly spaced cells from 0.3R_p to 1.9R_p, and an outer stretched grid with 128 cells between 1.9R_p and 6.0R_p where R_p is the (fixed) planet location. Concerning the vertical direction, the grid is refined around the midplane with 256 evenly-spaced cells between −0.4R_p and 0.4R_p, and two stretched grids symmetric around the midplane with 64 cells on both sides such that the total vertical extension spans between −6.0R_p and 6.0R_p. Around the planet location at R = R_p = 1, the radial, azimuthal and vertical resolutions correspond to 16 points per disk pressure scale height H(R) = hR. The total resolution is therefore N_R × N_ϕ × N_z = 640 × 2048 × 384. The aspect ratio h is considered constant, such that h = h₀ = 0.05.

2.2 Average conventions

For a quantity Q, we define the azimuthal average, the temporal average between T₁ and T₂, and the vertical integral between altitudes z₋ and z₊ respectively in Eqs. (1), (2) and (3): (1) (2) (3)

We make use of the notation 〈Q〉_t if we temporally average a quantity over the last 100 orbits of the simulation. We note the value of Q at the disk midplane. Finally, we define as (4)

which is the difference of Q between the disk’s upper layer at z_w and its lower layer at −z_w, with z_w ≃ 4H(R).

2.3 Planet

The planet is held on a fixed circular orbit at r_p = 1 in the disk midplane. Several planet-to-primary mass ratios q_p are chosen for the planet to explore the planet-disk-wind interaction: 3 × 10⁻³, 10⁻³, 3 × 10⁻⁴ and 3 × 10⁻⁵, which correspond respectively to M_p = 3M_j, M_j, M_s and 10M_⊕ for a solar-mass star. We make use of the planet’s Hill radius, defined as (5)

Accretion on the planet is neglected, as well as planet migration. The mass of the planet is gradually increased over its first N₀ = 10 orbits until it reaches q_p, according to , with t in units of Ω⁻¹. The total gravitational potential Φ considered in the simulation is the following: (6)

We note that Φ_★ is the term due to the central star: (7)

Φ_p is the gravitational potential of the planet: (8)

where ϵ is a softening length that avoids divergence at the planet location. ϵ is set to 0.1H(r_p), which is slightly larger than the characteristic size of a grid cell at rp. Finally, we take into account the indirect term Φ_ind arising from the acceleration of the star by the planet: (9)

All simulations have been run for at least 200 orbits, which is a good compromise between the computational cost and the characteristic timescales of disk/planet interactions.

2.4 Gas

The setup that we use for the gas is similar to the one presented in Riols et al. (2020); Martel & Lesur (2022). We summarize here some of the main features of this setup.

2.4.1 Nonideal MHD

We solve the nonideal MHD equations, for the density ρ, the velocity field υ and the magnetic field B: (10) (11) (12)

where Φ is the total gravitational potential defined in Eqs. (6–9), P the gas pressure, J = ∇ × B is the current density, e_b the unit vector parallel to the magnetic field line, η_A the ambipolar diffusivity and η_O the Ohmic diffusivity. We neglect the contribution of the Hall effect, and the Ohmic resistivity is included only as a damping process close to the inner radial boundary (see Sect. 2.4.3).

Regarding ambipolar diffusion, we define the dimensionless ambipolar Elsässer number (13)

with the Alfvén speed and Φ_K the Keplerian angular frequency. To keep the model as simple as possible, we follow Lesur (202l) and prescribe an ambipolar diffusion profile instead of computing a ionization model which would make the computation prohibitive. In the simple prescription that we use here, the midplane Elsässer number is constant and equal to 1, similarly to full radiative models (e.g., Thi et al. 2019). At higher altitudes, A_m increases abruptly to mimic the ionization due to X-rays and far UVs at the disk surface. The diffusivity profile we use is (14)

with ϵ_id = 6.0 in all our simulations, which quantifies the thickness of the nonideal layer. We cap A_m in order to save some computational time. Note that we would actually expect the ionization by X-rays from the central star to increase when the gas is depleted in a planet gap, leading to a stronger coupling with magnetic fields, as reported in Kim & Turner (2020). As we prescribe A_m, its dependence on Σ is not taken into account here. For such , we expect the MRI to emerge at low amplitude and generate some weak turbulence (see Cui & Bai 2022, and Appendix B.2).

The disk is initially threaded by a large-scale vertical magnetic field given by the profile (15)

with β₀ the initial value of the plasma parameter β, defined as the ratio of the thermal pressure P over the magnetic pressure: (16)

Two β₀ are chosen in order to explore the influence of the magnetization on the interaction between the planet and the disk: 10⁴ and 10³.

Our numerical parameters are chosen to be close to more complete thermo-chemical models. For instance, a constant matches the thermo-chemical model of Thi et al. (2019) with in the midplane (see, e.g., the green region in the middle left panel of their Fig. 8). Moreover, Ohmic diffusion is negligible for R > 10 au as these models have Elsässer number ≫ 1 (see the top left panel of Thi et al. 2019 Fig. 8). Hence, our numerical models with ambipolar diffusion only closely resemble thermo-chemical models in the 10–100 au range that corresponds also to the regions that are typically unveiled by radio interferometry with ALMA. In particular, sequences of dark and bright rings along with other substructures are commonly detected in such regions at millimeter wavelengths, in the continuum emission and CO emission lines (see, e.g., the DSHARP Andrews et al. 2018 and MAPS Öberg et al. 2021 programs).

2.4.2 Initial disk properties

The initial gas volume density arises from the hydrostatic disk equilibrium: (17)

with the sound speed at disk midplane and Σ the surface density following (18)

which is coherent with some spectral lines observations (see Williams & McPartland 2016, for the disk HD 163296). We neglect self-gravity throughout this study, which is a priori valid for Class II objects.

The simulation domain can be divided into a cold dense disk at temperature and a low-density hot corona with temperature to account for the inefficient gas cooling via thermal accommodation on dust particles in the low-density hot corona (see Thi et al. 2019, e.g., their Appendix F). We model the transition between these two regions with the following smooth vertical profile: (19)

We use a fast ℬ-cooling prescription, with ℬ = 0.1Ω⁻¹, in order to quickly relax the temperature T(R, z) toward the prescribed profile T_eff(R, z). With such finite but short cooling timescale, we expect to be on the verge of the VSI-unstable regime (Manger et al. 2021). We discuss whether the VSI is present in our simulation in Appendix B.2. The pressure and density are related by P(R, z) = ρ(R, z)T(R, z). The cooling prescription is used to solve the energy equation, assuming an ideal equation of state for the gas, with an internal energy density , and an adiabatic index γ = 5/3. Regarding the initial gas azimuthal velocity, we use the following sub-Keplerian velocity profile (20)

with the Keplerian speed. Note that the above expression is not an equilibrium in the hot corona because of the larger temperature in this region. This is not a problem since the MHD wind, which is neither a dynamical equilibrium, is launched very rapidly once the simulation has started and clear the initial condition in the corona.

2.4.3 Boundary and internal conditions

At the inner radius, we copy the gas density and pressure, as well as B_ϕ and B_z from the innermost active zone to the ghost zone. We impose v_z = 0 and a Keplerian velocity for υ_ϕ. The boundary condition is reflecting (symmetric) for υ_R, such that the fluxes are null at the interface with the ghost zone. At the outer radius, we apply power-law expressions for all quantities to ghost zones, with radial power-law exponents of −1.5, −2.5, −0.5, −0.5, −0.5, −5/4 and −5/4 for ρ, P, υ_R, υ_ϕ, υ_z, B_ϕ and B_z, respectively to mimic self-similar models (Lesur 2021). The boundary conditions are periodic in the azimuthal direction, and we apply outflow conditions for all the quantities in the vertical direction.

To avoid too small time steps, we limit the Alfvén speed in the simulation domain to the value of the Keplerian speed at the inner radial boundary of the disk R_in = 0.3. This assumption leads to an overestimation of the density when the criterion υ_A > υ_K(R_in) is encountered, typically in the disk corona. We also add a threshold in density at ρ_MIN = 10⁻¹². Because the planet is kept fixed at R_p = 1 and far from the outer radial boundary of the disk R_out = 6, we do not apply an outer damping zone. However, when R < 1. 1R_in, we prescribe a simple damping zone for the radial and vertical components of the momentum (ρυ_R, ρυ_z). In this disk innermost region, we also add a band of resistivity such that η_O = H²(R_in)Ω_K(R_in) to avoid artificial magnetic coupling between the inner boundary condition and the computational domain.

2.5 Numerical protocol and parameter study

Before introducing the planet, we run two axisymmetric simulations for 100 orbits at R = 1, for two values of the initial midplane β₀ parameter. This allows the system to launch a MHD-driven wind and reach a quasi-steady-state at R = 1. The final snapshot of these axisymmetric planet-free simulations is then extended in the azimuthal direction and used as an initial conditions for the full 3D simulations. This instant marks the t = 0 of the simulations presented in this paper. We then introduce the planet as described in Sect. 2.3, which evolves thereafter during N_orb orbital periods. In practice, we performed several 3D simulations varying the planet mass M_p and the initial β₀-plasma parameter. In Table 1, we detail these parameters that we used in eight 3D MHD simulations. The names and duration of the various runs are indicated. We also show additional simulations that were used for initialization or comparison purpose, in particular the axisymmetric 2D simulations, a 3D planet-free simulation with β₀ = 10³ and a purely hydrodynamic inviscid run. Concerning the computational cost, 475 000 GPU hours have been granted for the simulations presented here, which corresponds to an estimated energy consumption of about 200 MWh, equivalent of an emission on the order of 20 in tons of CO₂.

2.6 Disk accretion and wind launching

We define here some important parameters that are classically used to describe disk accretion and wind launching. We first define the accretion rate as (21)

Using the equation of angular momentum conservation, Ṁ_acc is linked to the transport of angular momentum by radial and surface stresses, respectively noted 𝒲_Rϕ and 𝒲_zϕ (22)

with 𝒲_Xϕ = ρυ_X(υ_ϕ − υ_K) − B_XB_ϕ the ϕ component of the stress, and is the torque per unit of volume exerted by the planet on the gas, with Φ_p defined in Eq.(8). In the r.h.s. of Eq. (22), we have therefore four terms that can contribute to the total planet-free torque 𝒯, and therefore to Ṁ_acc: the radial Reynolds (ℛ_R) and Maxwell (ℳ_R) torques, and the surface Reynolds (ℛ_surf) and Maxwell (ℳ_surf) torques. Figure 1 shows the contribution of the absolute value of these four torques in the accretion, as a function of the distance to the star, for the case Mj-β₃, that is for the largest magnetization and the Jupiter-mass planet. The dominant term is M_surf (yellow curve), which means that the vertical extraction of angular momentum at the disk surface seems to be the main driver of accretion, in particular when β₀ = 10³. On the contrary, the surface Reynolds torque ℛ_surf (green curve) has little impact on the accretion mechanism. Concerning the radial torques (ℛ_R in blue and ℳ_R in red), they have a similar influence on the accretion, and can be comparable to ℳ_surf if combined. In summary, ℳ_R ≃ ℛ_R, ℳ_surf ≫ ℛ_surf and ℳ_surf + ℛ_surf ≳ ℳ_R + ℛ_R in general, except for strong negative gradients in the gas density, for example near the gaps’ inner edge. Although there is a variability in the radial profile of the four torques, the torques per unit of surface density are relatively constant radially on average. Due to the importance of the ℳ_surf term, we define the dimensionless parameter u as in Lesur (2021): (23)

which, when positive, gives an idea of the vertical evacuation of angular momentum from the disk by the magnetic torque at the disk surface. In practice with this definition, υ is positive when the wind extracts angular momentum from the disk. We note that for β₀ = 10⁴, that is at lower magnetization, assuming that ℛ_surf is still negligible, ℛ_R is actually the dominant term in the internal parts of the disk (R < 3, where ℛ_R/ℳ_surf < 10), whereas ℳ_surf takes over in the external parts of the disk. It is safe to neglect ℛ_surf as it is proportional to the density, which is small at the disk surface z = 4H We also define the ejection efficiency dimensionless parameter ξ as: (24)

which quantifies the fraction of the accreted material that is ejected. Another important quantity related to the radial transport efficiency is the α_v (Shakura & Sunyaev 1973) parameter: (25)

where and are respectively related to the Reynolds and Maxwell components of the radial stress. We also define a dimensionless parameter α_dw, similar to the one introduced in Tabone et al. (2022), and which we can divide in Reynolds () and Maxwell () contributions: (26)

We note that , using Eq. (23), and .

In Appendix B.2, we address the radial transport efficiency and the level of turbulence in the 3D planet-free run 3D-β₃ temporally averaged between 50 and 100 orbits. In particular, we compare α_v but also the turbulent component of with their corresponding values in the run 10Me-β₃. We show that the total radial transport of angular momentum is barely affected by the presence of a low-mass planet, except close to its coorbital region. We also show that the level of turbulence is similar in both cases near the coorbital region, but 2–3 times lower elsewhere for the low-mass planet case. This may be due to the presence of nonaxisymmetric spiral features that weaken the turbulence, as in Ziampras et al. (2023), but is beyond the scope of this study.

Fig. 1 Radial profile of the four torques ℛR (blue), ℳR (red), ℛsurf (green), ℳsurf (yellow) and their sum 𝒯 (black) normalized by , azimuthally averaged and temporally averaged over the last 100 orbits for the run Mj-β3. We consider here only the logarithm of the absolute value of these quantities.

3 Results

3.1 Gap opening

3.1.1 Overview

We consider the opening of a planet gap in a wind-launching disk. The disk surface density structures after 200 planet orbits are shown in the cartesian maps of Fig. 2, for the different planet masses (four rows), and β₀ parameters (two columns). The colormap indicates the logarithm of the gas surface density, with purple for low-density regions, and red for high-density regions. In these maps, we detect the presence of an annular planet-induced gap, whenever q_p > 3 × 10⁻⁴ (M_s). In terms of parameter, this is coherent with the classical condition for nonlinear effects in the flow whenever , although this condition also depends on the level of turbulence α_v (see the gap-opening criterion in Crida et al. 2006, henceforth CR6) and should be less true for highly magnetized disks (see the extended gap-opening criterion in Elbakyan et al. 2022). The gap is wider and deeper when the planet mass is higher, with at least a factor 10 between the 3 M_j and the M_s cases. For a given M_p > M_s, the gap is denser when β₀ is larger, with a factor on the order of 5. For all β₀, when M_p > M_j, the horseshoe region appears clearly asymmetrical, with a difference in density in front of and behind the planet in azimuth (see panel a. of the sketch presented in Fig. 3). As in Baruteau et al. (2011), this kind of asymmetry in the horseshoe region could generate a positive or negative contribution in the total gravitational torque exerted by the gas onto the planet due to a static corotation torque (see Sect. 3.3). For massive planets (M_p > M_j) and highly magnetized disk (β₀ = 10³), in addition to this azimuthal asymmetry, there is also a two-fold radial asymmetry, with the outer gap wider (see panel b of Fig. 3) and deeper (see panel c of Fig. 3) than the inner gap. In practice, the planet is therefore closer to the gap’s inner edge than the gap’s outer edge. After 400 orbits, for the case Mj-β₃, the inner gap can be three times narrower and three times denser than the outer gap. The sketch in Fig. 3 isolates these three asymmetries, labeled as azimuthal asymmetry, radial width asymmetry and radial depth asymmetry.

We are able to retrieve the inner and outer planet wakes in all cases, although the fixed colorbar that we choose does not allow to disentangle for the lowest planet mass case the wakes-induced perturbations from perturbations independent of the planet presence. For all planet masses and disk magnetizations, spontaneous zonal flows emerge in the simulations, associated with concentric density rings, following a “self-organization” process already described in a number of planet-free simulations (Béthune et al. 2017; Suriano et al. 2017, 2018, 2019; Riols et al. 2020; Cui & Bai 2021). Comparing the planet-free case 3D-β₃ and 10Me-β₃ in Appendix B.1 confirms that the self-organized density rings and gaps are present at high magnetization regardless of the presence of a planet. Nevertheless, we show that despite its low mass, a 10 M_⊕ planet is still able to impact its environment, especially regarding the distribution of density, the accumulation of poloidal magnetic field and the strength of the wind torque. The amplitude of gas surface density perturbations (Σ − Σ_th)/Σ_th are larger for smaller β₀. This is particularly visible in the top row of Fig. 2, with ~ 80% in amplitude for the case 10Me-β₃ and ~ 20% in amplitude for the case 10Me-β₄, compared to a simple theoretical power-law profile in Σ_th = Σ₀R^−0.5. These gaseous rings near R = 1 tend to be smoothed out when the planet is able to generate sufficient perturbations. The competition between self-organization rings and density gaps induced by massive planets can be seen in the right column of Fig. 2, with the low-density ring near the planet location and the high-density ring just outside the planet orbit. For the low-mass planet case (10Me-β₃), the gas perturbations are dominated by self-organization. For the intermediate mass case (Ms-β₃), the competition between self-organization and the planet tends to diminish the contrast in density. For larger planet masses (Mj-β₃, 3Mj-β₃), the strong gas perturbations near the planet location are dominated by planet-driven flows. Along with these density perturbations, we detect multiple pressure bumps that are shaped at the same time by the planet and self-organization. For example, in the run Mj-β₃, we obtain four pressure maxima outside the planet location after 400 orbits. The three outermost pressure maxima are probably linked to self-organization, as we retrieve them for all planet masses, and they are less pronounced at lower magnetization.

When the planet is sufficiently massive to carve a deep annular gap around its orbit, it will form a pressure maximum at the outer edge of its gap due to the deposition of the angular momentum flux carried away by the planet’s outer wake (Baruteau et al. 2014; Bae et al. 2016). Dust traps, and therefore dust rings, can be naturally obtained at pressure bumps (see, e.g., Pérez et al. 2019; Wafflard-Fernandez & Baruteau 2020), which is of prime interest as we commonly detect annular substructures with ALMA in radio (see, e.g., Huang et al. 2018; Guzmán et al. 2018; Pérez et al. 2019) and SPHERE in near-IR (see, e.g., Avenhaus et al. 2018). Post-processing our simulations with dust radiative transfer calculations, we would expect multiple bright and dark rings of emission, either due to the planet’s outer pressure bump (as in Nazari et al. 2019), or due to self-organization (as in Riols et al. 2020). Note that the emission rings induced by self-organization are expected to be deeper and more extended in 3D than in the 2D axisymmetric models of Riols et al. (2020). In any case, we have here two different mechanisms able to generate multiple annular substructures all at once, which could be valuable for comparison purpose (in particular concerning the kinematic signatures of such gaps in the gas).

We note that we also retrieve the formation of several small-scale vortices at the gap’s outer edge corresponding to vortensity minima, probably arising because of the Rossby-Wave Instability (RWI, Lovelace et al. 1999). Such vortices eventually merge in one large-scale vortex. This is particularly visible early in the run 3Mj-β₄. Then, the vortex gradually disappears and becomes axisymmetric. Such vortex is actually composed of numerous small-scale structures that resemble the ones obtained with elliptical instability in Lesur & Papaloizou (2009). Large-scale vortices are mostly transient in our nonideal MHD simulations, but some nonaxisymmetric structures can also appear quite late until the end of the simulation, for example at R ≃ 1.3 for the run 10Me-β₄.

In Fig. 2, there is evidence that an evacuated region is generated by the inner boundary condition, especially when β = 10³. This is probably due to the fact that magnetic field lines do not cross the inner edge of the grid, and therefore tend to efficiently accumulate in that region. This accumulation of magnetic field is accompanied by a depletion of matter (see Sect. 3.2.2), and tends to drift outward with the same mechanism that makes inner cavities (as presented in Martel & Lesur 2022) and planet gaps (as presented here, see panel b. of Fig. 3 and Sect. 3.2.3) drift outward. Although the mechanisms behind the radial asymmetries of the planet gap and the inner region are similar in their evolution, their origin is different, as they result respectively from the carving of a gap and an artifact in the inner boundary.

Fig. 2 Gas surface density maps for the eight MHD runs, after 200 planet orbits. The four rows show the different values of the planet mass. The two columns show the different values of the initial β0 parameter. The white circle marks the planet location.

Fig. 3 Sketch of the different asymmetries of the planet gap, for massive planets and highly magnetized disks. The yellow and black circles represent respectively the central star and the planet. The horseshoe region is represented in shades of gray, with high-density zones in light gray and low-density zones in dark gray. The dashed line corresponds to the planet orbit. The global shape of the gap results from the combination of three asymmetries: (a) Azimuthal asymmetry of the horseshoe region. (b) Radial width asymmetry, with the outer gap wider than the inner gap. (c) Radial depth asymmetry, with the outer gap deeper than the inner gap.

3.1.2 Temporal variability

Figure 2 gives a good idea of the gas surface density in the disk near the planet at a given time. In order to apprehend the temporal evolution of the gas distribution, we compute in Fig. 4 the space-time diagram of the logarithm of the azimuthally averaged 2 for the eight nonideal MHD runs. We also show the space-time diagram of the inviscid hydrodynamic simulation (3Mj-α₀) for comparison purpose. These diagrams confirm that a planet gap forms when M_p > M_s, and that the global shape in depth and width of this gap depends on the planet mass and the initial magnetization. For the low-mass planet cases, we notice the formation and evolution of self-organized rings and gaps, visible for example in the case 10Me-β₃ with the under-density that appears near R = 1.2.

The planet gap is denser when the initial magnetization increases (β₀ decreases), which is particularly visible when M_p = 3M_j (bottom group of Fig. 4). A similar trend was found in simulated magnetized cavities of transition disks (Martel & Lesur 2022, their Fig. 22), and can be explained mostly with the interdependence between the radial profiles of the accretion rate Ṁ_acc, β and Σ, by comparing their values in the outer disk (e.g.,β_disk ≃ β₀) and their counterpart in the outer gap (e.g.,β_gap). Firstly, Ṁ_acc is known to be a function of Σβ^−b (0 < b < 1, with in particular b = 0.78 in Lesur 2021). Secondly, when a gap forms, the magnetization in the gap self-regulates with β_gap ~ 10 (see bottom panel of Fig. C.1). As aresult, in the gap, Ṁ_acc ∝ Σ_gap. If we make the assumption of a quasi steady-state, which implies that Ṁ_acc in the gap matches Ṁ_acc in the disk, then .

Along with the increase in density, we detect episodes of activity in the gap when increasing the magnetization and for M_p > M_j. In particular, we detect episodes of partial filling of the gap due to bursts of accretion in the gap from the outer disk. We come back to this accretion behavior in Sect. 3.2.3. We also recover the radial asymmetries as suggested in Fig. 2 for the cases 3Mj-β₃ and Mj-β₃, whether it be in depth or in width. This is particularly visible for the Jupiter-mass planet after 400 orbits, with the outer gap slowly spreading outward. We discuss this gap dynamics in Sect. 3.2.3.

Looking again at the bottom group of Fig. 4 (e.g., for a 3 Jupiter mass planet), we evaluated how the gap density evolves with time, following the procedure of Fung et al. (2014) and Fung & Chiang (2016, see Appendix C for more details). Using this metrics, we notice that the gap reaches its minimum density earlier when the magnetization is higher, but in that case the gap-opening speed is also larger at least early on in the simulations. Qualitatively, this suggests that the accretion, and therefore mainly the wind torque (via the dominant term ℳ_surf), tends to help the gap opening until β in the gap reaches a threshold value (on the order of 1–10) which in turn makes the density reach a quasi-equilibrium value. When the initial magnetization is larger, the wind torque is also larger, helping even more the gap to be carved and to reach earlier its quasi-equilibrium state.

Finally, we observe a strong time-variability in the gap density in magnetized models, with a larger amplitude when the magnetization increases. This implies that the accretion varies in time, with stronger episodes as the magnetization increases (see snapshots in Fig. A.1 that illustrate such temporal variability).

3.1.3 Gap-opening criterion in a wind-launching disk

The goal of this paragraph is to obtain qualitatively a criterion for gap opening in a disk where accretion is mostly driven by MHD winds via vertical extraction of angular momentum. We start by stating that a gap may form if the time it takes for a parcel of gas being accreted at velocity υ_acc to cross the Hill sphere of the planet τ_acc = R_hill/υ_acc is longer than the time it takes for the planet to clear up the Hill sphere τ_p ≃ R_p/(R_hillΩ). This last term is obtained by saying that the speed at which a fluid element is evacuated from the Hill sphere by the planet is R_hill therefore the typical timescale to clean a full ring at the planet location R_p is approximately given by τ_p. Said differently, the characteristic emptying timescale of the Hill sphere is linked to half the libration time, which is given by as in Masset (2008, see their Sect. 4.2.1), where we replaced the half-width of the horseshoe region (used to estimate the characteristic emptying timescale of that region) by the Hill radius. scales directly with τ_p, hence the relevance of this quantity.

With these two timescales, the criterion for gap opening reads (27)

Using the equation right before Eq. (17) in Lesur (2021), we can link the accretion speed υ_acc to the normalized magnetic wind torque u. We get υ_acc = 2R_pΩhυ. The gap-opening criterion then becomes (28)

Using the expression for the Hill radius in Eq. (5), we eventually get (29)

This should also be combined to another criterion which states that the disk thickness should be smaller than the Hill radius (30)

Using Eqs. (29) and (30), we can eventually get a combined criterion for gap opening as in CR6, when accretion is dominated by MHD winds (31)

Note that we have kept the factor 3/4 in the original CR6 formulation for the disk thickness criterion. When accretion is dominated by viscosity, we can replace τ_acc by , with C a constant factor and rewrite the gap-opening criterion of Eq. (31). By considering the Reynolds number, Eq. (29) would become 3C/(q_pℛ) < 1 in the fully viscous case. We therefore retrieve the same scaling ∝ 1/q_pℛ as in CR6 criterion (see their Eq. (15)). Note that their factor 50 comes from a fitting procedure. Coming back to the MHD wind-driven accretion disk, we can estimate the critical to open a gap in a MHD wind-driven accretion disk, with h > 0 and υ > 0: (32)

In order to make this expression useful, one must estimate the wind torque coefficient υ. Here, we propose two approaches: first using the self-similar solutions from Lesur (2021) with , giving . In this case our gap-opening criterion becomes essentially a function of β₀ and h. The second approach is to measure υ in our simulations with a 10M_⊕ planet, in the radial interval [, ] and assuming that the torque obtained in this case is only barely affected by the planet.

We plot in Fig. 5 the gap-opening criterion obtained from Eq. (31) using these two approaches and assuming an aspect ratio h = h₀ = 0.05. The hatched area above the dashed line corresponds to Eq. (31) with the self-similar scaling for υ (Lesur 2021) while the pink box plots indicate the values critical estimated from υ measured in our simulations. The 8 colored dots correspond to the 8 (q_p, β₀) pairs in our nonideal MHD simulations, with the 3M_j, M_j, M_s and 10M_⊕ cases respectively in blue, red, yellow and green. We find that our expression correctly predicts gap formation where we find them (see Fig. 2). For a given planet mass, this plot confirms that it is harder for a planet to carve a gap at higher initial magnetization (decreasing β₀). More precisely, the minimum mass for a planet to open a gap lies between a mass of Saturn and a mass of Jupiter for β₀ = 10³, whereas it is inferior to a mass of Saturn for β₀ = 10⁴.

For this estimation of a gap-opening criterion in a wind-launching disk, it should be mentioned that we did not consider the self-organized structures, and neglected the mechanism of magnetic accumulation in the gaps. Because β decreases also in the gap, u tends to increase as described by the self-similar scaling law, and the carving of the gap will be easier. We mentioned this runaway gap opening in the last paragraph of Appendix C (see, e.g., the dashed and solid blue lines in the top panel of Fig. C.1).

Fig. 4 Space-time diagrams showing the evolution of the logarithm of the gas surface density Σ for the nine runs divided in four groups of different planet masses. The top plot of each group has the lowest magnetization, whereas the bottom plot has the highest. The radial extent of all the diagrams range from 0.5 to 1.5 code units, with the planet at R = 1. We pinpoint as ℰ1 and ℰ2 two episodes in the run Mj-β3 that will be discussed in Sect. 3.2.3.

Fig. 5 Gap-opening criterion with MHD winds displayed in the (qp,β0) domain. The dashed line indicates the minimum mass to open a gap from this criterion, using a scaling law between u and β0. Planets in the hatched area are able to open a gap. To obtain the pink box plots, we estimated υ in our low-mass planet simulations.

3.2 Case study: high planet mass and high magnetization runs

We focus now mainly on the simulations with M_p > M_j and β₀ = 10³.

3.2.1 Meridional flows: Accretion layers and sonic streamers

We discuss here the behavior of the flow in the vicinity of a planet in a fixed circular orbit. In the next paragraphs, we focus on the mass flux instead of the velocity, as mass flux brings a more precise view of the gas bulk motion. In Fig. 6, we show in background color the magnitude of the poloidal mass flux 〈ρυ_p〉_ϕ,t in log scale, averaged azimuthally and temporally over the last 100 orbits for the run 3Mj-β₄. The texture in LIC (line integral convolution) and the white arrows indicate the direction of the radial (〈ρυ_R〉_ϕ,t) and vertical (〈ρυ_z〉_ϕ,t) components of this mass flux. Following the nomenclature in Fung & Chiang (2016), the black arrows help to schematically classify the different flows in the (R–z) plane in three categories.

• Planet-driven flows (dashed black arrows) are probably linked to the planet’s repulsive Lindblad torques.

• Upper layer accretion flows (solid black arrows) are driven by the accretion of material from the surface layers of the outer disk to the planet poles and from the inner gap to the surface layers of the inner disk.

• Wind-driven flows (dotted black and white arrows) act to carry material away from the lower altitudes of the disk.

These flows eventually collide, driving the merged flow upward and resulting in a large-scale meridional circulation in the outer disk (as in Fung & Chiang 2016). We also observe a small-scale localized meridional circulation close to the planet, with bigger loops in the inner gap when R ∈ [0.9, 1] and |z| ∈ [0.0, 0.1] than in the outer gap when R ∈ [1, 1.1] and |z| ∈ [0.0, 0.02]. This asymmetry is due to the two accretion layers in R ∈ [1.0, 1.3] and |z| < 0.1. On the one hand they carry material from the outer disk onto the poles of the planet, and on the other hand they pinch and confine the localized outer meridional circulation even closer to the planet by counteracting the planet-driven flows. If we increase the magnetization, the mass flux is globally increased in the accretion layers compared to the case β₀ = 10⁴. This behavior is the result of an enhanced ℳ_surf component (magnetic torque at the disk’s surface) and to a lesser extent to an enhanced radial torque (ℛ_R + ℳ_R) leading overall to a more efficient accretion. In addition, we notice that the accretion layers are more collimated toward the midplane, thus carrying denser material compared to the lower magnetization case. This behavior prevents even more the planet-driven flows to counterbalance accretion. In that high magnetization case, accretion onto the planet – as well as planet growth – should therefore be more efficient because of this denser inflow coming from lower altitudes toward the planet poles, as observed in Gressel et al. (2013).

When M_p > M_j and β₀ = 10³, in addition to the structures observed in Fig. 6, we note the presence of a streamer at z < 0 which accretes even more material toward the inner disk. In Fig. 7, we show that this top-down nonsymmetrical flow (for z ∈ [−0.1, 0]) has nearly sonic velocities (on the order of 60% of the sound speed) in the gap. Several other studies have pointed out that sonic accretion was expected in strongly magnetized disks, as is the case in our planet gaps (Combet & Ferreira 2008; Martel & Lesur 2022). This is also coherent with the fact that the mass flux 〈ρυ_p〉_ϕ,t is approximately radially constant (see Fig. 6), whether it be in the outer disk, the planet gap or the inner disk, in spite of a deep depletion of gas in the gap. Because the gap density is low, the velocity has necessarily to increase to keep a constant mass flux, reaching a nearly sonic speed in the gap and a β much smaller than β₀ (reaching locally values close to unity, see e.g., the bottom panel in Fig. C.1). Note that this approach is time-averaged and the large-scale gas motion is actually extremely variable in time, with complex redistribution of matter in the gap and evacuation in the wind (see Appendix A).

In conclusion, we retrieve in our MHD simulations a large-scale meridional circulation of the gas, as in purely hydrody-namic 3D simulations (Szulágyi et al. 2014, 2022; Morbidelli et al. 2014; Fung & Chiang 2016; Teague et al. 2019), though less clearly symmetrical between outer disk region and inner disk region due to the different dynamical processes at stake. The gas in the outer disk flows from the upper accretion layers into the gap and eventually falls toward the midplane. The planet is able to maintain a sustainable gap via its tidal torque, which tends to empty the planet’s coorbital region. Therefore, part of the accretion flow is advected inward in the inner gap, whereas the rest is pushed back in the outer disk. At the gap’s outer edge, the gas still at low altitude is then progressively driven upward to reach again the surface accretion layers. Note that this meridional circulation is not a closed loop. On the one hand, the accretion occurs from the disk’s surface layers, and inevitably transports gas material inward, regardless of the planet presence. On the other hand, part of the gas is evacuated vertically in the wind at the gap’s outer edge during its descent toward the midplane. We see in the next sections that this gas dynamics in the gap comes from the distribution and dynamics of magnetic field lines.

Fig. 6 Meridional flows around the planet for the run 3Mj-β4, averaged in azimuth and over the last 100 orbits. The background color represents the logarithm of the poloidal mass flux. The white streamlines and the LIC correspond to the radial and vertical components of the poloidal mass flux. Dashed black arrows, solid black arrows and dotted black and white arrows show respectively the planet-driven flows, the upper-layer accretion flows and the wind-driven flows. The upper and lower surfaces ±zw ≃ ±4H(R) are indicated by the black dotted lines. The planet is indicated by the white circle.

3.2.2 Magnetic field

We focus here on the large-scale magnetic field threading the disk and its long term behavior. We consider first the temporal evolution of the vertical component of the azimuthally averaged magnetic field in the midplane . Figure 8 shows the space-time diagram of this quantity for the runs Mj-β₄ (top panel) and Mj-β₃ (bottom panel), to be compared with the corresponding space-time diagrams of log〈Σ〉_ϕ in Fig. 4. In both cases, during the first orbits, several small-scale structures are visible, each one corresponding to the self-organization of vertical magnetic field (Riols & Lesur 2019) in the axisymmetric 2.5D simulation used for the initial condition. In 3D, they tend to quickly merge into more extended and diffuse structures, see for example at R ≃ 2 in both cases. We note that this large zonal field is associated to a dark ring of matter (i.e., a gap induced by self-organization). The anticorrelation between density and magnetic fluctuations in gaps formed in planet-free disks have been reported in various works (see, e.g., Béthune et al. 2017; Suriano et al. 2018; Riols & Lesur 2019).

At high magnetization (bottom panel of Fig. 8), as the planet is growing, the gap is formed and reaches its minimum density (on average) after ≃50 planet orbits, as indicated in Sect. 3.1.2. After the planet gap has opened, it becomes a prime area of magnetic field accumulation between 75 and 200 orbits, with reaching its highest value. If we focus on log〈Σ〉_ϕ, this strong magnetic accumulation is concomitant to a strong episode of gas depletion and a fast outward drift of the outer gap (episode noted ℰ₁ in Sect. 3.2.3). This accumulation at the gap center has been observed for low-mass planets in various ideal (Papaloizou et al. 2004; Zhu et al. 2013; Carballido et al. 2017) and non-ideal (Keith & Wardle 2015) MHD simulations. Subsequently, after 200 orbits, the magnetic field is divided in two regions of accumulation, one at the transition between the inner disk and the inner gap (near R ≃ 0.85), and the other one at the transition between the outer gap and the outer disk (near R ≃ 1.3). Carballido et al. (2017) obtained as well for their highest planet masses (a few M_j) a double accumulation for the vertical magnetic field (as well as α_v), larger close to the gap edges than at the gap center. In addition, they argued that the radial band of large is smaller than the band of low-density (i.e., the density gap), which is also what we retrieve in our simulations: ≃0.45 code units for at 400 orbits instead of ≃0.55 code units for 〈Σ〉_ϕ, corresponding to a difference of 1 au for a planet orbiting at 10 au. Concerning the temporal variability of magnetic field accumulation, the two accumulation regions in Fig. 8 do not exhibit a constant in time, with episodes of strong magnetic field accumulation (e.g., between 310 and 350 orbits) alternating with episodes of weaker accumulation (e.g., between 285 and 305 orbits). This behavior may be coupled to the sporadic accretion of matter illustrated in Appendix A.

At lower magnetization (top panel of Fig. 8), the behavior seems a little bit more complex. We first witness a phase of accumulation at the gap center (< 100 orbits), though less efficient than in the case of stronger magnetization (see maximum value of the colorbars between the two panels). This episode is then followed by an oscillation of the magnetic field accumulation between the inner gap edge and the outer gap edge, which could be due to a complex behavior of matter in the planet’s horseshoe region.

Coming back to the two-fold accumulation of the magnetic field, we represent in Fig. 9 the distribution and topology of the time-average poloidal magnetic field and the disk surface density for the run Mj-β₃. We retrieve the double accumulation of the poloidal magnetic field at the transitions between the inner disk and the inner gap and between the outer gap and the outer disk, as presented in Carballido et al. (2017). The magnetic field lines are nearly vertical in the outer disk. We however detect the presence of kinks in the 〈B_p〉_ϕ,300–400 lines, as described in Lesur (2021); Martel & Lesur (2022). Such structures trace the two disk’s upper accretion layers that bend magnetic field lines as material is carried inward in the gap. If we now focus on the gap, the magnetic field lines are pinched in a plane parallel to the midplane but such that z < 0. This localization actually matches the sonic accretion streamer, indicating here again that poloidal magnetic fields lines are bent by the streamer.

Magnetic field accumulation in the gap is so strong that it could be responsible for sporadic protoplanetary jets. We demonstrate this in Fig. 10 with an azimuthal cut of log (Σ) near the location of the planet, at ≃114 orbits for the run 3Mj-β₄. This zoom near the planet shows several interesting structures: the gap, the density maxima on both gap edges, a circumplanetary disk (CPD) and an outflow from the planet poles. We can characterize the CPD by fitting a power law to the azimuthal velocity υ_ϕ − υ_K in the two regions close to R = 1 in the Hill sphere. We obtain 0.038(R − 1)^−0.459 and 0.034(R − 1)^−0.453 respectively for the outer CPD and the inner CPD, while with the planet mass we have , confirming that the CPD is rotating at the expected Keplerian velocity around the planet. Turning now to the meridional flows, we find that, in addition to the circulation presented in Sect. 3.2.1, we get an outflow from the planet poles (r ≃ 1, |z| < 0.25). At higher altitudes, the outflow is slightly bent in the disk wind direction where the outflow material eventually falls back toward the planet. We note however that this outflow is only occasionally present, as it is not visible in time-average plots. Outflows from CPDs have already been observed in 3D nonisothermal MHD simulations (Gressel et al. 2013) where it was suggested that they could be magnetocentrifugally driven. For the sake of conciseness, we however postpone the detailed investigation of these outflows to a future study.

Overall, these results indicate that the poloidal magnetic field is very dynamical and strongly affected by the planet, with a significant accumulation in the gap. This in turn impacts the strength of magnetic torques and of the accretion flow, which we explore next.

Fig. 7 2D map of the Mach number, defined as , in the (R-z) plane for the run Mj-β3, averaged in azimuth and over the last 100 orbits. The background color represents the logarithm of the Mach number. The black streamlines and the LIC correspond to the radial and vertical components of the poloidal Mach number. The nearly sonic top-down nonsymmetrical streamer is visible in white, in the gap and under the planet indicated by the white circle.

Fig. 8 Space-time diagram of the vertical component of the magnetic field in the midplane , for the runs Mj-β4 (top panel) and Mj-β3 (bottom panel). The horizontal dashed lines indicate the location of the planet. The two episodes ℰ1 and ℰ2 are indicated in white.

3.2.3 Accretion rate, wind torque, and gap asymmetries

In this section, we focus on the case Mj-β₃. We look at the link between the accretion rate Ṁ_acc defined in Eq. (21), the specific wind torques via ℳ_surf and the υ parameter defined in Eq. (23), and the radial asymmetries of the planet gap. In particular, we focus on two specific moments in the simulation: the first one, noted 61, corresponds to the strong episode of gas depletion between 165 and 200 orbits, the second one, ℰ₂, corresponds to the last 100 orbits of the simulation, after the magnetic field has been divided in two regions of accumulation. During ℰ₁, the gap in the density lies approximately between R_p – 2R_hill and R_p + 3R_hill, whereas during ℰ₂, the outer gap edge is further out in the disk, near R_p + 5R_hill. In the graphs of this section, we highlight these three characteristic radii with vertical red dashed lines.

Figure 11 displays the radial profiles of Ṁ_acc (top panel) and the ejection efficiency ξ (bottom panel) for ℰ₁ (dashed line) and ℰ₂ (solid line). For these two episodes, the curves of Ṁ_acc and ξ are red in the gap and black beyond the gap. As a general remark, we find that the accretion is stronger on average in the gap than in its vicinity (e.g., at R ≃ 0.7 and R ≃ 1.3–1.5), and even stronger at higher magnetization. In addition, the ejection efficiency ξ is very small in the gap² (bottom panel of Fig. 11), indicating that almost all of the matter that enters the outer gap eventually leaves via the inner gap (in other words, the gas “leak” due to the outflow is negligible). The step in Ṁ_acc and the fact that ξ ≪ 1 in the gap are crucial to interpret the radial width asymmetry of the gap: on average, the gap progressively erodes the outer disk, while material tends to accumulate at the inner gap edge. The outward drift of the gap is therefore mainly due to a slight mismatch of the accretion rate in the gap and in the outer/inner disks. Such phenomenon has already been reported in MHD-driven disks, with the outward drift of the cavity in transition disks (Martel & Lesur 2022). If we compare the two episodes, Ṁ_acc is larger in the gap during ℰ₁ than during ℰ₂. Actually, the accretion rate reaches its highest value in the gap during ℰ₁, which is related to the saturation of the B_z everywhere in the gap during this epoch (see Fig. 8). Moreover, the step of Ṁ_acc in the gap and out of the gap is bigger during ℰ₁, with a steeper gradient (see in particular between R = R_p + 3R_hill and R ≃ 1.35 for ℰ₁ and between R = R_p + 5R_hill and R ≃ 1.5 for ℰ₂). We therefore expect a faster outward drift of the outer gap during ℰ₁ than during ℰ₂, which is indeed what we observe in the space-time diagram of log〈Σ〉_ϕ in the sixth panel of Fig. 4.

Furthermore, Ṁ_acc is approximately constant in the gap on average for both episodes, which is key for the interpretation of the gap asymmetry. If we consider the torques exerted on the gas in the gap, we have two main contributions. On one hand, the tidal torque 𝒫 is due to the presence of the planet and deposits angular momentum in the flow in the form of spiral wakes. These wakes, in turn, transport angular momentum in the radial fluid torques ℛ_R + ℳ_R. We therefore gather these three torques into a planet-induced torque 𝒯_R = ℛ_R + ℳ_R + 𝒫. On the other hand, the magnetic wind torque ℳ_surf is the main driver of accretion (see Fig. 1). In Fig. 12, we show the radial profiles in the gap of 𝒯_R/Ω_k (blue curve), ℳ_surf/Ω_k (yellow curve), and the total torque exerted on the gas (black curve) as expressed in Eq. (22), averaged over the last 100 orbits (episode ℰ₂). We note that the total torque (𝒯 + 𝒫)/Ω_k is close to the accretion rate Ṁ_acc/4π displayed in the top panel of Fig. 11, with a relative difference of ≃22% on average in the gap. In particular, we retrieve that the total torque is approximately constant in the gap, which means that the torques 𝒯_r/Ω_k and ℳ_surf/Ω_k that contribute to the accretion somehow counterbalance each other. We also retrieve a small step in the total torque at the gap’s outer edge, necessary to understand the gradual erosion of the outer disk by the gap. Regarding 𝒯_R, we can divide it in 𝒯_R– and 𝒯_R+, respectively for the planet-related torques in the inner gap and the outer gap. 𝒯_R+ is negative and tends to push material outward from the planet, with angular momentum that is deposited by the planet’s outer wake. However, 𝒯_R– is quite small and positive, and therefore seems inefficient at extracting angular momentum via the planet’s inner wake.

We can interpret the redistribution of angular momentum in the gap as follows. Near the gap outer edge, the planet-induced torque deposits an excess of angular momentum (𝒯_R+ is negative and minimum, see label a in Fig. 12) that is efficiently transferred to the wind and eventually vertically evacuated at the disk surface (ℳ_surf positive and maximum, see label b in Fig. 12). Near the gap inner edge, the planet does not extract efficiently angular momentum from the gas (𝒯_R– positive but close to zero), and the wind torque is also negligible at this location (ℳ_surf ≳ 0, see label c in Fig. 12). The end result is an apparent aborted wind in the inner gap regions, and an enhanced wind in the outer gap region.

Due to the predominance of the wind torque in the accretion, we focus in Fig. 13 on the evolution in the gap of the υ parameter, which is a dimensionless measure of ℳ_surf (Eq. (23)). This plot shows the 2D map of υ in (R-ϕ) coordinates for ℰ₂, as well as its azimuthally averaged version (solid black line). We also show the azimuthally averaged υ for ℰ₁ (dashed black line). During both episodes, υ is larger in the outer gap than in the inner gap, with an asymmetry of the angular momentum extraction Δυ ≃ 0.2 in both cases. In the (R-ϕ) plane, we detect a clear azimuthal asymmetry of υ, with a larger value at ϕ > 0. The wind torque seems affected by the outer planet’s spiral wake, with an enhanced υ along the wake toward the inner gap near the CPD, from R = 1.25 to R = 0.95. This increased υ is also accompanied by a localized enhancement of the magnetization near the outer wake, which would suggest that magnetic field accumulates efficiently from the outer wake to the CPD. This scenario of magnetic accumulation in the CPD seems to be confirmed by the launching of a collimated protoplanetary jet (see Sect. 3.2.2).

Near the gap’s inner edge, the wind torque is particularly small, and can even be negative close to the planet’s azimuth. This region of negative is represented in Fig. 13 with the dark blue area near R = 0.8, and can be related to label c in Fig. 12 at that same radial location. At this location, the wind is actually deflected due to the planet-induced deflection of the poloidal magnetic field lines at the disk surface toward the midplane, leading to a rearranging of (see in Fig. 9 the LIC texture and magnetic field lines in black that are nearly tangential to the surface layers at ±z_w ≃ ±4H(R), also near R = 0.8). We have checked that this deflection is all the more pronounced as the planet is massive. We can therefore have a strong magnetic accumulation both near the inner and outer gap edges as indicated in Figs. 8 and 9, and yet have an apparent inefficient (efficient) wind with a small (large) fraction of mass flux escaping from the inner (outer) gap edge, as suggested in Fig. 6 (see in particular the wind-driven flows pinpointed by the dotted black and white arrows).

We would like to draw attention to the fact that the 1D estimation of the velocity structure of the wind (and more generally the gas flow) can be misleading, as it actually results from a complex 3D dynamics. In a theoretical point of view, the radial evolution of 1D or 2D quantities like ℳ_surf or could suggest that there is no wind near the gap’s inner edge (the wind torque is almost zero), while the 3D dynamics indicates in reality that this behavior is due to the local deflection of the wind downward and then upward (with an N-shaped structure of the wind for the most massive planets, see e.g., the intersection of the upper and lower surfaces ±z_w ≃ ±4H(R) with the white streamlines in Fig. 6). Similarly in an observational point of view, for such a wind we would expect a predominant role of the emitting surface on the kinematic signature of the wind. At high altitude, for example for an optically thick tracer such as CO, one would expect to measure an outflow at all radii. At lower altitude (but still far enough from the disk bulk), for example for an optically thinner tracer such that ¹³CO, one would expect to measure an outflow signature above the outer gap and an inflow signature above the inner gap due to the deflection of the wind in this region. At even lower altitude, for example for an even optically thinner observational tracer like C¹⁸O, one would this time expect to be sensitive to the nearly sonic accretion flows, from the outer disk’s surface layers to the planet poles and from the inner gap to the inner disk’s surface layers.

To sum up, one of the key point here is the connection between the magnetic field accumulation, the wind torque and the planet-related torques. If Bp accumulates at the gap edges, we would expect a strong wind torque at these locations. Because the outer planet wake deposit angular momentum near the outer gap edge and because magnetic field lines are deflected by the planet at a certain altitude above the inner gap edge, our estimate of the wind torque is larger in the outer gap than in the inner gap. We next explore how the erosion of the outer disk by the gap affects the effective torque exerted by a magnetized disk onto massive planets.

Fig. 9 Poloidal magnetic field in the gap for the run Mj-β3, averaged in azimuth and over the last 100 orbits (episode ℰ2). The background color represents the logarithm of the poloidal magnetic field. The black streamlines and the LIC correspond to the radial and vertical components of the poloidal magnetic field. The white curves represents the radial profile of log(Σ), also averaged in azimuth and over the last 100 orbits. The red cross and solid black line pinpoint the location of the planet. Vertical black dashed lines delineate the gap edges at R ≃ Rp − 2Rhill and R ≃ Rp + 5Rhill. The upper and lower surfaces ±zw ≃ ±4H(R) are indicated by the black dotted lines.

Fig. 10 Azimuthal cut of the log (Σ), at the azimuth of the planet ϕp, at ≃ 114 orbits for 3Mj-β4. A protoplanetary jet is clearly visible, launching from the planet poles.

Fig. 11 Radial profile of the accretion rate Ṁacc (top panel) and the ejection parameter ξ (bottom panel), averaged azimuthally and temporally over the last 100 orbits (solid line) in the run Mj-β3 for ℰ2 and over 25 orbits during the episode ℰ1 of gas depletion and fast outward drift of the outer gap (dashed line). We removed the contribution close to the planet and its CPD in the evaluation of Ṁacc and ξ.

Fig. 12 Radial profile of the torques exerted by the wind (ℳsurf, yellow curve) and the sum of the radial torques and planet torques (ℛR, ℳr and 𝒫, blue curve), averaged over the episode ℰ2 in the run Mj-β3. The black curve displays the sum of all torques enumerated in Eq. (22). We removed the contribution close to the planet and its CPD in the evaluation of the planet and radial torques. We add three labels to identify different angular momentum redistribution phenomena that we believe are essential in apprehending the gap asymmetries. (a) Angular momentum deposition by the outer planet wake at the gap’s outer edge. (b) Angular momentum evacuation in the wind. (c) Inefficient angular momentum extraction by the inner planet wake at the gap’s inner edge, and small value of the wind torque.

Fig. 13 2D map of υ, i.e., the magnetic torque per unit of surface density, at the disk surface in (R-ϕ) coordinates, averaged over the episode ℰ2 in the run Mj-β3. We highlight in dark blue the region where the magnetic torque is negative, which corresponds here to a slight deposition of angular momentum. The solid and dashed black lines show respectively υ during ℰ2 and ℰ1.

3.3 Gravitational torques

We discuss here the evaluation of the gravitational torque Γ exerted by the gas onto the different planets in all 3D simulations, and focus more specifically on high-mass planets. Such torque is defined in IDEFIX as (33)

with the torque per unit of volume defined in Sect. 2.6 and evaluated for the cell i, j, k of volume . The estimate of this quantity is valuable for disk/planet interactions modeling, because it has a direct impact on how the planet would migrate if its orbital evolution were taken into account. In this case indeed, numerically, the value of the instantaneous torque is used to update the orbital parameters of the planet. Figure 14 displays the temporal evolution of Γ for our 8 nonideal MHD simulations and the inviscid hydrodynamic simulation 3Mj-α₀, with a moving average of 10 orbits and normalized by . Each color corresponds to a planet mass M_p, and the line-style indicates the value of β₀.

We first notice that Γ is negative for the three most massive planets, with a smaller |Γ| for larger planet masses. It suggests a slower inward migration for more massive planets. We observe indeed a factor ≃10 for |Γ| between the cases 3M_j and M_s. This is in agreement with the models of interactions between a gap-opening planet and a viscous disk: the amount of gas inside this gap is larger if the planet is not massive enough to empty its coor-bital region (see, e.g., Baruteau & Masset 2013). This excess of gas boosts the so-called Lindblad torque, which is an essential (usually) negative component in the total torque. A more negative torque for a less massive planet is indeed what we observe in Fig. 14.

Second, at fixed M_p (except for the least massive planet 10M_⊕), we show that β₀ seems to have little influence on the value of the torque, at least on the first 200 orbits. Still, it would seem that the absolute value of the torque ends up being slightly larger when the magnetization increases (i.e., when β₀ decreases), which is especially visible between ≃100 and 200 orbits. It means that the migration is faster when the magnetization is higher. It is consistent, because a larger magnetization leads to a denser gap (see in particular Sect. 3.1.2 and Fig. 4), and therefore an amplified Lindblad torque and a faster migration.

However, for the longer run Mj-β₃, we expect an additional mechanism to occur. Because the outer gap widens and drifts outward (i.e., the planet is closer to the inner gap edge), we expect the external Lindblad torque to weaken and decrease compared to the internal Lindblad torque, leading to a less and less negative total torque. We checked that the total outer torque |Γ₊| decreases more than the total inner torque |Γ₋|. We also noticed that the fluctuations of |Γ₊| and |Γ₋| are at first equivalent until |Γ₋| overcomes |Γ₊| near ≃280 orbits. At that time, the inner torque starts to strongly fluctuate, with fluctuations ≃3 times larger with respect to the outer torque. These fluctuations may be due to the formation of a thin density ring (corresponding to a pressure maximum) confined between the planet gap and a cavity that develops in the inner grid edge, probably linked to an artifact in the inner boundary (see end of Sect. 3.1.1). This ring is unstable to the RWI as a vortex develops, and becomes particularly strong after ≃330 orbits. In practice, the frequency of the planet torque fluctuations may be linked to the periodic transit of this vortex in the vicinity of the inner planet wake. Due to the weakening of the external torque, we expect an increasingly slower inward migration, and even a potential reverse migration, as shown in Fig. 14 by the sign reversal of Γ at 270 orbits. This phenomenon could also be obtained by the action of a positive component of the corotation torque exerted by the MHD wind, as proposed by some studies which have added hydrodynamic prescriptions of the impact of these winds on the migration of massive planets (Kimmig et al. 2020). The main difference here is that the outward migration they observe is due to a dynamical corotation torque linked to the inward accretion flow and leading to a type-III runaway migration. Even without mentioning the dynamical corotation torque, it could be challenging to study migration of such planets by 3D nonideal MHD simulations, because the dynamical evolution of the gap is complex and the characteristic timescale for migration reversal can be quite long (≳ 300 orbits for the run Mj-β₃).

Third, we notice that the total torque exerted by the gas on the least massive planet (10M_⊕) seems stochastic during the simulation (for the two magnetization cases), which makes it difficult to guess the direction of migration (i.e., the response of the planet to the total torque). This is similar to what has been observed by Nelson & Papaloizou (2004); Uribe et al. (2011), with strong fluctuations of the total torque experienced by low-mass planets on an orbital timescale, oscillating between positive and negative values. They suggested that the orbital evolution of such planet would proceed as a random walk. For the disk models without and with a low-mass planet, we show in Appendix B.2 that there is a low level of turbulence in terms of angular momentum transport (i.e., accretion is mostly wind-driven). However, we also show in Appendix B.3 that this turbulence plays a nonnegligible role in the stochastic nature of the planet torques. This relation between the turbulence intensity and the planet torque dispersion is not linear, because the level of turbulence increases faster than the torque dispersion. As a first experiment, we activated planet migration after the first ≃50 orbits of the run 10Me-β₄, for an initial disk density of Σ₀ = 1.25 × 10⁻⁴ in code units, which corresponds to ≥ 10 g cm⁻² at 10 au for a solar-mass star. The migration was found to be slow and directed outward during at least 10 planet orbits, at a speed of ≃3.5 × 10⁻³ au per orbit. Outward planet migration in wind-driven accretion disks could thus help the survival of planets at large radii that are often invoked to interpret substructures (Zhang et al. 2018; Isella et al. 2018; Lodato et al. 2019; Baruteau et al. 2019) and gas kinematics (Pinte et al. 2018, 2019, 2020; Izquierdo et al. 2021, 2022) in ALMA disks.

When there is an episode of strong accretion, Kimmig et al. (2020) expect a mass excess in front of the planet because the gas performs an inward U-turn and then is accreted in the inner disk before it can perform an outward U-turn. In other words, it is easier for the gas to cross the horseshoe region in front of the planet than behind due to the inward accretion flow (see their Fig. 11). This strong azimuthal asymmetry of the horseshoe U-turns (there are more inward U-turns than outward U-turns) would be associated with a positive dynamical corotation torque, and therefore an outward migration. However, the azimuthal asymmetry in our nonideal MHD runs is not the same depending on time and on the simulation, with sometimes an under-density in front, sometimes behind the planet in azimuth. This could be due to the joint influence of magnetic field accumulation and accretion sporadicity, either leading to a complex redistribution of matter inside the horseshoe region or an evacuation in the wind (see Appendix A). Besides, by looking at the streamlines in Fig. 15 for the run Mj-β₃ averaged over the last 100 orbits, we remark that the horseshoe region shrinks to a region of limited azimuthal extent, as a result of the fast radial accretion flow driven by the MHD wind. This region is centered around a libra-tion point closer to the planet in azimuth than the Lagrange point L5. This is consistent with Pepliński et al. (2008a,b) who observe such “tadpole-like” region around generalized Lagrange points G4 and G5 if the migration is respectively inward and outward. For slowly migrating planets, G4 and G5 correspond to L4 and L5 and the horseshoe region is fairly symmetric on the 2π angle of the corotation region. On the contrary, if planet migration is fast, the libration points move azimuthally closer to the planet as the horseshoe region shrinks. In the planet’s reference frame, because in our case the libration point is G5#L5, the inward accretion flow that crosses the planet gap on average over the last 100 orbits is equivalent to a fast outward migration, regarding the shape of the horseshoe region. It would therefore be worthwhile to determine how the migration of Jupiter-mass planets with such inward nearly-sonic accretion streamer inside the gap differs from the classical picture of type-II migration in viscous disks. In this latter case, we indeed expect a slow migration proportional to the disk viscosity, with a negligible role of gap-crossing flows (Robert et al. 2018).

As mentioned in Sect. 3.1.1, for the most massive planets (M_p ≥ M_j) and the low-magnetization case (β₀ = 10⁴), a strong vortex appears at the gap’s outer edge, before decaying after ≃ 100 orbits. It would be interesting to compare planet migration with the two-phase vortex-driven migration presented in Lega et al. (2021, 2022). Note that in all our simulations, small-scale nonaxisymmetric structures seem to develop everywhere in the disk and could interact with the planet during its migration. Note also that the torque calculated here is the one exerted by the gas on a fixed planet, which does not take into account the coupling between Γ and the migration. This link between migration and torque operates in what is called the dynamical corotation torque, which tends to amplify (positive feedback) or attenuate (negative feedback) the migration according to the physical properties of the disk and the coorbital region (Masset 2008; Paardekooper 2014). The next step is therefore to confirm the first results presented in Fig. 14 by letting the planet evolve radially according to Γ, and thus to study the impact of MHD winds on planetary migration. Even if the gaps seem to open in less than 100 orbits for the most massive planets whatever the magnetization (see Sect. 3.1.2), the strong variability in the accretion and the secular evolution of the gap (outer gap’s outward drift) can change the planet response: for the case Mj-β₃, the torque seems to reach a steady-state positive value solely after ≃300 planet orbits, which is challenging in terms of computational costs. It also raises the question of the typical timescale at which to activate planet migration.

Fig. 14 Temporal evolution of the gravitational torques exerted by the gas onto the planets for the nine runs. The torques are normalized by . The different colors show different planet masses: 10M⊕ (green), Ms (yellow), Mj (red) and 3Mj (blue). The different line-styles indicate the initial magnetization: β0 = +∞ (inviscid hydrodynamic run, dotted line), β0 = 104 (solid line) and β0 = 103 (dashed line). The transparent variable curves show the instantaneous torques, on top of which are displayed the torques with a moving average of 10 orbits.

Fig. 15 Gas surface density log (Σ), averaged over the last 100 orbits, for Mj-β3. The streamlines in the midplane are shown in black.

4 Summary and discussion

4.1 Summary

In this study, we performed with the IDEFIX code 3D nonideal MHD simulations of a planet in a fixed circular orbit, interacting with a wind-launching protoplanetary disk. We focused in particular on the ability of massive planets to open a gap in the presence of a vertical magnetic field. One of the challenges of this work is to achieve a balance between the large-scale study of a disk with wind and dynamical processes localized around a planet, as well as to isolate the multiple phenomena that emerge all at once in the simulations.

We find that a gap opening occurs, with density gaps being induced by self-organization everywhere in the disk, especially at higher magnetization, as well as with them being induced by the planet if the planet mass is high enough (Sect. 3.1). When comparing a gap opening in our MHD simulations, we notice that the planet gap is denser in a counter-intuitive way when the initial magnetization increases (i.e., when β₀ decreases), which we interpret as a consequence of the strong dependence of the mass accretion rate with the disk magnetization. We derived a gap-opening criterion when accretion is dominated by MHD winds, which depends on the planet-to-primary mass ratio q_p, the aspect ratio h, and the normalized magnetic wind torque υ (Sect. 3.1.3). In addition, the combination of upper-layer accretion flows, planet-driven flows, and wind-driven flows drives a large-scale gas motion that takes the form of a meridional circulation, especially in the outer disk (Sect. 3.2.1). The accretion of material comes from the surface layers of the outer disk to the inner disk. More specifically, two accreting surface layers fall directly from the outer disk toward the planet poles. For this study, we discarded the accretion onto the planet, but this should have a strong impact on the planet’s internal and orbital evolutions, as presented in Gressel et al. (2013).

The magnetic field tends to efficiently accumulate in the planet gap, which is coherent with Zhu et al. (2013), Keith & Wardle (2015), and Carballido et al. (2017). Depending on the initial magnetization and the planet mass, it is found to accumulate firstly at the gap center, and later on at both gap edges. Such magnetic accumulation could be at the origin of complex effects, such as the launching of a protoplanetary jet, magnetic reconnection, and magnetic braking. We observe a well-known anticorrelation between the density and magnetic field, particularly visible in the run Mj-β₃ with an episode (episode ℰ₁) of strong gas depletion in the gap associated with an efficient magnetic field accumulation. The radial extent of the density gap is larger than the magnetic accumulation region, as observed in other works (Sect. 3.2.2).

When β₀ = 10³, that is when the amplitude of the magnetic field is large in the gap, we have an efficient angular momentum removal by the wind torque, corresponding to a high equivalent ≃ 5 in the gap for the run Mj-β₃ averaged over the last 100 orbits (episode ℰ₂). We note that the radial stresses that act on the gas are also larger on average when the magnetization increases, leading to an equivalent Shakura & Sunyaev α_v ≲ 10⁻² in the cases β₀ = 10⁴, and α_v > 10⁻² when β₀ = 10³ in most of the disk, reaching ≃1 for the deepest gaps.

The torque exerted by the planet onto the gas (planet-induced torque) mainly deposits angular momentum near the gap’s outer edge. On the other hand, the magnetic wind torque is almost always extracting angular momentum from the disk and has therefore almost the same sign on both sides of the planet. The only exception is close to the gap’s inner edge, and all the more so when the planet mass is larger. In that case, the outflow and its poloidal magnetic field lines are deflected by the planet’s gravitational field at the disk surface, leading to a decrease in the wind torque which can even become negative. We note that this deflection could be interpreted as a local inflow signature (this is a purely geometrical effect as deflected streamlines eventually escape the disk). For low-mass planets, the planet-induced torque is relatively small, as the planet hardly generates a gap. The wind torque is therefore dominant and approximately constant in the gap, regardless of the initial magnetization. For high-mass planets, the wind torque is quenched in the inner gap due to this effect of magnetic field lines’ deflection, while planet-induced and wind torques have opposite signs in the outer gap. This forces a symmetry breaking between the inner and outer gaps. The end result is an outflow that is seemingly quenched from the inner gap because of the bending of magnetic lines. On the contrary, an enhanced wind is emitted from the outer gap, fed by the positive planet torque. The resulting asymmetry of the wind torque in the gap is somehow at the origin of the step in the accretion rate which in turn leads to a gradual erosion of the outer disk, only visible in our high planet mass (M_p ≥ M_j) and high magnetization (β₀ = 10³) runs. For weaker magnetizations (β₀ = 10⁴), the gap’s outer drift is small, if not negligible, on the small timescales considered here (200 orbits).

The complex redistribution of angular momentum in the gap is related to both the planet torque and the wind torque and plays a key role in the accretion process. In particular, we find that the accretion rate Ṁ_acc is nearly constant and comparable that of the disk in the gap. Furthermore, at high magnetization (β₀ = 10³), the enhanced wind torque induces a top-down non-symmetric, nearly sonic accretion streamer in the planet gap, as proposed by Frank et al. (2014). We note that the accretion in the gap is sporadic, leading to azimuthal asymmetries of the horseshoe region that could generate an unsaturated corotation torque (Ogihara et al. 2017; Kimmig et al. 2020). The temporal variability of accretion and magnetic field accumulation is also echoed in the temporal variability of the vertical redistribution of matter in and around the gap. For example, gas is sometimes evacuated directly into the wind from the outer disk, or, only a few orbits later, wraps poloidally around the planet (Appendix A). Although the accretion rate does not evolve much radially in the disk, it is on average slightly larger in the gap than at the gap edges. Such a mismatch in the radial profile of Ṁ_acc leads to strong radial asymmetries of the gap both in width and in depth: there is a mass loss (accumulation) in the outer (inner) gap, which is slightly emptier (denser). More specifically, the outer gap is wider and emptier at higher magnetization due to a stronger wind torque.

Regarding tidal torques exerted on massive planets in magnetized disks, the radial asymmetries of the planet gap gradually lead to a slightly positive gravitational torque exerted by the gas onto the planet. This means that the inward migration of a massive planet in such a disk could be slowed down, then stopped, and even reversed. For low-mass planets, the gravitational torque is stochastic because of turbulence oscillating between positive and negative values. It would seem, however, that the migration of such planets is slow and directed outward (Sect. 3.3). This is to be confirmed by enabling planet migration in future numerical simulations.

4.2 Comparison to works using prescribed wind torques

We can estimate α_v and α_dw in our simulations and compare them with Fig. 4 from Elbakyan et al. (2022), which shows the dependence of the minimum planet mass needed to open a gap as a function of α_v and α_dw. In our cases where β₀ = 10⁴, with the exact definitions of α_v and α_dw in Eqs. (4) and (6) by Tabone et al. (2022, neglecting ℛ_surf and averaged in time over the last 100 orbits), we find outside the planet gap that α_v = 10⁻³ − 10⁻² and α_dw = 4 × 10⁻². This corresponds to a minimum pessimistic (concerning their model) planet mass needed to open a gap of ≃ 1M_j, but our run Ms-β₄ is able to produce a gap. Likewise, when β₀ = 10³ and outside the planet gap, α_v = 3 × 10⁻³ − 6 × 10⁻² and α_dw = 6 × 10⁻² − 0.9 (also averaged in time over the last 100 orbits). When β₀ = 10³, the Maxwell stresses dominate their Reynolds counterparts, both radially and at the disk surface. Looking at Fig. 4 of Elbakyan et al. (2022), this corresponds to a minimum (very) pessimistic planet mass needed to open a gap of ≃2M_j, but here again, the run Mj-β₃ is able to produce a gap. One major difference with Elbakyan et al. (2022) is that α_v and α_dw strongly depend on the magnetic field. Because the distribution of magnetic field relies on self-organized gaps and planet gaps, α_v and α_dw also depend on the localization in the disk. One aspect that could be improved when prescribing the wind torque is to use a α_dw that is not constant in time (accretion bursts) and in radius (especially in the gap, with values 10–100 times larger). We would also need a way to model the accumulation of the magnetic field in the gaps, as it is important to understand the gaps’ radial asymmetries. It is important to note that we retrieved a correct order of magnitude for the critical planet-to-primary mass ratio to open a gap with the criterion we derived in Eq. (31).

In Kimmig et al. (2020), they consider the migration of massive planets in 2D, with a wind torque prescribed via a constant lever arm λ and constant mass-loss parameter b. This latter quantity is related to the vertical mass flux . However, if we compute from our 3D simulations and use it to define an equivalent mass loss, we find that b is not constant and varies radially. In particular, for (Mj-β₃), b reaches a few 10⁻³ in the disk, slightly increasing radially, but an order of magnitude higher in the gap during the episode ℰ₂ (i.e., during the last 100 orbits), reaching 1.75 × 10⁻². Considering Fig. 7 of Kimmig et al. (2020, or the two bottom panels of their Fig. 5), it could have a strong influence on the Jupiter-mass planet migration, as in one case we would have almost no migration (b = 10⁻³) and in the other case a fast runaway outward migration (b ≃ 10⁻²). The lever arm λ is more complicated to characterize as there is some inflow at the gap edges, leading to sign reversal of the mass outflow rate. However, in the gap, λ is larger than 2.5, reaching ≃7.5 in some gaps. Kimmig et al. (2020) chose a constant lever arm λ = 2.25, which is coherent with the values we obtain in the disk but underestimated in the planet-induced gaps.

Lega et al. (2022) show that when the planet-induced torque was larger than the wind torque, matter should pile up at the gap outer edge leading to a fast inward migration and the blocking of the accretion flow. In our simulations, we never observed the blocking of the accretion flow. On the contrary, we find that a large positive planet-induced torque is usually compensated for by a stronger wind torque (Fig. 12) so as to keep a constant accretion rate through the gap. This discrepancy is probably due to the fact that the wind responds dynamically to the planet-induced torque, but also to the redistribution of the magnetic field that changes the efficiency of magnetic wind launching. The overall conclusion is that a constant wind torque, which is often prescribed in effective 2D or 3D hydrodynamic models, is too simplistic to capture the dynamics of planet-disk-wind interactions.

Finally, as observed in McNally et al. (2020), our runs involving low-mass planets (M_p = 10M_⊕) confirm that surface accreting layers and the midplane seem decoupled regarding accretion, even if we detect highly variable gravitational torques Γ exerted on such planets. We would therefore need to activate planetary migration in order to assess the impact of wind-driven accretion on the orbital evolution of low-mass planets.

Acknowledgements

The authors acknowledge support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (Grant agreement no. 815559 (MHDiscs)). This work was granted access to the HPC resources of IDRIS under the allocation A0100402231 made by GENCI. The 3D simulations were performed with the version v0.8.0 of IDEFIX, commit bbffc9390353ab917728734e77543f69041eb003. Some of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. The data presented in this work were processed and plotted with Python via various libraries, in particular numpy, matplotlib and scipy, but also the following personal projects under development: nonos (https://github.com/volodia99/nonos) to analyze results from IDEFIX simulations, cblind (https://github.com/volodia99/cblind) for most of the colormaps and lick (https://github.com/volodia99/lick) for the visualizations in line integral convolution. G.W.F. wishes also to thank Clément Baruteau for fruitful discussions.

Appendix A Temporal variability of the flow

Fig. A.1 Temporal variability of the azimuthally averaged mass flux after 119.75 (top panel), 121.25 (middle panel) and 129.5 (bottom panel) orbits for the run 3Mj-β3, with specific features represented by the black arrows: a. South-North U-turn of material at the transition between the outer disk and the gap. b. South-North U-turn of material at the transition between the gap and the inner disk. The material either rotates in the poloidal plane around the planet, or is evacuated in the upper layer’s wind c. Top-down nonsymmetrical accretion streamer from the outer disk to the inner disk.

The aim of this appendix is to illustrate the temporal variability of the accretion streamer and of matter transport from the outer disk to the inner disk via the planet gap. In Figure A.1, we show three snapshots between 119.75 and 129.5 orbits of the magnitude of the poloidal mass flux 〈ρυ_p〉_ϕ in log scale, averaged azimuthally for the run 3Mj-β₃. As in Figure 6, the texture in LIC and the white arrows indicate the direction of the radial (〈ρυ_R〉_ϕ) and vertical (〈ρυ_z〉_ϕ) components of the mass flux. Black arrows illustrate schematically streamlines of interest for the gas. The mass flux is extremely stochastic, with a fast temporal variability (inferior to the orbital timescale). Panel c. displays how the material sometimes crosses the gap from the outer disk to the inner disk via an accretion streamer at z < 0 unperturbed by the disk morphology, as described in Section 3.2.1. On the contrary, the flow is sometimes affected by the presence of the gap and is forced to perform sporadic South-North U-turns at the transition with the gap edges. Such material is either evacuated in the wind or redistributed in the outer disk and/or in the gap. In panel a., the U-turns carry material in the outer disk from the lower accretion layer to an upper wind-launching layer, and thus partially prevent the material from entering the planet gap. In panel b., the fast accreting material crosses R = 1, under the planet, from the outer gap to the inner gap. Then, the South-North U-turns bring a fraction of this material directly to the wind, while the rest seems to accumulate in the partially refilled gap, rotating in a poloidal way around the planet radial location (note that the graphs are azimuthally averaged maps, not azimuthal cuts). These U-turns are beyond the scope of our study, but could result from the accumulation of magnetic field at the gap edges.

Appendix B Planet-free disk model, turbulence, and stochasticity

Appendix B.1 Planet-free disk model

In this paragraph, we develop how the protoplanetary disk model behaves in the absence of the planets for an initial β₀ = 10³, and compare with the run 10Me-β₃ in order to determine if the structures observed in that simulation are mainly due to the planet or not. If we look at the structures in the disk, they are similar without and with a low-mass planet, at least concerning the distributions of density and poloidal magnetic field, with in particular the emergence of self-organized structures (see also Section 3.1.1 and, e.g., the top-right panel of Figure 2). We can however highlight the fact that the presence of such planet can create rings of underdensity locally on both sides of its radial location (see the first panel in Figure B.1, near R = 0.85 and R = 1.3) that accumulate the magnetic field more efficiently (see the two peaks of ≃ 15% and ≃ 40% in the second panel in Figure B.1 at those same locations). At the same time, this accumulation of magnetic field compared to a run without planet slightly enhances the wind torque by ≃ 30% (see the third panel in Figure B.1). Note also that the accretion via the radial mass flux is also increased by ≃ 20% in the run 10Me-β₃ around the planet between R = 0.7 and R = 1.3, although quite variable in that region (see the fourth panel in Figure B.1). The step in Ṁ_acc between R = 1.2 and R = 1.4 seems to indicate that the structures observed are drifting outward during the episode considered here.

Fig. B.1 Radial profiles of different quantities averaged between 50 and 100 orbits, compared in percent between a disk model without and with a 10 Earth mass planet: gas surface density Σ in blue (first panel), amplitude of the poloidal magnetic field Bp in green (second panel), wind torque at the disk surface ℳsurf in yellow (third panel), radial mass flux Ṁacc in red (fourth panel). The location of the planet is indicated by the vertical black dashed line. Note that such disk model evolves dynamically, and the structures that are created at a given time will evolve radially and in amplitude, making the comparison difficult.

Appendix B.2 Turbulence

Our nonideal MHD disk models are susceptible to be MRI and possibly VSI unstable. Applying the VSI cooling timescale criterion of Lin & Youdin (2015) to our parameters (γ = 5/3, q = 1, h = 0.05, with –q being the power-law exponent of the radial temperature profile near the midplane (R)), we need (B.1)

to get “vigorous” VSI. Hence, our models are likely to have severely quenched (or even suppressed) VSI signatures since we always use Ωℬ = 0.1.

Fig. B.2 Space-time diagram in log scale of the vertical component of the mass flux in the midplane and azimuthally averaged 〈ρυz〉ϕ,z=0, for the low-mass planet and low magnetization case 10Me-β4.

In Figure B.2, we show the space-time diagram of the vertical component of the mass flux at the midplane in our 3D run with a low-mass planet, averaged in ϕ across the disk. When the VSI is present, one expects to see obvious signatures of large-scale corrugation modes (e.g., fig. 3 in Svanberg et al. 2022). Our diagram does not exhibit any of these patterns, suggesting that the VSI does not seem to develop in the simulation, which is coherent with the criterion of Eq. (B.1). We also checked the nondevelopment of the VSI at higher magnetization. In particular, we note the absence of large-scale coherent vertical motions most of the time during the simulation. We therefore expect the MRI to be the main driver of turbulence in our models, but not to be the dominant process for angular momentum transport compared to MHD winds since . Other studies observe a low level of turbulence due to MRI for such nonideal MHD simulations with similar parameters (Cui & Bai 2022).

We show in Figure B.3 how the presence (in green) or absence (in gray) of a low-mass planet impacts the radial profiles of the total α_v (solid lines, as defined in Eq. 25) and (dash-dotted lines). This last quantity is linked to the turbulent component of the Reynolds stress, defined as , the second term of this expression corresponding to the laminar component of the Reynolds stress. Note that the Reynolds stress is defined from the deviation from Keplerian rotation, but is not supposed to be vanishing when averaged in time (see Eq. 20 in Balbus & Papaloizou 1999). Hence this deviation may be split into a laminar term (i.e., nonzero when averaged in time) and a turbulent term (i.e., zero when averaged in time). We plot the turbulent stress in Figure B.3 (dash-dotted lines). We find that there is indeed turbulence in our disk model, both with and without a planet, but at a low level (a few ) and therefore not sufficient to explain alone angular momentum transport. This confirms the predominance of the ℳ_surr term in the angular momentum budget, as presented in Figure 1.

Fig. B.3 Radial profiles of the logarithm of the total αv (solid lines) and (dash-dotted lines) related to the turbulent Reynolds stress, with (green curves) and without (gray curves) a 10M⊕ planet, for an initial β0 of 103 and averaged between 50 and 100 orbits at R = 1.

Besides, is even not dominant in the total radial transport as stated earlier, (it is about ≃ 5 times smaller than the total radial stress). In particular, for the run 10Me-β₄, we checked that for β₀ = 10⁴, which is similar to the value obtained in Cui & Bai (2022) for their run with instantaneous cooling and A_m = 1 constant within ±3.5H. Far from the planet position (R > 2), the total radial transport of angular momentum does not seem to be impacted by the planet, with the overlap of the green and gray solid lines. Close to R = 1, the total stress is slightly larger (see, e.g., near R = 1.25) in the simulation with planet, which is due to the slight carving of underdensity rings on both sides of the planet compared to a planet-free case. This underdensity is deeper in the outer part (≃ 20%) than in the inner part (10%) compared to the run without planet (see top panel of Figure B.1), which is consistent with a α_v larger in the outer part than in the inner part. Note, as indicated in Appendix B.1, that these underdensities are also associated to a slightly more efficient wind torque at these locations for the simulation with planet (29% and 37% for ℳ_surf respectively in the outer part and the inner part). If now we try to determine the impact of the low-mass planet on turbulence and , we observe that the level of turbulence is globally lower in the simulation with planet than in the simulation without planet by a factor 2 – 3, except close to the planet location (between R = 0.8 and R = 1.3) where the level of turbulence is similar in both cases.

Appendix B.3 Stochasticity

In this paragraph, we aim at making this link between the level of turbulence in the disk with and without a low-mass planet at R = 1 and the standard deviation of the probability density function of the planet torque Γ, and in particular its stochastic nature.

Fig. B.4 Frequency distribution of the torques exerted on the least massive planets plotted in green in Figure 14, e.g., for the runs 10Me-β3 and 10Me-β4. These two histograms show the stochasticity of the torques experienced by the 10 Earth mass planets, via a Gaussian distribution with standard deviation of σ3 ≃ 4 × 10−2 and σ4 ≃ 2 × 10−2 for an initial β0 of 103 and 104 respectively.

To do so, we plot in Figure B.4 the frequency distribution of the planet torques Γ for the runs 10Me-β₃ and 10Me-β₄ that were displayed in green in Figure 14. We sample these planet torques every twentieth of an orbit during the simulations. We then compute histograms of the values taken by the torques in each bin with a bin width of Δ(Γ/q) = 2.8 × 10⁻³, with a procedure similar to the one presented in (Nelson 2005). Fitting these distributions with Gaussian profiles, we obtain a standard deviation of σ₃ ≃ 4 × 10⁻² and σ₄ ≃ 2 × 10⁻² for an initial β₀ of 10³ and 10⁴ respectively. On the other hand, if we quantify the

Reynolds stresses, and more precisely in both runs 10Me-β₃ and 10Me-β₄, we realize that decreasing the β₀ by a factor 10 increases at R = 1 by a factor ≃ 5, whereas the estimated standard deviation of the probability density function of the planet torque increases only by a factor ≃ 2. Therefore, the dispersion of the planet torque does not scale directly with .

As a summary, when quantifying the link between the turbulence and the stochastic torques exerted on low-mass planets, we find that the turbulence has certainly a non-negligible role (the stochasticity seems to increase when the magnetic field strength increases) but not a proportional one, because increases faster than the torque dispersion.

Appendix C Density and magnetization in a gap

We focus in this appendix on the opening of a gap for a 3 Jupiter mass planet, and especially the evolution of the density and magnetization inside the gap.

Fig. C.1 Temporal evolution of the gas surface density (top panel) and the magnetization (bottom panel), in a gap carved by a 3 Jupiter mass planet, for three initial magnetizations (103, 104, and a viscous case).

To do so, we evaluate in the top panel of Figure C.1 the temporal evolution of the azimuthally averaged gap density, radially averaging the annulus spanning from R_p – 2R_hill to R_p + 2R_hill and excised from −2R_hill/R_p to 2R_hill/R_p in order to remove the contribution of the planet and its circumplanetary disk (CPD). R_hill is the planet’s Hill radius and is defined in Eq. (5). This procedure follows the evaluation of the gap’s surface density Σ_gap presented in Fung et al. (2014) and Fung & Chiang (2016). We also evaluate in the bottom panel of Figure C.1 the temporal evolution of the azimuthally averaged β in the gap, considering its minimum in the radial interval [R_p – R_hill, R_p + R_hill] and also excised from −2R_hill/R_p to 2R_hill/R_p. We define here several quantities that can be useful to understand and analyze qualitatively the process of gap opening:

• is the quasi-equilibrium value for the density in the gap, that is such that the mean value of Σ_gap does not evolve much with time.

• is the variability of the density in the gap once it has reached .

• is the instantaneous gap-opening speed, from the initial density at the planet location to .

• is the time needed by Σ_gap to reach .

Using this metrics and focusing on the 3M_j planet cases, we notice in the top panel of Figure C.1 that increases with the magnetization. More precisely, Σ_gap reaches ≃ 3 × 10⁻³ of its initial value Σ₀ (and still decreasing) for 3Mj-α₀ after 200 planet orbits, whereas it reaches Σ₀ for 3Mj-β₄ after planet orbits and Σ₀ for 3Mj-β₃ after planet orbits. This is consistent with our interpretation based on the interplay between the accretion rate Ṁ_acc, β and Σ and presented in the second paragraph of Section 3.1.2. also increases with the magnetization, and could be linked to stronger time-variability in the gap when the initial magnetization increases, as mentioned in the last paragraph of Section 3.1.2. We do not focus on these first two quantities here, but rather on υ_gap and . On the one hand, the gap-opening speed is larger when β₀ d_fecreases, at least early in the simulations and before (see, e.g., between 20 and 90 orbits). It would suggest that a torque linked to the ac_fcretion is actually helping to form the gap. On the other hand, decreases when β₀ decreases, which means that the gap reaches its minimum density earlier when the magnetization is higher, respectively after 200, 90 and 50 planet orbits for 3Mj-α₀, 3Mj-β₄ and 3Mj-β₃ respectively.

In order to reconcile these two processes, we can argue that the magnetic torque linked to the vertical extraction of angular momentum helps the planet-related torque to open a gap, and even more if β₀ decreases. The deeper the gap is, the more magnetic field accumulates in it, which leads to an increase of the wind torque that helps even more to open the gap. We thus witness a runaway gap opening, particularly visible in the concave segment of the curves just before (see the β₀ = 10⁴ and 10³ cases in the top panel of Figure C.1). At that moment, when β_gap reaches a threshold on the order of 1 – 10, a quasi steady state is achieved in the gap-opening process. This happens earlier when the magnetization is initially lower. It is confirmed by looking at the bottom panel of Figure C.1, which shows that β_gap reaches a value < 10 near 50 and 90 orbits respectively when β₀ = 10³ and 10⁴.

References

All Tables

All Figures

	Fig. 1 Radial profile of the four torques ℛR (blue), ℳR (red), ℛsurf (green), ℳsurf (yellow) and their sum 𝒯 (black) normalized by , azimuthally averaged and temporally averaged over the last 100 orbits for the run Mj-β3. We consider here only the logarithm of the absolute value of these quantities.
In the text

	Fig. 2 Gas surface density maps for the eight MHD runs, after 200 planet orbits. The four rows show the different values of the planet mass. The two columns show the different values of the initial β0 parameter. The white circle marks the planet location.
In the text

Fig. 3 Sketch of the different asymmetries of the planet gap, for massive planets and highly magnetized disks. The yellow and black circles represent respectively the central star and the planet. The horseshoe region is represented in shades of gray, with high-density zones in light gray and low-density zones in dark gray. The dashed line corresponds to the planet orbit. The global shape of the gap results from the combination of three asymmetries: (a) Azimuthal asymmetry of the horseshoe region. (b) Radial width asymmetry, with the outer gap wider than the inner gap. (c) Radial depth asymmetry, with the outer gap deeper than the inner gap.

In the text

Fig. 4 Space-time diagrams showing the evolution of the logarithm of the gas surface density Σ for the nine runs divided in four groups of different planet masses. The top plot of each group has the lowest magnetization, whereas the bottom plot has the highest. The radial extent of all the diagrams range from 0.5 to 1.5 code units, with the planet at R = 1. We pinpoint as ℰ1 and ℰ2 two episodes in the run Mj-β3 that will be discussed in Sect. 3.2.3.

In the text

	Fig. 5 Gap-opening criterion with MHD winds displayed in the (qp,β0) domain. The dashed line indicates the minimum mass to open a gap from this criterion, using a scaling law between u and β0. Planets in the hatched area are able to open a gap. To obtain the pink box plots, we estimated υ in our low-mass planet simulations.
In the text

Fig. 6 Meridional flows around the planet for the run 3Mj-β4, averaged in azimuth and over the last 100 orbits. The background color represents the logarithm of the poloidal mass flux. The white streamlines and the LIC correspond to the radial and vertical components of the poloidal mass flux. Dashed black arrows, solid black arrows and dotted black and white arrows show respectively the planet-driven flows, the upper-layer accretion flows and the wind-driven flows. The upper and lower surfaces ±zw ≃ ±4H(R) are indicated by the black dotted lines. The planet is indicated by the white circle.

In the text

Fig. 7 2D map of the Mach number, defined as , in the (R-z) plane for the run Mj-β3, averaged in azimuth and over the last 100 orbits. The background color represents the logarithm of the Mach number. The black streamlines and the LIC correspond to the radial and vertical components of the poloidal Mach number. The nearly sonic top-down nonsymmetrical streamer is visible in white, in the gap and under the planet indicated by the white circle.

In the text

	Fig. 8 Space-time diagram of the vertical component of the magnetic field in the midplane , for the runs Mj-β4 (top panel) and Mj-β3 (bottom panel). The horizontal dashed lines indicate the location of the planet. The two episodes ℰ1 and ℰ2 are indicated in white.
In the text

Fig. 9 Poloidal magnetic field in the gap for the run Mj-β3, averaged in azimuth and over the last 100 orbits (episode ℰ2). The background color represents the logarithm of the poloidal magnetic field. The black streamlines and the LIC correspond to the radial and vertical components of the poloidal magnetic field. The white curves represents the radial profile of log(Σ), also averaged in azimuth and over the last 100 orbits. The red cross and solid black line pinpoint the location of the planet. Vertical black dashed lines delineate the gap edges at R ≃ Rp − 2Rhill and R ≃ Rp + 5Rhill. The upper and lower surfaces ±zw ≃ ±4H(R) are indicated by the black dotted lines.

In the text

	Fig. 10 Azimuthal cut of the log (Σ), at the azimuth of the planet ϕp, at ≃ 114 orbits for 3Mj-β4. A protoplanetary jet is clearly visible, launching from the planet poles.
In the text

Fig. 11 Radial profile of the accretion rate Ṁacc (top panel) and the ejection parameter ξ (bottom panel), averaged azimuthally and temporally over the last 100 orbits (solid line) in the run Mj-β3 for ℰ2 and over 25 orbits during the episode ℰ1 of gas depletion and fast outward drift of the outer gap (dashed line). We removed the contribution close to the planet and its CPD in the evaluation of Ṁacc and ξ.

In the text

Fig. 12 Radial profile of the torques exerted by the wind (ℳsurf, yellow curve) and the sum of the radial torques and planet torques (ℛR, ℳr and 𝒫, blue curve), averaged over the episode ℰ2 in the run Mj-β3. The black curve displays the sum of all torques enumerated in Eq. (22). We removed the contribution close to the planet and its CPD in the evaluation of the planet and radial torques. We add three labels to identify different angular momentum redistribution phenomena that we believe are essential in apprehending the gap asymmetries. (a) Angular momentum deposition by the outer planet wake at the gap’s outer edge. (b) Angular momentum evacuation in the wind. (c) Inefficient angular momentum extraction by the inner planet wake at the gap’s inner edge, and small value of the wind torque.

In the text

	Fig. 13 2D map of υ, i.e., the magnetic torque per unit of surface density, at the disk surface in (R-ϕ) coordinates, averaged over the episode ℰ2 in the run Mj-β3. We highlight in dark blue the region where the magnetic torque is negative, which corresponds here to a slight deposition of angular momentum. The solid and dashed black lines show respectively υ during ℰ2 and ℰ1.
In the text

Fig. 14 Temporal evolution of the gravitational torques exerted by the gas onto the planets for the nine runs. The torques are normalized by . The different colors show different planet masses: 10M⊕ (green), Ms (yellow), Mj (red) and 3Mj (blue). The different line-styles indicate the initial magnetization: β0 = +∞ (inviscid hydrodynamic run, dotted line), β0 = 104 (solid line) and β0 = 103 (dashed line). The transparent variable curves show the instantaneous torques, on top of which are displayed the torques with a moving average of 10 orbits.

In the text

	Fig. 15 Gas surface density log (Σ), averaged over the last 100 orbits, for Mj-β3. The streamlines in the midplane are shown in black.
In the text

Fig. A.1 Temporal variability of the azimuthally averaged mass flux after 119.75 (top panel), 121.25 (middle panel) and 129.5 (bottom panel) orbits for the run 3Mj-β3, with specific features represented by the black arrows: a. South-North U-turn of material at the transition between the outer disk and the gap. b. South-North U-turn of material at the transition between the gap and the inner disk. The material either rotates in the poloidal plane around the planet, or is evacuated in the upper layer’s wind c. Top-down nonsymmetrical accretion streamer from the outer disk to the inner disk.

In the text

Fig. B.1 Radial profiles of different quantities averaged between 50 and 100 orbits, compared in percent between a disk model without and with a 10 Earth mass planet: gas surface density Σ in blue (first panel), amplitude of the poloidal magnetic field Bp in green (second panel), wind torque at the disk surface ℳsurf in yellow (third panel), radial mass flux Ṁacc in red (fourth panel). The location of the planet is indicated by the vertical black dashed line. Note that such disk model evolves dynamically, and the structures that are created at a given time will evolve radially and in amplitude, making the comparison difficult.

In the text

	Fig. B.2 Space-time diagram in log scale of the vertical component of the mass flux in the midplane and azimuthally averaged 〈ρυz〉ϕ,z=0, for the low-mass planet and low magnetization case 10Me-β4.
In the text

	Fig. B.3 Radial profiles of the logarithm of the total αv (solid lines) and (dash-dotted lines) related to the turbulent Reynolds stress, with (green curves) and without (gray curves) a 10M⊕ planet, for an initial β0 of 103 and averaged between 50 and 100 orbits at R = 1.
In the text

	Fig. B.4 Frequency distribution of the torques exerted on the least massive planets plotted in green in Figure 14, e.g., for the runs 10Me-β3 and 10Me-β4. These two histograms show the stochasticity of the torques experienced by the 10 Earth mass planets, via a Gaussian distribution with standard deviation of σ3 ≃ 4 × 10−2 and σ4 ≃ 2 × 10−2 for an initial β0 of 103 and 104 respectively.
In the text

	Fig. C.1 Temporal evolution of the gas surface density (top panel) and the magnetization (bottom panel), in a gap carved by a 3 Jupiter mass planet, for three initial magnetizations (103, 104, and a viscous case).
In the text