Three-dimensional Global Simulations of Type-II Planet-disk Interaction with a Magnetized Disk Wind: I. Magnetic Flux Concentration and Gap Properties

We use Athena++, a grid-based, higher-order Godunov MHD code (Stone et al., 2020) to solve the following fluid equations in conservative form

where $\rho,\Phi$ are the gas density and gravitational potential, $v$ and $B$ are vectors of gas velocity and magnetic field, $\mathsf{M}=\rho\mbox{\boldmath$v$}\mbox{\boldmath$v$}-\mbox{\boldmath$BB$}+P^{*}\mathsf{I}$ is the stress tensor where $\mathsf{I}$ is the rank two unit tensor and $P^{*}=P+B^{2}/2$ is the total pressure with $P$ being the thermal pressure, and $\mbox{\boldmath$E$}=\eta_{\mathrm{O}}\mbox{\boldmath$J$}+\eta_{\mathrm{A}}\mbox{\boldmath$J$}_{\perp}$ is the non-ideal electric field in the rest fluid frame, where $\eta_{\mathrm{O}}$ and $\eta_{\mathrm{A}}$ are the Ohmic and ambipolar diffusivities, respectively, and $\mbox{\boldmath$J$}=\nabla\times\mbox{\boldmath$B$}$ is the current density vector, $\mbox{\boldmath$J$}_{\perp}$ is the current density vector perpendicular to $B$, respectively. In the code, the units for $B$ is chosen such that the magnetic pressure is $B^{2}/2$ and hence magnetic permeability is 1.

We adopt a locally isothermal equation of state with $P=\rho c_{\mathrm{s}}^{2}$, where $c_{\mathrm{s}}$ is the isothermal sound speed to be specified in Section 2.3. In doing so, we are free from the buoyancy resonance (e.g., Bae et al., 2021) which allows us to focus on MHD effects, and this condition is likely satisfied in the outer PPDs (e.g. Lin & Youdin, 2015; Pfeil & Klahr, 2019). On the other hand, the system is subject to the vertical shear instability (e.g., Nelson et al., 2013), which can coexist with the MRI as a source of disk turbulence (Cui & Bai, 2022). In our simulation setup ($Am=1$ and 3, see Section 2.3), the MRI generally dominates.

We primarily focus on the outer regions of PPDs, where ambipolar diffusion (AD) is the dominant non-ideal MHD effect (Wardle, 2007; Bai, 2011a). The ambipolar diffusivity $\eta_{A}$ is characterized by the dimensionless Elsasser number defined as

where $v_{\mathrm{A}}=\sqrt{B^{2}/\rho}$ is the Alfvén speed, $\Omega_{\mathrm{K}}=\sqrt{GM_{*}R^{-3}}$ is the Keplerian angular velocity, $M_{*}$ is the central star mass, and $R=r\sin{\theta}$ is the cylindrical radius. Note that in the absence of abundant small grains, $\eta_{A}\propto B^{2}$, so that $Am$ is independent of field strength (Bai, 2011b). We further add some Ohmic resistivity $\eta_{O}$ near the inner region to stabilize the inner boundary (see Section 2.3).

Our simulations are conducted in a frame corotating with the planet, and describe the implementation of rotating frame appropriate source terms for angular momentum conservation in Appendix A. On the other hand, in the presentation of this paper, we consider $v_{\phi}$ as in lab-frame for convenience unless otherwise noted.

For 3D and 2D simulations, we solve these equations in the spherical polar coordinates ($r$, $\theta$, $\phi$) and cylindrical coordinate ($R$, $\phi$), respectively, centered at the central star. Note that we also use the 3D cylindrical coordinate ($R$, $\phi$, $z$) in analyzing the results of 3D simulations, where $z=r\cos{\theta}$.

The gravitational potential in the rotating frame is given as (e.g., Fung & Chiang, 2016)

where $q=M_{\mathrm{P}}/M_{*}$ is the planetary mass ratio to the central star where $M_{\mathrm{P}}$ is the planetary mass, $R=r\sin{\theta}$ is the cylindrical radius, $R_{\mathrm{P}}$ is the planet orbital radius, and $r_{\mathrm{s}}$ is the softening radius. The third term in Eq. (5) is the indirect potential due to fixing the star at the center of the simulation domain. We set $r_{\mathrm{s}}$ to be twice the cell size in $\phi$-direction at the planet. As a first study, we fix the planetary orbit with circular motion, while we also note that the gap structure can be affected by planet migration (e.g., Meru et al., 2019; Kanagawa et al., 2020) and orbital eccentricity (e.g., Sánchez-Salcedo et al., 2023).

The planet gravity is introduced after a time of $t_{\mathrm{P0}}$ (see the Table 1). To minimize the influence of sudden mass addition, we make the planet mass linearly increase to the final mass at the rate of

where

is the thermal mass, which is used as a criteria whether planet can open a density gap at its orbit, $H_{\mathrm{P}}$ is the disk scale height at the planetary orbit (see next subsection), $\tau_{\mathrm{P}}=2\pi/\Omega_{\mathrm{P}}$ is the planetary orbital period, and $\Omega_{\mathrm{P}}$ is the planetary angular velocity. This growth rate is comparable or longer than that in Fung et al. (2014). For reference, the planetary Hill radius is defined as

and hence for $M_{\mathrm{P}}=3M_{\mathrm{th}}$, $r_{H}=H_{P}$.

In our setup, the disk temperature is characterized by $c_{s}$, which we specify as a function of $R$ and is represented by the local disk scale height $H$ and disk aspect ratio $h$

