# Planet-disk-wind interaction: The magnetized fate of protoplanets

Gaylor Wafflard-Fernandez, Geoffroy Lesur

► To cite this version:

Gaylor Wafflard-Fernandez, Geoffroy Lesur. Planet-disk-wind interaction: The magnetized fate of protoplanets. Astronomy & Astrophysics - A&A, 2023, 677, pp.A70. 10.1051/0004-6361/202245305. hal-04961148

HAL Id: hal-04961148

https://hal.science/hal-04961148v1

Submitted on 21 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Planet-disk-wind interaction: The magnetized fate of protoplanets

Gaylor Wafflard-Fernandez<sup>D</sup>, and Geoffroy Lesur<sup>D</sup>

Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

e-mail: gaylor.wafflard@univ-grenoble-alpes.fr

Received 27 October 2022 / Accepted 18 June 2023

## ABSTRACT

**Context.** Models of a planet-disk interaction are mainly based on 2D and 3D viscous hydrodynamic simulations. In such models, accretion is classically prescribed by an $\alpha_v$ parameter which characterizes the turbulent radial transport of angular momentum in the disk. This accretion scenario has been questioned for a few years and an alternative paradigm has been proposed that involves the vertical transport of angular momentum by magneto-hydrodynamic (MHD) winds.

**Aims.** We revisit planet-disk interactions in the context of MHD wind-launching protoplanetary disks. In particular, we focus on the planet's ability to open a gap and produce meridional flows. The accretion, magnetic field, and wind torque in the gap are also explored, as well as the evaluation of the gravitational torque exerted by the disk onto the planet.

**Methods.** We carried out high-resolution 3D global nonideal MHD simulations of a gaseous disk threaded by a large-scale vertical magnetic field harboring a planet in a fixed circular orbit using the code IDEFIX, which is accelerated with graphics processing units. We considered various planet masses (10 Earth masses, 1 Saturn mass, 1 Jupiter mass, and 3 Jupiter masses for a solar-mass star) and disk magnetizations ($10^4$ and $10^3$ for the $\beta$-plasma parameter, defined as the ratio of the thermal pressure over the magnetic pressure).

**Results.** We find that a gap opening always occurs for sufficiently massive planets, typically on the order of a few Saturn masses for $\beta_0 = 10^3$, with deeper gaps when the planet mass increases and when the initial magnetization decreases. We propose an expression for the gap-opening criterion when accretion is dominated by MHD winds. We show that accretion is unsteady and comes from surface layers in the outer disk, bringing material directly toward the planet poles. A planet gap is a privileged region for the accumulation of a large-scale magnetic field, preferentially at the gap center or at the gap edges in some cases. This results in a fast accretion stream through the gap, which can become sonic at high magnetizations. The torque due to the MHD wind responds to the planet presence in a way that leads to a more intense wind in the outer gap compared to the inner gap. More precisely, for massive planets, the wind torque is enhanced as it is fed by the planet torque above the gap's outer edge, whereas the wind torque is seemingly diminished above the gap's inner edge due to the planet-induced deflection of magnetic field lines at the disk surface. This induces an asymmetric gap, both in depth and in width, that progressively erodes the outer gap edge, reducing the outer Lindblad torque and potentially reversing the migration direction of Jovian planets in magnetized disks after a few hundreds of orbits. For low-mass planets, we find strongly fluctuating gravitational torques that are mostly positive on average, indicating a stochastic outward migration.

**Conclusions.** The presence of MHD winds strongly affects planet-disk interaction, both in terms of flow kinematics and protoplanet migration. This work illustrates the tight dependence between the planet torque, the wind torque, and magnetic field transport that is required to get the correct dynamics of such systems. In particular, many of the predictions from “effective” models that use parameterized wind torques are not recovered (such as gap formation criteria, the migration direction, and speed) in our simulations.

**Key words.** accretion, accretion disks – protoplanetary disks – planet-disk interactions – magnetohydrodynamics (MHD) – methods: numerical

## 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^{-8} M_{\odot} \text{yr}^{-1}$ (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 $\nu$ with an $\alpha_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 α<sub>v</sub> lies in the range 10<sup>-4</sup>–10<sup>-2</sup>. Moreover, such hydrodynamic models are able to reproduce flows falling from the disk surface toward the disk midplane observed via <sup>12</sup>CO emission at the radial locations of dark rings, as in Teague et al. (2019). Such fast meridional flows (≈ 0.1 c<sub>s</sub>, with c<sub>s</sub> 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 α<sub>v</sub> < 10<sup>-4</sup>–10<sup>-3</sup> 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, differ-
ent from the viscous accretion model. In the MHD wind-driven
accretion model, angular momentum is evacuated from the sur-
face of the disk by a magnetized wind (Blandford & Payne 1982),
which induces a radial laminar transport of the gas mass follow-
ing the idea of Wardle & Koenigl (1993). This model reconciles
at once the high stellar accretion rate, the a priori low turbu-
lence, 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 com-
ing 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 phe-
nomenological models in 1D or 2D. It corresponds to the same
process that made it possible to carry out 2D hydrodynamic sim-
ulations with an α<sub>v</sub> 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 α<sub>v</sub>-like
dimensionless parameter, noted here as α<sub>dw</sub>. Such a prescription
is necessary in order to study long-term disk evolution or to per-
form 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 simpli-
fied 1D models that type-I migration can also be slowed down
and even reversed depending on the wind efficiency, with a pos-
itive 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 negligi-
ble 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<sub>j</sub> and M<sub>s</sub>). 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 mag-
netic 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, 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 pre-
scription, 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<fnref n="1" /> 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<sub>p</sub> to 1.9R<sub>p</sub>, and an outer stretched grid with 128 cells between 1.9R<sub>p</sub> and 6.0R<sub>p</sub> where R<sub>p</sub> is the (fixed) planet location. Concerning the vertical direction, the grid is refined around the midplane with 256 evenly-spaced cells between -0.4R<sub>p</sub> and 0.4R<sub>p</sub>, and two stretched grids symmetric around the midplane with 64 cells on both sides such that the total vertical extension spans between -6.0R<sub>p</sub> and 6.0R<sub>p</sub>. Around the planet location at R = R<sub>p</sub> = 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<sub>R</sub> × N<sub>φ</sub> × N<sub>z</sub> = 640 × 2048 × 384. The aspect ratio h is considered constant, such that h = h<sub>0</sub> = 0.05.

### 2.2. Average conventions

For a quantity Q, we define the azimuthal average, the temporal average between T<sub>1</sub> and T<sub>2</sub>, and the vertical integral between altitudes z<sub>-</sub> and z<sub>+</sub> respectively in Eqs. (1), (2) and (3):

$$ \langle Q \rangle_{\phi} = \frac{1}{2\pi} \int_{0}^{2\pi} Q d\phi, \qquad (1) $$

$$ \frac{\langle Q \rangle_{T_1 - T_2} = \frac{1}{T_2 - T_1}, \int_{T_1}^{T_2} Q dt}{\qquad (2)} $$

$$ \bar{Q} = \int_{z_{-}}^{z_{+}} Q dz. \qquad (3) $$

We make use of the notation $\langle Q \rangle_t$ if we temporally average a quantity over the last 100 orbits of the simulation. We note $\bar{Q}$ the value of Q at the disk midplane. Finally, we define as

$$ [Q]_{-}^{+} = Q|_{z_{w}} - Q|_{-z_{w}}, \qquad (4) $$

which is the difference of Q between the disk's upper layer at z<sub>w</sub> and its lower layer at -z<sub>w</sub>, with z<sub>w</sub> ≈ 4H(R).

