Local semi-analytic models of magnetic flux transport in protoplanetary discs

TL;DR

This paper develops semi-analytic models to study magnetic flux transport in protoplanetary discs, revealing how non-ideal MHD effects influence flux movement and the role of Hall drift in different regimes.

Contribution

It introduces two one-dimensional semi-analytic models for flux transport in discs, analyzing the impact of various non-ideal MHD effects and disc magnetizations.

Findings

Ohmic and ambipolar diffusivities cause outward flux transport.

Wind outflows promote inward flux transport, enhanced by Hall drift.

Hall effect can reduce or drive flux transport depending on conditions.

Abstract

The evolution of a large-scale poloidal magnetic field in an accretion disc is an important problem because it determines the launching of winds and the feasibility of the magnetorotational instability to generate turbulence or channel flows. Recent studies, both semi-analytical calculations and numerical simulations, have highlighted the crucial role non-ideal MHD effects (Ohmic resistivity, Hall drift and ambipolar diffusion), relevant in the protoplanetary disc context, might play in magnetic flux evolution in the disc. We investigate the flux transport in discs through the use of two one-dimensional semi-analytic models in the vertical direction, exploring regimes where different physical source terms and effects dominate. The governing equations for both models are derived by performing an asymptotic expansion in the limit of a thin disc, with the different regimes isolated through setting the relative order of the leading terms between variables. Flux transport rates and vertical structure profiles are calculated for a range of diffusivities and disc magnetisations. We found that Ohmic and ambipolar diffusivities drive radially outward flux transport with an outwardly inclined field. A wind outflow drives inward flux transport, which is significantly enhanced in the presence of Hall drift in the positive polarity case, $\eta_H (\boldsymbol{B}_z \cdot \boldsymbol{\Omega}) > 0$, an effect which has only been briefly noted before. Coupled only with outward inclination, the Hall effect reduces the flux transport given by a background Ohmic and/or ambipolar diffusivity, but drives no flux transport when it is the only non-ideal effect present.

Figures

Equations

$${\partial }_t\rho +\nabla \cdot (\rho \mathbf{u})=0,$$

$$\rho \left({\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\rho \nabla \Phi -\nabla p+\frac{1}{{\mu }_0}(\nabla \times \mathbf{B})\times \mathbf{B}+\nabla \cdot \mathrm{T},$$

$$\mathrm{T}=\rho \nu \left[\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\mathrm{T}}-\frac{2}{3}(\nabla \cdot \mathbf{u})\mathrm{I}\right],$$

$$p=c_s^2\rho .$$

$${\partial }_t\mathbf{B}=$$

$$\displaystyle\qquad\qquad+\eta_{A}\left[(\nabla\times\boldsymbol{B})\times\boldsymbol{b}\right]\times\boldsymbol{b}\bigg{)},$$

$$-r{\Omega }^2=-{\partial }_r{\Phi }_{\mathrm{m}},$$

$${\Phi }_{\mathrm{m}}(r)=-\frac{GM}{r}$$

$$\mathbf{v}=\mathbf{u}-r\Omega {\mathbf{e}}_{\varphi },$$

$$\xi =\frac{r}{ϵ},\zeta =\frac{z}{ϵ},t_1=ϵt,t_2={ϵ}^{2}t,$$

$${\partial }_r\mapsto {ϵ}^{-1}{\partial }_{\xi }+{\partial }_r,{\partial }_z\mapsto {ϵ}^{-1}{\partial }_{\zeta },$$

$${\partial }_t\mapsto {\partial }_t+ϵ{\partial }_{t_1}+{ϵ}^2{\partial }_{t_2}.$$

$$\rho ={\rho }_0+ϵ{\rho }_1+\cdots ,$$

$$p={ϵ}^2\left[p_0+ϵp_1+\cdots \right],$$

$$\mathbf{v}=ϵ\left[{\mathbf{v}}_0+ϵ{\mathbf{v}}_1+\cdots \right],$$

$$\mathbf{B}=ϵ\left[{\mathbf{B}}_0+ϵ{\mathbf{B}}_1+\cdots \right],$$

$${\eta }_O={ϵ}^2[{\eta }_{O0}+ϵ{\eta }_{O1}+\cdots ],$$

$${\eta }_H={ϵ}^2[{\eta }_{H0}+ϵ{\eta }_{H1}+\cdots ],$$

$${\eta }_A={ϵ}^2[{\eta }_{A0}+ϵ{\eta }_{A1}+\cdots ].$$

$$\Phi ={\Phi }_{\mathrm{m}}(r)+{ϵ}^2\Psi (r)\frac{1}{2}{\zeta }^2+O({ϵ}^4),$$

$$\Psi ={\Omega }^2$$

$${\partial }_{\tau }{\Sigma }_0=$$

$$\displaystyle\qquad\qquad\left[\partial_{r}\left(r^{2}\left\{\langle\rho_{0}v_{\phi 0}v_{r0}\rangle-\frac{1}{\mu_{0}}\langle B_{r0}B_{\phi 0}\rangle\right\}\right)\right]\bigg{\}},$$

$$\displaystyle\partial_{t_{1}}\langle B_{z0}\rangle+\frac{1}{r}\partial_{r}\bigg{(}r\bigg{[}\langle v_{r0}B_{z0}\rangle-\langle v_{z0}B_{r0}\rangle+\left\langle\eta_{O0}\partial_{\zeta}B_{r0}\right\rangle$$

$$\displaystyle\qquad+\bigg{\langle}\eta_{H0}\frac{1}{B_{0}}\bigg{\{}B_{r0}\partial_{\xi}B_{\phi 0}+B_{z0}\partial_{\zeta}B_{\phi 0}\bigg{\}}\bigg{\rangle}$$

$$\displaystyle\qquad+\bigg{\langle}\eta_{A0}\frac{1}{B_{0}^{2}}\bigg{\{}B_{r0}B_{\phi 0}\partial_{\zeta}B_{\phi 0}+\left(B_{r0}^{2}+B_{z0}^{2}\right)\left(\partial_{\zeta}B_{r0}-\partial_{\xi}B_{z0}\right)$$

$$\displaystyle\qquad\qquad+B_{\phi 0}B_{z0}\left(\partial_{\xi}B_{\phi 0}\right)\bigg{\}}\bigg{\rangle}\ \bigg{]}\ \bigg{)}=0.$$

$$⟨{\eta }_{H0}\left(\frac{B_{z0}}{B_0}\right){\partial }_{\zeta }{B}_{\varphi 0}⟩,$$

$$⟨{\eta }_{A0}\left(\frac{B_{r0}^2+B_{z0}^2}{B_0^2}\right){\partial }_{\zeta }{B}_{r0}⟩,$$

$$⟨{\eta }_{A0}\left(\frac{B_{r0}{B}_{\varphi 0}}{B_0^2}\right){\partial }_{\zeta }{B}_{\varphi 0}⟩,$$

$$c_s=1⟹v_a/c_s=v_a,$$

$$v_a=\frac{B}{\sqrt{{\mu }_0\rho }}=\sqrt{\frac{B_x^2+1}{\rho }},v_a(0)=\sqrt{\frac{1}{\rho (0)}}.$$

$$\mathbf{E}+\mathbf{v}\times \mathbf{B}={\eta }_O\mathbf{J}+{\eta }_H\mathbf{J}\times \mathbf{b}-{\eta }_A(\mathbf{J}\times \mathbf{b})\times \mathbf{b},$$

$$-2\rho \Omega v_y=J_y{B}_z,$$

$$\frac{1}{2}\rho \Omega v_x=-J_x{B}_z,$$

$$0=-\rho {\Omega }^{2}z-\frac{\mathrm{d}p}{\mathrm{d}z}+J_x{B}_y-J_y{B}_x.$$

Full text

††pubyear: 2019††pagerange: Local semi-analytic models of magnetic flux transport in protoplanetary discs–References

1 Introduction

Magnetic processes play a significant role in the evolution of accretion discs. Discs are often ionised to some degree (Alexander et al., 2004), making them susceptible to magnetic effects in the magnetohydrodynamic (MHD) approximation. Magnetic instabilities, such as the magnetorotational instability (MRI), generate turbulence which lead to inward mass accretion by acting as an effective viscosity transporting angular momentum outwards (Balbus & Hawley, 1991). The presence of a large-scale magnetic field can launch magnetocentrifugal winds (MCWs) which remove angular momentum vertically from the disc surface (Blandford & Payne, 1982), also resulting in mass accretion. Besides these two aforementioned mechanisms which have been subject to extensive studies, magnetic fields have also been theorised to affect disc evolution through other means such as magnetic braking (Krasnopolsky & Königl, 2002), fields linking different radii of the discs (Wardle, 2007), and many more.

A key parameter governing the significance and nature of these processes is the strength and geometry of the large-scale magnetic field threading the disc. The launching of a MCW requires the presence of a large-scale poloidal field above a threshold strength (Wardle & Koenigl, 1993), while the MRI operates most effectively in a ‘weak field’ regime (Fleming et al., 2000), and is quenched when the flux exceeds a certain value (Balbus & Hawley, 1998).

In protoplanetary discs, further complexities arise from the presence of non-ideal MHD effects (Ohmic resistivity, ambipolar diffusion and the Hall drift) from their weak ionisation level. These provide additional parameters affecting the launch criterion of the MCW (Wardle & Koenigl, 1993), and whether or not the MRI is stabilised (Wardle, 1999; Balbus & Terquem, 2001; Kunz & Balbus, 2004; Salmeron & Wardle, 2005). There is growing evidence that non-ideal MHD effects, when coupled with a large-scale magnetic field, can induce the formation of self-organised structures such as axisymmetric rings (Béthune et al., 2016, 2017; Suriano et al., 2018). Such features have been observed in many protoplanetary systems (ALMA Partnership et al., 2015) and are important for the study of planet formation.

So far, two possible origins for this large scale magnetic field have been proposed: it is either created locally by a dynamo process (Brandenburg et al., 1995; Hawley et al., 1996), or it arises from the concentration of the flux that is intrinsic to the accreting gas (Spruit et al., 1995; Lubow & Spruit, 1995). It is yet unclear whether the former can generate a significant magnetic field coherent over a scale comparable to the radius (Spruit, 2010).

The latter case, where magnetic flux of the accreting gas is concentrated as it accumulates in the inner regions of the disc, has been subject to increasing study. Semi-analytic work such as carried out by Lovelace et al. (2009) and Guilet & Ogilvie (2012, 2013) studied the competing effects of inward accretion due to advection with an accreting flow, and outward diffusion due to a turbulent resistivity, in what is known as the ‘advection-diffusion’ paradigm. The disc is assumed to be turbulent due to processes like the MRI. This, as mentioned before, leads to an effective viscosity, $\nu$, which drives accretion of the material towards the star, also advecting magnetic flux radially inwards. However, the same turbulence also mixes up fluid parcels and their magnetic flux, driving reconnection processes on small scales. This leads to an effective resistivity, $\eta$, which allows the magnetic field to diffuse, and in the presence of an outwardly bending global poloidal field, results in an outward transport of flux. The turbulent viscosity and diffusivity are linked via the effective magnetic Prandtl number $\mathcal{P}\equiv\nu/\eta$. For MHD turbulence, $\mathcal{P}$ is expected to be of order unity (Fromang & Stone, 2009; Guan & Gammie, 2009; Lesur & Longaretti, 2009).

