Hall Effect–Mediated Magnetic Flux Transport in Protoplanetary Disks

We use Athena++, a newly developed, grid-based, higher-order Godunov MHD code with constrained transport to conserve the divergence-free condition for magnetic fields (J. M. Stone et al. 2016, in preparation). It is the successor of the widely used Athena MHD code (Gardiner & Stone 2005, 2008; Stone et al. 2008), and is highly optimized in several aspects. In particular, it employs flexible grid spacings, allowing simulations to be performed over large dynamical ranges. Geometric source terms in curvilinear coordinate systems (e.g., cylindrical and spherical-polar coordinates) are carefully implemented, which ensures exact angular momentum conservation.

Using Athena++, we solve the standard MHD equations in conservation form

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{v}}{\boldsymbol{v}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{4\pi }+{{\mathsf{P}}}^{* }\right)=-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot \left[(E+{P}^{* }){\boldsymbol{v}}-\displaystyle \frac{{\boldsymbol{B}}({\boldsymbol{B}}\cdot {\boldsymbol{v}})}{4\pi }\right]=-{\rm{\Lambda }},\end{eqnarray} \tag{ 13 }$$

where P is gas pressure, ${P}^{* }=P+{B}^{2}/8\pi$ is total pressure, $E=P/(\gamma -1)+\rho {v}^{2}/2\,+\,{B}^{2}/8\pi$ is total energy density, γ is the adiabatic index, ${{\mathsf{P}}}^{* }\equiv {P}^{* }{\mathsf{I}}$ with ${\mathsf{I}}$ being the identity tensor, ${\rm{\Phi }}=-{GM}/r$ is the gravitational potential of the protostar, and Λ is the cooling rate. These equations are coupled with the induction Equation (4), which incorporates non-ideal MHD effects, to evolve the magnetic field. Also note that in the code, factors of $4\pi$ are absorbed into the definition of B so that magnetic permeability is $\mu =1$.

With these advantages, we perform 2D global MHD simulations of PPDs in spherical polar coordinates ($r\mbox{--}\theta$). The radial grid spans from r = 1 to 100 in code units with logarithmic grid spacing. The θ grid extends from the midplane all the way to near the poles (leaving only a 2° cone at each pole), with nonuniform grid spacing where ${\rm{\Delta }}\theta$ increases by a constant factor per grid cell from midplane to pole, with a contrasting factor of four between the midplane and the polar region. This allows us to properly resolve the disk, and in the meantime accommodate the MHD disk wind so that the simulation results are not affected by the outer boundary conditions. Note that in a spherical polar grid, θ increases from 0 to π from the north to the south pole. For notational convenience, we define $\delta \equiv \pi /2-\theta$, namely the elevation angle about the midplane.

For the initial condition, we adopt a self-similar disk density and temperature profiles in the form of $\rho ={\rho }_{0}{(r/{R}_{0})}^{-\alpha }\,f(\theta )$, $T=P/\rho ={T}_{0}{(r/{R}_{0})}^{-1}g(\theta )$. In code units, we set ${\rho }_{0}={T}_{0}={R}_{0}=1$. Furthermore, we take GM = 1 for the gravity of the central protostar. Once we specify the dimensionless function $g(\theta )$, which corresponds to the square of local disk aspect ratio ${(H/r)}^{2}$, the density profile can be solved assuming hydrostatic equilibrium. In the θ direction it yields

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}F}{d\,\mathrm{ln}\,\sin \,\theta }=\displaystyle \frac{{GM}}{{T}_{0}{R}_{0}}\displaystyle \frac{1}{g(\theta )}-(\alpha +1),\end{eqnarray} \tag{ 14 }$$

where we have defined $F(\theta )\equiv f(\theta )g(\theta )$. Hydrostatic equilibrium is possible when the right side is positive, which limits the maximum value of $g(\theta )$. Force balance in the radial direction determines v_ϕ:

$$\begin{eqnarray}&&{v}_{\phi }^{2}=\displaystyle \frac{{GM}-(\alpha +1){T}_{0}{R}_{0}g(\theta )}{r}.\end{eqnarray} \tag{ 15 }$$

Throughout this paper, we choose $\alpha =2$, corresponding to a surface density profile of ${\rm{\Sigma }}\propto {r}^{-1}$.