### 2.3. Planet

The planet is held on a fixed circular orbit at r<sub>p</sub> = 1 in the disk midplane. Several planet-to-primary mass ratios q<sub>p</sub> are chosen for the planet to explore the planet-disk-wind interaction: 3 × 10<sup>-3</sup>, 10<sup>-3</sup>, 3 × 10<sup>-4</sup> and 3 × 10<sup>-5</sup>, which correspond respectively to M<sub>p</sub> = 3M<sub>j</sub>, M<sub>j</sub>, M<sub>s</sub> and 10M<sub>⊕</sub> for a solar-mass star. We make use of the planet's Hill radius, defined as

$$ R_{\text{hill}} = (q_p/3)^{1/3} R_p. \qquad (5) $$

Accretion on the planet is neglected, as well as planet migration. The mass of the planet is gradually increased over its first N<sub>0</sub> = 10 orbits until it reaches q<sub>p</sub>, according to q<sub>p</sub>(t) = q<sub>p</sub> sin<sup>2</sup>(t/(4N<sub>0</sub>)), with t in units of Ω<sup>-1</sup>. The total gravitational potential Φ considered in the simulation is the following:

$$ \Phi = \Phi_{\star} + \Phi_{p} + \Phi_{\text{ind}}. \qquad (6) $$

We note that Φ<sub>★</sub> is the term due to the central star:

$$ \Phi_{\star}(\mathbf{r}) = -\frac{1}{|\mathbf{r}|}, \qquad (7) $$

Φ<sub>p</sub> is the gravitational potential of the planet:

$$ \Phi_{p}(\mathbf{r}) = -\frac{q_{p}}{\sqrt{|\mathbf{r} - \mathbf{r}_{p}|^{2} + \epsilon^{2}}}, \qquad (8) $$

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

$$ \Phi_{\text{ind}}(\mathbf{r}) = \frac{q_{p}(\mathbf{r} \cdot \mathbf{r}_{p})}{|\mathbf{r}_{p}|^{3}}. \qquad (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.

<footnote n="1">¹ https://github.com/idefix-code/idefix</footnote>
### 2.4.1. Nonideal MHD

We solve the nonideal MHD equations, for the density $\rho$, the velocity field $v$ and the magnetic field $B$:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho v) = 0, \qquad (10) $$

$$ \frac{\partial \rho v}{\partial t} + \nabla \cdot (\rho v \otimes v) = -\rho \nabla \Phi - \nabla P + J \times B, \qquad (11) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{B}) + \nabla \times [\eta_A (\mathbf{J} \times \mathbf{e}_b) \times \mathbf{e}_b - \eta_O \mathbf{J}], \qquad (12) $$

where $\Phi$ is the total gravitational potential defined in Eqs. (6–9), $P$ the gas pressure, $\mathbf{J} = \nabla \times \mathbf{B}$ is the current density, $\mathbf{e}_b$ the unit vector parallel to the magnetic field line, $\eta_A$ the ambipolar diffusivity and $\eta_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

$$ A_m = \frac{v_A^2}{\Omega_K \eta_A}, \qquad (13) $$

with $v_A = B/\sqrt{\rho}$ the Alfvén speed and $\Omega_K$ the Keplerian angular frequency. To keep the model as simple as possible, we follow Lesur (2021) 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 $\tilde{A}_m$ 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

$$ A_m = \max \left\{ \tilde{A}_m \exp \left[ \left( \frac{z}{\epsilon_{\text{id}} H(R)} \right)^4 \right], \frac{1}{10} \left( \frac{v_A}{c_s} \right)^2 \right\}, \qquad (14) $$