The first studies by Lubow et al. (1994) and Heyvaerts et al. (1996) concluded that outward diffusion would significantly outweigh any accretion for thin discs with significant field bending. These studies, however, neglected the vertical structure of the disc, which, as subsequently pointed out (Ogilvie & Livio, 2001; Rothstein & Lovelace, 2008; Guilet & Ogilvie, 2012), has a significant impact on the averaging procedure used in the calculations. These studies examined the impact of disc structure using quasi-steady state radially local models, and found that for weak magnetic fields, the contribution of the accretion flow to the flux transport velocity is much larger than the mass transport velocity. This is due to larger radial velocities away from the mid-plane, which barely affect mass flow but do affect flux transport. Further work by Guilet & Ogilvie (2014); Okuzumi et al. (2014) and Takeuchi & Okuzumi (2014) applied the results from these local models to the global evolution of large-scale magnetic fields, and verified that an equilibrium is reached and that the maximum attainable field strength is in agreement with the steady-state accretion rate observed in actual systems. Hence a possible solution to the too-efficient-diffusion problem of an inclined magnetic field is provided.

This flux transport paradigm, however, is more relevant to well ionised discs such as around black holes and compact stars, where turbulence can readily occur. Recent studies and global simulations have shown that PPDs are likely to be laminar in most parts, and therefore should not be modeled as turbulent (Béthune et al., 2017; Bai & Stone, 2017; Suriano et al., 2017, 2018). Current advection-diffusion models also do not account for the non-ideal effects of Hall drift and ambipolar diffusion, which are likely to be significant throughout much of the PPD. Recently numerical simulations have begun examining the flux transport problem globally in the non-ideal MHD protoplanetary disc context (Bai & Stone, 2017) and in the ideal MHD weak field regime (Zhu & Stone, 2018). However, they still suffer from current computational limitations of evolving a disc for long enough to determine its long term global flux evolution (G. Lesur, priv. comm.). A new model addressing these concerns and including all three non-ideal MHD contributions is required for modelling the flux transport problem for the PPD context, which is what we seek to address in this work.

In this paper, we investigate the flux transport problem in the PPD context through 1D radially local semi-analytic models similar to that of Guilet & Ogilvie (2012), but where non-ideal effects also have a significant role in the flux evolution (Fromang, 2013). However, as many of the same principles can be applied to accreting discs in general, it is hoped that the framework could be extended to the study of stellar mass black hole accretion and AGN discs as well.

The structure of this paper is as follows. In Section 2, we construct the multiscale aysmptotic approach which forms the basis for the two semi-analytic models used in this work. Section 3 discusses the shearing box equations vertical structure model, where we estimate flux transport rates due to an inclined $B$ field or a magnetic torque at the disc surface. Section 4 extends the Guilet & Ogilvie (2012) model (hereafter GO model) to include all three non-ideal effects. It also examines the flux transport driven by large scale radial gradients, in addition to that by a surface inclined field or magnetic torque. The implications and limitations of our results are discussed in Section 5, along with directions for future study. Finally, we summarise and conclude in Section 6.

2 A multiscale asymptotic approach

Guilet & Ogilvie (2012) were distinctive in using a multiscale asymptotic approach in deriving their 1D radially local vertical structure equations, allowing radial gradients of various quantities to be placed on the same footing as inclination and outflow in driving flux transport. This approach is often employed in the study of warped discs where processes vary at different spatial and time scales (Ogilvie, 1999). We follow the same multiscale asymptotic approach, but here derive a more general formalism. All our subsequent work is based on the equations derived using this formalism, and here we use the equations derived to also draw some general observations on the interplay between various source terms and the angular momentum transport and magnetic flux evolution.

2.1 Governing equations

The continuity equation is given by

[TABLE]

and the conservation of momentum by

[TABLE]

where $\rho$ is the density, ${\boldsymbol{u}}$ is the velocity, $\Phi$ is the gravitational potential, $p$ is the pressure, $\boldsymbol{B}$ is the magnetic field, and $\mathbf{T}$ is the viscous stress. The full viscous stress is given by

[TABLE]

where $\nu$ is the kinematic viscosity and $\boldsymbol{\mathrm{I}}$ is the unit tensor of second rank. We neglect self-gravity and do not consider a thermal energy equation. For thermodynamic closure, we assume an isothermal relation such that

[TABLE]

The full non-ideal MHD induction equation is given by (Balbus, 2011):

[TABLE]

where $\boldsymbol{b}=\boldsymbol{B}/\lvert\boldsymbol{B}\rvert$, and $\eta_{O},\eta_{H},\eta_{A}$ are the Ohmic, Hall and ambipolar diffusivities respectively. Ohmic and ambipolar diffusivities are always positive. The Hall coefficient, on the other hand, can have either sign, depending on whether the more massive charge carrier in the disc is positively ($\eta_{H}>0$) or negatively ($\eta_{H}<0$) charged. Ideal, Ohmic and/or ambipolar diffusive only discs are insensitive to the polarity of alignment of the large scale field with the rotation of the disc. The Hall effect, on the other hand, does depend on the polarity of the magnetic field. Reversing the sign of $(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})$ has the equivalent effect as reversing the sign of $\eta_{H}$ on the equations. In this paper, we define the positive (negative) polarity case as when $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>(<)~{}0$.

We consider the problem in cylindrical polar coordinates $(r,\phi,z)$, and assume an axisymmetric potential $\Phi(r,z)$ that is symmetric about the midplane $z=0$. The orbital angular velocity $\Omega(r)$ is defined via

[TABLE]

where

[TABLE]

is the midplane potential. We also define the residual velocity as

[TABLE]

which is the departure from Keplerian motion.

2.1.1 Asymptotic expansion

We perform an asymptotic expansion for the fluid variables, using an ordering scheme that is well understood for thin discs. The characteristic angular semi-thickness of the disc is given by the disc aspect ratio $H/r=\mathcal{O}(\epsilon)$, where $0<\epsilon\ll 1$ is a small dimensionless parameter. The sound speed to orbital velocity ratio in an isothermal disc is then of $\mathcal{O}(\epsilon)$.

For a magnetised disc, the ordering scheme of the magnetic field depends on the strength and orientation of the field (Ogilvie, 1997). We consider here a situation where the magnetic energy density is comparable to the thermal energy density, and set both the Alfvén speed and the adiabatic sound speed at $\mathcal{O}(\epsilon)$ compared with the orbital velocity. For specific field configurations, such as vertical field dominance in the Guilet & Ogilvie model described in Section 4, we might set some of the asymptotic terms to zero for that particular situation. We neglect the effects of a viscous stress, setting the viscosity to zero, giving us a laminar disc, as recent PPD simulations have shown them to be.

We assume axisymmetry and so neglect all $\partial_{\phi}$ terms. The internal structure of the thin disc and its slow evolution in time are resolved through rescaled spatial and time coordinates

[TABLE]

where $\eta$ is a scaled angular variable, $t_{1}$ is the flux evolution timescale and $\tau$ is the accretion timescale. To incorporate the multiscale nature of our approach, we separate the radial lengthscale into small variations of order $H$ within a small radial region of the disc ($\xi\sim\mathcal{O}(1)$) and global radial variations of order $r$ ($r\sim\mathcal{O}(1)$). Timescales are also separated into a fast dynamical timescale of order $\Omega^{-1}$ ($t\sim\mathcal{O}(1)$), an intermediate magnetic flux transport timescale ($t_{1}\sim\mathcal{O}(1)$), and a slow accretion timescale ($t_{2}\sim\mathcal{O}(1)$). The coordinate transformations are hence given by

[TABLE]

We then use the following expansion of fluid variables (each of the terms in the expansion is a function of $r,\xi,\zeta,t,t_{1},t_{2}$):

[TABLE]

We can see that $p/\rho\sim\mathcal{O}(\epsilon^{2})$, giving us $c_{s}\sim\mathcal{O}(\epsilon)$ as we would expect. We consider departures from Keplerian motion comparable to the sound speed, and a similar order for the Alfvén speed, hence set both ${\boldsymbol{v}}$ and $\boldsymbol{B}$ to be of $\mathcal{O}(\epsilon)$. Variations across the disc height due to diffusive effects are assumed to be at the dynamical timescale, hence $\eta_{O,H,A}\sim\mathcal{O}(\epsilon^{2}r^{2}/t)\sim\mathcal{O}(\epsilon^{2})$. Finally, within the thin disc, the gravitational potential has the Taylor expansion

[TABLE]

where

[TABLE]

for the point-mass potential we have chosen.

At leading order ($\mathcal{O}(1)$ for mass conservation and $\mathcal{O}(\epsilon)$ for the momentum equation), we recover the standard shearing box equations. The next order equations can be thought of as linear equations with source terms given by the solution of the leading order equations. The expansion of these equations at the first two leading orders can be found in Appendix A (see online supplementary materials).

2.1.2 Angular momentum and magnetic flux transport

Of particular interest to us are the leading terms for angular momentum and magnetic flux transport rates. Assuming that disc variables are periodic in $\xi$, we simplify the equations by averaging over the horizontal directions and over the fast dynamical timescale $t$. The leading order angular momentum transport equation is then recovered by combining the $\mathcal{O}(\epsilon^{2})$ $\phi$-component of the momentum equation with the $\mathcal{O}(\epsilon)$ mass equation, and integrating over the vertical extent with no outflow boundary condition to give

[TABLE]

where $\langle\rangle$ denote the spatial and temporal averagings and the vertical integration mentioned, and $\Sigma_{0}$ is the first order disc surface density. We recover the well known ‘diffusion equation’ for an accretion disc, that mass transport is facilitated through ‘viscous’ stresses arising from a Reynolds component associated with internal motions ($\langle\rho_{0}v_{\phi 0}v_{r0}\rangle$) and a Maxwell component associated with the horizontal magnetic field components ($-\langle B_{r0}B_{\phi 0}\rangle$).

The vertical flux evolution is calculated from the $z$ component of the induction equation and using the condition $\nabla\cdot\boldsymbol{B}=0$:

[TABLE]

The terms within the square brackets give us the vertical magnetic flux radial transport rate. The first term describes advection with the radial flow, while the second term shows us the effect from a vertical outflow from the disc. Standard accretion disc models assume the large scale magnetic field to have an open field geometry with the poloidal component bending outward from the star (Ogilvie, 1997) (see Figure 1). In this configuration, a vertical outflow would therefore drive radially inward accretion of the flux ($-\langle v_{z0}B_{r0}\rangle<0$). For the Ohmic term, we expect $\langle\partial_{\zeta}B_{r0}\rangle>0$ in the standard disc magnetic configuration. This means an overall effect of diffusing the magnetic field radially outwards, in agreement with other studies of the effect of Ohmic resistivity in the presence of an inclined field (Guilet & Ogilvie, 2012).