In this work, we consider modestly thin disks with a midplane aspect ratio ${H}_{\mathrm{mid}}/r=0.1$ within $\sim 2{H}_{\mathrm{mid}}$ about the disk midplane. It then increases smoothly toward disk surface and reaches $H/r=0.5$ near the pole.⁴ This allows the gas density to drop much more slowly with θ and hence alleviates the timestep constraint in the MHD wind zone. Physically, this is motivated by the fact that the wind zone is strongly heated by external UV and X-rays and becomes significantly hotter than the disk interior (e.g., Glassgold et al. 2004; Walsh et al. 2010). Part of the temperature profile $g(\theta )$ can be found in the left panels of Figure 3. We adopt an ideal gas equation of state with adiabatic index $\gamma =5/3$, which is appropriate for atomic gas in the wind zone, but the results are insensitive to the choice of γ. We seek for simple prescriptions of thermodynamics, achieved by a simple cooling prescription (the Λ term) that relaxes gas temperature to the initial value at the rate of local Keplerian frequency (based on spherical radius r). This prescription also avoids the development of hydrodynamic instabilities (e.g., Nelson et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Poloidal magnetic fields are initialized with vector potential generalized from Zanni et al. (2007)

$$\begin{eqnarray}&&{A}_{\phi }{(r,\theta )=\displaystyle \frac{2{B}_{z0}{R}_{0}}{3-\alpha }{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{-\tfrac{\alpha -1}{2}}{[1+(m\tan \theta )}^{-2}]}^{-\tfrac{5}{8}},\end{eqnarray} \tag{ 17 }$$

where $R\equiv r\sin \theta$, and m is a parameter that specifies the degree that poloidal fields bend, with $m\to \infty$ giving a pure vertical field. Poloidal field is given by ${\boldsymbol{B}}={\rm{\nabla }}\times ({A}_{\phi }\hat{\phi })$, so that in the midplane, ${\boldsymbol{B}}={B}_{z0}\hat{z}{(r/{R}_{0})}^{-(\alpha +1)/2}$, maintaining a constant ratio of gas to magnetic pressure, defined by plasma ${\beta }_{0}$. In practice, we choose m = 0.5, and we have tested that the results are insensitive to the choice of m. Fiducially, we choose ${\beta }_{0}={10}^{4}$, appropriate for the outer region of PPDs (Simon et al. 2013a; Bai 2015), but we also consider stronger and weaker fields in Section 7.

We have implemented all three non-ideal MHD terms in Athena++. In particular, Ohmic resistivity and AD are implemented using operator splitting, as in the original Athena code (Bai & Stone 2011) with super-timestepping (Simon et al. 2013b). We have tested the operator-split implementation of the Hall term following Bai (2014), which was shown to be marginally stable in Cartesian coordinates, and found that it becomes unstable in spherical coordinates. We thus adopt the non-operator-split implementation following Lesur et al. (2014), where the Hall term is incorporated to the Harten–Lax–van Leer (HLL) Riemann solver. While the HLL solver is very diffusive, recent shearing-box simulations using this method (Lesur et al. 2014; Simon et al. 2015) have yielded results consistent with Bai (2014, 2015), who used operator-split and the more accurate HLLD solver (Miyoshi & Kusano 2005).

In this paper, we focus on regions of PPDs where the Hall effect and AD are the dominant non-ideal MHD effects, and only include the two terms in our simulations. It typically corresponds to regions of intermediate radii ($r\sim 5\mbox{--}30$ au) (Wardle 2007; Bai 2011a). We set Am = 0.5 throughout the disk zone (within $\delta \sim \pm 2({H}_{\mathrm{mid}}/r)$), which is motivated from ionization chemistry calculations (Bai 2011a, 2011b) where Am is found to be of order unity in the outer disk and becomes smaller toward the inner disk. We take Am to be on the low side, which helps suppress the MRI, or at least to significantly reduce its growth rate so that the gas remains largely laminar during the simulations. We further set ${l}_{H}=2{H}_{\mathrm{mid}}$ in the inner edge of the disk midplane, which is taken from estimates in Bai (2015). The value of l_H elsewhere is simply determined from the relation ${\eta }_{H}/{\eta }_{A}\propto {(B/\rho )}^{-1}$, which leads to ${l}_{H}\propto {r}^{1/2}\sqrt{\rho }$ for constant Am. Beyond $\delta \sim \pm 2({H}_{\mathrm{mid}}/r)$, we smoothly reduce both Hall and AD diffusivities to zero, mimicking the fact that external far-UV ionization substantially increases the ionization fraction, which brings the gas to the ideal MHD regime (Perez-Becker & Chiang 2011).

Our fiducial simulations are performed with 560 × 216 grid cells in $r\times \theta$. With nonuniform grid spacing, we achieve a resolution of about 16 cells per ${H}_{\mathrm{mid}}$ in θ and 12 cells per ${H}_{\mathrm{mid}}$ in r around the disk midplane. For the fiducial runs, we also conduct simulations using twice the resolution for convergence study. As an initial effort, our simulations are 2D instead of 3D, which on the one hand substantially reduces the computational cost, and moreover, with the MRI largely suppressed or damped due to strong AD (Bai & Stone 2011), we expect our 2D simulations to capture the essential aspects of disk dynamics.

The outer radial boundary follows from standard outflow boundary prescriptions, where hydrodynamic variables are copied from the last grid zone assuming $\rho \propto {r}^{-2}$, ${v}_{\phi }\propto {r}^{-1/2}$, with v_r and v_θ unchanged except that we set v_r = 0 in case of inflow. At the inner radial boundary, hydrodynamic variables are fixed to initial state. We make this choice because the flow near the polar region is presumably originated from the part of the disk that is located within the inner radial boundary, whose dynamics is beyond the reach of the simulation. An outflow-type boundary condition prescription would violate causality, which can become unstable in the presence of magnetic fields and further interfere with the wind flow in the main computational domain. The fixed state boundary condition alleviates the causality issue, and as the gas flow in our simulations is largely laminar, it also guarantees stability. This allows us to pursue our study without being affected from the inner boundary. Magnetic variables in the inner/outer ghost zones are copied from the nearest grid zone assuming ${B}_{r}\propto {r}^{-2}$ and ${B}_{\phi }\propto {r}^{-1}$, with B_θ unchanged. Moreover, we smoothly reduce the Hall diffusivity to zero within about ${H}_{\mathrm{mid}}$ from the inner radial boundary to avoid dramatic flux transport due to the Hall drift. Reflection boundary conditions are applied in the θ boundaries.

We list all our simulation runs in Table 1. We will focus on our fiducial runs labeled "Fid±" in the main text, where the +/− signs correspond to simulations with poloidal field aligned/anti-aligned with disk rotation. For comparison, we also conduct a run "Fid0," where we turn off the Hall effect. Results from our high-resolution run are discussed in Appendix A to address numerical convergence. We further discuss the dependence of flux transport rate on the poloidal field strength in Section 7.

Table 1.  List of Simulation Runs

Run	Polarity	Resolution	${H}_{\mathrm{mid}}/r$	${\beta }_{0}$
Fid+	+	560 × 216	0.1	104
Fid0	No Hall	560 × 216	0.1	104
Fid−	−	560 × 216	0.1	104
Fid-hires+	+	1104 × 432	0.1	104
Fid-hires−	−	1104 × 432	0.1	104
B3+	+	560 × 216	0.1	103
B30	No Hall	560 × 216	0.1	103
B3−	−	560 × 216	0.1	103
B5+	+	560 × 216	0.1	105
B50	No Hall	560 × 216	0.1	105
B5−	−	560 × 216	0.1	105

Download table as:  ASCIITypeset image

We are most interested in the evolution in the inner part of our simulation domain where the Hall effect is dominant at the midplane. The simulations are run for about 400 rotations at the innermost disk radius ($2400{{\rm{\Omega }}}_{0}^{-1}$, where ${{\rm{\Omega }}}_{0}=1$ is the Keplerian frequency at the innermost orbit). This is much longer than the dynamical timescale within $r\lesssim 10\mbox{--}15$ so that the flow structure (e.g., disk winds and accretion flow) is approximately steady. On top of such flow structure, magnetic flux evolves on longer timescales. We call such a situation a quasi-steady state, and over the course we can measure magnetic flux evolution. Moreover, with weak magnetization, the rate of angular momentum transport is relatively slow, and within the duration of these simulations, the surface density remains largely unchanged.

The early evolution follows exactly from the expectations discussed in Section 2. In the aligned case, the Hall effect efficiently transports magnetic flux inward at the midplane, leading to some flux accumulation at the inner radial boundary. Beyond the midplane, however, the toroidal field gradient reverses, and hence magnetic flux is transported outward. The above two processes stretch the poloidal field into a highly radially elongated configuration, as seen in the second snapshot. As the radial field grows, shear produces stronger toroidal field, which in turn leads to faster flux transport, and further growth of the radial field. In fact, the runaway process described here is a global manifestation of the Hall-shear instability (Kunz 2008), previously discussed in local shearing-box simulations (Bai 2014; Lesur et al. 2014). Saturation of this instability is owing to additional dissipation, and here AD in our simulations, which we will discuss in more detail in Section 6.2.

Upon saturation, the system achieves a quasi-steady state, with magnetic field bending sharply across the disk midplane. This configuration allows dissipative processes (here AD), which generally tend to straighten field lines, to effectively transport magnetic flux outward in the midplane region. It is the competition between inward transport by the Hall effect and outward transport by AD that determines the overall direction of the transport at the midplane: looking at the third snapshot, the direction of transport is pointing outward.

In the anti-aligned case, the opposite occurs. We see that first, outward transport of magnetic flux takes place at the disk midplane as a result of the shear-produced toroidal field gradient, as discussed in Section 2. The Hall drift pushes the poloidal field into an unusual configuration, which first bends radially inward and then outward. Later on, poloidal field lines around the midplane straighten. This was discussed in Bai (2014), where the situation is exactly the opposite of the aligned case: horizontal components of the field are reduced toward zero instead of undergoing runaway amplification.

Upon achieving a quasi-steady state, poloidal field lines near the midplane region are largely vertical, and the toroidal field gradient across the midplane is greatly reduced but nonzero. With this field configuration, magnetic flux is transported outward almost entirely due to the Hall effect at the midplane, with negligible contribution from AD. Toward the disk surface, on the one hand, the Hall effect weakens compared with AD because of the density drop, and on the other hand, the concavely shaped field configuration is more favorable for AD to transport magnetic flux outward. Overall, we find that magnetic flux is systematically transported outward due to the Hall effect at the midplane and AD at the surface.

Without including the Hall term in run Fid0, we see that besides the production of a toroidal field from the radial field, there is no significant change of magnetic field configuration from the initial condition. The system reaches a quasi-steady state relatively quickly, with toroidal field generation compensated by ambipolar dissipation in the disk region and outward advection in the disk wind. Despite the significant difference in field configuration, we find that magnetic flux is slowly transported outward at a rate comparable to the Fid ± runs.

In Figure 3, we show the θ profiles of various diagnostic quantities of interest measured at spherical radius r = 8 and time $t=1200{{\rm{\Omega }}}_{0}^{-1}$ (corresponding to the third snapshot in Figure 2). The first (left) column of panels shows the density and temperature (expressed in H/R) profiles. The latter is well preserved from the initial profile due to our cooling prescription. The density profiles in the bulk disk are largely hydrostatic. Beyond $\delta \sim \pm 3{H}_{\mathrm{mid}}/R$, the density profiles in the aligned and Hall-free cases deviate from the anti-aligned case. This is related to the fact that magnetic pressure exceeds gas pressure beyond about $\delta \sim \pm 3{H}_{\mathrm{mid}}/R$ in the aligned and Hall-free cases, as shown in the second column of panels.

The strength of the Hall and AD terms are measured by their respective Elsässer numbers, shown in the third panels from left. The AD Elsässer number Am is independent of field strength, and hence the profiles of Am are the same for all cases. As we specified in Section 3, Am is fixed to 0.5 in the main disk body and it increases toward infinity starting from $2{H}_{\mathrm{mid}}$ above/below the midplane. The Hall Elsässer number depends on field strength, and compared with AD, the Hall term is more important in a weaker field. Overall, in both aligned and anti-aligned cases, the Hall term dominates within about $1-2{H}_{\mathrm{mid}}$ about the midplane.

The fourth and fifth panels show the profiles for the three components of magnetic fields and velocities. They are useful for interpreting angular momentum transport in the next subsection. Also note that the azimuthal gas velocity is always sub-Keplerian due to pressure support.

In all cases, MHD disk winds are launched. The MHD winds extract disk angular momentum vertically via a wind stress ${T}_{z\phi }\equiv -{B}_{z}{B}_{\phi }/4\pi$ exerted at the disk surface (wind base), leading to an accretion rate $\dot{M}\approx (8\pi R/{\rm{\Omega }}){T}_{z\phi }$. Approximately, we specify ${\delta }_{b}=\pm 0.25$ to be the location of the wind base⁵ , which roughly corresponds to the location where ${Am}\sim 1$ and the gas transitions from being dominated by non-ideal MHD to satisfying ideal MHD conditions (which permits efficient wind launching).

We see from Figure 3 that in the aligned case, with B_R and B_ϕ amplified via the Hall shear instability, we measure ${T}_{z\phi }\approx 8.5\times {10}^{-3}{P}_{\mathrm{mid}}$ at r = 8, where ${P}_{\mathrm{mid}}$ is midplane gas pressure. Without the Hall effect, B_ϕ is amplified to similar strength but not B_R, which remains very small. We measure ${T}_{z\phi }\approx 5.8\times {10}^{-3}{P}_{\mathrm{mid}}$, where it is smaller mainly because of smaller B_z in this run as a result of magnetic flux evolution. In the anti-aligned case, while the horizontal field in the midplane is reduced, B_ϕ grows steadily toward the surface, and we measure ${T}_{z\phi }\approx 6.3\times {10}^{-3}{P}_{\mathrm{mid}}$ at the wind base, which is only slightly smaller than the aligned case. These results are similar to those discussed in Bai (2014).

The aligned case also produces a relatively strong Maxwell stress ${T}_{R\phi }=-{B}_{R}{B}_{\phi }/4\pi$ around the disk midplane, leading to non-negligible radial transport of angular momentum (magnetic braking). To an order of magnitude, the resulting accretion rate is $\dot{M}\sim (2\pi /{\rm{\Omega }}){\int }_{-{\theta }_{b}}^{{\theta }_{b}}{T}_{R\phi }{dz}$, which is about a factor of H/R less efficient than wind-driven accretion (for similar stress levels). Defining $\alpha ={\int }_{-{\theta }_{b}}^{{\theta }_{b}}{T}_{R\phi }{dz}/{\int }_{-{\theta }_{b}}^{{\theta }_{b}}{Pdz}$, we find $\alpha \approx 0.013$. Comparing α with ${T}_{z\phi }/{P}_{\mathrm{mid}}$, we see that α is not sufficiently large to offset the R/H factor to dominate angular momentum transport, consistent with local studies (Bai 2014). In the anti-aligned case, on the other hand, due to the reduction of the horizontal field, radial transport is completely negligible.

The accretion flow associated with the wind can be directly seen in the rightmost panels of Figures 2 and 3. Accretion is driven by the torque associated with the vertical gradient of the wind stress, which is strongest when B_ϕ varies the fastest. Examining the fourth panels from the left in Figure 3, it becomes clear why the mass flux of the accretion flow is mostly concentrated in the midplane in the aligned and Hall-free cases, whereas it is more uniformly distributed with modest concentration around $z=2{H}_{\mathrm{mid}}$ in the anti-aligned case. The total mass accretion rates in all three cases, on the other hand, are comparable because of their similar ${T}_{z\phi }$ values.

The global distribution of mass flux is best viewed from the rightmost panels of Figure 2. The efficiency of wind-driven accretion is characterized by the ratio of the mass-loss rate ${\dot{M}}_{\mathrm{wind}}$ to the wind-driven accretion rate ${\dot{M}}_{\mathrm{acc}}$. It is closely related to the location of the Alfvén surface, at which poloidal flow velocity is equal to the poloidal Alfvén velocity ${v}_{{Ap}}={B}_{p}/\sqrt{4\pi \rho }$. For wind launched from radius R₀, following the field line, one can define the Alfvén radius R_A to be the cylindrical radius at the Alfvén surface. Assuming ideal MHD, axisymmetry, and steady state, we have (e.g., Spruit 1996; Bai et al. 2016)

$$\begin{eqnarray}&&\displaystyle \frac{d{\dot{M}}_{\mathrm{wind}}/d\mathrm{ln}R}{{\dot{M}}_{\mathrm{acc}}}=\displaystyle \frac{1}{2}\displaystyle \frac{1}{{R}_{A}^{2}/{R}_{0}^{2}-1}.\end{eqnarray} \tag{ 18 }$$

The Alfvén surface obtained in all our simulations is low, indicating a large fractional mass loss per unit mass of accretion (${R}_{A}/{R}_{0}\sim 2$). This is related to our adopted low level of magnetization (${\beta }_{0}={10}^{4}$), and in this case, the winds are largely driven by the magnetic pressure gradient (Bai et al. 2016).

We can also see from Figure 2 that in simulations with the Hall effect, the change of magnetic field configuration leads to a segregation of magnetic fluxes in our simulation domain between the inner boundary region and the main disk body. This is largely a numerical artifact because we do not cover disk regions within the inner boundary but still have their magnetic flux contained in the polar region. At $t=1200{{\rm{\Omega }}}_{0}^{-1}$, there is a lack of magnetic flux between $R\sim 1\mbox{--}6$, where both accretion and outflow mass fluxes tend to diminish, and the Alfvén surface is no longer well defined. We will not discuss this region any further. Beyond this region, there are also local concentrations/rarefactions of magnetic flux as a result of intrinsic flux evolution, leading to stronger/weaker local poloidal field (better seen in the anti-aligned case). They make the Alfvén surface move further/retreat, which is consistent with theoretical expectations (Bai et al. 2016).

We start from the Hall-free simulation Fid0 as a reference. Together with Figure 4, we further show in Figure 5 with a more detailed analysis on the contribution from individual terms to magnetic flux transport at spherical radius r = 8. In the left panel, we show the θ profiles of radial and vertical components (in cylindrical coordinates) of ambipolar drift and flow velocities. In the right panel, we show the profile of

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{B}}{{v}_{K}}\equiv \displaystyle \frac{{{ \mathcal E }}_{\phi }}{{v}_{K}{B}_{z}},\end{eqnarray} \tag{ 19 }$$

which is dimensionless and can be considered as the effective velocity of flux transport v_B normalized by the Keplerian velocity. Overall, upon reaching a quasi-steady state, magnetic flux is transported outward (${v}_{B}\gt 0$) at all heights at approximately the same rate. In other words, the system adjusts/relaxes itself (mainly in its magnetic field configuration) in such a way that the rate of transport at all heights converges toward a constant value. We now analyze the contribution from individual physical effects.

Zoom In Zoom Out Reset image size

6.1.1. Contribution from AD

6.1.2. Contribution from Fluid Advection

6.1.3. Overall Rate of Flux Transport

In a quasi-steady state, the rate of flux transport v_B slowly evolves with time and is weakly dependent on disk radius (see Appendix B). For reference, we quote the value of ${v}_{B}\approx 4\times {10}^{-3}\,{v}_{K}$ from our Fid0 run, which corresponds to the rate measured at r = 8 and t = 1200 based on Figure 5 averaged within $\delta =\pm 0.3$.

We note that this rate is consistent with a simple order-of-magnitude estimate discussed in Section 2:

$$\begin{eqnarray}&&{v}_{B}\sim {\eta }_{A}/H\approx \displaystyle \frac{2}{\mathrm{Am}}\displaystyle \frac{1}{\beta }\displaystyle \frac{{H}_{\mathrm{mid}}}{R}\,{v}_{K},\end{eqnarray} \tag{ 20 }$$

where the plasma β is the ratio of gas to total magnetic pressure. We have $\mathrm{Am}=0.5$, ${H}_{\mathrm{mid}}/R=0.1$ in run Fid0. The plasma β depends on field strength. Note that while the initial poloidal field has ${\beta }_{0}={10}^{4}$, in a quasi-steady state, the toroidal field dominates, and is of the order ∼10 times stronger (see the fourth column of Figure 3) near the midplane. It gives $\beta \sim {10}^{2}$, which together yields ${v}_{B}\sim 4\times {10}^{-3}\,{v}_{K}$.

We now consider the Hall simulation Fid+ with an aligned poloidal field, and in Figure 6, we show the contributions to magnetic flux transport from individual terms at r = 8 and t = 1200, with the Hall-drift term (Equation (6)) included. Again, in this quasi-steady state, magnetic flux is transported outward at approximately the same rate at all heights. The roles of fluid advection and AD are similar to those discussed in the Hall-free case. Below, we focus on the Hall effect–mediated flux transport.

Zoom In Zoom Out Reset image size

6.2.1. Hall Effect–Mediated Flux Transport

6.2.2. Overall Direction and Rate of Flux Transport

From Figure 6, we quote the rate of outward flux transport in run Fid+ to be ${v}_{B}\approx 4\times {10}^{-3}\,{v}_{K}$ measured at r = 8 and t = 1200 averaged within $\delta =\pm 0.3$. This is about the same as the Hall-free case, despite the dramatic influence of the Hall effect. We now ask, why is magnetic flux eventually transported outward, and at a rate that is comparable to the Hall-free case?

The direction of flux transport in the midplane region is determined by the competition between the Hall-effect-driven inward transport, and outward transport due to AD or other dissipative processes (e.g., resistivity). Based on our earlier discussion, we expect the rate of outward flux transport to increase as the radial field becomes more stretched.⁷ In this sense, by adjusting poloidal field configurations, flux transport in both directions with a wide range of rates can be accommodated. Therefore, we expect that the global rate of magnetic flux transport is not determined in the midplane region.

From Figure 6, we see that at vertical height $\delta \sim 0.1\mbox{--}0.3$ where the toroidal field gradient reverses, both the Hall effect and AD lead to outward flux transport yet no other mechanism can provide any significant inward transport. We thus conclude that it is this region that determines the overall rate of flux transport, while other regions (midplane and the wind zone) adjust their field configuration to achieve the same rate of transport in response.

To an order of magnitude, we may estimate the rate of outward transport owing to the Hall effect using Equation (10), but applied to the intermediate layer between the midplane and the wind zone ($\delta \sim \pm 0.2$ in our case):

$$\begin{eqnarray}&&{v}_{B}\sim {v}_{A}\displaystyle \frac{{l}_{H}}{H}=\sqrt{\displaystyle \frac{2}{\beta }}\displaystyle \frac{{l}_{H}}{{H}_{\mathrm{mid}}}\displaystyle \frac{{H}_{\mathrm{mid}}}{R}\,{v}_{K}.\end{eqnarray} \tag{ 21 }$$

At r = 8 and $\delta =\pm 0.2$, we can infer from Figure 3 that $\beta \sim 5$ and ${l}_{H}\sim 0.3H$. We then obtain ${v}_{B}\sim 2\times {10}^{-3}\,{v}_{K}$, which reasonably approximates the measured rate of flux transport.

Finally, we discuss the Hall simulation Fid− with anti-aligned poloidal field, and show in Figure 7 the contributions to magnetic flux transport from individual terms at r = 8 and t = 1200. With an anti-aligned poloidal field, we find that the mechanism of magnetic flux transport, as seen from the decomposition shown in the figures, is qualitatively different from the Fid0 and Fid+ cases.

Zoom In Zoom Out Reset image size

6.3.1. Contribution from Individual Terms

6.3.2. Overall Direction and Rate of Flux Transport

The conventional advection–diffusion framework has been developed further over the past decade, aiming to address the issue of rapid loss of magnetic flux in thin accretion disks expected from Lubow et al. (1994). The majority of the works have been constructed under the assumption that the accretion disk is turbulent as a result of the MRI, with the effect of the MRI turbulence represented by an effective viscosity and resistivity. Their ratio, called the magnetic Prandtl number, has been measured to be of order unity in full MRI turbulence (Fromang & Stone 2009; Guan & Gammie 2009; Lesur & Longaretti 2009), and it leads to the conventional wisdom that magnetic flux is transported outward for thin accretion disks. More recently, disk vertical structure and boundary conditions have been recognized as important elements in the flux transport process. Bisnovatyi-Kogan & Lovelace (2007), Rothstein & Lovelace (2008), and Bisnovatyi-Kogan & Lovelace (2012) postulate that the surface of the accretion disks can be largely laminar because it is magnetically dominated (e.g., Miller & Stone 2000) and the MRI is suppressed. They suggest that inward advection of magnetic flux can be easily achieved in the nonturbulent, conducting surface layer. Guilet & Ogilvie (2012, 2013) presented a family of asymptotic solutions in the thin-disk limit, assuming that disk magnetic field is largely vertical. They find that magnetic flux can be advected inward much faster than the mass flux due to the large radial velocities at the low-density disk surface, and this effect can substantially alleviate the issue of magnetic flux loss by outward diffusion. With prescriptions of disk resistivity/viscosity, the theories have also been applied to PPDs (Guilet & Ogilvie 2014; Okuzumi et al. 2014; Takeuchi & Okuzumi 2014), and reasonable steady-state distributions of magnetic flux have been obtained.

These updated theories may be better applicable to fully MRI turbulent disks (e.g., black hole accretion disks), although they remain to be tested against simulations. However, when applied to PPDs, these theories suffer from several major issues. First, PPDs are largely laminar with the MRI suppressed or significantly damped (Bai & Stone 2013b; Bai 2013; Simon et al. 2013a; Gressel et al. 2015). Second, while physical resistivity is present in PPDs, as considered in Guilet & Ogilvie (2013) and Okuzumi et al. (2014), it dominates only in very limited regions in PPDs. The Hall effect and AD, which govern almost the entire range of radii in PPDs, are missing. Third, in the presence of a poloidal field threading the disk, MHD disk wind launching appears to be inevitable (Suzuki & Inutsuka 2009; Suzuki et al. 2010; Bai & Stone 2013a, 2013b; Fromang et al. 2013; Lesur et al. 2013), which has been ignored in most of the conventional studies (with the exception of the Guilet & Ogilvie series, which partially incorporated its effect). Besides the wind-driven accretion process, the wind provides necessary surface boundary conditions that are not easily incorporated into the existing framework.

Our work advances our understanding of magnetic flux transport in PPDs over several major aspects. First, we have shown that the Hall effect and AD play distinct roles in magnetic flux transport that are dramatically different from resistivity. Because of the anisotropic nature of the Hall effect and AD, magnetic flux transport depends on the gradient of magnetic fields in all directions, with the most sensitive being the vertical gradient of the toroidal field. It implies that the process of magnetic flux transport is strongly coupled with the gas dynamics of the disk itself, leading to substantial complications. Second, our simulations have self-consistently incorporated the launching and propagation of MHD disk winds. On one hand, it provides realistic "boundary conditions" at disk surfaces, and on the other hand, it allows the wind dynamics to adjust itself to match the rate of flux transport demanded from the disk.

Our results are also in line with some of the updated conventional theories in that flux transport is mediated by different mechanisms at different heights in the disk. Therefore, a complete theory of magnetic flux transport in PPDs must properly take into account the disk vertical structure, and use realistic vertical profiles of magnetic diffusivities.

Magnetic flux is not only the controlling factor for PPD evolution, it also plays a fundamental role controlling the process of star and disk formation. Molecular clouds are known to be relatively strongly magnetized (Crutcher 2012), and the star and disk formation processes naturally inherit some magnetic flux from the parent cloud. The role of AD has been discussed extensively in the literature in the context of the "magnetic flux problem" of star formation during core collapse (see McKee & Ostriker 2007 and the references therein for a thorough review), where substantial magnetic flux must be lost. The Hall effect has generally been considered not to be important during core collapse, as long as there is no significant rotation to develop strong toroidal magnetic field (e.g., Kunz & Mouschovias 2010).

Formation of rotationally supported PPDs, on the other hand, unavoidably involves the winding up of poloidal fields into toroidal fields, leading to significant magnetic braking. Further loss of magnetic flux is necessary to enable disk formation (Mellon & Li 2008, and see Li et al. 2014 and the references therein for a review). In this context, all non-ideal MHD effects prove to be important (Krasnopolsky et al. 2011; Li et al. 2011; Tomida et al. 2013). In particular, it has recently been found that AD enables significant magnetic flux loss (Tomida et al. 2015), and further inclusion of the Hall effect leads to a bimodality on the initial disk size depending on the polarity of the background magnetic field with respect to the initial angular momentum vector (Tsukamoto et al. 2015; Wurster et al. 2016).

While we have focused on magnetic flux transport in PPDs, the same physics is applicable in the disk formation process. Tsukamoto et al. (2015) and Wurster et al. (2016) have found that a larger disk is formed when field polarity is anti-aligned, whereas aligned field polarity leads to smaller initial disk size, although the underlying physical reasons were not addressed. The fact that both polarities lead to disk formation agrees with our conclusion that magnetic flux is systematically transported outward. Our finding that the anti-aligned case loses flux faster than the aligned case implies that less magnetic flux would be preserved in the former case, and hence leads to weaker magnetic braking. Therefore, we expect larger/smaller disks to be formed in the anti-aligned/aligned cases, offering an explanation for the findings of Tsukamoto et al. (2015) and Wurster et al. (2016).

Our results serve as a first step toward a better understanding of magnetic flux transport in PPDs, which largely controls global disk evolution. We discuss below how our results should be interpreted for this purpose, as well as cautions that must be exercised.

We first note that the rate of transport obtained in this work is rather high. For ${v}_{B}/{v}_{K}=5\times {10}^{-3}$, the magnetic flux depletion timescale would be on the order of only ∼30 local orbits, which amounts to only ∼10³ years at the distance of 10 au! However, owing to various caveats to be summarized in the next subsection, especially the artificial prescriptions of non-ideal MHD diffusivities, the rate of flux transport measured in this work is unlikely to be realistic, and readers should not take the values too seriously. As already emphasized, we aim to clarify the physics before incorporating more realistic prescriptions, which are left for future works.

The trend that the rate of outward flux transport increases with increasing net field strength suggests that the strongly magnetized phase of PPDs, if present, is short-lived because of relatively rapid loss of magnetic flux. The bulk of the disk lifetime is likely associated with relatively weak disk magnetization. Because stronger magnetic flux leads to a higher accretion rate, our results further imply that the rapid phase of disk evolution with high accretion rate is brief, and accretion rate decreases with time. As the loss of magnetic flux slows down over time, we expect the deceleration of the accretion rate to slow down as well. Although they are still too premature to be incorporated into global disk evolution models, our results also suggest that the assumption that magnetic flux is conserved in recent global disk evolution models is unlikely to be valid, while models with decreasing magnetic flux are more appropriate (see, e.g., Armitage et al. 2013; Bai 2016).

At this point, the speculations above are qualitative. With the current tools available, it has become feasible to conduct realistic simulations of PPDs that incorporate more realistic prescriptions, and we expect the physics learned from this work to be greatly beneficial for future explorations. While we have ignored Ohmic resistivity in this study, making it more applicable to the outer regions of PPDs ($\gtrsim 10$ au), it is the inner disk ($\lesssim 15$ au) that is expected to be almost fully laminar (Bai 2013, 2014). A further complication in the inner disk is that the vertical structure can become asymmetric (e.g., Bai & Stone 2013b; Gressel et al. 2015). With Ohmic resistivity dominating the midplane region and suppressing the current, the strong current layer typically lies at a few $\sim {H}_{\mathrm{mid}}$ offset from the midplane. It may be present on only one side of the midplane, leading to an asymmetric current distribution. When the Hall effect is turned on, it was found in local shearing-box simulations that maintaining a physical wind geometry with poloidal field lines bending away from the star can hardly be achieved (Bai 2014). Based on what we have found, such asymmetric current distribution would lead to asymmetric magnetic flux transport due to the Hall drift. It is unclear whether a quasi-steady state is possible in such an asymmetric configuration, which is an intriguing question for future investigations.

As a first study, our simulations are subject to several caveats and limitations, some of which were already mentioned earlier in the text and are explored further in Appendix A. We briefly summarize them below.

A major limitation of the present work is the assumption of axisymmetry. While it is likely a valid approximation in the inner disk ($r\lesssim 10\mbox{--}15\,\mathrm{au}$; Bai 2013, 2014), extending it to the outer disk can be problematic because the MRI can develop in both midplane (weak) and surface (strong), as has been previously studied in local shearing-box simulations (Perez-Becker & Chiang 2011; Simon et al. 2013a, 2013b; Bai 2015). Additional contributions from turbulence introduce further complications and uncertainties that need to be addressed using full 3D simulations.

Another limitation arises from the use of the HLL solver. We have shown in Appendix A that the rate of flux transport converges with resolution, implying that the system is able to self-adjust to compensate for the excess numerical dissipation from the HLL solver. However, it also implies that in reality, the strong current layer would be much thinner, which raises concerns on its stability. Corrugation of the strong current layer has already been observed in some of our simulations, which may eventually destroy the strong current layer. A similar phenomenon has also been observed in 3D shearing-box simulations of Bai (2015), where the midplane strong current layer corrugates but without destroying itself. In general, properly capturing this dynamical behavior would require full 3D simulations, as well as using less diffusive solvers (at limited resolution). Therefore, further improvements on the Hall MHD algorithm would be highly desirable.

For clarity, we have chosen simplified disk models (including thermodynamics) and prescribed the non-ideal MHD diffusion coefficients. The prescriptions are motivated from yet do not fully reflect realistic disk conditions. In particular, the non-ideal MHD diffusivities are mainly determined by the disk ionization level, and the vertical extent where they dominate largely depends on the penetration depth of external far-UV radiation. None of these processes are included. For instance, around 5 au, the disk would be thinner with ${H}_{\mathrm{mid}}/r\sim 0.05$ instead of 0.1, and non-ideal MHD dominates up to $\sim 4{H}_{\mathrm{mid}}$ for typical far-UV penetration depth, instead of $\sim 2.5{H}_{\mathrm{mid}}$. Therefore, we caution on the interpretation of the measured rates of flux transport. Investigations are underway to use more realistic prescriptions that incorporate ionization chemistry, far-UV irradiation, flared disk geometry, etc. (see also initial results from Béthune et al. 2016).