In this work, we fix the aspect ratio of the bulk disk $h_{\mathrm{disk}}\equiv H_{\rm disk}/R$ to be $0.1$ (constant, and without causing ambiguity, we drop the subscript “disk"). Then, the vertical density profile in hydrostatic equilibrium is given by (e.g., Gressel et al., 2015)

where we use $r_{0}=1$, $\rho_{0}=1$, and $q_{\rho}=2.25$ for the fiducial runs (see Table 1). We also set a density floor value of $10^{-8}(r/r_{0})^{-q_{\rho}}$ in the code units. The initial surface density $\Sigma_{0}\approx\sqrt{2\pi}H_{\rm disk}\rho_{\mathrm{ini}}\propto r^{-1.25}$. This is also adopted in 2D simulations. To achieve force balance, the $\phi$-direction velocity is initially set to

in the corotating frame with the planet. The velocities in other directions $v_{r}$ and $v_{\theta}$ are initialized with random values up to $2.5\%$ of the local sound speed.

At high altitudes, the disk atmosphere is tenuous but is expected to be hot due to the UV and X-ray heating (e.g., Glassgold et al., 2004). To mimic such a hot atmosphere, we allow the disk aspect ratio $h$ to transition from $h_{\rm disk}$ to some higher value $h_{\rm atm}$ in the following functional form:

and we choose $h_{\mathrm{atm}}=0.4$ and $z_{\mathrm{t}}=4H$ is the transition height. Similar approaches have been adopted in our earlier works (e.g., Bai & Stone, 2017; Cui & Bai, 2021). We note that although hydrostatic equilibrium is not initially maintained at the disk upper layer, the subsequent wind launching process will quickly override any imbalance, and the outcome is largely independent of the detailed treatment here.

The poloidal magnetic field is initialized as

where the vector potential is given in Zanni et al. (2007); Bai & Stone (2017) as

where $B_{z0}$ is the initial poloidal field at the midplane of the $r$-inner boundary ($r=r_{0}=1$) and $m=1$ is the geometry parameter. We set $B_{z0}=\sqrt{2\beta_{0}\rho_{0}c_{\mathrm{s}}^{2}}$ so that the entire disk is magnetized with the same plasma beta ($\beta_{0}$, ratio of initial gas pressure to magnetic pressure), and we choose $\beta_{0}=10^{4}$. Hence, the initial magnetic flux at the midplane is $B_{z0\mathrm{,mid}}=\sqrt{2\beta_{0}\rho_{\mathrm{ini}}hv_{\mathrm{K}}}$, where $v_{\mathrm{K}}=R\Omega_{\mathrm{K}}$ is the Keplerian velocity.

Following earlier works (Bai, 2011a, b; Simon et al., 2013; Bai & Stone, 2017; Riols & Lesur, 2019; Cui & Bai, 2021), we set the ambipolar Elsasser number to be constant of order of unity over the bulk disk column, while the value of $Am$ transitions to much higher values in the disk atmosphere to mimic higher ionization due to the penetration of far-UV radiation (Perez-Becker & Chiang, 2011). Since the physics are similar to the heating of atmosphere, we use the same functional form as Eq. (12) for the transition in $Am$ value. We set the $Am_{\mathrm{disk}}=1$ or 3 and $Am_{\mathrm{atm}}=100$ in the bulk disk and the atmosphere, respectively. The midplane $Am$ of 1 or 3 leads to a moderately strong MRI turbulence that overwhelms the potential influence of the VSI (Cui & Bai, 2022). Although gap formation would modify the ionization profile and hence the $Am$ values in the vicinity of the planet orbit, we note that the previous studies on transitional disks infer that $Am$ is still of order unity even inside the deep density cavity (Wang & Goodman, 2017; Martel & Lesur, 2022). As the first study of planet-disk interaction with both MRI turbulence and MHD disk wind, this simplified treatment allows us to focus on purely MHD effects, and we defer more self-consistent treatment of ionization chemistry to a future study. It is beyond the scope of this work to explore broader ranges of $Am$ values and/or more self-consistent chemistry (e.g. Bai, 2017). In this study, we treat $Am_{\mathrm{disk}}=3$ as the fiducial, which is on the high side but is also in line with the realization that cosmic-ray traveling along magnetic fields can enhance the ionization rate in the disk outer region (Fujii & Kimura, 2022). Additionally, following § 2.3 in Cui & Bai (2021), we apply Ohmic resistivity in the inner region ($r<2$) to help stabilize the inner boundary.

For the main 3D simulations, the simulation domain spans from $r=1$ to $100$ in code units with logarithmic grid spacing, and the $\theta$ and $\phi$ domains cover full $\pi$ and $2\pi$, respectively, with the polar boundary condition (Zhu & Stone, 2018). The planet orbit is fixed at $R_{p}=6$ in code units¹¹1It is known that global MHD simulations of disk winds tend to be affected by inner boundary conditions (e.g. Cui & Bai, 2021). Despite of rendering the simulations more expensive (per planet orbit), we choose a relatively large $R_{p}$ in order to leave sufficient dynamical range to minimize this influence., and in the corotating frame, its coordinate is $(r,\theta,\phi)=(6,\pi/2,0)$. We note that the large radial extent relative to the planetary orbit and the $\theta$ extent to both poles are essential to properly accommodate the MHD disk wind. Similar to Cui & Bai (2021), the $\theta$ grid size increases from the midplane to the pole by a constant factor, and the pole to midplane grid size ratio is four. We use three levels of static mesh refinement on top of the root grid, which has 112, 48, and 96 grid cells in $r$-, $\theta$-, and $\phi$-directions, respectively. The finest grid ranges from $2.5<r<12$, $1.5<\theta<1.64$, and full 2$\pi$ in phi-direction, respectively. In the $\theta$-direction, the finest grid cell resolves the disk scale height by $\sim 26$ cells. At the planet orbit ($R=6$), the cell size is $(\Delta r,r\Delta\theta,r\Delta\phi)=(0.030,0.023,0.049)$ in code units.

At both radial boundaries, $\rho$ and $T$ are copied from the nearest grid cell to the ghost zones with extrapolating following the initial radial gradient. At the outer boundary, the $v_{r}$ and $v_{\theta}$ are copied from the nearest cell to the ghost cell, while for $v_{r}<0$, $v_{r}$ is forced to zero to prevent inflow. We also copy $v_{\phi}$ with extrapolation according to $r^{-1/2}$. At the inner boundary, both $v_{r}$ and $v_{\theta}$ are fixed to zero, and $v_{\phi}$ is fixed as initial, which was found useful to stabilize the inner boundary region (see also Bai & Stone, 2017). Magnetic field in the ghost zone is copied from the nearest grid cell assuming $B_{r}\propto r^{-2}$ and $B_{\phi}\propto r^{-1}$, with $B_{\theta}$ unchanged. Besides, the electric field in $\phi$-direction is forced to zero at the radial inner boundary to anchor the initial magnetic flux in the inner boundary.

In the 2D simulations, with cylindrical coordinates $(R,\phi)$, no magnetic field is included, and we include the standard $\alpha$ viscosity $\nu=\alpha c_{\mathrm{s}}H$. Other settings are same as the 3D simulation described above.

The main parameters for the simulations are listed in Table 1, where $t_{\mathrm{P0}}$, $t_{\mathrm{Pf}}$, and $t_{\mathrm{e}}$ are the times when the planet is introduced, when the planet reaches the final mass, and when the simulation stops, respectively. In this paper, we consider three different planet masses, with $M_{\mathrm{P}}=1M_{\mathrm{th}}$, $3M_{\mathrm{th}}$ and $5M_{\mathrm{th}}$, respectively, and take the case with $M_{\mathrm{P}}=3M_{\mathrm{th}}$ and $Am_{0}=3$ as fiducial. Our presentation of simulation results mostly focus on the $Am=3$ simulations, and we mainly summarize the results from the $Am=1$ simulations in Appendix D which are consistent with our fiducial runs. Given that we anticipate spontaneous formation of ring-like substructures even without planets (Cui & Bai, 2021), we choose planet masses that are expected to be largely enough to compete with the planet-free gap-formation mechanisms.

Before we present simulation results, we introduce and define certain concepts and quantities helpful for diagnostics. The simulation domain is divided into the disk and wind regions which coincide with the temperature and $Am$ transition (Bai & Stone, 2013; Gressel et al., 2015). Although the characteristic height for the transition is set to $z=4H$, the temperature and $Am$ transitions start at a bit lower height as following Eq. 12. Judging from the vertical profile of the velocities (see also figure 10), we define the disk zone as $|z|<3.5H$. Hence, the disk surface density is defined as

Note that the disk boundary is constant at $\theta=\theta_{\mathrm{t^{\prime}}}\equiv\arctan(3.5h)$ because the disk aspect ratio is assumed to be constant.

To illustrate the evolution of poloidal magnetic flux, we use linearly equally spaced contour for the poloidal magnetic flux function defined as

We also define the mean vertical magnetic field in the midplane region

which is more convenient when considering, e.g., correlations with disk surface density profiles. Since the poloidal magnetic field is nearly vertical at the midplane but bends radially outward at high altitudes (see figure 11), we average only within $1H$ rather than within the whole disk ($3.5H$).

The disk structure including the planet-induced gap is governed by angular momentum transport, which is more naturally analyzed in spherical-polar coordinates, and the equation for angular momentum conservation can be expressed as

where $j=r\sin{\theta}v_{\phi}$ is the specific angular momentum in the $\phi$-direction (in the lab frame), and $T^{\mathrm{Max}}_{ij}=-B_{i}B_{j}$ is the $ij$th component of the Maxwell stress tensor.

We can define the standard disk integral as

where $\theta_{\mathrm{d,u}}=\pi/2\pm\theta_{\mathrm{t^{\prime}}}$ (note $\theta_{\mathrm{t^{\prime}}}=0.35$ in our case). In particular, we can define the (spherical-shell-integrated) surface density and accretion rate as

We can further define the disk average of a quantity $A$, denoted by angle bracket as

In our analysis, unless otherwise noted, we further average the data over $130<t/\tau_{\mathrm{P}}<140$ in time, and $\pm 0.25~{}H$ in $r$ to enhance statistics, while excluding the region where $|\mbox{\boldmath$r$}-\mbox{\boldmath$R_{\mathrm{P}}$}|<r_{\mathrm{H}}$. In particular, we define

and local deviations from these averaged quantities can be denoted by $\delta(\rho v_{r})$ and $\delta j$.

By integrating over $\phi$- and $\theta$-directions, and by utilizing the continuity equation, Eq. (18) can be cast into

where

In the above, $\gamma_{\mathrm{P}}$ and $\gamma_{\mathrm{wind}}$ are the torque densities due to the planet gravity and MHD wind, $J_{\mathrm{wave}}$ and $J_{\mathrm{Max}}$ are the angular momentum flux due to the Reynolds and Maxwell stresses, respectively.

In addition, we define the cumulative torque $\Gamma$, which has a same dimension as the angular momentum flux, as

Note that $-\{\Gamma_{\mathrm{P}}(r_{\mathrm{out}})-\Gamma_{\mathrm{P}}(r_{\mathrm{in}})\}$ is the torque acting on the planet as the backreaction ($r_{\rm in}$, $r_{\rm out}$ are the radius of the inner and outer radial boundary), so-called the planet migration torque, where $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ are the radius of the inner and outer boundary of the disk.

We have already defined $T^{\rm Max}$, and the Reynolds stress is defined as

Note that the Reynolds stress has a contribution from the MHD-driven flows (e.g., the MRI turbulence), as well as the planetary density wave/shock, although it is not straightforward to separate them. The remaining term $\gamma_{A}$ is the torque density associated with bulk accretion, which should be understood as $\dot{M}_{\mathrm{acc}}$ being the response to the action of the other four terms.

The stresses are related to the Shakura-Sunyaev $\alpha$ values (Shakura & Sunyaev, 1973), which are given, when averaged over the disk, as (e.g., Balbus & Hawley, 1998)

and their sum gives the total $\alpha$ value.

There are a few points worth further clarification. First, while $\gamma_{\mathrm{wave}}$ comprises of the Reynolds stress, it is largely dominated by the planetary density wave (here the spiral shock) at least in the vicinity of the planet. Following Kanagawa et al. (2017), we refer to this term as the “wave" term representing the angular momentum flux carried by the density wave/shock. Second, while the MRI turbulence produces both the Reynolds and Maxwell stresses, the total stress and hence $\alpha$ is dominated by the Maxwell stress by a factor of a few (e.g., Bai & Stone, 2011). Thus, we consider $J_{\mathrm{Max}}$ as the proxy for total viscous flux of angular momentum (see also § 3.4). Third, we analyze angular momentum transport in spherical polar coordinates, which is more rigorous than in cylindrical coordinates, despite of being less optimal for wind diagnostics since the direction of the wind deviates substantially from the $\theta$-direction. Here, the local wind mass-loss rate can be defined as

Finally, we note that in deriving Equation (23), we have already accounted for the angular momentum loss associated with bulk mass loss $\dot{M}_{\mathrm{wind}}\langle j\rangle$, and hence in $\gamma_{\mathrm{wind}}$, the first term only accounts for the deviation $\delta j$ instead of the total $j$. Also note that we have omitted a term of $\iint ds\partial(\rho\delta j)/(\partial t)$ from Eq. (23), which is expected to be small at least in the quasi-steady state that we have achieved.

Similarly to the above diagnostics in 3D spherical coordinate, the torque in the 2D viscous simulations is categorized into $\gamma_{\mathrm{P}}$, $\gamma_{\mathrm{A}}$, $J_{\mathrm{A}}$, $J_{\mathrm{wave}}$, and $J_{\mathrm{vis}}$. Here, there is no $\gamma_{\mathrm{wind}}$, and the $\gamma_{\mathrm{Max}}$ is replaced by the viscous torque density defined as (e.g., Kanagawa et al., 2015)

The direct contribution from the planet arises from both the $\gamma_{\mathrm{P}}$ term of planet gravity, together with the $J_{\mathrm{wave}}$ term representing wave propagation. The net torque deposition rate can be approximately expressed as (e.g., Kanagawa et al., 2017)

which is associated with the dissipation of the density waves and directly leads to planet-induced substructure formation.

Without external forces, fluid parcels in the planet vicinity experience horseshoe orbits over the radial extent of $\sim r_{\mathrm{H}}$ (e.g., Masset et al., 2006; Paardekooper & Papaloizou, 2009), and this region is thus called the corotation region or horseshoe region. For convenience, we just refer to $|R-R_{\mathrm{P}}|<r_{\mathrm{H}}$ as the corotation region in this work.

Figure 1 shows the major MHD quantities in our fiducial case (Mt3Am3). The panels from (a) to (d) are the snapshot at $t=3$, 30, 100, and $140~{}\tau_{\mathrm{P}}$, respectively. Since we introduce the planet at $t=3.4~{}\tau_{\mathrm{P}}$, there is no planet and no significant structure in panel (a). Note that this snapshot is common with the three runs with $Am=3$, and the density feature in this snapshot is part of the initial relaxation process. The disk is subject to both the MHD-disk wind and the MRI turbulence. While the MRI is not yet well developed at the planetary orbit, the disk wind is established and starts to drive the disk accretion.

The planet reaches to the final mass of $3~{}M_{\mathrm{th}}$ at $t=28.4~{}\tau_{\mathrm{P}}$. In the panel (b) of $t=30~{}\tau_{\mathrm{P}}$, the planet excites the density wave (spiral shock) and starts to open a density gap around its orbit. Interestingly, we see that poloidal magnetic flux gets concentrated in the gap region, which we identify as a major characteristic of planet-disk interaction in the presence of MHD winds. We note that this phenomenon was also identified in earlier local shearing-box simulations (Zhu et al., 2013; Carballido et al., 2016, 2017), nevertheless, these simulations were unstratified thus do not launch disk winds, and they were conducted under ideal MHD conditions. Also note that, at this time, the MRI turbulence is fully developed at the planetary orbit.

In panel (c) at $t=100~{}\tau_{\mathrm{P}}$, the planetary gap gets wider and deeper. In addition, a planet-free gap appears inside the planetary orbit. This phenomenon has been found in previous simulations of the MRI turbulence and/or disk winds without a planet (Bai & Stone, 2014; Bai, 2015; Béthune et al., 2017; Suriano et al., 2018, 2019; Riols & Lesur, 2019; Riols et al., 2020; Cui & Bai, 2021), where magnetic flux spontaneously concentrate into quasi-axisymmetric flux sheets, leading to the formation of ring-like substructures. This process is expected to be stochastic, and applies for a wide range of $Am$ values (Cui & Bai, 2021). Both the planetary and the planet-free gaps get wider over time, and by the time of $t=140~{}\tau_{\mathrm{P}}$, they start to overlap with each other, namely the density between the two gaps is getting depleted over time relative to the initial (white), and the two gaps are on a merging trajectory (see Figure 3 in the next subsection). Due to substantial computational cost, we terminate our simulation at this time. Also note that in practice, it will likely take $10^{3}$ planetary orbits for the gap profile to reach a steady state as found in pure hydrodynamic simulations (e.g., Fung & Chiang, 2016).

To further illustrate the gap opening process and magnetic flux concentration, we show in figure 2 the snapshots of $\Sigma$ (upper) and $B_{z\mathrm{,mid}}$ (bottom) at $t=10$, $30$, $100$ and $140~{}\tau_{\mathrm{P}}$. The planetary gap ($R=R_{\mathrm{P}}(=6)$) is clearly associated with the spiral density waves/shocks. For the vertical magnetic flux threading the disk, we see that in the inner planet-free gap, magnetic flux largely concentrates in the gap region, consistent with previous findings (Cui & Bai, 2021). On the other hand, magnetic flux distribution around the planet-induced gap is clearly asymmetric, akin to the spiral density waves/shocks. Moreover, a patch of negative magnetic flux is present at the trailing side of the spiral shock near the planet. We will discuss the magnetic field concentration in § 3.3. Its effects on the disk turbulence and gap opening will be discussed in § 3.4 and § 5, respectively.

To highlight the density waves/shocks, we further show in the middle panels of Figure 2 (and also in Figure 6) the azimuthal surface density variations $\Sigma/\left\langle\Sigma\right\rangle_{\phi}$. Besides the primary planet density waves/shocks already visible in the top panels, we can identify a secondary spiral interior of the planet orbit, especially at early time of $t=30\tau_{\mathrm{P}}$. Secondary spirals have been found in previous hydrodynamic simulations of planet-disk interaction in inviscid/low-viscosity disks (Dong et al., 2015; Zhu et al., 2015; Bae et al., 2017), and is a result of constructive interference of a set of spiral modes Bae & Zhu (2018). They can also open additional gaps in low viscosity ($\alpha\lesssim 10^{-3}$) disks (Bae et al., 2017). While they might contribute to the initial formation of the planet-free gap in our simulations, we note that the effective viscosity in our simulations is higher, and that the amplitudes the secondary spirals in our simulations are relatively low ($\lesssim 20\%$). We also observe that at later time, the secondary spiral appears to break apart towards disk interior upon encountering the planet-free gap, and quantifying their role in our simulation is thus not straightforward. We thus refrain from further discussion of secondary spirals, and regardless of the onset of the planet-free gap, they are clearly sustained by magnetic flux concentration in our simulations.

In this subsection, we study the structure of the planet-induced density gap in detail. Firstly, we show in the top panels of figure 3 the surface density profiles of the three runs Mt1Am3, Mt3Am3, and Mt5Am3 at different snapshots, normalized by the initial surface density profile. As discussed before, two gaps are developed. The one peaking at $R\sim R_{\mathrm{P}}$ is induced by planet, and the other one at $R\sim 3$ is induced by magnetic flux concentration. As the two gaps are relatively close, they get overlapped with time, and in the case of run Mt5Am3, they get fully merged at $t=140~{}\tau_{\mathrm{P}}$.

To quantify the gap properties, we define the gap depth as the ratio of the local minimum to the initial surface density profile, and gap width being the full width at half of the initial surface density, normalized by the planet orbital radius. Figure 4 shows the temporal evolution of the depth (black solid) and width (red solid) of the planet-induced gap in the three runs. The gap width gradually get wider with time but shows a jump at $t\sim 110$ and $\sim 90\tau_{\mathrm{P}}$ for run Mt3Am3 and Mt5Am3, respectively. It corresponds to when the $\Sigma$ between the planet-induced and planet-free gaps gets smaller than $\Sigma_{0}/2$. After the jump, the width in Mt3Am3 almost stays constant while the radial gap structure of $\Sigma$ is still evolving. On the contrary, the width in Mt5Am3 gets slightly narrower due to the smoothing after the two gaps get fully merged (see the purple and red lines in figure 3). In run Mt1Am3, the jump has not yet occurred, but $\Sigma$ between two gaps is getting lower towards the end of the simulation.

The width and depth of planet-induced gaps have been widely studied in previous works (e.g., Crida et al., 2006; Duffell & MacFadyen, 2013; Fung et al., 2014; Duffell, 2015; Kanagawa et al., 2015, 2016; Ginzburg & Sari, 2018; Duffell, 2020), almost exclusively under the framework of viscous disks. To compare our results from these studies, in figure 4, we overplot the depth (grey) and width (pink) in the viscous disks. Hence, for comparison, we also compare our results with the semi-analytical fit for the gap depth and width in $\alpha$-viscous disks by Kanagawa et al. (2016)

where

$\Sigma_{\mathrm{th}}=\Sigma_{0}/2$ is the threshold surface density, and $h_{\mathrm{P}}$ is the disk aspect ratio at the planetary orbit. In our runs, the turbulent viscosity is effectively evaluated as $\alpha=6\times 10^{-3}$ far from the plane ($R\sim 20$; see § 3.4). By substituting this $\alpha$ and $h_{\mathrm{P}}=0.1$ into Eqs. (37)-(40), the gap depth and width in the viscous accretion disk are predicted to be $(\Sigma_{\mathrm{min}}/\Sigma_{0},\Delta_{\mathrm{gap}}(\Sigma_{0}/2)/R_{\mathrm{P}})=(6.0\times 10^{-1},0.26)$, $(1.4\times 10^{-1},0.45)$, and $(5.6\times 10^{-2},0.59)$ for $M_{\mathrm{P}}=1$, 3, and $5M_{\mathrm{th}}$, respectively. They are shown as the dotted lines in figure 4. At the end of our simulations, the gap profile is wider and deeper than the prediction of Eqs. (37) and (38) by factors of (4, 1.7), (6, 1.5), and (10, 1.3) for the three runs (here the width are measured just before the jump for the latter two runs Mt3Am3 and Mt5Am3 for fair comparison). Note that our viscous simulations of Mt1$\alpha$6, Mt3$\alpha$6, and Mt5$\alpha$6 have not yet reach to the full steady state by the end of simulation and hence the depth is shallower than the prediction.

As we shall see in Section 3.4, the $\alpha$ values in our simulations are not spatially constant, and become much higher in the gap region due to magnetic flux concentration. If we adopt $\alpha$ in the gap region (averaged within $|R-R_{\mathrm{P}}|<r_{\mathrm{H}}$) for runs Mt1Am3, Mt3Am3, and Mt5AM, their values measured at $t=140~{}\tau_{\mathrm{P}}$ are found to be $3\times 10^{-2}$, $2\times 10^{-1}$, and $9\times 10^{-1}$, respectively. Applying the same semi-analytical formula, we find the predicted gaps would be much shallower and narrower (thus deviate much more substantially from simulation results), as shown in dashed lines in figure 4.

In addition, we compare the gap depth and width with those in the inviscid disk (pale grey and pale pink in figure 4). Even comparing with the inviscid disk, the gap depth is deeper by a factor of 3, 5, and 10, respectively. On the contrary, interestingly, the gap width evolves similarly in the wind and inviscid disks, at least before the jump due to gap merging occurs.

Overall, we conclude that for planet mass $\gtrsim$ thermal mass, the planet-induced gap in the MHD-wind-driven disk is much deeper and modestly wider than that in a viscous disk. Rather, in terms of the gap width, the windy disk is more similar to the inviscid disk, while the gap depth is still much deeper in the windy disk. This result is due to a combination of the nature of wind-driven accretion and magnetic flux concentration, as will be detailed in Section 5.

We have already mentioned that our simulations form a planet-free gap due to magnetic flux concentration at $R\sim 3$, together with the planet-induced gap, which also shows magnetic flux concentration. To better comprehend how these structures come into being, we show in figure 5 the temporal evolution of the disk surface density $\Sigma$ and mean vertical magnetic field at the midplane $B_{z\mathrm{,mid}}$, for all three runs with $Am=3$.

We start by discussing the gap around the planet orbit. Magnetic flux concentration in the planet gap is one of the most significant new findings in our study. It is the main reason that leads to the formation of wider and deeper gaps, as we analyze in more detail in Section 5. How magnetic flux spontaneously concentrates to the gap at first place, on the other hand, is a much more complex question, which generally involves the complex interplay between the gas dynamics (leading to flux advection), and non-ideal MHD and (turbulent) dissipative processes (leading to flux diffusion), yet both of which are dependent upon the magnetic flux configuration itself that drives wind and turbulence. Despite numerous attempts (e.g., Lubow et al., 1994; Guilet & Ogilvie, 2013; Okuzumi et al., 2014; Riols et al., 2020), most of these works rely on phenomenological assumptions of advection/diffusion that lack dynamical feedback, while simulation results also appear to depend on implementation details (e.g., Bai & Stone, 2017; Lesur, 2021). Clearly, flux concentration in planet-induced gap reflects the back-reaction of gravitational influence of the planet on magnetic flux transport, which further complicates the picture. Here, we focus our discussion mainly at phenomenological level. Further analysis will be presented in Appendix B, where we show that flux concentration in planet-induced gaps firstly occurs at high altitudes and is likely associated with vertical flows above the planet Hill sphere.

There are several noticeable features in magnetic flux concentration around the planet gap. First, this process requires the presence of net poloidal flux. By contrast, earlier simulations with toroidal magnetic flux (but zero net poloidal flux) showed that toroidal magnetic flux gets expelled from the gap region (Nelson & Papaloizou, 2003; Winters et al., 2003; Papaloizou et al., 2004). Second, the radial profile of $B_{z\mathrm{,mid}}$ is not necessarily centrally peaked, but can exhibit double peaks at higher planet mass, as seen in figures 3 and  5. We note that similar findings were reported in local simulations by Carballido et al. (2016) despite of being in the ideal MHD regime without stratification. This double-peak structure is related to the asymmetric distribution of $B_{z}$ in the $R-\phi$ plane as seen in figure 2 for run Mt3Am3, as well as shown in figure 6 for the other two runs. We also find that the double-peak profile merge into a single peak at higher altitudes, as can also be seen in figures 1 and 11. Third, the level of magnetic flux concentration is relatively strong, with $B_{z\mathrm{,mid}}$ reaching a factor of $\sim 2\text{--}4$ the background level within the $R_{H}$ from the planetary orbit. This requires gathering magnetic flux from a wide range of neighboring regions, especially as the corotation region (and the gap) becomes wider for more massive planets.

The planet-free magnetic flux concentration has been reported in various MHD simulations of PPDs without planet (e.g., Bai & Stone, 2014; Bai, 2015; Béthune et al., 2017; Suriano et al., 2018, 2019; Riols & Lesur, 2019; Riols et al., 2020; Cui & Bai, 2021). In the presence of the MRI turbulence, the location and evolution of magnetic flux concentration is stochastic and generally unpredictable. We clarify that the location of the planet-free gap in our three fiducial runs are largely the same because they are restarted from the same initial condition. The specific location of $R\approx 3$ is likely the outcome of the stochasticity from our initial condition.²²2We cannot rule out that the planet influences the formation of the planet-free gap due to secondary spirals (Bae et al., 2017) and subsequent magnetic flux concentration, especially as the planet is inserted relatively early. On the other hand, the $R=3$ gap is not present in our simulation with $Am=1$ presented in Appendix D, where the planet is inserted at much later time (and note that the “gap” further inward is too close to the inner boundary to be trustworthy). Also note that there are additional flux concentration beyond the radial range of figure 5 at larger radii. The planet-free $B_{z\mathrm{,mid}}$ concentration moves slightly outward in the Mt1Am3 and Mt5Am3 cases, but stays at the initial location in the Mt3Am3 case. We anticipate such migration to be less predictable. On the other hand, the amount of magnetic flux trapped at $R\sim 3$ appears to anti-correlate with the planet mass. Being in the vicinity of the planet orbit, reflecting that more flux gets “attracted" to the more prominent planet gap for more massive planets.

We further look at the spatial distribution of the magnetic flux in figures 2 and 6, in the $R-\phi$ plane. In the planet-free gap, the distribution of $B_{z\mathrm{,mid}}$ is largely axisymmetric. However, in the planet gap, the distribution is clearly asymmetric, and magnetic flux concentration traces the spiral density shocks in the leading and trailing side in the vicinity of the planet. The positive concentration of $B_{z\mathrm{,mid}}$ has two arms starting from the planet; the inner (outer) arm azimuthally extends within (to the outer side) the corotation region ($|R-R_{\mathrm{P}}|<r_{\mathrm{H}}$) towards the leading (tailing) side. On the other hand, the tailing side of the gap has a localized region with negative $B_{z\mathrm{,mid}}$. As a result, the two separated $B_{z\mathrm{,mid}}$ peaks in the radial profile is a continuous distribution in the $r$-$\phi$ plane. We also observe that the double peaks in run Mt5Am3 turn to a single peak at $t\gtrsim 120~{}\tau_{\mathrm{P}}$ as the gap merges with the planet-free gap and gather more flux from there.

It is well known that the strength of the MRI turbulence, together with the Maxwell and Reynolds stresses, directly correlates with the amount of net vertical magnetic flux threading the disk (e.g., Hawley et al., 1995; Bai & Stone, 2011). The spontaneous flux concentration thus implies that the level of turbulence in the disk is inhomogeneous.

Figure 7 shows the radial distribution of the stresses for the three runs with $Am=3$. As the $\alpha$ values correspond to stresses normalized by pressure, radial variations in surface density (and hence pressure) yields additional variations in $\alpha$ values. To facilitates the analysis, we show the $\alpha$ values normalized by both the initial pressure profile (top) and the concurrent pressure profile (bottom). Without the planet, and/or sufficiently far away from the planet, we see that $\alpha^{\mathrm{Max}}$ is on average greater than $\alpha^{\mathrm{Rey}}$ by a factor of several, and we quote $\alpha^{\mathrm{Max}}\sim 6\times 10^{-3}$ for $Am=3$ as a reference value. Below, we focus on the influence of the planet on the $\alpha$ profiles.

The radial profile of $\alpha^{\rm Max}$ (and hence $\alpha$) largely reflects magnetic flux concentration, which peaks at the radii where magnetic flux concentrates at $R=R_{\mathrm{P}}$ as well as in the planet-free gaps. With both density depletion and magnetic flux concentration, effective magnetization is greatly enhanced in the gap region, leading to a very large $\alpha\gtrsim 0.1$. This value increases with increasing planet mass, reflecting deeper density gap. Similar behavior was found in previous unstratified sharing-box simulations with net poloidal magnetic flux (Zhu et al., 2013; Carballido et al., 2016, 2017). These works also found that the radial profile of the Maxwell stress (not normalized by pressure) is largely flat across the gap. In other words, the larger $\alpha^{\mathrm{Max}}$ is mainly due to the lower density within the gap. In our simulation, as shown in the upper panels of figure 7, the Maxwell stress in fact increases from the inner to the outer sides of the gap. This is mainly because poloidal fields are inclined radially outwards at higher altitudes accompanying wind launching, which greatly enhances $\alpha^{\mathrm{Max}}$ at outer radii. This fact will be further discussed in Section 5.3. On the other hand, we observe that for the Reynolds stress, $\alpha^{\mathrm{Rey}}$ becomes negative around the planet orbit, as well as in the planet-free gap, which can be partly attributed to the strong accretion flow (see Section 4.1) in the midplane region in addition to the density waves/shocks.³³3The Reynolds stress has major contributions from the density wave/shock exhibiting strong azimuthal variations. We find that another major contribution arises from the highly non-uniform distribution of the accretion flow concentrated in the midplane region with large accretion velocities (see Section 4.1) that yields a negative $\delta(\rho v_{r})$. Given our definition of the Reynolds stress (see Equation 31), the midplane region also tends to share a small excess of (positive) $\delta j$ thanks to the $\sin\theta$ factor in $j$. The combined effect thus promotes a negative Reynolds stress in the midplane.

To further look into turbulence and stresses, we show in figure 8 the vertical profiles of the velocity fluctuations (upper) and $\alpha_{\theta}\equiv T_{r\phi}/P$ (lower) at $r=3$ (left), $6$ (center) and $10$ (right), for all three runs with $Am=3$, where the location $r=10$ is chosen to be well outside of the planet-induced density gap. There are several noticeable features. First, the value of $\alpha^{\mathrm{Max}}_{\theta}$ increases from the midplane to disk surface, consistent with standard expectations (e.g., Simon et al., 2013; Bai, 2015), and it largely dominates the total stress at all heights and both radial locations. Second, the vertical profile of the Reynolds stress can be oscillatory, with amplitudes comparable to the Maxwell stress. This can be related to the periodic passing of the density shocks, though its net contribution to angular momentum transport is reduced through vertical averaging. At highest altitude, the Reynolds stress is dominated by the wind. Third, there are strong velocity fluctuations. At $r=10$, the level of velocity fluctuations in all cases are similar and are stronger than the planet-free case reported in Cui & Bai (2021) (see their figure 11), by a factor of $\gtrsim 3$, and the dominant component is $v_{r}$. We anticipate that the main contribution is from the spiral shocks rather than the MRI turbulence, given that the physical conditions remains similar to the planet-free case at that radius. The gap region is characterized by even stronger velocity fluctuations, with $\delta v$ increasing towards higher planet mass, and reaching $\sim c_{s}$ for $M_{p}\gtrsim 3M_{t}$. In fact, the flow structure in the planet orbital region is very complex, which we address further in the next section. The planet-free gap at $r=3$ (for $M_{\mathrm{P}}\leq 3M_{\mathrm{th}}$) also has strong velocity fluctuations with a positive correlation to planetary mass, which also indicates that the planet-induced density waves/shocks, including the secondaries, can enhance turbulence in these strongly magnetized planet-free gaps not far from the planet orbit. As the planet-free gap already merges with the planetary gap in the $M_{\mathrm{P}}=5M_{\mathrm{th}}$ run, the magnetization at $r=3$ is weaker, thus showing weaker velocity fluctuations there.

Figure 9 shows the radial profiles of accretion rate and rotational velocity (density weighted over the bulk disk). Without the planet, given the radial density gradient and $h=0.1$, we anticipate $|v_{\phi}-v_{\mathrm{K}}|$ to be $0.98~{}v_{\mathrm{K}}$. With the planet, we see that the disk rotating velocity is super- and sub-Keplerian outside and inside the planet, respectively. The deviation from the local Keplerian velocity increases with the planetary mass. The maximum and minimum $v_{\phi}/v_{\mathrm{K}}$ are measured as $(1.02,0.95)$, $(1.06,0.92)$, and $(1.07,0.89)$ for runs Mt1Am3, Mt3Am3, and Mt5Am3, respectively. We also note the sub-/super-Keplerian modulation due to the planet-free gap. The variation in $v_{\phi}$ is at a similar amplitude as the planet-induced gap. This suggests that such variations are insufficient to distinguish between gap formation mechanisms.

The accretion rate profile is not flat in our simulations, implying that a steady state is not reached. This is not a major concern for us because firstly, in wind-driven accretion, the accretion rate profile is largely set by magnetic flux distribution (e.g., Bai et al., 2016; Lesur, 2021), which is not known a priori⁴⁴4Our initial radial profile of surface density is relatively steep, leading to relatively large amount of magnetic flux in the inner disk, and hence higher accretion rate in the inner than outer disk regions.. In addition, long-term disk evolution is subject to the not-so-predictable nature of magnetic flux evolution, which we are not in a position to address in this work. Therefore, we mainly focus on the main stage of gap formation, with more detailed analysis in Section 5. We speculate the gap profile may undergo some long-term secular evolution, though this is beyond the scope of this work.

Focusing on the planet orbital region, the mass accretion rate is on the order of $1$–$2\times 10^{-3}~{}[2\pi R_{\mathrm{P}}^{2}\Sigma_{\mathrm{P}}\tau_{\mathrm{P}}^{-1}]$, where $\Sigma_{\mathrm{P}}=\Sigma_{0}(R_{\mathrm{P}})$, and remains at a similar level as the accretion rates at neighboring radii. With the surface density depleted by a factor of ten to hundreds (see figure 4), the mean accretion velocity increases by a similar factor to cancel out density depletion. While there are accretion rate bumps in the planet orbital region, they are related to the exclusion of the planetary Hill sphere. Further incorporating the Hill sphere indicates roughly constant accretion flow across the planetary gap region (dotted black lines in figure 9). Note that without a mass sink at the planet, there is relatively large mass accumulation within the circumplanetary region, giving substantial mass weighting in the measurements, as also reflected in the $v_{\phi}$ profiles. We will study the dynamics of this region with further mesh refinement and sink prescriptions in a follow up work.

Next, we look into more details of the global flow structure and its relation to the magnetic field configuration. Figure 10 shows the vertical profiles of mean velocity and magnetic fields at the planetary gap (r=6), at some further distance away from the gap ($r=10$), and at the planet-free gap ($r=3$), and figure 11 shows the azimuthally-averaged configurations of the flow and magnetic field properties. At the gaps with magnetic flux concentration ($r=3$ and 6), we see that a strong inflow appears at the midplane (run Mt3Am3 and Mt5Am3) or in the disk upper layer (run Mt1Am3) to achieve the mass accretion rate. For wind-driven accretion, the accretion mass flux is directly proportional to the vertical gradient of toroidal field, and is the strongest when $B_{\phi}$ changes sign (e.g., Bai & Stone, 2013). This can be clearly seen by comparing the 2nd and 3rd rows in figure 11, as well as in figure 10. We note that the location where $B_{\phi}$ flips in the presence of the MRI turbulence is generally stochastic and evolves over time (Cui & Bai, 2021), whereas in runs Mt3Am3 and Mt5Am3, the strongly magnetized nature in the gap region (with plasma $\beta\lesssim 1$) makes toroidal field to flip in the midplane, confining the accretion flow there (see also Martel & Lesur, 2022). In this case, the accretion velocity becomes on the order of the sound speed.

The segregation of magnetic flux due to magnetic flux concentration into the planet gap also strongly affect the disk wind properties. It is generally expected that the wind is launched from the disk surface around the height (here $|z|\sim 3.5H$) where the flow approaches the ideal MHD regime (Bai & Stone, 2013; Gressel et al., 2015), above which various conservation laws along poloidal field lines apply. As we see from figure 11, this is still largely the case, except that the outflow from the planetary gap region can be launched from lower height at $z\sim 2$–$3H$. This is likely because of the strongly magnetized nature there, where the Lorentz force starts to dominate at lower heights.

The outflow launched from the planetary gap is asymmetric. As poloidal magnetic fields bend radially outward from the bulk disk, the magnetic flux threading the disk gap region, when approaching the wind base, is already outside of the planet radius. Therefore, most of the wind mass flux arises from the outer gap edge, as can be clearly identified in figure 11. This means that mass loss and angular momentum extraction from the wind is strongly concentrated on the outer side of the gap.

Thanks to the asymmetry, the wind from the planetary gap region is launched with super-Keplerian rotation. Together with the strong magnetic flux concentration, the wind efficiency is strongly enhanced (Bai et al., 2016). We also show the location of the Alfvén surface in the top panel of figure 11, where the poloidal wind speed matches the poloidal Alfvén velocity. Along a poloidal field line, the ratio of the Alfvén radii to wind launching radii ($R_{\mathrm{A}}/R_{0w}$, in cylindrical coordinates) characterizes the efficiency of angular momentum transport. In standard form (e.g., Ferreira & Pelletier, 1995; Bai et al., 2016), one expects

where the left shows the ratio of wind mass loss rate (per logarithmic radius) to the wind-driven accretion rate.

Quantitatively analyzing the wind properties based on the above is, however, no longer straightforward in our simulations. This is for several reasons. We first notice from figure 9 that the local wind mass loss rate, when measured through the mass flux normal to the surface of “wind base" at $z=0.35H$ (which is along the $\hat{\theta}$ direction), can become negative. This is also evident from figure 11, where we see the negative mass flux is mainly in regions where magnetic flux is sparse. It reflects the limitation of analyzing wind in spherical coordinates when the disk is not very thin, as the wind streamlines can be nearly radial, especially in regions with very little magnetic flux. Second, the wind launching radius $R_{\mathrm{0w}}$ may not be well characterized, especially in the planet gap region due to its strong magnetization as mentioned earlier. Lastly, we also show the location of the Alfvén surface in the top panels of figure 11. This surface is meaningful mainly in regions above the wind base as it is derived from the conservation laws under ideal MHD. We see that this surface is only well-defined in regions originating from strong magnetic flux concentration. In regions with magnetic flux deficit, this surface is below the wind base and becomes invalid (not shown).

Here, we seek for qualitative understanding of the outflow properties. The most prominent feature is that the wind is mostly launched from regions with magnetic flux concentration. This can be clearly identified from the streamlines in the bottom panel of figure 11, and is also accompanied by super-Keplerian $v_{\phi}$ in the wind zone as expected. This applies to both the planet gap and the planet-free gaps around $R=3$ in runs Mt1Am3 and Mt3Am3. The wind still shows substantial mass loading, as can be seen in the mass loss rate measured slightly outward of the gap region (as field lines bend radially outward). We may still apply Equation (41) in the wind region for a rough examination of wind properties. From the top panels of figure 11 which mark the Alfvén surface, we may approximately estimate the ratio of $R_{A}/R_{\mathrm{0w}}$ to be $\lesssim 1.3$, $\lesssim 1.3$, and $\lesssim 1.4$ for field lines originated from the planet gaps in runs with $M_{P}=1,3,5M_{\mathrm{th}}$, respectively. The small lever arm is consistent with the high value of $\mathrm{d}\dot{M}_{\rm wind}/\mathrm{d}\ln R/\dot{M}_{\mathrm{acc}}$ seen in figure 9.

There are essentially no poloidal magnetic field lines originate between the gaps, and consequently, no outgoing flow streamlines, in line with the wind mass loss profile shown in figure 9. Instead, wind launched from smaller radii essentially flows over such ring-like regions (at $R\sim 4\text{--}6$), nearly along the radial direction. In contrast to the strong accretion flow associated with the wind, we notice from figure 11 that the flow structures inward and outward of the planetary gaps in the bulk disk are more stochastic, lacking the characteristic wind-driven accretion flows associated with the vertical gradient of $B_{\phi}$, but instead exhibit circulation patterns in the meridional plane (which are different from the expected meridional flow seen in hydrodynamic simulations, see next subsection).

Previous 3D hydrodynamic simulations found that massive planets trigger meridional circulation in the gap region: the flow gets pushed away from the planet at the midplane, lifts up, then approaches the planet at the high altitudes, and vertically falls onto the planet’s orbital region (Morbidelli et al., 2014; Szulágyi et al., 2014; Fung & Chiang, 2016). There is also observational evidence for such meridional patterns(Teague et al., 2019; Yu et al., 2021).

In the windy disk, however, we do not see obvious signs of this meridional circulation. When averaged over $\phi$, we have already discussed that the flow pattern around the planet is primarily governed by MHD forces, which drives strong outflow, downflow, and accretion flow, etc. We do notice that in the case of run Mt3Am3 and Mt5Am3, there are flows from the gap outer edge towards the planet orbital region at up to $\pm H$, which might disguise as meridional flow if that region happens to be the emission surface. However, at upper altitude at the gap outer edge, gas consistently fly away in the disk outflow. Also note that there is higher-altitude ($z\sim\pm 2H$) flows towards the planet orbital region from the inner side of the gap, which are also similar to meridional flow, though they are not present at lower altitudes.

On the other hand, we identify that such flow pattern occurs only around the planet itself, namely, $\phi\sim 0$ (rather than the entire orbital/gap region). At the $\phi=0$ slice where the planet locates, we find that the meridional flow extends to no more than $\pm 2H$, where strong outflow still takes over at higher altitudes ($z\gtrsim 2H$). As this pattern is more closely related to circumplanetary disks, we will defer for more detailed study in a follow-up paper.

Typically, in a viscous disk, the fluid in the planetary co-orbital region executes a horseshoe orbit with U-turn at the planet vicinity. Since the fluid exchanges angular momentum with the planet at this turn, this horseshoe region is crucial for understanding planet migration (e.g., Ward, 1991; Masset, 2001; Paardekooper & Papaloizou, 2009).

In the windy disk, however, this circulation in the corotation region is prevented by the fast radial accretion flow discussed in § 4.1. Figure 12 shows the horizontal velocity vector in the planet-corotating frame (arrow) and radial velocity (color in the arrows) in $z$-integrated $R$-$\phi$ plane. The velocity vector within the corotation region clearly directs radially inward (upper), while the counterpart in 2D-viscous runs shows the horseshoe circulation (lower). This can be understood by comparing the timescales of the horseshoe libration and Hill-region passage. The libration timescale can be, neglecting the pressure gradient, estimated as (Paardekooper & Papaloizou, 2009)

where $\Delta R$ is the half-width of the horseshoe region. Since $\Delta R\approx r_{\mathrm{H}}$ when $M_{\mathrm{P}}\gtrsim M_{\mathrm{th}}$ (Masset et al., 2006; Paardekooper & Papaloizou, 2009), $\tau_{\mathrm{lib}}\sim 19~{}(M_{\mathrm{P}}/M_{\mathrm{th}})^{-1/3}~{}\tau_{\mathrm{P}}$ in our disk model. On the other hand, the mean radial velocity within the corotation region is measured as to be $v_{R}\sim-0.002$ (Mt1Am3), $-0.01$ (Mt3Am3), and $-0.06~{}v_{\mathrm{K}}$ (Mt5A3), resulting in the timescale passing the corotation region to be $\tau_{\mathrm{pass}}=16$, $3.2$, and $0.53~{}\tau_{\mathrm{P}}$, all of which are shorter than $\tau_{\mathrm{lib}}$. Namely, the fluid passes the corotation region before being able to finish the horseshoe orbit.

In addition, we note that there is outward radial flow at the planetary tailing side ($-\pi/4<\phi<0$), which coincides with regions of negative $B_{z\mathrm{,mid}}$ (see § 3.3, and the black contours of figure 12). We find that in this localized region, $B_{\phi}$ remains relatively smooth, therefore, a negative $B_{z\mathrm{,mid}}$ leads to the reversal of the Lorentz force, and hence this radially-outward gas flow. We highlight that the presence of this localized radial outflow is a robust consequence of the negative $B_{z\mathrm{,mid}}$, which itself is a robust phenomenon that we speculate to be associated with magnetic flux concentration. As a result, the flow in the trailing side of the gap exhibits certain localized circulation between $-\pi/2<\phi<0$. These patterns differ dramatically from the viscous counterpart that clearly show the horseshoe orbits, and could have important implications to planet migration (see Section 6).

We start from a brief introduction of the torque balance in the viscous accretion disk of the Mt3$\alpha$6 run. The torques in this viscous disk is shown in the left panel of figure 13. Due to the spiral density waves/shocks, the planet gravity torque is positive and negative at outward and inward of the planet, respectively, indicating that $\gamma_{\mathrm{P}}$ drives gap opening. Close to the planet ($5\lesssim R\lesssim 7$), $\gamma_{\mathrm{wave}}$ has the opposite sign and a similar amplitude to $\gamma_{\mathrm{P}}$, leading to torque balance. The total planet torque $\gamma_{\mathrm{P}}+\gamma_{\mathrm{wave}}$ is deposited at further distances (e.g., Takeuchi et al., 1996; Goodman & Rafikov, 2001; Crida et al., 2006; Kanagawa et al., 2015, 2017), in this case at about $3\lesssim R\lesssim 5$ and $7\lesssim R\lesssim 12$. The deposited planetary torque is then balanced by the viscous torque $\gamma_{\mathrm{vis}}$ to fill in the gap. It is the balance of these torques that determine the gap structure (e.g., Lin et al., 1993). Note that the torque density associated with accretion $\gamma_{\mathrm{A}}$ is negative here, namely, the disk gas flows outward. This is because of the relatively steep surface density slope adopted here. Nevertheless, this component is negligible in the torque balance in the planetary gap and hardly affects the gap structure.

From now on, we focus on our fiducial run (Mt3Am3), where the relevant torques are shown in the right panel of figure 13. Here, the Maxwell stress ($\gamma_{\mathrm{Max}}$) plays the role of viscosity ($\gamma_{\mathrm{vis}}$). In addition, the MHD-driven wind torque ($\gamma_{\mathrm{wind}}$) is added to the analysis.

Firstly, we see that $\gamma_{\mathrm{P}}$ and $\gamma_{\mathrm{wave}}$ share similar shapes compared to the viscous case, but have smaller amplitudes especially in the vicinity of the planet orbit. Therefore, the physics of gap opening by the planet remains similar. As the torque density is proportional to the surface density of the gas, the smaller amplitudes is primarily due to the fact that the gap in windy disks are deeper and wider. As a result, the overall torque balance is largely determined by the MHD torques and the accretion flow, which are strongly affected by magnetic flux concentration.

The wind torque $\gamma_{\mathrm{wind}}$ mainly dominates in regions outward of the planet gap. As has been discussed in Section 4.2, this is the result of field lines bent radially outward at higher altitudes: concentrated poloidal magnetic fields in the gap region emanate from the disk surface from the outer side of the gap (see figure 11). Within the gap, the wind torque diminishes as there is a deficit of magnetic flux originating from inward of the planet orbit. We also see the additional bump in the wind torque at $r\sim 4$, which results from magnetic flux concentration in the planet-free gap at $r\sim 3$.

The radial profile of the torque density resulting from the Maxwell stress is more complex, which is a result of the radial gradient of the angular momentum flux. We can see from the top panels of figure 7, that although the $\alpha^{\mathrm{Max}}$ profiles peaks within the gap, the unnormalized Maxwell stress itself peaks at outer radii. This leads to a negative torque density in the gap region, and a positive torque density beyond the gap. Although $\mathrm{d}/\mathrm{d}r(J_{\mathrm{Max}})$ is the counterpart of $\mathrm{d}/\mathrm{d}r(J_{\mathrm{vis}})$, it plays a very different role because the highly inhomogeneous radial distribution of magnetic flux. In particular, in the gap region, the net angular momentum flux is nearly entirely from this component, which drives gas accretion across the gap, leading to gap depletion. Outward of the gap, it is mainly the radial gradient of this $J_{\mathrm{Max}}$ that balances the wind torque.

Piecing together, the aforementioned torques drive the accretion flow through the gap. We note that the torque density contributed from accretion $\gamma_{\mathrm{A}}$ varies substantially across the gap region, in contrast to the relatively smooth variation in the accretion rate (see figure 9). This can be accounted for by the modulation in disk rotation rate (see Equation (24)). For instance, $\mathrm{d}\langle j\rangle/\mathrm{d}r$ is steeper in the gap region than the Keplerian case, leading to enhanced $\gamma_{\mathrm{A}}$ in the gap.

The comparison between the left and right panels in figure 13 infers why the gap is deeper and wider in the windy disk than in the viscous disk. In the viscous disk, the gap is carved by the planet-induced torques ($\gamma_{\mathrm{P}}+\gamma_{\mathrm{wave}}$). Since the planet-induced torque is proportional to density while the viscous torque scales with the density gradient, there is a limit on gap depth once the gap becomes deep enough. In contrast, in the windy disk, gap opening is primarily due to MHD torques boosted by magnetic flux concentration at the planetary orbit. Since the $\mathrm{d}J_{\mathrm{Max}}/\mathrm{d}r$ is insensitive to gas density, it keeps depleting the gap regardless the gap depth. We will further show in the next subsection that the role of $J_{\mathrm{Max}}$ and $J_{\mathrm{wind}}$ are in fact the two sides of the same coin, thus this gap depletion can be equally understood as angular momentum removal from the wind.

As discussed in the previous subsection, the local torque balance in the windy disk is greatly different from that in the viscous case, primarily due to the MHD torques. However, the inhomogeneity in the torque profiles adds to the complexity against a clear understanding of the underlying physics. In this subsection, we aim to provide a more transparent picture of the torque balance, focusing on the MHD torques.

In figure 14, we show the cumulative torque from the two components of the MHD torques, corresponding to the radial angular momentum flux $J_{\mathrm{Max}}$, and the radially-integrated (negative) wind torque $\Gamma_{\rm wind}$. We similarly also compute the cumulative (negative) torques from other components (planetary gravity $\Gamma_{P}$ and the wave angular momentum flux $J_{\mathrm{wave}}$). Note that they share the same dimensions, and positive/negative slopes correspond to negative/positive torque density. Here, the cumulative torques $\Gamma$ are set to zero at the planetary orbit $r=6$.

The radial profile of $J_{\mathrm{Max}}$ peaks at about the outer gap edge as a result of magnetic flux concentration, similar to that shown and discussed in the top panels of figure 7. As its value falls off beyond the outer gap edge, the wind torque picks up, as a result of winds emanating from included field lines emanating the outer gap edge. In fact, the increase of the wind torque largely cancels the decrease of the radial angular momentum flux $J_{\mathrm{Max}}$. To see this, we show in black lines the total cumulative MHD torque $J_{\mathrm{MHD}}\equiv J_{\mathrm{Max}}+\Gamma_{\rm wind}$. We find that beyond the planetary gap, $J_{\mathrm{MHD}}$ approximately follows a $r^{1/2}$ scaling, which extends all the way to larger radii. This is consistent with a standard wind torque balancing steady state accretion, which is proportional to $\dot{M}_{\mathrm{acc}}\mathrm{d}/\mathrm{d}r(r^{2}\Omega_{\mathrm{K}})\propto\sqrt{r}$. At the planetary gap region, there is little contribution from $\gamma_{\mathrm{wind}}$ as we discussed in the previous subsection, the MHD torque is dominated by the rising side of $J_{\mathrm{Max}}$, with a steeper slope.

Inward of the planetary gap, we see that the slope of $J_{\mathrm{MHD}}$ first remains similar to that outward of the planetary gap. At around $r=3$ where the planet-free gap is located, there is another bump in $J_{\mathrm{Max}}$ as a result of magnetic flux concentration there. Again, the fall-off of the bump is compensated from the rise in $\Gamma_{\rm wind}$. Overall, $J_{\mathrm{MHD}}$ develops another steeper slope in the planet-free gap region from the rising side of $J_{\mathrm{Max}}$, similar to the planetary gap case.

To summarize the analysis above, we find that despite of the inhomogeneities in individual components of the MHD torque, the combined effect of radial and vertical angular momentum is in fact quite simple: away from the gap region, the total MHD torque behaves as a standard wind torque, whereas in the gap region, the wind torque gets enhanced due to magnetic flux concentration. To be more quantitative, when we normalize the torque density (i.e., the slope in figure 14) by $R_{\mathrm{P}}^{2}\Sigma_{\rm P}v_{\mathrm{P}}\Omega_{\mathrm{P}}$, we find the slope in $J_{\mathrm{MHD}}$ is about $2.3\times 10^{-3}$, as opposed to $5.0\times 10^{-4}$ based on the $r^{1/2}$ fit outward of the gap region.⁵⁵5Note that the $J$ in the figure 14 is normalized by $R_{\mathrm{P}}^{3}\Sigma_{\mathrm{P}}v_{\mathrm{P}},\Omega_{\mathrm{P}}$, which differs by $R_{\mathrm{P}}=6$ from the dimension of torque-density we use here. There is a factor of $\sim 5$ difference, which is largely owing to magnetic flux concentration. Also note that as the $\dot{M}_{\mathrm{acc}}$ is largely constant across the gap, this enhancement of torque density is mainly compensated by the steeper $j$ gradient within the gap (see figure 9 and discussion in the previous subsection).

To better understand why the stresses that correspond to radial and vertical transport of angular momentum are so smoothly connected that leads to the simple picture above, we show in figure 15 the vectors of the angular momentum flux from the poloidal Maxwell stress $-{\boldsymbol{B}}_{p}B_{\phi}$. We see that in regions with magnetic flux concentration, the angular momentum flows along the poloidal magnetic field lines. This is characteristic of angular momentum extraction by disk winds. However, the standard diagnostics decomposes this flux into the radial ($-B_{r}B_{\phi}$) and vertical ($-B_{\theta}B_{\phi}$) components, treating them separately as viscous and wind stresses. This approach artificially separates a single piece of physics into two disparate components, causing certain degree of conceptual confusion. In this figure, we also see that, outward of gap regions, angular momentum transport is largely radially outward in the bulk disk (despite of the relatively small angular momentum flux), conforming to the standard scenario of viscous transport. In the meantime, there is additional wind transport that originate from the disk surface. Therefore, the overall situation in such regions is similar to the planet-free case (e.g. Cui & Bai, 2021).

In this subsection, we further examine the torque balance in the other two runs Mt1Am3 and Mt5Am3. The radial profiles of the torque densities are shown in figure 24 in Appendix D, which are qualitatively similar to that of the fiducial run. Here, we mainly focus on the MHD torques, showing the resulting cumulative torques in figure 17.

The role of the MHD torque on gap formation is largely the same for three planetary masses in our simulations, namely, there is an enhanced torque in the gap region, and the torque becomes similar to a standard wind torque beyond the gap. This is likely the reason why the gap profile is similar to that in the inviscid disk (see figure 16). Interestingly, the enhanced slope of $J_{\mathrm{MHD}}$ in the gap region is quantitatively similar in different planetary masses ($3.8$, $3.8$, and $4.0\times$ $10^{-4}R_{\mathrm{P}}^{3}\Sigma_{0}(R_{\mathrm{P}})v_{\mathrm{P}}\Omega_{\mathrm{P}}$ for 1, 3, and 5 $M_{\mathrm{th}}$, respectively). Correspondingly, the negative torque acting on the gap center is same for the three cases within $5\%$. This is consistent with the fact that the $B_{z\mathrm{,mid}}$ has a similar magnitude around the planet in different mass cases (see figure 3). Moreover, figures 14 and 17 show that the radial range where the MHD torque is enhanced approximately scales with the Hill length $r_{\mathrm{H}}$, consistent with radial range where $B_{z\mathrm{,mid}}$ is enhanced through magnetic flux concentration (see figure 3).

For run Mt5Am3, we note that the radial range with enhanced MHD torques extends further outward ($\sim 2r_{\mathrm{H}}$) relative to the other two runs ($\sim r_{\mathrm{H}}$ for runs Mt1Am3 and Mt3Am3). However, we find that such difference appears only after the planet-free gap merges to the planetary gap ($t\gtrsim 110~{}\tau_{\mathrm{P}}$; see § 3.2). After the merger, the planetary gap gathers more poloidal magnetic fluxes, which enhances both the $B_{z\mathrm{,mid}}$ magnitude and its radial extent (see the red line of the bottom right panel in figure 3). Before the merger, both the slope and radial extent are more similar to the other runs.