The Hall and ambipolar terms are more complicated. We further simplify the problem by assuming that only large scale radial gradients and vertical variations are present (in other words we set $\partial_{\xi}=0$). Then the Hall term,

[TABLE]

depends on the vertical gradient of the azimuthal field. Its direction is also dependent on the sign of the Hall coefficient and the vertical magnetic field, reversing if one of them changes. Its effect is strongest when there is a strong vertical field threading the disc. From simulations (Bai & Stone, 2017; Béthune et al., 2017), $\partial_{\zeta}B_{\phi 0}<0$ near the midplane and $>0$ at the disc surface, hence we would expect the Hall flux transport to switch direction as we go up from the disc midplane. For $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>(<)~{}0$, we expect this transport to be radially inward (outward) near the midplane. This is indeed what is found in recent Hall global simulations (Bai & Stone, 2017).

The ambipolar term on the other hand has both an ‘Ohm-like’ and a ‘Hall-like’ component. The ‘Ohm-like’ component,

[TABLE]

is outwardly directed assuming the standard configuration, and is strongest and equivalent in form to Ohmic diffusion when $\boldsymbol{B}$ is purely poloidal. The ‘Hall-like’ component,

[TABLE]

depends on the correlation between the horizontal components of the magnetic field. Normally, we expect the sign of $B_{\phi}$ to be opposite to that of $B_{r}$ because of Keplerian shear, and $B_{r0}B_{\phi 0}<0$. Hence the ‘Hall-like’ ambipolar term would act in the opposite (same) direction as the Hall term for $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>(<)~{}0$.

It is worth noting that Section 2 (and Figure 1) of Bai & Stone (2017) also provides an analytic explanation on why the Hall term drives magnetic flux transport when a surface toroidal field accompanied by a wind is present, with Equation 9 of Bai & Stone (2017) essentially the same as our Equation (20). Although our analysis is more rigorous, the essence of the Hall-mediated magnetic flux transport is the same.

3 A shearing box (radially local) 1D vertical structure model

In our first investigation, we developed a semi-analytic radially local 1D vertical model that includes all non-ideal MHD effects for the flux transport problem. Our first model is essentially a local shearing box, where we solve the leading order equations, and assume variations only in the vertical direction. This assumption is motivated by global non-ideal MHD simulations such as those of Bai & Stone (2017), where quasi-steady states were achieved and the 1D vertical disc profiles were measured. The goal of our approach is to gain insight into the interpretation of their results.

We begin by examining the case with no MCW present by restricting the range of field inclinations explored to below the wind launching condition, but then mimic the presence of a vertical outflow by setting a non-zero azimuthal field, hence magnetic torque, at the boundary. In our model, flux transport is driven by the interplay between the bending of the poloidal field and the diffusivities present.

3.1 Leading order equations

We use the thin disc approximation ($J_{z}=0$, see Equation (25)), and consider a Keplerian rotating frame. We want to find the quasi-steady state equilibrium profiles, so the disc is set to be steady on the dynamical timescale (i.e. $\partial_{t}=0$). As before, we assume an isothermal and inviscid flow. For now, there is also no vertical outflow, so we set $v_{z}=0$.

Without loss of generality, we use units such that $B_{z},\Omega,c_{s},\mu_{0}=1$, with $B_{z}$ set to constant in this model. The vertical coordinate $z$ is then in units of the isothermal scaleheight, $H=c_{s}/\Omega$, while diffusivities are given in units of $H^{2}\Omega$. The relation between the Alfvén speed and the sound speed (a measure of the field strength) is then defined by our choice of the midplane ($z=0$) density, $\rho(0)$:

[TABLE]

[TABLE]

$\rho(0)\gg 1$ therefore corresponds to a weak magnetic field and $\rho(0)\ll 1$ to a strong one. When $\rho(0)=1$, $c_{s}=v_{a}$ at the midplane. For an isothermal disc profile, $\rho(z)$ decreases as $z$ increases, so we expect magnetic effects to become stronger compared to hydrodynamic effects as we move away from the midplane.

The local equations are identical to the shearing box equations, and correspond to the leading order equations of the multiscale asymptotic approach. We consider only vertical variations, setting $\partial_{\xi}=0$. To incorporate non-ideal effects, we use the modified electric field in the rest frame of the multi-component fluid $({\boldsymbol{E}}+{\boldsymbol{v}}\times\boldsymbol{B})$, with ${\boldsymbol{E}}$ being the electric field in the shearing-box frame:

[TABLE]

where ${\boldsymbol{J}}=\nabla\times\boldsymbol{B}$ is the current density, ${\boldsymbol{v}}$ throughout this section is the fluid velocity (rather than the Keplerian-shear subtracted velocity in Section 2), and the other variables are as before. This is another way of writing the modified induction equation, since $\partial_{t}\boldsymbol{B}=-\nabla\times{\boldsymbol{E}}$. In our scheme here, the positive (negative) polarity configuration for the Hall effect is achieved when $\eta_{H}>(<)~{}0$.

We then have the following governing equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $x,y,z$ denoting the radial, azimuthal and vertical directions respectively, and we have removed the ordering subscripts. Equations (26) to (28) are the three components of the equation of motion, while the horizontal components of the induction equation are given by Equations (29) and (30). These terms on the RHS of Equations (26) and (27) indicate that the departure from Keplerian motion are driven by the Maxwell terms. We can also see that the density structure is determined by the effects of gravity, pressure balance, and magnetic compression, which are represented by the first, second and final pair of terms respectively of Equation (28).

3.2 Flux transport

The flux transport rate is given by the azimuthal electric field $E_{y}$, which is constrained to be a constant with varying disc height in our model. We can see that it measures the flux transport rate by examining the vertical component of the induction equation in cylindrical geometry:

[TABLE]

In strict steady state, the left hand side of Equation (31) is zero. But since $B_{z}$, and equivalently the poloidal magnetic flux function $\Psi(r,z)$ (see Sections 2 & 3.3 in Guilet & Ogilvie (2012) for a justification using that approach), can evolve on the long accretion time-scale $\tau_{a}\equiv r/\lvert v_{r}\rvert$ as the poloidal field drifts radially through the disc, we can define a radial flux transport rate to be given by

[TABLE]

(Ogilvie & Livio, 2001; Königl et al., 2010; Guilet & Ogilvie, 2012). This is found in our model as an eigenvalue of the solution. $E_{y}>(<)~{}0$ would hence imply a radially outward (inward) transport.

3.3 Boundary conditions and numerical method

We use the standard assumption that the disc is antisymmetric about the midplane in $B_{x}$ and $B_{y}$, and symmetric in $v_{x}$ and $v_{y}$. The value of $B_{x}(z\to\infty)$ sets the inclination of the poloidal field, which is one of the parameters we vary. In a real disc, the inclination would be determined by the global magnetic field geometry. For the case where there is no wind launch, we expect no external magnetic torque to be acting on the disc, and set $B_{y}=0$ at $z\to\infty$. In the case where we mimic the presence of a wind, we set $B_{y}$ to be non-zero at $z\to\infty$, since there is now an external magnetic torque acting on the disc removing angular momentum vertically. For an actual outflow, we would expect $B_{y}(\infty)<0$. The value of $\rho(z=0)$ is another parameter we vary and sets the field strength. In our units, the midplane density $\rho_{0}$ is equivalent to $\beta_{0}/2$, where $\beta=p_{\text{thermal}}/p_{\text{magnetic}}$ is the plasma beta. A large value of $\rho_{0}$ would therefore indicate a weak midplane magnetic field.

We recast Equations (26)$-$(30) into a fourth order ordinary differential system with variables $\rho$, $B_{x}$, $B_{y}$ and $E_{x}$ (see Appendix B in the online supplementary materials for details), where $E_{y}$ becomes an eigenvalue of the problem, and solve them using the shooting method. We integrate upwards from the midplane with values determined by the boundary conditions and guessed midplane values for $E_{x}$ and $E_{y}$. The ODE system is then integrated using fourth order Runge-Kutta up to a sufficient height to mimic the solution as $z\to\infty$ (we set $z_{\text{end}}=10H$ in our case, where $H$ is the isothermal scale height. Extending $z_{\text{end}}$ beyond this value is found to yield no significant variation in our results). There, boundary conditions at infinity are applied. This is followed by Newton-Raphson iterations to adjust the guessed solution until the variables converge to the desired upper boundary conditions.

3.3.1 Diffusivity profiles

We used two different diffusivity profiles in our calculations. The first profile, ‘Uniform’, assumes uniform $\eta_{O}$, $\eta_{H}$ and $\eta_{A}$ across all scale heights, and is useful for determining the general behaviour of the disc due to each diffusivity. The second scheme, ‘CstIon’, tries to capture PPD vertical diffusivity profiles more accurately by assuming only constant ionisation across the vertical extent, leading to $\eta_{O}=\mathrm{constant}$, $\eta_{H}\propto B/\rho$ and $\eta_{A}\propto B^{2}/\rho^{2}$ (Wardle, 2007; Balbus, 2011). We also impose a rapid transition to the ideal MHD regime for $z>2$ to model the effect of photoionisation of the disc surface from the star. In all the plots that follow, the labelled values of $\eta$ are true for the whole disc in the ‘Uniform’ case, whereas for the ‘CstIon’ case, they are the values of the midplane diffusivities.

3.4 Marginal stability analysis

According to Ogilvie (1998), steady state solutions where field lines bend several times as they pass through the disc indicate their instability to the MRI in the ideal MHD regime. Although it yet remains to be rigorously proven, we follow Ogilvie & Livio (2001) and Guilet & Ogilvie (2012) in assuming that this to be true when non-ideal effects are included as well. The multiple bending corresponds to the “channel mode” of the MRI for a vertical field. The lowest order “$n=1$” mode, which determines the overall marginal stability, has opposite symmetry in the disc to the fluid variables. There, $\Delta B_{x}$ and $\Delta B_{y}$, where $\Delta$ denotes a perturbation, would be symmetric about the midplane, while $\Delta u_{x}$ and $\Delta u_{y}$ would be antisymmetric. Hence we can compute the level of magnetisation (in our case given by $\rho_{0}=\beta_{0}/2$) that leads to a marginally stable state to the MRI by solving for a set of linearised equations of small perturbations on top of the nonlinear equations, with the appropriate symmetry conditions.

This was done for a range of diffusivity coefficients using the shooting method described in Section 3.3. Midplane values of $\Delta B_{y}$ and $\rho$ were guessed along with $E_{x}$ and $E_{y}$ as before ($\Delta E_{x}$ is calculated as an output, while $\Delta E_{y}$ is zero due to the symmetry of the perturbation). The midplane value of $\Delta B_{x}$ determines the amplitude of the marginal mode, and was arbitrarily set to $10$. Boundary conditions for the unperturbed variables were the same as before. For the perturbed variables, midplane values are derived from the symmetry conditions ($\Delta u_{x}=\Delta u_{y}=0$), while the upper boundaries are set to enforce the fixed inclination ($\Delta B_{x}=0$) and magnetic torque ($\Delta B_{y}=0$) conditions. Threshold field strengths are then computed, which mark the field strengths below which the disc would be unstable to the MRI, and where we would expect multiple bending in the vertical structure solutions.