with $\epsilon_{\text{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 $\Sigma$ is not taken into account here. For such $\tilde{A}_m$, 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

$$ B_z(R) = h_0 \sqrt{\frac{2}{\beta_0}} R^{-5/4}, \qquad (15) $$

with $\beta_0$ the initial value of the plasma parameter $\beta$, defined as the ratio of the thermal pressure $P$ over the magnetic pressure:

$$ \beta = \frac{2P}{B_z^2}. \qquad (16) $$

Two $\beta_0$ are chosen in order to explore the influence of the magnetization on the interaction between the planet and the disk: $10^4$ and $10^3$.

Our numerical parameters are chosen to be close to more complete thermo-chemical models. For instance, a constant $\tilde{A}_m \sim 1$ matches the thermo-chemical model of Thi et al. (2019) with $\tilde{\beta} = 10^4$ 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 $\gg 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:

$$ \rho(R, z) = \frac{\Sigma(R)}{\sqrt{2\pi}H(R)} \exp \left[ \frac{1}{\tilde{c}_s^2(R)} \left( \frac{1}{\sqrt{R^2 + z^2}} - \frac{1}{R} \right) \right], \qquad (17) $$

with $\tilde{c}_s = h_0/\sqrt{R}$ the sound speed at disk midplane and $\Sigma$ the surface density following

$$ \Sigma(R) = \Sigma_0 R^{-0.5}, \qquad (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 $\tilde{T}(R) = h_0^2/R$ and a low-density hot corona with temperature $T_{\text{cor}}(R) = 36 \tilde{T}(R)$ 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:

$$ T_{\text{eff}}(R, z) = \frac{1}{2} (T_{\text{cor}}(R) + \tilde{T}(R)) + \frac{1}{2} (T_{\text{cor}}(R) - \tilde{T}(R)) \tanh \left( 6 \ln \left| \frac{z}{\epsilon_{\text{id}} H(R)} \right| \right). \qquad (19) $$

We use a fast $\mathcal{B}$-cooling prescription, with $\mathcal{B} = 0.1\Omega^{-1}$, in order to quickly relax the temperature $T(R, z)$ toward the prescribed profile $T_{\text{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) = \rho(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 $e(R, z) = \frac{P(R, z)}{\gamma - 1}$, and an adiabatic index $\gamma = 5/3$. Regarding the initial gas azimuthal velocity, we use the following sub-Keplerian velocity profile

$$ v_{\phi}(R, z) = v_{\kappa} \left( \frac{R}{\sqrt{R^2 + z^2}} - 1.5h_0^2 \right), \qquad (20) $$
with $v_K = R\Omega_K = 1/\sqrt{R}$ 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_\phi$ and $B_z$ from the innermost active zone to the ghost zone. We impose $v_z = 0$ and a Keplerian velocity for $v_\phi$. The boundary condition is reflecting (symmetric) for $v_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 $\rho, P, v_R, v_\phi, v_z, B_\phi$ 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 $v_A > v_K(R_{in})$ is encountered, typically in the disk corona. We also add a threshold in density at $\rho_{MIN} = 10^{-12}$. 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 ($\rho v_R, \rho v_z$). In this disk innermost region, we also add a band of resistivity such that $\eta_O = H^2(R_{in})\Omega_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 $\beta_0$ 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 $\beta_0$-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 $\beta_0 = 10^3$ 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$.

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

<table><caption>**Table 1.** Simulations, characterized by their name, the planet mass, and the $\beta$ parameter.</caption><thead><tr><th>Name</th><th>M<sub>p</sub></th><th>β</th><th>N<sub>orb</sub></th></tr></thead><tbody><tr><td>3Mj-β<sub>3</sub></td><td>3M<sub>j</sub></td><td>10<sup>3</sup></td><td>200</td></tr><tr><td>3Mj-β<sub>4</sub></td><td>3M<sub>j</sub></td><td>10<sup>4</sup></td><td>200</td></tr><tr><td>Mj-β<sub>3</sub></td><td>M<sub>j</sub></td><td>10<sup>3</sup></td><td>400</td></tr><tr><td>Mj-β<sub>4</sub></td><td>M<sub>j</sub></td><td>10<sup>4</sup></td><td>200</td></tr><tr><td>Ms-β<sub>3</sub></td><td>M<sub>s</sub></td><td>10<sup>3</sup></td><td>200</td></tr><tr><td>Ms-β<sub>4</sub></td><td>M<sub>s</sub></td><td>10<sup>4</sup></td><td>400</td></tr><tr><td>10Me-β<sub>3</sub></td><td>10M<sub>⊕</sub></td><td>10<sup>3</sup></td><td>200</td></tr><tr><td>10Me-β<sub>4</sub></td><td>10M<sub>⊕</sub></td><td>10<sup>4</sup></td><td>200</td></tr><tr><td>3Mj-α<sub>0</sub></td><td>3M<sub>j</sub></td><td>∞</td><td>200</td></tr><tr><td>3D-β<sub>3</sub></td><td>0</td><td>10<sup>3</sup></td><td>100</td></tr><tr><td>axi-β<sub>3</sub></td><td>–</td><td>10<sup>3</sup></td><td>200</td></tr><tr><td>axi-β<sub>4</sub></td><td>–</td><td>10<sup>4</sup></td><td>200</td></tr></tbody></table>

**Notes.** Each simulation was run during $N_{orb}$ orbital periods.

define the accretion rate as

$$ \dot{M}_{\text{acc}} = -2\pi R \overline{\langle \rho v_R \rangle}_{\phi}, \quad (21) $$

Using the equation of angular momentum conservation, $\dot{M}_{acc}$ is linked to the transport of angular momentum by radial and surface stresses, respectively noted $W_{R\phi}$ and $W_{z\phi}$:

$$ \begin{aligned} \frac{\dot{M}_{\text{acc}}\Omega_K}{4\pi} &= \frac{1}{R}\partial_R(R^2\overline{\langle W_{R\phi} \rangle}_{\phi}) + R[\overline{\langle W_{z\phi} \rangle}_{\phi}]^+ - \overline{\langle \Gamma_{p \to d} \rangle}_{\phi}, \\ &= \underbrace{R_R + M_R + R_{\text{surf}} + M_{\text{surf}}}_{\mathcal{T}} + \mathcal{P}, \end{aligned} \quad (22) $$

with $W_{X\phi} = \rho v_X(v_\phi - v_K) - B_X B_\phi$ the $\phi$ component of the stress, and $\mathcal{P} = -\overline{\langle \Gamma_{p \to d} \rangle}_{\phi}$. $\Gamma_{p \to d} = -\rho (\mathbf{r}_p \times \nabla \Phi_p) \cdot \mathbf{e}_z$ is the torque per unit of volume exerted by the planet on the gas, with $\Phi_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 $\mathcal{T}$, and therefore to $\dot{M}_{acc}$: the radial Reynolds ($R_R$) and Maxwell ($M_R$) torques, and the surface Reynolds ($R_{surf}$) and Maxwell ($M_{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-\beta_3$, 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 $\beta_0 = 10^3$. On the contrary, the surface Reynolds torque $R_{surf}$ (green curve) has little impact on the accretion mechanism. Concerning the radial torques ($R_R$ in blue and $M_R$ in red), they have a similar influence on the accretion, and can be comparable to $M_{surf}$ if combined. In summary, $M_R \approx R_R$, $M_{surf} \gg R_{surf}$ and $M_{surf} + R_{surf} \gtrsim M_R + 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 $M_{surf}$ term, we define the dimensionless parameter $\nu$ as in Lesur (2021):

$$ \nu = - \frac{[\langle B_{\phi} B_{z} \rangle_{\phi}]^{+}}{\langle \Sigma \rangle_{\phi} \Omega_{K}^{2} H} = \frac{M_{\text{surf}}}{h_{0} v_{K}^{2} \langle \Sigma \rangle_{\phi}}, \quad (23) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Radial profile of the four torques $R_R$ (blue), $M_R$ (red), $R_{surf}$ (green), $M_{surf}$ (yellow) and their sum $\mathcal{T}$ (black) normalized by $h_0 v_K^2 \langle\Sigma\rangle_{\phi,t}$, azimuthally averaged and temporally averaged over the last 100 orbits for the run $M_j-\beta_3$. We consider here only the logarithm of the absolute value of these quantities.</figcaption></figure>

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, $v$ is positive when the wind extracts angular momentum from the disk. We note that for $\beta_0 = 10^4$, that is at lower magnetization, assuming that $R_{surf}$ is still negligible, $R_R$ is actually the dominant term in the internal parts of the disk ($R < 3$, where $R_R/M_{surf} \le 10$), whereas $M_{surf}$ takes over in the external parts of the disk. It is safe to neglect $R_{surf}$ here as it is proportional to the density, which is small at the disk surface $z = 4H$. We also define the ejection efficiency dimensionless parameter $\xi$ as:

$$ \xi = \frac{2\pi R^2 [\langle \rho v_z \rangle_{\phi}]_{-}^{+}}{\dot{M}_{\text{acc}}}, \qquad (24) $$

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

$$ \alpha_v = \frac{\overline{\langle W_{R\phi} \rangle_{\phi}}}{\overline{\langle P \rangle_{\phi}}} = \alpha_v^R + \alpha_v^M. \qquad (25) $$

where $\alpha_v^R$ and $\alpha_v^M$ are respectively related to the Reynolds and Maxwell components of the radial stress. We also define a dimensionless parameter $\alpha_{dw}$, similar to the one introduced in Tabone et al. (2022), and which we can divide in Reynolds ($\alpha_{dw}^R$) and Maxwell ($\alpha_{dw}^M$) contributions:

$$ \alpha_{dw} = \frac{R \left[ \langle W_{z\phi} \rangle_{\phi} \right]_{-}^{+}}{\overline{\langle P \rangle_{\phi}}} = \alpha_{dw}^{R} + \alpha_{dw}^{M}. \qquad (26) $$

We note that $v \approx h_0 \alpha_{dw}^M \approx h_0 \alpha_{dw}$, using Eq. (23), $\overline{\langle P \rangle}_{\phi} \approx \langle \Sigma \rangle_{\phi} \tilde{c}_s^2$ and $\alpha_{dw}^R \ll \alpha_{dw}^M$.

In Appendix B.2, we address the radial transport efficiency and the level of turbulence in the 3D planet-free run 3D-$\beta_3$, temporally averaged between 50 and 100 orbits. In particular, we compare $\alpha_v$ but also the turbulent component of $\alpha_v^R$ with their corresponding values in the run $10Me-\beta_3$. 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.

## 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 $\beta_0$ 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 \ge 3 \times 10^{-4}$ ($M_s$). In terms of $q_p/h_p^3$ parameter, this is coherent with the classical condition for non-linear effects in the flow whenever $q_p/h_p^3 \le 1$, although this condition also depends on the level of turbulence $\alpha_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 $3M_j$ and the $M_s$ cases. For a given $M_p \ge M_s$, the gap is denser when $\beta_0$ is larger, with a factor on the order of 5. For all $\beta_0$, when $M_p \ge 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 \ge M_j$) and highly magnetized disk ($\beta_0 = 10^3$), 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 $M_j-\beta_3$, 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”
<figure><img src="image_3.png" /><figcaption>**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 $\beta_0$ parameter. The white circle marks the planet location.</figcaption></figure>

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-$\beta_3$ and 10Me-$\beta_3$ 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_\oplus$ 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 $(\Sigma - \Sigma_{\text{th}})/\Sigma_{\text{th}}$ are larger for smaller $\beta_0$. This

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

is particularly visible in the top row of Fig. 2, with ~80% in amplitude for the case 10Me-$\beta_3$ and ~20% in amplitude for the case 10Me-$\beta_4$, compared to a simple theoretical power-law profile in $\Sigma_{\text{th}} = \Sigma_0 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-$\beta_3$), the gas perturbations are dominated by self-organization. For the intermediate mass case (Ms-$\beta_3$), the competition between self-organization and the planet tends to diminish the contrast in density. For larger planet masses (Mj-$\beta_3$, 3Mj-$\beta_3$), 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-$\beta_3$, 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 $3M_j\beta_3$. 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 \approx 1.3$ for the run $10Me\beta_4$.

In Fig. 2, there is evidence that an evacuated region is generated by the inner boundary condition, especially when $\beta = 10^3$. 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.

### 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 $\Sigma$ for the eight nonideal MHD runs. We also show the space-time diagram of the inviscid hydrodynamic simulation ($3M_j\alpha_0$) for comparison purpose. These diagrams confirm that a planet gap forms when $M_p \ge 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\beta_3$ with the under-density that appears near $R = 1.2$.

The planet gap is denser when the initial magnetization increases ($\beta_0$ 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 inter-dependence between the radial profiles of the accretion rate $\dot{M}_{acc}$, $\beta$ and $\Sigma$, by comparing their values in the outer disk (e.g., $\beta_{disk} \approx \beta_0$) and their counterpart in the outer gap (e.g., $\beta_{gap}$). Firstly, $\dot{M}_{acc}$ is known to be a function of $\Sigma\beta^{-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 $\beta_{gap} \sim 10$ (see bottom panel of Fig. C.1). As a result, in the gap, $\dot{M}_{acc} \propto \Sigma_{gap}$. If we make the assumption of a quasi steady-state, which implies that $\dot{M}_{acc}$ in the gap matches $\dot{M}_{acc}$ in the disk, then $\Sigma_{gap}/\Sigma_{disk} \propto \beta_{disk}^{-b} \approx \beta_0^{-b}$.

Along with the increase in density, we detect episodes of activity in the gap when increasing the magnetization and for $M_p \ge 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 $3M_j\beta_3$ and $M_j\beta_3$, 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 $M_{surf}$), tends to help the gap opening until $\beta$ 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 $v_{acc}$ to cross the Hill sphere of the planet $\tau_{acc} = R_{hill}/v_{acc}$ is longer than the time it takes for the planet to clear up the Hill sphere $\tau_p \approx R_p/(R_{hill}\Omega)$. 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}\Omega$, therefore the typical timescale to clean a full ring at the planet location $R_p$ is approximately given by $\tau_p$. Said differently, the characteristic emptying timescale of the Hill sphere is linked to half the libration time, which is given by $\tau_{hill}^{\text{lib}}/2 = 2\pi R_p/(3/2)\Omega R_{hill}$ 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. $\tau_{hill}^{\text{lib}}$ scales directly with $\tau_p$, hence the relevance of this quantity.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Space-time diagrams showing the evolution of the logarithm of the gas surface density $\Sigma$ 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 $\mathcal{E}_1$ and $\mathcal{E}_2$ two episodes in the run $M_j-\beta_3$ that will be discussed in Sect. 3.2.3.</figcaption></figure>

With these two timescales, the criterion for gap opening reads

$$ \tau_p/\tau_{\text{acc}} < 1. \qquad (27) $$

Using the equation right before Eq. (17) in Lesur (2021), we can link the accretion speed $v_{\text{acc}}$ to the normalized magnetic wind torque $\nu$. We get $v_{\text{acc}} = 2R_p\Omega h\nu$. The gap-opening criterion then becomes

$$ \frac{2R^2 h \nu}{R_{\text{hill}}^2} < 1. \qquad (28) $$

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

$$ 2h\nu(q_p/3)^{-2/3} < 1. \qquad (29) $$

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

$$ h(q_p/3)^{-1/3} < 1. \qquad (30) $$

Using Eqs. (29) and (30), we can eventually get a combined criterion for gap opening as in CR6, when accretion is dominated
<figure><img src="image_6.png" /><figcaption>Fig. 5. Gap-opening criterion with MHD winds displayed in the ($q_p, \beta_0$) domain. The dashed line indicates the minimum mass to open a gap from this criterion, using a scaling law between $v$ and $\beta_0$. Planets in the hatched area are able to open a gap. To obtain the pink box plots, we estimated $v$ in our low-mass planet simulations.</figcaption></figure>

by MHD winds

$$ \frac{3h}{4(q_p/3)^{1/3}} + \frac{2hv}{(q_p/3)^{2/3}} < 1. \qquad (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 $\tau_{acc}$ by $\tau_v = R_{hill}^2/(C\nu)$, with $C$ a constant factor and rewrite the gap-opening criterion of Eq. (31). By considering $\mathcal{R} = R_p^2\Omega/\nu$ the Reynolds number, Eq. (29) would become $3C/(q_p\mathcal{R}) < 1$ in the fully viscous case. We therefore retrieve the same scaling $\propto 1/q_p\mathcal{R}$ 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 $q_p = q_p^C$ to open a gap in a MHD wind-driven accretion disk, with $h > 0$ and $v > 0$:

$$ q_p^C = \frac{192(hv)^3}{\left(\sqrt{\frac{9h^2}{16} + 8hv} - \frac{3h}{4}\right)^3}. \qquad (32) $$

In order to make this expression useful, one must estimate the wind torque coefficient $v$. Here, we propose two approaches: first using the self-similar solutions from Lesur (2021) with $\tilde{A}_m = 1$, giving $v = 3\beta_0^{-0.78}$. In this case our gap-opening criterion becomes essentially a function of $\beta_0$ and $h$. The second approach is to measure $v$ in our simulations with a $10M_\oplus$ planet, in the radial interval $[R_p - 2R_{hill}^{\text{sup}}, R_p + 2R_{hill}^{\text{sup}}]$ 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 = 0.05$. The hatched area above the dashed line corresponds to Eq. (31) with the self-similar scaling for $v$ (Lesur 2021) while the pink box plots indicate the values critical $q_p^C$ estimated from $v$ measured in our simulations. The 8 colored dots correspond to the 8 ($q_p, \beta_0$) pairs in our nonideal MHD simulations, with the $3M_j, M_j, M_s$ and $10M_\oplus$ 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 $\beta_0$). More precisely, the minimum mass for a planet to open a gap lies between a mass of Saturn and a mass of Jupiter for $\beta_0 = 10^3$, whereas it is inferior to a mass of Saturn for $\beta_0 = 10^4$.

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 $\beta$ decreases also in the gap, $v$ 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).

## 3.2. Case study: high planet mass and high magnetization runs

We focus now mainly on the simulations with $M_p \ge M_j$ and $\beta_0 = 10^3$.

### 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 $\langle\rho v_p\rangle_{\phi,t}$ in log scale, averaged azimuthally and temporally over the last 100 orbits for the run $3M_j-\beta_4$. The texture in LIC (line integral convolution) and the white arrows indicate the direction of the radial $(\langle\rho v_R\rangle_{\phi,t})$ and vertical $(\langle\rho v_z\rangle_{\phi,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 \in [0.9, 1]$ and $|z| \in [0.0, 0.1]$ than in the outer gap when $R \in [1, 1.1]$ and $|z| \in [0.0, 0.02]$. This asymmetry is due to the two accretion layers in $R \in [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 $\beta_0 = 10^4$. This behavior is the result of an enhanced $M_{\text{surf}}$ component (magnetic torque at the disk's surface) and to a lesser extent to an enhanced radial torque ($\mathcal{R}_R + M_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
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Meridional flows around the planet for the run 3*Mj*-*β*₄, 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 ±z<sub>w</sub> ≈ ±4H(R) are indicated by the black dotted lines. The planet is indicated by the white circle.</figcaption></figure>

planet – as well as planet growth – should therefore be more effi-
cient because of this denser inflow coming from lower altitudes
toward the planet poles, as observed in Gressel et al. (2013).

When $M_p \ge M_j$ and $\beta_0 = 10^3$, 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 \in [-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 $\langle \rho v_p \rangle_{\phi,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 $\beta$ much smaller than $\beta_0$ (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 lay-
ers 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 verti-
cally 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.

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 ⟨B̃z⟩_φ. 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
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** 2D map of the Mach number, defined as $\sqrt{\langle\rho v_R\rangle_{\phi,t}^2 + \langle\rho v_z\rangle_{\phi,t}^2 / \langle\rho\rangle_{\phi,t} / \langle c_s\rangle_{\phi,t}}$, in the $(R-z)$ plane for the run $Mj-\beta_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.</figcaption></figure>

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 $\langle \tilde{B}_z \rangle_\phi$ reaching its highest value. If we focus on $\log\langle\Sigma\rangle_\phi$, this strong magnetic accumulation is concomitant to a strong episode of gas depletion and a fast outward drift of the outer gap (episode noted $\mathcal{E}_1$ in Sect. 3.2.3). This $\langle \tilde{B}_z \rangle_\phi$ 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 \approx 0.85$), and the other one at the transition between the outer gap and the outer disk (near $R \approx 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 $\alpha_v$), larger close to the gap edges than at the gap center. In addition, they argued that the radial band of large $\langle \tilde{B}_z \rangle_\phi$ 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 $\langle \tilde{B}_z \rangle_\phi$ at 400 orbits instead of ≈0.55 code units for $\langle \Sigma \rangle_\phi$, 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 $\langle \tilde{B}_z \rangle_\phi$ in time, with episodes of strong magnetic field accumulation (e.g., between 310 and 350 orbits) alternating with episodes of weaker $\langle \tilde{B}_z \rangle_\phi$ 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 $\langle \tilde{B}_z \rangle_\phi$ 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-\beta_3$. 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 $\langle B_p \rangle_{\phi,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
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Space-time diagram of the vertical component of the magnetic field in the midplane $\langle \tilde{B}_z \rangle_\phi$, for the runs $Mj-\beta_4$ (top panel) and $Mj-\beta_3$ (bottom panel). The horizontal dashed lines indicate the location of the planet. The two episodes $E_1$ and $E_2$ are indicated in white.</figcaption></figure>

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-\beta_4$. 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 $v_\phi - v_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 $\sqrt{q_p} \approx 0.032$, 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 \approx 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.

### 3.2.3. Accretion rate, wind torque, and gap asymmetries

In this section, we focus on the case $Mj-\beta_3$. We look at the link between the accretion rate $\dot{M}_{acc}$ defined in Eq. (21), the specific
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Poloidal magnetic field in the gap for the run *Mj*-β₃, averaged in azimuth and over the last 100 orbits (episode *E*₂). 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* ≈ *R*<sub>p</sub> − 2*R*<sub>hill</sub> and *R* ≈ *R*<sub>p</sub> + 5*R*<sub>hill</sub>. The upper and lower surfaces ±*z*<sub>w</sub> ≈ ±4*H*(*R*) are indicated by the black dotted lines.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Azimuthal cut of the log (Σ), at the azimuth of the planet φ<sub>p</sub>, at ≈114 orbits for 3*Mj*-β<sub>4</sub>. A protoplanetary jet is clearly visible, launching from the planet poles.</figcaption></figure>

wind torques via *M*<sub>surf</sub> and the *v* 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 *E*₁, corresponds to the strong episode of gas depletion between 165 and 200 orbits, the second one, *E*₂, corresponds

to the last 100 orbits of the simulation, after the magnetic field has been divided in two regions of accumulation. During *E*₁, the gap in the density lies approximately between *R*<sub>p</sub> − 2*R*<sub>hill</sub> and *R*<sub>p</sub> + 3*R*<sub>hill</sub>, whereas during *E*₂, the outer gap edge is further out in the disk, near *R*<sub>p</sub> + 5*R*<sub>hill</sub>. In the graphs of this section, we highlight these three characteristic radii with vertical red dashed lines.

Figure 11 displays the radial profiles of $\dot{M}_{\text{acc}}$ (top panel) and the ejection efficiency $\xi$ (bottom panel) for $\mathcal{E}_1$ (dashed line) and $\mathcal{E}_2$ (solid line). For these two episodes, the curves of $\dot{M}_{\text{acc}}$ and $\xi$ 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 \approx 0.7$ and $R \approx 1.3-1.5$), and even stronger at higher magnetization. In addition, the ejection efficiency $\xi$ is very small in the gap<sup>2</sup> (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 $\dot{M}_{\text{acc}}$ and the fact that $\xi \ll 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, $\dot{M}_{\text{acc}}$ is larger in the gap during $\mathcal{E}_1$ than during $\mathcal{E}_2$. Actually,

<footnote n="2">² We note that beyond the gap on both sides, the ejection in the wind is particularly efficient ($\xi > 1$) where the accretion rate is minimum.</footnote>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Radial profile of the accretion rate $\dot{M}_{\text{acc}}$ (top panel) and the ejection parameter $\xi$ (bottom panel), averaged azimuthally and temporally over the last 100 orbits (solid line) in the run $M_j-\beta_3 \in \mathcal{E}_2$ and over 25 orbits during the episode $\mathcal{E}_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 $\dot{M}_{\text{acc}}$ and $\xi$.</figcaption></figure>

the accretion rate reaches its highest value in the gap during $\mathcal{E}_1$, which is related to the saturation of the $B_z$ everywhere in the gap during this epoch (see Fig. 8). Moreover, the step of $\dot{M}_{\text{acc}}$ in the gap and out of the gap is bigger during $\mathcal{E}_1$, with a steeper gradient (see in particular between $R = R_p + 3R_{\text{hill}}$ and $R \approx 1.35$ for $\mathcal{E}_1$ and between $R = R_p + 5R_{\text{hill}}$ and $R \approx 1.5$ for $\mathcal{E}_2$). We therefore expect a faster outward drift of the outer gap during $\mathcal{E}_1$ than during $\mathcal{E}_2$, which is indeed what we observe in the space-time diagram of $\log\langle\Sigma\rangle_{\phi}$ in the sixth panel of Fig. 4.

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Radial profile of the torques exerted by the wind ($M_{\text{surf}}$, yellow curve) and the sum of the radial torques and planet torques ($R_R$, $M_R$ and $\mathcal{P}$, blue curve), averaged over the episode $\mathcal{E}_2$ in the run $M_j-\beta_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.</figcaption></figure>

Furthermore, $\dot{M}_{\text{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 $\mathcal{P}$ 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 + M_R$. We therefore gather these three torques into a planet-induced torque $\mathcal{T}_R = R_R + M_R + \mathcal{P}$. On the other hand, the magnetic wind torque $M_{\text{surf}}$ is the main driver of accretion (see Fig. 1). In Fig. 12, we show the radial profiles in the gap of $\mathcal{T}_R/\Omega_K$ (blue curve), $M_{\text{surf}}/\Omega_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 $\mathcal{E}_2$). We note that the total torque $(\mathcal{T} + \mathcal{P})/\Omega_K$ is close to the accretion rate $\dot{M}_{\text{acc}}/4\pi$ 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 $\mathcal{T}_R/\Omega_K$ and $M_{\text{surf}}/\Omega_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 $\mathcal{T}_R$, we can divide it in $\mathcal{T}_{R-}$ and $\mathcal{T}_{R+}$, respectively for the planet-related torques in the inner gap and the outer gap. $\mathcal{T}_{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, $\mathcal{T}_{R-}$ is quite small and positive,
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** 2D map of $v$, i.e., the magnetic torque per unit of surface density, at the disk surface in $(R-\phi)$ coordinates, averaged over the episode $E_2$ in the run $M_j-\beta_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 $v$ during $E_2$ and $E_1$.</figcaption></figure>

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 ($\mathcal{T}_{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 ($\mathcal{M}_{\text{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 ($\mathcal{T}_{R-}$ positive but close to zero), and the wind torque is also negligible at this location ($\mathcal{M}_{\text{surf}} \gtrsim 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 $v$ parameter, which is a dimensionless measure of $\mathcal{M}_{\text{surf}}$ (Eq. (23)). This plot shows the 2D map of $v$ in $(R-\phi)$ coordinates for $E_2$, as well as its azimuthally averaged version (solid black line). We also show the azimuthally averaged $v$ for $E_1$ (dashed black line). During both episodes, $v$ is larger in the outer gap than in the inner gap, with an asymmetry of the angular momentum extraction $\Delta v \approx 0.2$ in both cases. In the $(R-\phi)$ plane, we detect a clear azimuthal asymmetry of $v$, with a larger value at $\phi > 0$. The wind torque seems affected by the outer planet's spiral wake, with an enhanced $v$ along the wake toward the inner gap near the CPD, from $R = 1.25$ to $R = 0.95$. This increased $v$ 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 $[\langle B_\phi B_z \rangle_t]^+_{-}$ 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 $[\langle B_\phi B_z \rangle_t]^+_{-}$ (see in Fig. 9 the LIC texture and magnetic field lines in black that are nearly tangential to the surface layers at $\pm z_w \approx \pm 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 $\mathcal{M}_{\text{surf}}$ or $[\langle B_\phi B_z \rangle_t]^+_{-}$ 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 $\pm z_w \approx \pm 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 $^{13}$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 $^{18}$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 $B_p$ 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.

### 3.3. Gravitational torques

We discuss here the evaluation of the gravitational torque $\Gamma$ 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

$$ \Gamma = - \sum_{i,j,k}^{N_R,N_\phi,N_\zeta} \Gamma_{p \to d}^{i,j,k} V_{\text{cell}}^{i,j,k} \quad (33) $$

with $\Gamma_{p \to d}^{i,j,k}$ the torque per unit of volume defined in Sect. 2.6 and evaluated for the cell $i, j, k$ of volume $V_{\text{cell}}^{i,j,k}$. 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
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Temporal evolution of the gravitational torques exerted by the gas onto the planets for the nine runs. The torques are normalized by $\Gamma_0 = q_p^2 \Sigma_0 R_p^4 \Omega_K (R_p)^2 h_0^{-2}$. The different colors show different planet masses: $10 M_\oplus$ (green), $M_s$ (yellow), $M_j$ (red) and $3 M_j$ (blue). The different line-styles indicate the initial magnetization: $\beta_0 = +\infty$ (inviscid hydrodynamic run, dotted line), $\beta_0 = 10^4$ (solid line) and $\beta_0 = 10^3$ (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.</figcaption></figure>

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 $\Gamma$ for our 8 nonideal MHD simulations and the inviscid hydrodynamic simulation $3M_j-\alpha_0$, with a moving average of 10 orbits and normalized by $\Gamma_0 = q_p^2 \Sigma_0 R_p^4 \Omega_K (R_p)^2 h_0^{-2}$. Each color corresponds to a planet mass $M_p$, and the line-style indicates the value of $\beta_0$.

We first notice that $\Gamma$ is negative for the three most massive planets, with a smaller $|\Gamma|$ for larger planet masses. It suggests a slower inward migration for more massive planets. We observe indeed a factor $\approx 10$ for $|\Gamma|$ 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 coorbital 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 $10 M_\oplus$), we show that $\beta_0$ 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 $\beta_0$ decreases), which is especially visible between $\approx 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 $M_j-\beta_3$, 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 $|\Gamma_+|$ decreases more than the total inner torque $|\Gamma_-|$. We also noticed that the fluctuations of $|\Gamma_+|$ and $|\Gamma_-|$ are at first equivalent until $|\Gamma_-|$ overcomes $|\Gamma_+$ near $\approx 280$ orbits. At that time, the inner torque starts to strongly fluctuate, with fluctuations $\approx 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 $\approx 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 $\Gamma$ 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
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Gas surface density log (Σ), averaged over the last 100 orbits, for *Mj*-*β*₃. The streamlines in the midplane are shown in black.</figcaption></figure>

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 (10*M*<sub>⊕</sub>) 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 10*Me*-*β*₄, 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 libration 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*<sub>p</sub> ≥ *M*<sub>i</sub>) 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.
# 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 $β_0$ 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 $v$ (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-β_3$ with an episode (episode $E_1$) 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 $β_0 = 10^3$, 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 $α_{dw} ≈ α_{dw}^M ≈ 5$ in the gap for the run $Mj-β_3$ averaged over the last 100 orbits (episode $E_2$). 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 \lesssim 10^{-2}$ in the cases $β_0 = 10^4$, and $α_v > 10^{-2}$ when $β_0 = 10^3$ 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 \ge M_j$) and high magnetization ($β_0 = 10^3$) runs. For weaker magnetizations ($β_0 = 10^4$), 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 $\dot{M}_{acc}$ is nearly constant and comparable that of the disk in the gap. Furthermore, at high magnetization ($β_0 = 10^3$), 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 horse-shoe 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 $\dot{M}_{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 $\alpha_v$ and $\alpha_{dw}$. In our cases where $\beta_0 = 10^4$, with the exact definitions of $\alpha_v$ and $\alpha_{dw}$ in Eqs. (4) and (6) by Tabone et al. (2022), neglecting $R_{surf}$ and averaged in time over the last 100 orbits), we find outside the planet gap that $\alpha_v \simeq 10^{-3} - 10^{-2}$ and $\alpha_{dw} \simeq 4 \times 10^{-2}$. This corresponds to a minimum pessimistic (concerning their model) planet mass needed to open a gap of $\simeq 1M_j$, but our run $Ms-\beta_4$ is able to produce a gap. Likewise, when $\beta_0 = 10^3$ and outside the planet gap, $\alpha_v \simeq 3 \times 10^{-3} - 6 \times 10^{-2}$ and $\alpha_{dw} \simeq 6 \times 10^{-2} - 0.9$ (also averaged in time over the last 100 orbits). When $\beta_0 = 10^3$, 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 $\simeq 2M_j$, but here again, the run $M_j-\beta_3$ is able to produce a gap. One major difference with Elbakyan et al. (2022) is that $\alpha_v$ and $\alpha_{dw}$ strongly depend on the magnetic field. Because the distribution of magnetic field relies on self-organized gaps and planet gaps, $\alpha_v$ and $\alpha_{dw}$ also depend on the localization in the disk. One aspect that could be improved when prescribing the wind torque is to use a $\alpha_{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 $q_p^C$ 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 $\lambda$ and constant mass-loss parameter $b$. This latter quantity is related to the vertical mass flux $\dot{\Sigma}_{wind} = [\langle\rho v_z\rangle_{\phi}]_+^+$. However, if we compute $\dot{\Sigma}_{wind}$ 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 $(M_j-\beta_3)$, $b$ reaches a few $10^{-3}$ in the disk, slightly increasing radially, but an order of magnitude higher in the gap during the episode $E_2$ (i.e., during the last 100 orbits), reaching $1.75 \times 10^{-2}$. 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^{-3}$) and in the other case a fast runaway outward migration ($b \simeq 10^{-2}$). The lever arm $\lambda$ 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, $\lambda$ is larger than 2.5, reaching $\simeq 7.5$ in some gaps. Kimmig et al. (2020) chose a constant lever arm $\lambda = 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_\oplus$) confirm that surface accreting layers and the midplane seem decoupled regarding accretion, even if we detect highly variable gravitational torques $\Gamma$ 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, blind (https://github.com/volodia99/blind) for most of the column maps 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.

## References

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41

Aoyama, Y., & Bai, X.-N. 2023, ApJ, 946, 5

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

Bae, J., Zhu, Z., & Hartmann, L. 2016, ApJ, 819, 134

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650

Baruteau, C., & Masset, F. 2013, in Lecture Notes in Physics, 861 (Berlin: Springer Verlag), eds. J. Souchay, S. Mathis, & T. Tokieda, 201

Baruteau, C., Fromang, S., Nelson, R. P., & Masset, F. 2011, A&A, 533, A84

Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 667

Baruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304

Béthune, W., & Latter, H. 2022, A&A, 663, A138

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Carballido, A., Matthews, L. S., & Hyde, T. W. 2017, MNRAS, 472, 3277

Combet, C., & Ferreira, J. 2008, A&A, 479, 481

Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587

Cui, C., & Bai, X.-N. 2021, MNRAS, 507, 1106

Cui, C., & Bai, X.-N. 2022, MNRAS, 516, 4660

de Valon, A., Dougados, C., Cabrit, S., et al. 2020, A&A, 634, L12

Elbakyan, V., Wu, Y., Nayakshin, S., & Rosotti, G. 2022, MNRAS, 515, 3113

Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 451

Fung, J., & Chiang, E. 2016, ApJ, 832, 105

Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88

Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509

Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59

Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097

Guilet, J., & Ogilvie, G. I. 2013, MNRAS, 430, 822

Guilet, J., & Ogilvie, G. I. 2014, MNRAS, 441, 852

Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L48

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49

Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021, A&A, 650, A179

Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2

Keith, S. L., & Wardle, M. 2015, MNRAS, 451, 1104

Kim, S. Y., & Turner, N. J. 2020, ApJ, 889, 159
Kimmig, C. N., Dullemond, C. P., & Kley, W. 2020, A&A, 633, A4

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Lega, E., Nelson, R. P., Morbidelli, A., et al. 2021, A&A, 646, A166

Lega, E., Morbidelli, A., Nelson, R. P., et al. 2022, A&A, 658, A32

Lesur, G. R. J. 2021, A&A, 650, A35

Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Lesur, G., Ecolano, B., Flock, M., et al. 2022, arXiv e-prints, [arXiv:2203.09821]

Lesur, G. R. J., Baghdadi, S., Wafflard-Fernandez, G., et al. 2023, A&A, 677, A9

Leung, P. K. C., & Ogilvie, G. I. 2019, MNRAS, 487, 5155

Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17

Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, MNRAS, 486, 453

Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120

Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805

Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402

Martel, É., & Lesur, G. 2022, A&A, 667, A17

Masset, F. S. 2008, in EAS Publications Series, 29, eds. M. J. Goupil, & J. P. Zahn, 165

Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494

McNally, C. P., Nelson, R. P., Paardekooper, S.-J., Benítez-Llambay, P., & Gressel, O. 2020, MNRAS, 493, 4382

Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2014, J. Comput. Phys., 257, 594

Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266

Nazari, P., Booth, R. A., Clarke, C. J., et al. 2019, MNRAS, 485, 5914

Nelson, R. P. 2005, A&A, 443, 1067

Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 584, A1

Ogihara, M., Kokubo, E., Suzuki, T. K., Morbidelli, A., & Crida, A. 2017, A&A, 608, A74

Paardekooper, S. J. 2014, MNRAS, 444, 2031

Papaloizou, J. C. B., Nelson, R. P., & Snellgrove, M. D. 2004, MNRAS, 350, 829

Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, ASP Conf. Ser., 534, 567

Pepliński, A., Artymowicz, P., & Mellema, G. 2008a, MNRAS, 387, 1063

Pepliński, A., Artymowicz, P., & Mellema, G. 2008b, MNRAS, 386, 179

Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15

Perez-Becker, D., & Chiang, E. 2011a, ApJ, 735, 8

Perez-Becker, D., & Chiang, E. 2011b, ApJ, 727, 2

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13

Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109

Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9

Rabago, I., & Zhu, Z. 2021, MNRAS, 502, 5325

Riols, A., & Lesur, G. 2019, A&A, 625, A108

Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95

Robert, C. M. T., Crida, A., Lega, E., Méheut, H., & Morbidelli, A. 2018, A&A, 617, A98

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2017, MNRAS, 468, 3850

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., Suzuki, T. K., & Shang, H. 2019, MNRAS, 484, 107

Svanberg, E., Cui, C., & Latter, H. N. 2022, MNRAS, 514, 4581

Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65

Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1

Tabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022, MNRAS, 512, 2290

Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378

Thi, W. F., Lesur, G., Woitke, P., et al. 2019, A&A, 632, A44

Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 85

Venuti, L., Bouvier, J., Flaccomio, E., et al. 2014, A&A, 570, A82

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

Wafflard-Fernandez, G., & Baruteau, C. 2020, MNRAS, 493, 5892

Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218

Williams, J. P., & McPartland, C. 2016, ApJ, 830, 32

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47

Zhu, Z., Stone, J. M., & Rafikov, R. R. 2013, ApJ, 768, 143

Ziampras, A., Kley, W., & Nelson, R. P. 2023, A&A, 670, A135
Appendix A: Temporal variability of the flow

<figure><img src="image_18.png" /><figcaption>**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-β₃, 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.</figcaption></figure>
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 $\langle\rho v_p\rangle_\phi$ in log scale, averaged azimuthally for the run $3M_j\beta_3$. As in Figure 6, the texture in LIC and the white arrows indicate the direction of the radial $(\langle\rho v_R\rangle_\phi)$ and vertical $(\langle\rho v_z\rangle_\phi)$ 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 $\beta_0 = 10^3$, and compare with the run $10Me\beta_3$ 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 $\approx 15\%$ and $\approx 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 $\approx 30\%$ (see the third panel in Figure B.1). Note also that the accretion via the radial mass flux is also increased by $\approx 20\%$ in the run $10Me\beta_3$ 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 $\dot{M}_{acc}$ between $R = 1.2$ and $R = 1.4$ seems to indicate that the structures observed are drifting outward during the episode considered here.

<figure><img src="image_19.png" /><figcaption>**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 $\Sigma$ in blue (first panel), amplitude of the poloidal magnetic field $B_p$ in green (second panel), wind torque at the disk surface $M_{surf}$ in yellow (third panel), radial mass flux $\dot{M}_{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.</figcaption></figure>

### 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 ($\gamma = 5/3$, $q = 1$, $h = 0.05$, with $-q$ being the power-law exponent of the radial temperature profile near the midplane $\tilde{T}(R)$), we need

$$ \Omega\mathcal{B} < 0.075 \qquad (B.1) $$

to get "vigorous" VSI. Hence, our models are likely to have severely quenched (or even suppressed) VSI signatures since we always use $\Omega\mathcal{B} = 0.1$.
<figure><img src="image_20.png" /><figcaption>**Fig. B.2.** Space-time diagram in log scale of the vertical component of the mass flux in the midplane and azimuthally averaged $\langle\rho v_z\rangle_{\phi,z=0}$, for the low-mass planet and low magnetization case $10\text{Me}-\beta_4$.</figcaption></figure>

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 $\phi$ 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 $\tilde{A}_m = 1$. 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 $\alpha_\nu$ (solid lines, as defined in Eq. 25) and $\alpha_T^R$ (dash-dotted lines). This last quantity is linked to the turbulent component of the Reynolds stress, defined as $\alpha_T^R = \alpha_\nu^R - \langle\rho\rangle_{\phi,t} \langle v_R \rangle_{\phi,t} \langle(v_\phi - v_K)\rangle_{\phi,t}$, 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 $\alpha_T^R$ 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

<figure><img src="image_21.png" /><figcaption>**Fig. B.3.** Radial profiles of the logarithm of the total $\alpha_\nu$ (solid lines) and $\alpha_T^R$ (dash-dotted lines) related to the turbulent Reynolds stress, with (green curves) and without (gray curves) a $10\tilde{M}_\oplus$ planet, for an initial $\beta_0$ of $10^3$ and averaged between 50 and 100 orbits at $R = 1$.</figcaption></figure>

at a low level (a few $\alpha_T^R \sim 10^{-3}$) and therefore not sufficient to explain alone angular momentum transport. This confirms the predominance of the $M_{\text{surf}}$ term in the angular momentum budget, as presented in Figure 1.

Besides, $\alpha_T^R$ is even not dominant in the total radial transport as stated earlier, (it is about $\approx 5$ times smaller than the total radial stress). In particular, for the run $10\text{Me}-\beta_4$, we checked that $\alpha_T^R \approx 7 \times 10^{-4}$ for $\beta_0 = 10^4$, 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 $\alpha_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 $M_{\text{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 $\alpha_T^R$, 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 $\Gamma$, and in particular its stochastic nature.

<figure><img src="image_22.png" /><figcaption>**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-$\beta_3$ and 10Me-$\beta_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 $\sigma_3 \approx 4 \times 10^{-2}$ and $\sigma_4 \approx 2 \times 10^{-2}$ for an initial $\beta_0$ of $10^3$ and $10^4$ respectively.</figcaption></figure>

To do so, we plot in Figure B.4 the frequency distribution of the planet torques $\Gamma$ for the runs 10Me-$\beta_3$ and 10Me-$\beta_4$ 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 $\Delta(\Gamma/q) = 2.8 \times 10^{-3}$, with a procedure similar to the one presented in (Nelson 2005). Fitting these distributions with Gaussian profiles, we obtain a standard deviation of $\sigma_3 \approx 4 \times 10^{-2}$ and $\sigma_4 \approx 2 \times 10^{-2}$ for an initial $\beta_0$ of $10^3$ and $10^4$ respectively. On the other hand, if we quantify the

Reynolds stresses, and more precisely $\alpha_T^R$ in both runs 10Me-$\beta_3$ and 10Me-$\beta_4$, we realize that decreasing the $\beta_0$ by a factor 10 increases $\alpha_T^R$ 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 $\alpha_T^R$.

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 $\alpha_T^R$ 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.

<figure><img src="image_23.png" /><figcaption>**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 ($10^3$, $10^4$, and a viscous case).</figcaption></figure>

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_{\text{hill}}$ to $R_p + 2R_{\text{hill}}$ and excised from $-2R_{\text{hill}}/R_p$ to $2R_{\text{hill}}/R_p$ in order to remove the contribution of the planet and its circumplanetary disk (CPD). $R_{\text{hill}}$ is the planet's Hill radius and is defined in Eq. (5). This procedure follows the evaluation of the gap's surface density $\Sigma_{\text{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 $\beta$ in the gap, considering its minimum in the radial interval $[R_p - R_{\text{hill}}, R_p + R_{\text{hill}}]$ and also excised from $-2R_{\text{hill}}/R_p$ to $2R_{\text{hill}}/R_p$. We define here several quantities that can be useful to understand and analyze qualitatively the process of gap opening:
- $\Sigma_{\text{gap}}^{\text{eq}}$ is the quasi-equilibrium value for the density in the gap, that is such that the mean value of $\Sigma_{\text{gap}}$ does not evolve much with time.

- $\delta\Sigma_{\text{gap}}^{\text{eq}}$ is the variability of the density in the gap once it has reached $\Sigma_{\text{gap}}^{\text{eq}}$.

- $v_{\text{gap}} = \left| \frac{d\Sigma_{\text{gap}}}{dt} \right|$ is the instantaneous gap-opening speed, from the initial density at the planet location to $\Sigma_{\text{gap}}^{\text{eq}}$.

- $t_{\text{gap}}^{\text{eq}}$ is the time needed by $\Sigma_{\text{gap}}$ to reach $\Sigma_{\text{gap}}^{\text{eq}}$.

Using these metrics and focusing on the $3M_1$ planet cases, we notice in the top panel of Figure C.1 that $\Sigma_{\text{gap}}^{\text{eq}}$ increases with the magnetization. More precisely, $\Sigma_{\text{gap}}$ reaches $\approx 3 \times 10^{-3}$ of its initial value $\Sigma_0$ (and still decreasing) for $3M_j-\alpha_0$ after 200 planet orbits, whereas it reaches $\Sigma_{\text{gap}}^{\text{eq}} \approx 10^{-2} \Sigma_0$ for $3M_j-\beta_4$ after $t_{\text{gap}}^{\text{eq}} \approx 90$ planet orbits and $\Sigma_{\text{gap}}^{\text{eq}} \approx 2 \times 10^{-2} \Sigma_0$ for $3M_j-\beta_3$ after $t_{\text{gap}}^{\text{eq}} \approx 50$ planet orbits. This is consistent with our interpretation based on the interplay between the accretion rate $M_{\text{acc}}$, $\beta$ and $\Sigma$ and presented in the second paragraph of Section 3.1.2. $\delta\Sigma_{\text{gap}}^{\text{eq}}$ 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 $v_{\text{gap}}$ and $t_{\text{gap}}^{\text{eq}}$. On the one hand, the gap-opening speed is larger when $\beta_0$ decreases, at least early in the simulations and before $\Sigma_{\text{gap}} \approx \Sigma_{\text{gap}}^{\text{eq}}$ (see, e.g., between 20 and 90 orbits). It would suggest that a torque linked to the accretion is actually helping to form the gap. On the other hand, $t_{\text{gap}}^{\text{eq}}$ decreases when $\beta_0$ 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 $3M_j-\alpha_0$, $3M_j-\beta_4$ and $3M_j-\beta_3$ 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 $\beta_0$ 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 $\Sigma_{\text{gap}} \approx \Sigma_{\text{gap}}^{\text{eq}}$ (see the $\beta_0 = 10^4$ and $10^3$ cases in the top panel of Figure C.1). At that moment, when $\beta_{\text{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 $\beta_{\text{gap}}$ reaches a value < 10 near 50 and 90 orbits respectively when $\beta_0 = 10^3$ and $10^4$.