3.5 Disc vertical structure profiles

We found disc vertical structure profiles to be divided into stable and unstable solutions (see Figure 2) when the field is aligned with the rotation. Weak field solutions with small non-ideal contributions show multiple bending of the poloidal field lines, indicating an unstable configuration due to the MRI or other instabilities (Ogilvie, 1998). Strong field solutions or those with large non-ideal contributions, on the other hand, share the same general shape with a single bend and are stabilised by the strong field and/or diffusive effects present. When the field is anti-aligned with the rotation and the Hall effect is present, almost all solutions have multiple bends, and the solver often failed to converge. This is interpreted as indicating that there are no stable solutions in the anti-aligned case with the Hall effect in our 1D equilibrium model.

Qualitatively, we find that the ‘Uniform’ and ‘CstIon’ diffusivity profiles give us similarly shaped vertical structures for the same set of midplane diffusivities. Differences lie in the magnitude of the horizontal $\boldsymbol{B}$ field, which may contrast by up to an order of magnitude (with those under ‘Uniform’ being larger than those under ‘CstIon’) due to the ‘CstIon’ profile allowing $\eta_{H}$ and $\eta_{A}$ to take significantly larger values away from the midplane as $\rho$ decreases. Another difference is that the ‘CstIon’ solutions have less smooth features in their vertical profiles at $z\sim 3-4$, when the diffusivity profile rapidly drops off to the ideal MHD regime, than their ‘Uniform’ counterparts. However, overall, their trends in the stability of solutions and dependence on diffusivities and field strength are the same. In all the figures that follow, we plot the solutions for the more realistic ‘CstIon’ diffusivity profile.

3.6 Threshold magnetisation for marginal stability and variation with diffusive effects

We conducted the marginal stability analysis described in Section 3.4 to find the transition field strength between stable and unstable configurations. Some representative plots of our results are shown in Figure 3. As diffusive effects tend towards zero, the threshold midplane $\rho_{0}$ for marginal stability tends towards the value of $0.4474$, giving us a corresponding $\beta_{0}\approx 0.89$. This agrees with the value found in other studies111In Paris & Ogilvie (2018), a value of $B_{z^{\prime}}=0.94$ is found, where $B_{z^{\prime}}=B_{z}/\sqrt{\mu_{0}\Sigma\Omega c_{s}}$ and $\Sigma=\int^{\infty}_{\infty}\rho\mathrm{d}z$. For an isothermal disc with small magnetic compression, $\rho(z)$ follows a Gaussian profile and $\rho_{0}/\Sigma=1/\sqrt{2\pi}$, giving us $\beta=p/p_{\text{mag}}\approx 0.894$. for the largest value of $B_{z}$ allowing the MRI to operate in the ideal MHD regime (Gammie & Balbus, 1994; Paris & Ogilvie, 2018).

When $\eta_{H}>0$ (postive polarity configuration), all three non-ideal effects help stabilise the disc - the threshold midplane $\rho_{0}$ increases (i.e. critical field strength decreases) with increasing strength of each non-ideal effect. On the other hand, for negative polarity (which in our case is when $\eta_{H}<0$), increasing the magnitude of the Hall diffusivity decreases the threshold $\rho_{0}$ (critical field strength increases), showing that the Hall term has the effect of destabilising the disc222The ”stability” here refers to the threshold magnetization for marginal stability, or in terms of local analysis, to the critical wavenumber for instability. However, the Hall term also has the effect of enhancing the maximum growth rate for $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>0$ (e.g., Figure 6 of Wardle & Salmeron (2012)), the phenomenon called the Hall-shear instability (Kunz, 2008). For this reason, some authors, e.g., Wardle & Salmeron (2012), describe the Hall effect for $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>0$ as ”destabilizing”, which is true for our disc model as well if the disc is already in the unstable regime, but irrelevant if our disc is in the regime stable to the MRI..

Multiple branches are observed when we varied the strength of ambipolar diffusivity under certain backgrounds. These were also observed as we varied the initial guessed value of $\rho_{0}$. We interpret these branches to correspond to the excitation of different modes of instability, with higher modes giving multiple bends in the disc. Marginal stability should therefore be given by the branch with the lowest $\rho_{0}$ values giving us the lowest order mode, where only one bend occurs.

3.7 Flux transport rate and variation with diffusive effects

The variation of $E_{y}$ with parameters was calculated by solving for a particular set of diffusivities, then iterating to either higher/ lower diffusivities using the previous solution. This procedure allows us to identify different branches (if any) of solutions, which occur when the disc is in an unstable configuration.

We studied most extensively the positive polarity case when $\eta_{H}>0$. When $\eta_{H}<0$ (negative polarity case), the solver failed to converge for the majority of cases. This may be indicative of the destabilising effect a negative polarity coupled with $\eta_{H}$ has, as found by Balbus & Terquem (2001).

3.7.1 Unstable solutions and bifurcations

In all cases, weak field solutions with low background non-ideal effects show multiple branches when the direction of iteration is varied, indicative of unstable configurations. From the marginal stability analysis, we expect the point of bifurcation at the largest value of the varied diffusivity to correspond to an equivalent point there, or at least occur at a lower diffusivity than the critical one calculated. Currently our study is not fully conclusive on this, as not all branches have been exhaustively identified using our method. However, we do find that all points of bifurcation occur at lower diffusivities than the critical points calculated. On the other hand, strong field solutions with large diffusivities appear to be stable, and solutions calculated by referencing in both directions agree.

3.7.2 Inclination driven flux transport

First, we explore the flux transport driven by inclination alone by setting the inclination through the value of $B_{x}(\infty)$. We eliminate the effect of outflow driven flux transport by setting $B_{y}(\infty)$ to zero.

In unstable configurations (top row of Figure 4), we find that the value of $E_{y}$ flips between positive and negative, with asymptotes in the diffusivity space where $E_{y}\to\pm\infty$. This may be indicative of collective effects such as the presence of MRI channel modes.

In stable configurations (second row of Figure 4), both Ohmic resistivity and ambipolar diffusion result in significant roughly linear increases in $E_{y}$, and facilitate radially outward diffusion ($E_{y}>0$) of the $\boldsymbol{B}$ field when $B_{x}(\infty)>0$ (i.e. field bends outwards) as in Figure 2. This agrees with the prediction from our simple analysis in Section 2.1.2 for the Ohmic case. For the ambipolar case, we expect the variation of $E_{y}$ with $\eta_{A}$ to follow the same pattern as that with $\eta_{O}$ as long as the ‘Ohm-like’ term in Equation (21) is significantly larger than the ‘Hall-like’ term in Equation (22). This condition is indeed met for the profiles considered, as both $B_{x}$ and $B_{y}$ throughout the disc are significantly smaller than $1$, leading the disc to be largely dominated by the poloidal field component.

The Hall effect gives no flux transport ($E_{y}=0$) when it is the only non-ideal effect present, but reduces the magnitude of $E_{y}$ when Ohmic and/or ambipolar diffusion are present, tending to a limit as $\eta_{H}$ becomes very large. When the Hall effect is the only non-ideal term present, $B_{y}=0$ throughout the disc, and there is a solution with a purely poloidal field, so that the current and Hall drift are purely in the azimuthal direction.

3.7.3 Outflow driven flux transport

Next, we explore the flux transport driven by outflow alone through setting a non-zero negative value for $B_{y}(\infty)$. This provides a crude way of mimicking the effect of angular momentum lost vertically in an outflow. We eliminate the effect of inclination driven flux transport by setting $B_{x}(\infty)$ to zero.

Again, in the unstable configurations (third row of Figure 4), we find that the value of $E_{y}$ flips between positive and negative, with asymptotes in the diffusivity space where $E_{y}\to\pm\infty$ indicative of collective effects such as the presence of MRI channel modes.

In stable configurations (bottom two rows of Figure 4), there is radially inward flux transport even when all diffusive effects disappear. This is due to flux being dragged in with the accreting gas as angular momentum is lost vertically through the disc. When $\lvert B_{y}(\infty)\rvert<1$ (fourth row of Figure 4), both Ohmic resistivity and ambipolar diffusion lead to roughly linear but small increases in the magnitude of the flux transport, keeping the same radially inward direction. This is expected as the “Ohm-like” component of the ambipolar term would dominate since the disc field is largely poloidal in nature. The Hall effect, on the other hand, contributes significantly to radially inward flux transport when coupled with the outflow, with roughly linear increases between $\lvert E_{y}\rvert$ and $\eta_{H}$. When $\lvert B_{y}(\infty)\rvert>1$ (bottom row of Figure 4), variations with Ohmic and Hall coefficients remain the same as before, but the trend with ambipolar diffusivity is reversed. Increasing the ambipolar effect now leads to a marked decrease in the inward flux transport, tending to a limiting value at high ambipolar diffusivities. This agrees with the picture described in Equation (5) when the “Hall-like” component (22) of the ambipolar term dominates because of the large poloidal field now present in the disc. This component operates in the opposite direction to the Hall effect in the aligned polarity configuration, hence leads to a dampening of the inward flux transport when present.

3.8 Incompressible limit - constant density profile

To help us understand and interpret the results, we examined solutions of the same set of equations in the incompressible limit for two density profiles: (i) a constant density profile with prescribed disc height, and (ii) a Gaussian density profile that takes into account the hydrostatic balance in the disc, but ignores the effect of magnetic compression. Here, we present the results of (i), which admit approximate analytic solutions, and how they inform the relations between flux transport, diffusive effects and inclination/outflow.

3.8.1 Approximate analytic solutions

We solve the same system of Equations (26)$-$(30) and set $\rho=\mathrm{constant}$. The permeability of free space, $\mu_{0}$ is set to 1. We assume a Keplerian disc and set the orbital parameter $q=-\partial\ln{\Omega}/\partial\ln{r}$ to $3/2$. For simplicity and the possibility of admitting analytic solutions, we also assume that $\eta_{O}=\mathrm{constant}$, and constant ionisation fractions in the disc. This means $\tilde{\eta}_{H},\tilde{\eta}_{A}=\mathrm{constant}$, where they are field-independent coefficients given by

[TABLE]

[TABLE]

This leads to the following form of the modified Ohm’s law including the non-ideal effects:

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

The ODE system is linear in the Ohmic and Hall only cases, and can be reduced to give us the following wave equation:

[TABLE]

with $k^{2}$ given by

[TABLE]

where $v_{az}=B_{z}/\sqrt{\rho}$ is the vertical component of the Alfvén velocity.

If we expect the disc to satisfy the symmetry and boundary conditions described in Section 3.3, then the relevant solution is

[TABLE]

where the surfaces of the disc are at $z=\pm H$, and above and below the disc we have $B_{x}=\pm B_{x}(\infty)$ respectively.

The other variables can be deduced from $B_{x}$ and the bonudary conditions. If $B_{y}=\pm B_{y}(\infty)$ at the upper and lower boundaries (where $B_{y}(\infty)$ can be zero), we have

[TABLE]

The ambipolar contribution, on the other hand, has both a linear and a nonlinear part. The linear part behaves as an additional Ohmic resistivity. The nonlinear part is of order $B_{\text{horizontal}}^{2}$, which is small if the surface value of $B_{\text{horizontal}}$ is also small. We limit our analysis here to the cases when $B_{y}(\infty)<1$, and therefore only the “Ohm-like” component of the ambipolar term contributes significantly to the flux transport. The system can hence be roughly solved as the same set of linear equations with a modified Ohmic resistivity $\eta_{O,\text{mod}}=\eta_{O}+\tilde{\eta}_{A}B_{z}^{2}$ when ambipolar diffusion is also included.

3.8.2 Relation to Balbus &

Terquem (2001)

Equation (40), which gives us the wavenumber of the field solution, is mathematically identical to the wavenumber formula (Equation (46)) of a marginally stable Ohm and Hall only MRI mode derived in Balbus & Terquem (2001). The problem they examined, however, is different from the one investigated here. They determined the local stability of a differentially rotating disc threaded by a weak vertical field and were looking at plane wave disturbances of the form $\exp{(ikz-i\omega t)}$. The Boussinesq limit was used as it corresponds to fluid displacements within a local approximation in $z$, rendering vertical structure unimportant. Our work, on the other hand, solves for the disc vertical structure, but deliberately makes the physically unrealistic uniform density approximation to simplify the equations. Our solution therefore must obey the boundary conditions determined by field inclination and symmetry criteria, which are not required in the problem examined by Balbus & Terquem (2001).

3.8.3 Marginal stability and variation with magnetisation and diffusivities

As in Section 3.4, we examine the case of a disc at marginal stability to the MRI. In our approximate analytic model, because the equations are already linear, we solve the same set of equations for the perturbed fluid variables. As before, these perturbations have opposite symmetry about the midplane to the fluid variables, which in our case means that $\delta B_{x}\propto\cos{kz}$. The boundary condition is now $\delta B_{x}=0$ and $\delta B_{y}=0$, to meet the fixed inclination and outflow conditions. We expect the marginally stable mode to allow just one bending of the field through the disc, hence we require the corresponding critical wavenumber to satisfy $k_{\text{crit}}H=\pi/2$. The condition for marginal stability then becomes:

[TABLE]

We can recast this in terms of a dimensionless magnetisation parameter:

[TABLE]

where $\Sigma=2\rho H$ is the disc surface density, giving us

[TABLE]

Magnetisations that give us a stable solution can be determined when the left hand side of Equation (45) is larger than the right hand side. We rearrange the terms in the equation to give us the following quartic, with the inequality giving us the condition for a stable solution as:

[TABLE]

where we have defined the following dimensionless parameters:

[TABLE]

[TABLE]

[TABLE]

The shape of $f(\mu)$ vs. $\mu$ (see Figure 5) is determined by the Hall diffusivity only, while Ohmic diffusivity stabilises the entire profile by adding a constant term to $f(\mu)$. At high magnetisations, the disc is always stabilised whether the field is aligned with the rotation or not. In the intermediate region, a large Hall parameter extends the region for instability when the polarity is negative, but has a stabilising effect when the polarity is positive. This agrees with our previous result that the Hall effect is responsible for destabilising the disc in the anti-aligned case, while it is stabilising in the aligned case.

3.8.4 Flux transport

For flux transport, the constant density model gives us this simple relation:

[TABLE]

with the radial flux transport velocity given by $E_{y}/B_{z}$.

We can see that flux transport is directly proportional to the inclination and Ohmic diffusivity, which is also the case in the numerical solutions of the compressible model. When no diffusive effects are present, the relation suggests that the presence of an outflow would still drive inward flux accretion, with a large magnetisation (small $\rho$) giving the fastest rate. This is indeed what we observed in the compressible model, and can be interpreted to be a result of the loss of angular momentum driving inward gas accretion, which in turn advects the flux frozen into the gas along with it. A larger magnetisation would mean a higher degree of flux freezing, hence advection with the accretion flow.

The flux transport rate modification due to the Hall effect coupled with inclination is not explained here, but the constant density model predicts that the Hall effect coupled with a non-zero $B_{y}$ at infinity would also lead to a non-zero flux transport. A positve $\eta_{H}$ coupled with the negative $B_{y}^{+}$ expected for a wind would therefore lead to a negative $E_{y}$ proportional to $\eta_{H}$ and $B_{y}^{+}$, signifying a radially inward accretion of flux. This is indeed what we find in the solutions of Section 3.7.3. Such an effect has been noted before in passing in Hall-wind shearing box simulations (Bai, 2014) and also in the global simulations of Bai & Stone (2017). It is also noteworthy that the $\eta_{H}$ term in our flux transport equation is essentially of the same form as Equation 10 of (Bai & Stone, 2017), although our analytic model is more rigorous and would therefore provide a more accurate estimate.

3.9 Incompressible limit - Gaussian density profile

For our second incompressible limit, we imposed a Gaussian density profile that takes into account the balance between pressure and gravity in an isothermal disc, but ignores the effect of magnetic compression from the horizontal field components. Solutions are calculated numerically using the same method as the compressible case, but with the density profile fixed. Comparing the solutions of this model to that of the compressible model allows us to determine whether magnetic compression has a significant effect on the flux transport observed.

The Gaussian density profile imposed is the following:

[TABLE]

where $\rho_{0}$ is the midplane density value.

We found that solutions from the Gaussian incompressible model are very similar both in form and magnitude to those of the compressible model. The only differences arise where bifurcations in the $E_{y}$-diffusivity space occur in the compressible model. There, the Gaussian incompressible model yields vertical asymptotes on either side of the point of divergence instead. This can be explained by the fact that as the instability encounters a point of resonance, magnetic compression can provide the required change in density structure to dampen its effect, hence allowing the solution branches to diverge rather than tend towards a vertical asymptote. Therefore we can conclude that magnetic compression has no significant effect on the flux transport, in the parameter space we have explored.

4 The Guilet & Ogilvie approach

Our second semi-analytic model follows the scalings used in Guilet & Ogilvie (2012) (hereafter GO1). This scheme considers the case when the vertical field $B_{z}$ dominates over $B_{r}$ and $B_{\phi}$ by a factor of $\epsilon^{-1}$. The motivation for this is to put the effects from large scale radial gradients (the $\partial_{r}$ terms) on the same footing as inclination and outflow. Again, instead of using the residual velocity ${\boldsymbol{v}}={\boldsymbol{u}}-r\Omega\hat{\mathbf{e}}_{\phi}$, we use the full velocity ${\boldsymbol{u}}$ (we do not relabel it as ${\boldsymbol{v}}$ here). Compared to the multiscale approach, we make the following adjustments:

[TABLE]

The next order terms are assumed to be dominant, and the higher order terms are neglected. We extend the analysis of GO1 to derive the case for a laminar disc with all non-ideal diffusivities present.

4.1 Leading order equations

The leading order radial component of the equation of motion is

[TABLE]

which describes the centrifugal force balance against the inward gravitational pull. For a Keplerian disc this gives us $\Omega_{0}=(GM/r^{3})^{1/2}$.

The vertical component of the equation of motion at leading order gives us the hydrostatic equilibrium between the vertical gravitational force and the vertical pressure gradient:

[TABLE]

For an isothermal disc, we have the familiar solution

[TABLE]

where $H_{0}=c_{s}/\Omega_{0}$ is the isothermal scaleheight and $\Sigma_{0}(r,\tau)$ is the surface density.

The horizontal components of the equation of motion and the induction equation come at the next order, due to the terms that we set to zero in Equation (52). They are identical to the equations in GO1 except for the addition of the Hall and ambipolar terms:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Under these scalings, we can see that Ohmic and ambipolar diffusivities have the same effect on the disc dynamics. The Hall term comes in at $\pi/2$ phase difference with the Ohmic and ambipolar terms in the induction equation. This reflects the nature of the Hall term being the magnetic field cross multiplied once with the Ohmic term, while the ambipolar term is cross multiplied twice with the Ohmic term hence acts in the same direction at leading order.

This analysis yields four linear equations for the unknowns $u_{r1}$, $u_{\phi 1}$, $B_{r1}$ and $B_{\phi 1}$. The linearity is a result of the assumption that we are examining the case of small deviations from orbital motion and a vertical magnetic field.

4.2 Non-dimensionalisation

We follow the same approach as GO1 in non-dimensionalising our equations. We relegate the definitions to Appendix C (see online supplementary materials), but note that $\tilde{\rho},u_{r},u_{\phi},b_{r},b_{\phi},\tilde{\eta}_{O},\tilde{\eta}_{H},\tilde{\eta}_{A},\zeta$ are the non-dimensionalised forms of $\rho,u_{r1},u_{\phi 1},B_{r1},B_{\phi 1},{\eta}_{O0},{\eta}_{H0},{\eta}_{A0},\zeta$ respectively.

We assume a point-mass potential and circular Keplerian orbital motion at leading order. An isothermal equation of state is also employed for simplicity.

The density profile in dimensionless form is

[TABLE]

and the differential equations (56)$-$(59) are rewritten as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A key thing to note here is that the LHS of Equations (61)-(64) forms the differential system, while the RHS are source terms that drive the advection and diffusion of flux. These equations are in many ways similar to those we obtained in Section 3.1. The first two terms on the LHS of Equation (61) correspond to the two terms in Equation (26), with the only difference being the addition of large scale radial gradient and Keplerian source terms on the RHS. Equation (62) is identical to Equation (27), while Equations (63) and (64) are linearised versions of Equations (29) and (30) respectively.

4.3 Boundary conditions

We use the same boundary conditions as determined in Section 3.2 of GO1, which come from analysing the same expected symmetry of the solutions as Section 3.3, with $v_{r1}$ and $v_{\phi 1}$ being even functions of $\zeta$, while $B_{r1}$ and $B_{\phi 1}$ are odd. The following quantities can then be determined to vanish exponentially fast as $\zeta\to\pm\infty$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The symmetry of the solutions also constrain the following mid-plane ($\zeta=0$) values to be the following:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A non-vanishing $b_{\phi s}$, the value of $b_{\phi}$ at $\zeta\to\infty$, is again included for mimicking the effect of angular momentum removal by a MCW or magnetic braking due to interaction with an external medium. As noted in GO1, the boundary conditions are homogeneous, except for the linear source terms proportional to $D_{B}$, $b_{rs}$ and $b_{\phi s}$ in Equations (71) and (72).

4.4 Form of the solution for a laminar inviscid disc

The linearity of the equations means that the general solution is a linear combination of the solution vectors corresponding to each source term, which appear either on the right hand side of equations (61)$-$(64), or as a non-vanishing boundary condition at infinity in the form of $b_{rs}$ and $b_{\phi s}$. We can thus write the general solution $\boldsymbol{X}=\{u_{r},u_{\phi},b_{r},b_{\phi}\}$ as:

[TABLE]

where $\boldsymbol{X}_{\mathrm{DH}}$ is the solution vector corresponding to the source term proportional to $D_{H}$ and so on. $\boldsymbol{X}_{K}$ corresponds to the solution vector when $D_{H},D_{\nu\Sigma},D_{B},b_{rs},b_{\phi s}=0$, where the source terms arise only from the radial derivatives of the leading order Keplerian, gravitational and geometric terms. Following GO1, we also define

[TABLE]

as the solution with hydrodynamic (‘Hydro’) source terms $D_{H}=1$ (corresponding to a disc with a constant aspect ratio $H/r$) and $D_{\nu\Sigma}=0$ (for a steady accretion flow far from the inner boundary). Under these definitions, we have

[TABLE]

The solution depends in a non-linear way only on the parameters $\beta_{0}$, $\eta_{O}$, $\eta_{H}$ and $\eta_{A}$ (we drop the $\tilde{}$ in the text from this point onwards but refer to the non-dimensionalised diffusivities). Since we consider a laminar inviscid disc, we set $\alpha$ to zero, and can neglect the $\mathcal{P}$ dependence. $\eta_{O}$ and $\eta_{A}$ have the same effect at the order we are considering, so we only need to examine one of them, which we choose to be $\eta_{O}$. For each triplet of values of the three parameters $\beta_{0}$, $\eta_{O}$ and $\eta_{H}$, one needs to compute the six solution vectors for each of the terms on the right hand side of Equation (77). The general solution is then given by a linear combination of these solution vectors with the appropriate coefficients.

4.5 Numerical method

To solve the system of ODEs, we employ the same pseudo-spectral method of GO1 using a decomposition on a basis of Whittaker cardinal functions. These functions are well suited to problems on an infinite interval (Boyd, 2001; Latter et al., 2010). The use of Whittaker cardinal functions implicitly imposes the condition that the variables have to vanish exponentially fast at infinity. Following GO1, we replace $b_{r}$ and $b_{\phi}$ with the following variables:

[TABLE]

[TABLE]

to satisfy this condition. We can then see that these new variables do vanish exponentially fast at infinity from the boundary conditions given in equations (71) and (72).

4.6 Flux transport for uniform diffusivities

Flux transport is again calculated by integrating the radial component of the induction equation:

[TABLE]

where

[TABLE]

is the dimensionless magnetic flux transport velocity, same as the one defined in equation (32).

For simplicity, we used constant diffusivity profiles in our calculations. We present here only results from the positive polarity case of $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>0$, as most solutions for the negative polarity case failed to converge, indicative of unstable configurations. Representative plots of $u_{\Psi}$ variation with diffusivities can be found in Figure 6. For all our solutions, we isolate a particular source term and set its coefficient to $1$. The case for $b_{\phi s}=1$ corresponds to a torque that spins the disc up, leading to decretion rather than the normal accretion. To examine the case for $b_{\phi s}=-1$, which mimics the effect of a MCW, we simply need to reverse the sign of $u_{\Psi}$ when interpreting the plots.

4.6.1 Variation with Ohmic resistivity

Characteristic plots of $u_{\Psi}$ variation with $\eta_{O}$ are displayed in the left and right columns of Figure 6. We found that qualitatively, the ‘Hydro’ and $D_{B}$ source terms share the same trend, while the $b_{rs}$ and $b_{\phi s}$ source terms follow a different trend. As the midplane field strength is decreased ($\beta_{0}$ increased), solutions change from a single smooth curve, to multiple curves separated by asymptotes in $\eta_{O}$. The multiple curves signify the onset of instabilities, and this result is in agreement with the picture given by our model in Section 3, with instabilities setting in below a critical field strength. The same is observed for when $\eta_{O}$ increases. Beyond a threshold diffusivity, the solution is stable and has only one branch.

In the stable configurations, both ‘Hydro’ and $D_{B}$ source terms contribute to radially outward flux transport ($u_{\Psi}>0$). Away from the asymptote, which appears to indicate a region of instability, the flux transport velocity increases fractionally ($<10\%$) as $\eta_{O}$ is increased.

On the other hand, for both $b_{rs}$ and $b_{\phi s}$ source terms, an increase in $\eta_{O}$ leads to a similar order of magnitude increase in $u_{\Psi}$, which is postive in the stable region. This confirms both the picture in our previous model (see Section 3.7) that a positive inclination of the poloidal field away from the star when coupled with diffusivity drives outward flux transport, and also that the rate at which it does so correspond to roughly linear increases. Remembering to reverse the sign for an actual wind, the results for $b_{\phi s}$ show us that a wind, coupled with resistivity, leads to an accretion of flux ($u_{\Psi}<0$), again in agreement with our results from the shearing box model.

4.6.2 Variation with Hall diffusivity

When the Hall term was the only non-ideal MHD effect present, the solver did not converge. Analytically, we can determine that we would have iso-rotation and no flux transport. We therefore examined the cases where a background Ohmic resistivity was also present.

$u_{\Psi}$ variation with $\eta_{H}$ are found to share the same trend between the ‘Hydro’, $D_{B}$ and $b_{rs}$ source terms, while those for the $b_{\phi s}$ source term follow a different trend. Again, as the midplane field strength is decreased, solutions go from a single smooth curve, to multiple curves separated by asymptotes in $\eta_{H}$. The same is observed for when $\eta_{H}$ increases. Beyond a threshold diffusivity, the solution is stable and only has one branch.

In the stable configurations, ‘Hydro’, $D_{B}$ and $b_{rs}$ source terms contribute to a radially outward flux transport ($u_{\Psi}>0$). Away from the final asymptote that marks the region of instability, the flux transport velocity decreases to zero as $\eta_{H}$ is increased. This is similar to the picture in our previous shearing box model that indicates a decrease in the inclination driven flux transport rate with an increase in the Hall coefficient, with the only difference being that the nonlinear Hall term in our previous model leads only to a minor correction of the flux transport rate, whereas the linearised Hall term here reduces it to zero. On the other hand, for the $b_{\phi s}$ source term, an increase in $\eta_{H}$ leads to a similar order of magnitude increase in $u_{\Psi}$. $u_{\Psi}$ is positive in the stable region. Therefore for $b_{\phi s}=-1$, mimicking the effect of a MCW, increasing the Hall contribution leads to a significant accretion of flux ($u_{\Psi}<0$), which is again in agreement with our results from the shearing box model.

4.7 Analytic models

To help us interpret the trends in flux transport, we developed three approximate analytic models following the same procedure as Section 4.1 and Appendix A of GO1. Full mathematical details of how these three models were calculated can be found in Appendix D (see online supplementary materials), but here we give a brief qualitative description of the procedure and assumptions used to derive these models.

We split the disc into two regions: a passive field region with weak (passive) magnetic field ($\beta\gg 1$) around the midplane where hydrodynamic effects dominate over magnetic effects; then further up the disc, we have a region with strong magnetic field ($\beta\ll 1$) where magnetic effects are dominant and the field is approximately force-free. The transition point between the regions is where the magnetic pressure is equal to the thermal pressure, and is given by:

[TABLE]

We can see from this that the lower the field strength, the larger the value of $\zeta_{B}$, and the better the approximation of the midplane region as under a ‘passive’ field.

The models are constructed by first calculating the general forms of the approximate analytic profiles in each of the two regions, constrained by the given midplane and disc surface ($\zeta\to\infty$) boundary conditions. We assume that the Lorentz force is negligible in the passive field region hence ignore the effect of the magnetic field in the velocity profile there. In the force-free region, nothing can compensate the Lorentz force, and the fluid is frozen into the magnetic field lines, with iso-rotation being enforced as in the ideal MHD case. The two regions are then connected across the transition region, which is assumed to be infinitesimally thin about $\zeta_{B}$. When connecting the regions, we assumed continuity of the magnetic field components across the transition, and integrated the horizontal components of the induction equation over the transition region to find the appropriate boundary conditions.

The first two of the models represent different diffusivity regimes: the first assumes only the presence of a constant Ohmic diffusivity, while the second addresses the situation when the disc is dominated by Hall drift but with a small Ohmic contribution. Our third model improves on the first model by using better boundary conditions that more appropriately address the presence of the intermediate transition region, and is derived using the description outlined in Appendix A of GO1.

4.7.1 Source terms and disc vertical profiles

Figures 7 and 8 show the vertical structure profiles in the case of a disc with only Ohmic diffusion and a Hall dominated disc with small Ohmic contribution respectively. The red lines are from the simple two-zone models that use the boundary condition of continuous magnetic fields, while the purple lines in the Ohm only case are from the improved two-zone model that take into account the intermediate region, and modifies the jump condition in the magnetic fields at the transition point. They both provide good qualitative descriptions of the solution, with the improved model matching very well to the actual numerical solution.

By analysing the mathematical forms of the analytic solutions (see Appendix D of the online supplementary materials), we can deduce how the various source terms affect the shape of the disc vertical profile in the different variables. The radial magnetic gradient $D_{B}$ provides the background gradient in $\zeta$ for the profile of $b_{r}$, while $b_{rs}$ sets the limit of $b_{r}(\infty)$ when $D_{B}$ is absent, defining the surface inclination of the poloidal field. Similarly, $b_{\phi s}$ sets the value of $b_{\phi}(\infty)$, and hence the magnetic torque acting on the disc. Ohmic resistivity causes bends to occur in the passive field region for the $b_{\phi}$ profile, while the Hall drift, which operates at $\pi/2$ phase to Ohmic diffusion, causes the bends to happen in $b_{r}$ instead. $u_{r}$ tends towards the flux transport value as $\zeta\to\infty$, but is largely zero in the passive field region. $u_{\phi}$ is similarly very small in the disc, but increases drastically in the force-free region as a result of iso-rotation with the magnetic fields.

4.7.2 Flux transport rates

Here, we analyse only the flux transport driven by inclination and outflow, while we leave those due to large scale radial gradients to a future investigation. Analytically, we find that the flux transport is given by these simple expressions:

[TABLE]

for the Ohm only case,

[TABLE]

for the Hall dominated case, and

[TABLE]

for the Ohm only case with improved boundary conditions accounting for the transition region.

These expressions do not involve the radial advection velocity in the ideal MHD limit like second term in Equation (50) because of the separation of the disc into two zones, with the flux transport value determined by the expression in the passive field region where the only non-diffusive contribution, $u_{r}$, is approximated to be zero. In the absence of any diffusivity, it can be easily shown that the original equations admit an analytic solution with uniform $u_{r}$ throughout the disc, and the flux transport is simply due to advection via the wind torque, with $u_{\Psi}=u_{r}=4b_{\phi s}/\beta_{0}$. Hence our simple two-zone analytic models should only be used in the high diffusivity regime, as is common in PPDs.

As in the uniform density shearing box model (see Section 3.8), we find Ohmic diffusivity to be linked with inclination, and Hall diffusivity with outflow, in driving flux transport. Again, it is worth noting that the $\eta_{H}$ term in Equation (86) is essentially of the same form as Equation 10 of (Bai & Stone, 2017). In the modified analytic model where jump conditions in the magnetic fields are calculated by taking the transition region into account, we have a modified flux transport rate which also couples Ohmic resistivity with outflow.

Figure 9 shows how these analytic predictions fit the actual flux transport rates calculated. We first look at the inclination cases which are the top two rows of the figure. Here the red lines are from the simple two zone models, while the orange lines in the Ohm only cases are from the third model that includes the modified boundary conditions. We see that the analytic models offer good descriptions in the limit of weak field ($\beta_{0}$ large) and high diffusivity. This is because under these conditions, the passive field region is extended and better matched by the assumptions used in the analytic model. Particularly, the asymptotes of the modified analytic model appears to be able to predict the transition point to instability, which occurs when the denominator of Equation (87) is equal to zero. It would be worth investigating in the future how the region of instability is influenced by the balance between diffusive effects (from $\eta_{O}$) and the field strength (from $\zeta_{B}$).

Next, we turn to the outflow cases, which are displayed in the bottom two rows of Figure 9. The orange lines in the top two Ohm only plots are from the modified analytic model, while the red lines in the bottom two Hall dominated panels are from the simple analytic model. We can see that again they predict the outflow driven flux transport fairly well in the high diffusivity, weak field ($\beta_{0}$ large) regime. The simple two zone models are unable to describe the presence of the asymptote and transition to the unstable regime, but the model with modified boundary condition matches the first asymptote in the Ohm-only case in the weak field limit.

The relative success of the analytic models so far (though with limitations, such as in the outflow case) in matching the numerical solutions suggest that it is also possible to understand the flux transport driven by other source terms similarly, and we will be looking at constructing these models and examining their interpretations in the future.

5 Discussion

5.1 Comparing the shearing box and Guilet & Ogilvie models

Both the shearing box and the GO models share the same qualitative trends of how diffusivities, coupled with inclination and outflow, can be effective in driving flux transport. In particular, approximate analytic solutions for both models under simplified schemes show how Ohmic resistivity (and the ‘Ohm-like’ term of ambipolar diffusivity) is coupled with inclination to drive radially outward flux transport, while Hall drift with $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>0$ is coupled with outflow in facilitating radially inward flux transport. The similarity between the expressions in Equations (50), (85) and (86) suggest that we can identify $\zeta_{B}$, the height at which magnetic pressure equals thermal pressure, as the relevant value for the disc height $H$ to be used in the constant density shearing box model of Section 3.8. Both models also show similar trends and behaviour in how the stability of solutions vary with diffusivity and magnetisation values. The large similarities between the results of the two models may be due to the fact that horizontal magnetic field components are largely relatively small compared with the vertical field for the parameter space explored, hence the effects due to the nonlinearity of the Hall and ambipolar terms are less significant overall.

5.2 Comparison with current global simulations

While outward inclination coupled with Ohmic resistivity has long been known to facilitate outward flux transport, the possibility of Hall drift coupled with an outflow driving significant flux transport has only been briefly noted (Bai, 2014) in the past, and never extensively investigated. Only recently have global simulations (Béthune et al., 2017; Bai & Stone, 2017; Suriano et al., 2017, 2018) begun exploring the flux transport problem in laminar PPDs where all three non-ideal effects are accounted for. In Bai & Stone (2017), the mechanism governing flux evolution can be described as a competition between the Hall effect, ambipolar diffusion, and the MCW-driven accretion. In the positive (negative) polarity configuration, the Hall effect transports flux (out)inwards rapidly in the midplane regions and (in)outwards slowly in the disc upper layers. Ambipolar diffusion, on the other hand, always transports flux outward, whereas MCW driven accretion advects flux inwards. In the aligned case, there is a cancellation effect between the inward flux transport due to Hall and the MCW and outward due to ambipolar, while all three effects work to transport flux outward in the negative polarity case. Our calculations also found a similar picture in that Hall drift reduces the effect of outward flux transport due to ambipolar/Ohmic diffusion in the positive polarity case. However, it is unclear from the results of Bai & Stone (2017) whether there is any specific coupling between Hall drift and the wind that further enhances flux transport inwards. Our models also agree with Béthune et al. (2017) and Bai & Stone (2017) that the polarity of the magnetic field is a significant parameter on disc dynamics and flux transport when Hall drift is present. However, we are not able to confirm the flux transport trends in the negative polarity cases due to the breakdown of stable solutions in our models. Our model restricts the overall flux transport rate in a radially local region to be constant across all scale-heights, which is indeed what is found in global simulations (Bai & Stone, 2017; Bai, 2017).

A significant proportion of the parameter space explored yielded unstable solutions characteristic of MRI channel modes, and it would be interesting to see in future studies if they might relate to flux transport mechanisms reported in global simulations that are cyclic in nature, such as those in Suriano et al. (2017, 2018, 2019).

5.3 Implications, model limitations and future directions

One important implication of our results is the need to further investigate the interplay between Hall drift and a wind outflow in facilitating flux transport. Bai (2016), which attempted to present a framework of global PPD evolution incorporating the latest advances in PPD research, found disc evolution to be largely dominated by wind-driven processes, while viscous spreading is suppressed. He noted that the timescale of disc evolution is largely governed by the global flux evolution in the disc, and admitted that an understanding of this and in particular of how the Hall effect affects magnetic flux evolution is still lacking. To bridge this gap, a more detailed way of modelling the wind outflow in the presence of the Hall effect than our simple presciption of a non-zero surface $B_{\phi}$ would be required, where disc boundary conditions are formally matched to the Blandford & Payne solutions, and where we allow for a vertical outflow $v_{z}$. The feasibility of such approach, tailored for the flux transport problem, will be examined in future work.

The question of how to model asymmetrical solutions semi-analytically should also be looked into. Disc solutions with vertical symmetry breaking have been observed in many recent global disc simulations (Béthune et al., 2016, 2017; Bai, 2017). They have drastic contrasts in their overall accretion and flux transport rates compared to their standard symmetric counterparts. It would be worth investigating the effect of symmetry breaking or having the opposite symmetry to the one prescribed in the current work using our semi-analytic approach in future work. Currently, switching the even-odd symmetry of the variables would yield the trivial result of zero net flux transport in our local models, as the transport from the bottom and top parts of the disc will cancel each other out. Modelling the flux transport meaningfully in such context therefore would require careful thought.

Future investigations should also seek to implement more realistic diffusivity profiles, such as those used in these simulations (Béthune et al., 2016, 2017; Bai, 2017). The flux transport mechanism reported in Suriano et al. (2018) is heavily dependent on the variation on the ambipolar diffusivity in the disc, and involves a time-dependent periodic flux concentration processes.

Given the theoretical work nature of the present work, shearing box simulations would be needed to verify the equilibrium solutions found by our semi-analytic models. The advantage of running a simulation is that it gives us insight into how the disc develops over time before reaching equilibrium, and it would be interesting to see whether any structures with periodic cycles are formed in the unstable regimes, and how they may inform us about the flux transport mechanisms reported in global simulations.

Finally, it would be good to combine these local models into a global flux evolution framework, such as done by Guilet & Ogilvie (2014); Okuzumi et al. (2014) and Takeuchi & Okuzumi (2014). The main challenge resides in the fact that Hall/wind-driven flux transport requires knowledge of the disc surface toroidal field, while the previous models only required the calculation of the poloidal field structure which is more readily determined through the Biot-Savart law. Careful modelling of the global wind solution is required, and will be looked into in a future investigation. In the absence of a magnetic wind, the Hall effect (when $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>0$) reduces the radially outward flux transport rate induced by Ohmic and/or ambipolar diffusion coupled with an inclined field. Hence we should expect the PPDs where the Hall effect is present, field is aligned with rotation, and where no outflow is launched, to settle into a more highly magnetised steady-state solution than those previously calculated.

6 Summary and conclusions

In this work, we have developed a formalism to compute the radially local effects of non-ideal diffusivities coupled with inclination, outflow and large scale radial gradients on disc dynamics through the use of a semi-analytic multiscale asymptotic approach. We investigated the flux transport due to inclination and outflow first at the shearing box order, and then at the Guilet & Ogilvie order (Guilet & Ogilvie, 2012), which also allows the computation of the additional contributions from large scale radial gradients. Our findings from both models are qualitatively similar, and we examined the trends of inclination and outflow driven flux transport and their variation with the three diffusivities. Using approximate analytic models, we gained insights into how diffusivities are coupled with the various other parameters in facilitating flux transport, and derived simple relations for estimating this.

Our main findings include:

Stable disc configurations arise from having a strong field (low $\beta_{0}$) and high diffusivity values, while weak field and low diffusivities give rise to unstable configurations characteristic of MRI channel modes. 2. 2.

In the positive polarity case, where $\eta_{H}(\boldsymbol{B}_{z}\cdot\boldsymbol{\Omega})>0$, all diffusivities are stabilising, while the Hall effect becomes destabilising in the anti-aligned case. 3. 3.

Outward inclination of the poloidal field coupled with Ohmic and ambipolar diffusion both drive radially outward flux transport which increases roughly linearly with diffusivity, while the Hall effect coupled with outward inclination reduces the flux transport rate but does not reverse its direction. 4. 4.

When the toroidal field in the disc is small compared to the overall field, outflow driven flux transport is not significantly affected by the presence of Ohmic and ambipolar diffusion, but is significantly enhanced when the Hall effect is present. In the high $\beta_{0}$, high diffusitivity limit common to PPDs, this scales roughly linearly with $\eta_{H}$ and the surface toroidal field strength. All outflow driven flux transport is radially inward, as would be expected from advection due to the accretion flow caused by vertical removal of angular momentum in the outflow. When the toroidal field in the disc is large, the flux transport behaviour with Ohmic and Hall coefficients remain the same, while it is reversed in the ambipolar case, as the "Hall-like" component of the ambipolar term dominates over the "Ohm-like" component, and acts in the opposite direction to the Hall term.

This work represents an initial effort toward modelling flux transport in PPDs incorporating all three non-ideal effects, inclination of the large scale field, outflow and the presence of large-scale radial gradients. At present, we have focused on the local transport in the disc, with a range of parameters characterising the contribution from each effect, and assumed a quasi-steady equilibrium state. Future work would need to address the time-dependent aspect of disc evolution, and also consider the flux transport globally under more realistic physical parameters for better comparison and understanding to the results of present and future PPD simulations.

Acknowledgements

We would like to thank the annonymous referee for a detailed report which improved the presentation of the manuscript. PKCL would like to thank Jérôme Guilet for useful discussions, and for letting him have a copy of the code used for the 2012 Guilet & Ogilvie paper. He would also like to thank the Croucher Foundation and the Cambridge Commonwealth, European & International Trust for their generous support in funding his PhD studentship through a Cambridge Croucher International Scholarship.

Appendix A Multiscale approach equations at different orders

A.1 Mass and momentum equations

The mass and momentum equations, assuming axisymmetry, are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

and the current is given by

[TABLE]

We introduce the following operators:

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

Also

[TABLE]

with

[TABLE]

[TABLE]

Assuming the potential of Equation (7), the leading-order equations are then

[TABLE]

[TABLE]

[TABLE]

[TABLE]

At the next order, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

A.2 Induction equation

We assume axisymmetry throughout this analysis. The full induction equation is given by

[TABLE]

where $\boldsymbol{b}=\boldsymbol{B}/\lvert\boldsymbol{B}\rvert$, and $\eta_{O},\eta_{H},\eta_{A}$ are the Ohmic, Hall and ambipolar diffusivities respectively.

We define

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

Note that this is of the same form as $\mathrm{D}$ with the exception of $\partial_{t},\partial_{\tau}$ terms.

At zeroth order, all the terms cancel out.

The order $\epsilon$ equations are:

[TABLE]

[TABLE]

[TABLE]

The next order $\epsilon^{2}$ equations are:

[TABLE]

[TABLE]

[TABLE]

We also have the zero divergence condition

[TABLE]

which when written in cylindrical coordinates, give:

[TABLE]

The resulting equations in increasing order are therefore:

[TABLE]

[TABLE]

Appendix B Equations for numerical solver in the 1D radially local vertical structure model

For the numerical solver, we recast Equations (26) to (30) to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We define the diffusion constants as the following:

[TABLE]

[TABLE]

such that

[TABLE]

Then the currents are given by

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Guilet & Ogilvie approach details

C.1 Non dimensionalisation

We used the following definitions:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A dimensionless vertical spatial coordinate is also defined as

[TABLE]

This $\zeta$ is different from the rescaled dimensional variable $\zeta$ defined earlier.

We can then define dimensionless laminar diffusivities:

[TABLE]

[TABLE]

[TABLE]

Appendix D Mathematical derivation of analytic models for analysing GO model

D.1 Analytic model Ohmic diffusion only

This section follows the same procedures as Section 4.1 of Guilet & Ogilvie (2012). We assume constant Ohmic diffusivity throughout the disc.

D.1.1 Passive field regime

$\beta\to\infty$ and the equations become:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The azimuthal velocity can then be written as:

[TABLE]

[TABLE]

[TABLE]

From (155), and the boundary condition $b_{r}(0)=0$, we can deduce that:

[TABLE]

We can then find $b_{\phi}$ by substituting (157) and (160) into (156), giving us:

[TABLE]

This, with the boundary condition $b_{\phi}(0)=0$, has the solution

[TABLE]

The flux transport in this region is given by Equation (82):

[TABLE]

D.1.2 Force-free regime

$\beta\ll 1$, nothing can compensate the Lorentz force resulting in force-free magnetic fields with current parallel to the field lines. The equations become:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For our case where the field is almost vertical, the radial and azimuthal currents vanish. Diffusive effects no longer have any importance, and:

[TABLE]

[TABLE]

With our boundary conditions at infinity, we find that throughout the force-free region:

[TABLE]

[TABLE]

where $\pm$ stands for the sign of $\zeta$.

For the azimuthal velocity, substituting (171) and (170) to (167):

[TABLE]

where $u^{\prime}_{\phi 0}$ is an integration constant to be determined.

From Equation (82), the radial velocity is equal to the flux transport velocity and is a constant:

[TABLE]

D.1.3 Two-zone model

The height at which the magnetic pressure is equal to the thermal pressure is given by:

[TABLE]

We assume that for $\zeta<\zeta_{B}$ we have a ‘passive field’, and for $\zeta>\zeta_{B}$, we have a ‘force-free field’. We neglect the thickness of the transition where $\beta$ is of order unity, and determine proper boundary conditions at $\zeta=\zeta_{B}$. Four conditions are needed to constrain the four unknowns.

Two boundary conditions can be obtained from the analysis of the induction equation. Using (163), (160) and (162) at $\zeta=0$, we obtain:

[TABLE]

The azimuthal component of the induction is more complicated because $b_{r}$ acts as a source term. We integrate between two heights $\zeta_{1}$ and $\zeta_{2}$ to get:

[TABLE]

We can use this relation between $\zeta_{B}^{-}$ and $\zeta_{B}^{+}$ to connect the two regions. We neglect the width of the intermediate region, so the RHS of the equation vanishes. This is equivalent to the condition that horizontal electric field components are continuous across the boundary. Including the force-free conditions that $\partial_{\zeta}b_{r}(\zeta_{B}^{+})=D_{B}$ and $\partial_{\zeta}b_{\phi}(\zeta_{B}^{+})=0$, we obtain:

[TABLE]

The other two boundary conditions come from assuming that the magnetic field does not vary significantly at the transition between the two regions:

[TABLE]

[TABLE]

The radial magnetic field profile in the passive-field region is then given by:

[TABLE]

While the azimuthal magnetic field profile in the passive-field region is:

[TABLE]

where

[TABLE]

[TABLE]

The azimuthal velocity field in the passive-field region is given by (157), while in the force-free field region we use (177) to find:

[TABLE]

[TABLE]

To calculate the flux transport, we evaluate the gradients of the horizontal magnetic fields in the passive-field region:

[TABLE]

[TABLE]

This gives us the flux transport velocity as:

[TABLE]

This shows us that the Ohm term can drive flux transport given the presence of an inclination in the poloidal field to the vertical. A large scale radial magnetic gradient also contributes to flux transport.

D.2 Approximate analytic model when $\eta_{H}\gg\eta_{O}$

We follow the same procedure as Appendix D.1 but with $\eta_{H}$ included.

D.2.1 Purely hydrodynamic, $\beta\to\infty$ [Passive field]

The equations become:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From (189), we can write

[TABLE]

[TABLE]

[TABLE]

Assuming constant diffusivities, we have:

(191):

[TABLE]

(192):

[TABLE]

We find $B_{r}$ via (197)$\times\eta_{H}$ $+$ (196) $\times\eta_{O}$, which after some rearrangement gives us:

[TABLE]

This can be solved as a second order ODE using standard techniques:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Applying midplane boundary conditions, we obtain $B_{1}=0$. Hence,

[TABLE]

[TABLE]

Note that when $\eta_{H}=0$, $\partial_{\zeta}^{2}b_{r}=0$, and

[TABLE]

where $b_{r1}$ is a constant to be determined.

We find $B_{\phi}$ via (196)$\times\eta_{H}$ $-$ (197)$\times\eta_{O}$, which after rearranging gives us:

[TABLE]

When $\eta_{O}\neq 0$, $\eta_{H}=0$:

[TABLE]

[TABLE]

When $\eta_{O},\eta_{H}\neq 0$:

[TABLE]

As $\eta_{H}\to 0$, $\kappa\to 0$ but the limit is such that $b_{\phi}\to\pm\infty$.

The flux transport in this region is given by Equation (82):

[TABLE]

D.2.2 Force-free magnetic field

By the same arguments outlined in Appendix D.1.2, we have:

[TABLE]

[TABLE]

With our boundary conditions at infinity, we find that throughout the force-free region:

[TABLE]

[TABLE]

where $\pm$ stands for the sign of $\zeta$.

The absence of a current means that the magnetic field cannot diffuse. Velocity is dtermined by the fact that the fluid is frozen into the magnetic field lines.

For the azimuthal velocity, substituting (214) and (213) to Equation (64) and integrating over $\zeta$:

[TABLE]

The radial velocity is equal to the flux transport velocity and is a constant:

[TABLE]

D.2.3 Two-zone model

The height at which the magnetic pressure is equal to the thermal pressure is given by:

[TABLE]

We make the same assumptions and follow through the argument of Appendix D.1.3 to determine the boundary conditions at $\zeta=\zeta_{B}$ to connect the two zones.

Two boundary conditions can be obtained from the analysis of the induction equation. Using (210), (203) and (209) at $\zeta=0$, we obtain:

[TABLE]

The azimuthal component of the induction is more complicated because $b_{r}$ acts as a source term. We integrate between two heights $\zeta_{1}$ and $\zeta_{2}$ to get:

[TABLE]

We use this relation between $\zeta_{B}^{-}$ and $\zeta_{B}^{+}$ to connect the two regions and neglect the width of the intermediate region, so the RHS of the equation vanishes. The force-free condition gives us $\partial_{\zeta}b_{r}(\zeta_{B}^{+})=D_{B}$ and $\partial_{\zeta}b_{\phi}(\zeta_{B}^{+})=0$. Hence,

[TABLE]

The two other conditions come from the two components of the equation of motion. Assuming that the magnetic field does not vary significantly at the transition between the two regions:

[TABLE]

[TABLE]

The radial magnetic field profile in the passive-field region is given by:

[TABLE]

where

[TABLE]

The azimuthal magnetic field profile in the passive-field region is:

[TABLE]

The azimuthal velocity field in the passive-field region is given by (193), while in the force-free field region it is:

[TABLE]

[TABLE]

The gradients of the horizontal magnetic fields in the passive-field region are:

[TABLE]

[TABLE]

The flux transport is then explicitly calculated as:

[TABLE]

D.3 Ohmic diffusion only incorporating effects from the intermediate transition region

Here, we assume intermediate region to the thin, and concentrate only on the $D_{B}$, $b_{rs}$ and $b_{\phi s}$ source terms (neglecting Keplerian and Hydrodynamic source terms).

The density profile around $\zeta_{B}$ is approximated by:

[TABLE]

where $x\equiv\zeta-\zeta_{B}$. The intermediate region where $\beta\sim 1$ has a typical thickness of $x\sim 1/\zeta_{B}$, and our assumption of a thin transition is valid if $\zeta_{B}^{2}\gg 1$.

The differential system is then rewritten as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We integrate between two heights $\zeta_{1}$ and $\zeta_{2}$ to get:

[TABLE]

As before, we use this relation between $\zeta_{B}^{-}$ and $\zeta_{B}^{+}$ to connect the two regions, neglecting the width of the intermediate region, so the RHS of the equation vanishes. The force-free condition is also employed where $\partial_{\zeta}b_{r}(\zeta_{B}^{+})=D_{B}$ and $\partial_{\zeta}b_{\phi}(\zeta_{B}^{+})=0$. Hence,

[TABLE]

and

[TABLE]

Thus,

[TABLE]

[TABLE]

where $\Delta u_{\phi}=u_{\phi}(\zeta_{B}^{+})-u_{\phi}(\zeta_{B}^{-})=\eta_{O}\partial_{\zeta}b_{\phi}(\zeta_{B}^{-})$.

We substitute these relations into (232) and (233):

[TABLE]

[TABLE]

We can treat them as simultaneous equations in $\partial_{x}b_{r}-D_{B}$ and $\partial_{x}b_{\phi}$:

[TABLE]

[TABLE]

We find $\partial_{x}b_{r}$ via (244)$\times\eta_{O}-$(243)$\times 4\mathrm{e}^{\zeta_{B}x}$, which after rearranging gives:

[TABLE]

We find $\partial_{x}b_{\phi}$ via (243)$\times\eta_{O}-$(244)$\times\mathrm{e}^{\zeta_{B}x}$, which after rearranging gives:

[TABLE]

Integrating over $x$, we compute the jump in $b_{r}$ across the intermediate region if $\mathrm{e}^{\zeta_{B}x}\gg 1$:

[TABLE]

The jump condition is thus:

[TABLE]

Doing the same for $b_{\phi}$, we have:

[TABLE]

The jump condition is thus:

[TABLE]

These new boundary conditions can be injected into the two zone model such that

[TABLE]

[TABLE]

Hence for $b_{r}$, we have

[TABLE]

After some algebra, we obtain

[TABLE]

For $b_{\phi}$, we obtain

[TABLE]

Under these new conditions, the flux transport is:

[TABLE]

[TABLE]