# Non-ideal Magnetohydrodynamic Instabilities in Protoplanetary Disks:
Vertical Modes and Reflection Asymmetry

Lile Wang The Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Department of Astronomy, School of Physics,
Peking University, Beijing 100871, China  Sheng Xu The Kavli Institute for
Astronomy and Astrophysics, Peking University, Beijing 100871, China
Department of Astronomy, School of Physics, Peking University, Beijing 100871,
China  Zhenyu Wang Institute of Plasma Physics, Chinese Academy of Sciences,
Hefei 230031, China  Department of Astrophysical Sciences, Princeton
University, Princeton NJ 08540, USA  Min Fang Purple Mountain Observatory,
Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210023, China
University of Science and Technology of China, Hefei 230026, China  Jeremy
Goodman Department of Astrophysical Sciences, Princeton University, Princeton
NJ 08540, USA Lile Wang<EMAIL_ADDRESS>Zhenyu Wang<EMAIL_ADDRESS>
###### Abstract

Magnetized disk winds and wind-driven accretion are an essential and
intensively studied dispersion mechanism of protoplanetary disks. However, the
stability of these mechanisms has yet to be adequately examined. This paper
employs semi-analytic linear perturbation theories under non-ideal
magnetohydrodynamics, focusing on disk models whose magnetic diffusivities
vary by a few orders of magnitude from the disk midplane to its surface.
Linear modes are distinguished by their symmetry with respect to the midplane.
These modes have qualitatively different growth rates: symmetric modes almost
always decay, while at least one anti-symmetric mode always has a positive
growth rate. This growth rate decreases faster than the Keplerian angular
velocity with cylindrical radius $R$ in the disk and scales as $R^{-2}$ in the
fiducial disk model. The growth of anti-symmetric modes breaks the reflection
symmetry across the disk equatorial plane, and may occur even in the absence
of the Hall effect. In the disk regions where fully developed anti-symmetric
modes occur, accretion flows appear only on one side of the disk, while disk
winds occur only on the other. This may explain the asymmetry of some observed
protoplanetary-disk outflows.

Magnetohydrodynamics (1964), Accretion (14), Stellar accretion disks (1579),
Protoplanetary disks (1300), Exoplanet formation (492)

##  1 Introduction

Protoplanetary disks (PPDs), which serve as the cradles for planet formation,
undergo a lifespan of approximately $10^{6}-10^{7}~{}\mathrm{yr}$. They follow
three distinct processes to disperse, namely: (1) creation of planets, (2)
accretion onto the central protostar, and (3) outflowing in winds. The latter
two processes directly compete with the first, impeding the time and mass
available to create planets. Compared to photoevaporative winds that are
mostly unable to drive accretion (Owen et al., 2012; Wang & Goodman, 2017a),
magnetized winds of PPD exert torque on the disks. Magnetic fields create the
linkage between disk winds and accretions and have recently been identified as
a crucial factor in PPDs (Bai & Stone, 2013; Bai, 2013; Bai et al., 2016; Bai,
2017; Wang et al., 2019). In contrast, other prospective mechanisms have
inadequate efficiency in viscous turbulent accretion via the magnetorotational
instability (MRI; e.g. Balbus & Hawley 1998) or any other hydrodynamic
instabilities under PPD conditions (e.g., Bai & Stone 2013; Bai 2013; Simon et
al. 2013b, a; see Turner et al. 2014 for a review).

In recent years, several studies have been conducted to investigate the
magnetized wind-driven accretion process in protoplanetary disks. In one such
study (Bai, 2017), global simulations were performed using 2.5-dimensional
axisymmetric full MHD simulations with non-ideal MHD effects. Magnetic
diffusivities were evaluated through a pre-calculated interpolation table, and
thermodynamics were calculated via a simple relaxation-time recipe that
accounted for temperature dependence on spatial location. Wang et al. (2019)
adopted non-equilibrium thermochemical networks that co-evolved with non-ideal
MHD to consistently determine all magnetic diffusivities and thermodynamic
properties in real time. Building on these works, Nemer et al. (2020)
performed calculations to predict observational evidence of magnetized disk
winds and the resulting predictions were subsequently confirmed in recent
observations. High-resolution observations of the [O i] $6300~{}\mathrm{\AA}$
emission line were conducted by Fang et al. (2023b), who confirmed that
magnetized disk winds, rather than photoevaporative winds, were essential in
explaining the spatially resolved emission line features. These studies
highlight the importance of non-ideal MHD effects and the role of magnetic
fields in understanding the complex dynamics of protoplanetary disks and their
associated winds. There have been, nevertheless, no direct measurements of
disk magnetic fields. Attempts up to now only yield the upper limits of field
strengths in a handful of disks (e.g., TW Hya, Vlemmings et al. 2019; AS 209,
Harrison et al. 2021).

Reflection symmetry over the equatorial plane has been commonly assumed in
research on magnetized winds in protoplanetary disks, mainly concentrating on
their general features, such as launching mechanisms, observables, and
kinematics. Until recently, asymmetries in magnetic winds have been largely
overlooked, with only a few observations indicating their existence, such as
in the case of HH 30 (Watson & Stapelfeldt, 2004). However, with advancements
in observational techniques, more observations are expected to emerge in
various wavelength bands. Understanding these asymmetries is crucial for
future observations that depict detailed characteristics of magnetized disk
winds. Bai & Stone (2017) discussed the Hall effect and attributed the
breaking of reflection symmetry to it. Nevertheless, some recent simulation
efforts discovered that the asymmetry may still emerge when the Hall effect is
intentionally turned off while growth rates seem to continue. This paper will
delve into the symmetry of wind-driven accreting protoplanetary disk systems,
paying attention to their physical origins involving magnetohydrodynamics and
thermochemistry and their potential observables.

This paper is structured as follows. §2 offers a detailed description of the
physical models and equations used to analyze the given system. Simplified
versions of the model are applied to idealized systems for qualitative
discussions in §3. The mathematical models are subsequently applied to the
numerical model cited in Nemer et al. (2020) to study its asymmetric
instabilities, elaborated in §4. Furthermore, in §5, we explore the parameter
space to understand how different physical assumptions impact the
instabilities studied. Finally, §6 provides discussions on possible
generalizations and observable tests and presents a summary of the principal
findings and conclusions. Some details of analytic derivations are presented
in the Appendices.

##  2 Method

The mathematical model of this work is based on the physical picture of wind-
driven accretion in PPDs. Unless otherwise specified, this paper uses
cylindrical coordinates ${R,\varphi,z}$ for describing everything associated
with the geometries, where the $z$-axis is aligned with the PPD’s axis of
rotation.

###  2.1 Non-ideal MHD for wind-driven accretion

This work adopts non-ideal MHD equations to describe the accretion and winds
of magnetized PPDs,

$\begin{split}&\partial_{t}\rho+\nabla\cdot(\rho\mathbf{v})=0,\\\
&\rho\partial_{t}\mathbf{v}+\rho\mathbf{v}\cdot\nabla\mathbf{v}=\frac{\mathbf{J}}{c}\times\mathbf{B}-\rho\nabla\Phi-\nabla
p,\\\ &\nabla\cdot\mathbf{B}=0,\
\nabla\times\mathbf{B}=\frac{4\pi}{c}\mathbf{J},\\\
&\partial_{t}\mathbf{B}=-c\nabla\times\mathbf{E},\
\nabla\cdot\mathbf{E}=0.\end{split}$ (1)

Here we use $\rho$ for the mass density, $\mathbf{v}$ the gas velocity ($c$
the vacuum speed of light), $p$ the gas pressure, $\Phi$ the gravitational
potential, $\mathbf{J}$ the electric current density, and $\mathbf{B}$ and
$\mathbf{E}$ the magnetic and electric fields. For non-ideal MHD systems with
finite conductivities, the $\mathbf{E}$ field is evaluated by the electric
field in the local fluid rest frame $\mathbf{E}^{\prime}$,

$\begin{split}\mathbf{E}^{\prime}&=\mathbf{E}+\frac{\mathbf{v}}{c}\times\mathbf{B}\\\
&=\frac{4\pi}{c^{2}}\left[\eta_{\mathrm{O}}\mathbf{J}+\eta_{\mathrm{H}}\frac{\mathbf{J}\times\mathbf{B}}{|\mathbf{B}|}+\eta_{\mathrm{A}}\frac{\mathbf{B}\times(\mathbf{J}\times\mathbf{B})}{|\mathbf{B}|^{2}}\right],\end{split}$
(2)

where $\eta_{\mathrm{O},\mathrm{H},\mathrm{A}}$ are used to denote the Ohmic
($\eta_{\mathrm{O}}$), Hall ($\eta_{\mathrm{H}}$), and ambipolar
($\eta_{\mathrm{A}}$) components of magnetic diffusivities (see also Wardle,
2007; Xu & Bai, 2016).

While previous studies like Wang et al. (2019) have already determined
diffusivity profiles through self-consistent calculations based on non-
equilibrium thermochemistry for magnetized wind-driven PPD simulations, this
paper aims to explore the potential impact of hypothetical diffusivity
profiles on the dynamical stability of magnetized PPDs. This approach
complements previous studies and contributes to the ongoing efforts to develop
more accurate and comprehensive models of the disks. Using hypothetical
diffusivity profiles allows for the exploration of a range of possible
scenarios and the identification of patterns or effects that may not have been
previously considered, deepening our understanding of the physics of PPDs and
informing future research efforts.

####  2.1.1 Radially local problems

In order to explore the breaking of reflection symmetry, it is necessary to
have steady-state, symmetric solutions as a foundation for perturbation
theories. Axisymmetric models local in cylindrical radius $R$ are used to
discuss reflection symmetry properly. In most radial ranges of typical PPDs,
the vertical density scale height $h$ is much smaller than $R$, which means
that the derivatives of physical variables with respect to $R$ are much
smaller than those with respect to $z$. It is assumed that the disk gas
temperature $T$ is only a function of $R$. This assumption is a reasonable
approximation within PPDs where $T$ is primarily controlled by thermal
accommodation with dust grains, which are in tight thermal equilibrium with
the radiation emanating from the central star (Chiang & Goldreich, 1997). It
may not hold very well above the disk surfaces, such as inside magnetized
winds (Wang et al., 2019), but our analysis focuses primarily on the behaviors
below disk surfaces. For more explicit discussions, we construct isothermal
models dependent only on $z$, with most radial derivatives $\partial_{R}$
considered negligible unless otherwise specified. Because of the axisymmetry,
all azimuthal derivatives ($\partial_{\varphi}$) vanish.

To better regularize and generalize the analyses in this paper, conversions to
dimensionless variables are necessary. For the radially local isothermal
models we discuss, we introduce the following conversions:

$\begin{split}&\zeta\equiv\frac{z}{h},\ \tau\equiv\Omega_{\rm\,K}t,\
\partial_{z}\rightarrow\frac{\mu_{k}}{R}\partial_{\zeta},\
\partial_{t}\rightarrow\Omega_{\rm\,K}\partial_{\tau};\\\
&\mu_{k}\equiv\frac{\Omega_{\rm\,K}R}{c_{s}},\
\mu_{\varphi}\equiv\frac{v_{\varphi}}{c_{s}}-\mu_{k},\
\mu_{R}\equiv\frac{v_{R}}{c_{s}},\ \mu_{z}\equiv\frac{v_{z}}{c_{s}};\\\
&\varrho\equiv\frac{\rho}{\rho_{0}},\ \beta_{0}\equiv=\frac{8\pi
c_{s}^{2}\rho_{0}}{B_{z0}^{2}},\ b_{i}\equiv\frac{B_{i}}{B_{z0}},\\\
&\varepsilon_{i}\equiv\frac{cE_{i}}{c_{s}B_{z0}},\ j_{i}\equiv\frac{4\pi
RJ_{i}}{cB_{z0}\mu_{k}},\ i\in\\{R,\varphi,z\\}.\end{split}$ (3)

Here $\Omega_{\rm\,K}$ is the angular velocity of the Keplerian orbital
motion, $\rho_{0}$ is the mid-plane mass density, $B_{z0}\equiv B_{z}|_{z=0}$
is the $z$-component of the magnetic fields at the equatorial plane, and the
isothermal sound speed satisfies $c_{s}^{2}=p/\rho$. The scale height is
adopted as $h=c_{s}/\Omega_{\rm\,K}$. The Mach number of Keplerian velocity
$\mu_{k}$ compares the vertical sound crossing timescales to the orbital
timescales. Its value can be estimated by,

$\begin{split}\mu_{k}\simeq&29\times\left(\frac{M_{*}}{M_{\odot}}\right)^{1/2}\left(\frac{T}{300~{}{\rm\,K}}\right)^{-1/2}\\\
&\quad\times\left(\frac{R}{\mathrm{AU}}\right)^{-1/2}\left(\frac{\langle
m_{\rm mol}\rangle}{2.35m_{p}}\right)^{1/2},\end{split}$ (4)

in which $M_{*}$ is the stellar mass, $m_{p}$ is the proton mass, and $\langle
m_{\rm mol}\rangle$ is the mean molecular mass. It can be easily verified that
almost all radial derivatives in the dimensionless forms are multiplied by
$\mu_{k}^{-1}$, meaning that their effects on the equations are suppressed by
more than one order of magnitude. The most important exception is
$\partial_{R}E_{z}$, in which the term $\partial_{R}(v_{\varphi}B_{R})$
reflects the radial velocity shear, and leads to the $(-3b_{R}/2)$ term in the
$\partial_{\tau}b_{\varphi}$ expression without the $\mu_{k}^{-1}$ suppression
(see also Wardle & Königl, 1993). One direct consequence of this approximation
is that the solenoidal condition of magnetic fields
$\nabla\cdot\mathbf{B}=R^{-1}\partial_{R}(RB_{R})+\partial_{z}B_{z}=0$ reduces
to $\partial_{z}B_{z}=0$, and hence $b_{z}=1$ always holds. Simulations of
magnetized PPD winds have confirmed this assumption with high confidence (e.g.
Bai, 2017).

Under such transforms, the MHD differential eqs. (1) and (2) are reduced into
dimensionless differential equations for the axisymmetric, radially local,
isothermal system,

$\begin{split}&\partial_{\tau}\mu_{z}+\partial_{\zeta}\ln\varrho+\mu_{z}\partial_{\zeta}\mu_{z}=\frac{2}{\beta_{0}\varrho}(j_{R}b_{\varphi}-j_{\varphi}b_{R})-\zeta,\\\
&\partial_{\tau}\mu_{R}+\mu_{z}\partial_{\zeta}\mu_{R}=\frac{2}{\beta_{0}\varrho}(j_{\varphi}b_{z}-j_{z}b_{\varphi})+2\mu_{\varphi}+g_{R},\\\
&\partial_{\tau}\mu_{\varphi}+\mu_{z}\partial_{\zeta}\mu_{\varphi}=\frac{2}{\beta_{0}\varrho}(j_{z}b_{R}-j_{R}b_{z})-\frac{\mu_{R}}{2},\\\
&\partial_{\tau}b_{\varphi}=-\partial_{\zeta}\varepsilon_{R}-\frac{3}{2}b_{R},\quad\partial_{\tau}b_{R}=\partial_{\zeta}\varepsilon_{\varphi}.\\\
\end{split}$ (5)

Here $g_{R}$ is a correction term for reducing radial gravitational force at
high altitudes (see Appendix A). The dimensionless electric current densities
$j$ and fluid frame electric fields $\varepsilon^{\prime}$ are related to the
reduced fields $b$ and the rest frame electric fields $\varepsilon$ by,

$\begin{split}&j_{R}=-\partial_{\zeta}b_{\varphi},\
j_{\varphi}=\partial_{\zeta}b_{R},\\\
&\varepsilon^{\prime}_{R}=\varepsilon_{R}+(\mu_{\varphi}+\mu_{k})b_{z}-\mu_{z}b_{\varphi},\\\
&\varepsilon^{\prime}_{\varphi}=\varepsilon_{\varphi}+\mu_{z}b_{R}-\mu_{R}b_{z},\
\varepsilon^{\prime}_{z}\simeq 0.\end{split}$ (6)

The value of $j_{z}$ is generally considered tiny, yet in some conditions, its
contribution is not negligible; Appendix A also estimates the $j_{z}$ values.
Similar to the dimensional case, the dimensionless $\varepsilon^{\prime}$ and
$j$ are related by the magnetic diffusivities,

$\begin{split}&\begin{bmatrix}\varepsilon^{\prime}_{R}\\\
\varepsilon_{\varphi}^{\prime}\end{bmatrix}=\begin{bmatrix}\alpha_{RR}&\alpha_{R\varphi}\\\
\alpha_{\varphi R}&\alpha_{\varphi\varphi}\end{bmatrix}\begin{bmatrix}j_{R}\\\
j_{\varphi}\end{bmatrix},\\\
&\alpha_{RR}\equiv\alpha_{\mathrm{O}}+\frac{\alpha_{\mathrm{A}}}{b^{2}}(b_{\varphi}^{2}+1),\
\alpha_{R\varphi}\equiv\frac{\alpha_{\mathrm{H}}}{b}-\frac{\alpha_{\mathrm{A}}}{b^{2}}b_{R}b_{\varphi},\
\\\ &\alpha_{\varphi
R}\equiv-\frac{\alpha_{\mathrm{H}}}{b}-\frac{\alpha_{\mathrm{A}}}{b^{2}}b_{R}b_{\varphi},\
\alpha_{\varphi\varphi}\equiv\alpha_{\mathrm{O}}+\frac{\alpha_{\mathrm{A}}}{b^{2}}(b_{R}^{2}+1)\
;\\\
&\alpha_{i}\equiv\frac{\eta_{i}}{hc_{s}}\equiv\frac{2b^{2}}{\beta_{0}\varrho}\Lambda_{i}^{-1},\
i\in\\{\mathrm{O},\mathrm{H},\mathrm{A}\\}.\end{split}$ (7)

Here $b^{2}\equiv b_{R}^{2}+b_{\varphi}^{2}+1$, and we also introduce
$\Lambda_{i}\equiv B^{2}/(4\pi\rho\eta_{i}\Omega_{\rm\,K})$ for the Elsasser
numbers of the three components. It is already known that the
$\eta_{\mathrm{O}}$, $\eta_{\mathrm{H}}/|\mathbf{B}|$, and
$\eta_{\mathrm{A}}/|\mathbf{B}|^{2}$ are mostly independent of $\mathbf{B}$
(e.g. Xu & Bai, 2016). This fact implies that we should adopt
$\alpha_{\mathrm{O}}$, $\alpha_{\mathrm{H}}/b$, and
$\alpha_{\mathrm{A}}/b^{2}$ for the diffusivity profiles in practice.

####  2.1.2 Steady states

To solve eqs. (5) for the steady states, it is necessary to set
$\partial_{\tau}\rightarrow 0$. Additionally, a few extra approximations can
be made to simplify the problem further. In the accretion layer, vertical gas
motion $\mu_{z}$ is significantly smaller than the horizontal motion $\mu_{R}$
and $\mu_{\varphi}$. The smallness of $\mu_{z}$ makes the steady-state
versions of eqs. (5) stiff via the operator $\mu_{z}\partial_{z}$; yet this
operator is unimportant for vertically smooth solutions. Therefore, setting
$\mu_{z}$ to zero for the system is practical. The transport of magnetic
fluxes mainly affects $\varepsilon_{\varphi}$, which is related to radial
transport, that is, $\partial_{t}B_{z}$. For a radially local solution, it is
safe to assume that $\varepsilon_{\varphi}$ is zero, allowing us to
concentrate on the vertical modes. Eq. (5) are simplified by these assumptions
into the ordinary differential equations (ODEs),

$\begin{split}&\partial_{\zeta}\ln\varrho=\dfrac{2}{\beta_{0}\varrho}(j_{R}b_{\varphi}-j_{\varphi}b_{R})-\zeta,\
\partial_{\zeta}\varepsilon_{R}=-\dfrac{3}{2}b_{R},\\\
&\partial_{\zeta}b_{R}=j_{\varphi},\
\partial_{\zeta}b_{\varphi}=-j_{R},\end{split}$ (8)

where the dimensionless current densities $j_{R,\varphi}$ and the subsequent
velocity components are determined by solving the combination of eqs. (6) and
(7) (see also Appendix B).

The symmetry required for the steady state solutions leads to
$b_{R}=b_{\varphi}=0$ at $\zeta=0$, and $\varrho|_{\zeta=0}=1$ by definition.
The free parameter to be determined is
$\varepsilon_{R0}\equiv\varepsilon_{R}|_{\zeta=0}$, whose value is determined
by matching physical parameters of wind solutions constructed with ideal MHD.
Given proper profiles of diffusivities, eqs. (8) are integrated to an altitude
at which the diffusivity is sufficiently low, characterized by a critical
dimensionless Ohmic diffusivity $\alpha_{\mathrm{O},c}$. In practice, we find
that $\alpha_{\mathrm{O},c}=10^{-4}$ appears to be a good cutoff, around which
small variations of $\alpha_{\mathrm{O}c}$ will not change the results
significantly. At this altitude, an ideal MHD wind solution is generated using
the Bai et al. (2016) scheme, whose Bernoulli parameter should match the value
given by the steady-state “disk” solution described here. We refer the readers
to Appendix C for more details.

###  2.2 Linear perturbations

The stability of steady-state solutions §2.1.2 should be analyzed by
perturbing eqs. (5). For any relevant dependent variables $x$, we decompose it
into the form $x\rightarrow x+\delta x$, where the $x$ now stands for the
steady-state solution, and $\delta x$ for the perturbation. This decomposition
is applied to eqs. (5), (6) and (7). The perturbations of variables are
preserved only up to the first order. Similar to §2.1.2, the perturbations of
vertical velocities also vanish ($\delta\mu_{z}\rightarrow 0$) as the focus of
our stability analyses is located below the disk surface.

To analyze the time evolution of the perturbations, we assume that all
perturbation terms have the same time dependence, $\delta
x\propto\mathrm{e}^{\nu\tau}$, where $\nu$ is the dimensionless growth rate.
The sign of $\nu$ determines whether a perturbative mode is stable ($\nu\leq
0$) or unstable ($\nu>0$), and the absolute value of $\nu$ indicates how fast
the mode grows or decays. The momentum equations for the perturbations then
read,

$\nu\delta\mu_{R}-2\delta\mu_{\varphi}=\frac{2\delta
j_{\varphi}}{\beta_{0}\varrho},\
\frac{\delta\mu_{R}}{2}+\nu\delta\mu_{\varphi}=-\frac{2\delta
j_{R}}{\beta_{0}\varrho}.$ (9)

The derivatives of the perturbed velocities are then,

$\begin{split}&\partial_{\zeta}\delta\mu_{R}=\zeta\delta\mu_{R}-\left[\frac{2}{\beta_{0}\varrho(1+\nu^{2})}\right](2\partial_{\zeta}\delta
j_{R}-\nu\partial_{\zeta}\delta j_{\varphi}),\\\
&\partial_{\zeta}\delta\mu_{\varphi}=\zeta\delta\mu_{\varphi}-\left[\frac{2}{\beta_{0}\varrho(1+\nu^{2})}\right]\left(\nu\partial_{\zeta}\delta
j_{R}+\frac{\partial_{\zeta}\delta j_{\varphi}}{2}\right).\end{split}$ (10)

To further simplify the analyses, we make two extra assumptions: (1) that
$\delta(\ln\varrho)=0$ (since $\mu_{z}\rightarrow 0$ and $\delta\mu_{z}=0$) as
the linear perturbations are insufficient to reshape the vertical density
profile, and (2) that the matrix elements $\alpha_{ij}$
($i,j\in\\{R,\varphi\\}$) in eqs. (7) are not susceptible to the perturbations
either. Taking the perturbations into the field equations in eqs. (5), (6) and
(7), one obtains the governing equations of the perturbations, the matrix
form,

$\partial_{\zeta}[(A+\xi V_{1})\delta b]=V_{0}\delta b,$ (11)

where $\delta b\equiv[\delta b_{R},\delta b_{\varphi},\delta j_{R},\delta
j_{\varphi}]^{T}$ is the vector of perturbation variables, and the other
involved matrices are defined as,

$\begin{split}&A\equiv\begin{bmatrix}&&\alpha_{\varphi
R}-2\xi&\alpha_{\varphi\varphi}\\\ &&-\alpha_{RR}&-\alpha_{R\varphi}-\xi/2\\\
&-1&&\\\ 1&&&\end{bmatrix},\\\ &V_{0}\equiv\begin{bmatrix}\nu&&&\\\
3/2&\nu&&\\\ &&1&\\\ &&&1\end{bmatrix},\quad
V_{1}\equiv\begin{bmatrix}&&&\nu\\\ &&-\nu&\\\ &0&&\\\
0&&&\end{bmatrix},\end{split}$ (12)

in which $\xi\equiv 2/[\beta_{0}\varrho(1+\nu^{2})]$, and we have used the
approximation $\partial_{\zeta}\xi=\zeta\xi$ which holds inside the concerned
regions (below the wind base). Eqs. (11) shall be integrated from $\zeta=0$ to
the wind bases with eigenvalue $\nu$. Two types of symmetries are possible for
$b_{R,\varphi}$:

  * •

Symmetry: $\delta b_{i}=0$ but $\partial_{\zeta}\delta b_{i}\neq 0$ at
$\zeta=0$;

  * •

Anti-symmetry: $\delta b_{i}\neq 0$ but $\partial_{\zeta}\delta b_{i}=0$ at
$\zeta=0$.

Note that, similar to §2.1, the definitions of “symmetry” and “anti-symmetry”
are for the field line morphologies, viz. the integral curves of $\delta b$
(see Figure 1). Because $B_{z}$ should not change its sign over the equatorial
plane, the symmetry in $b_{R,\varphi}$ leads to the anti-symmetry in field
line morphologies and vice versa.

Figure 1: Schematic illustration (not to scale) of symmetric (upper panel) and
anti-symmetric (lower panel) perturbation modes. Note that the symmetries are
indicated in terms of field line geometries (which should refer to
$\int\mathrm{d}\zeta\;\delta b$), not the values of $\delta b$ (§2.2).

For symmetric modes, we set $\delta b_{R0}=\delta b_{\varphi 0}=0$, $\delta
j_{R0}=1$, but $\delta j_{\varphi 0}$ as the other free parameter (the extra
subscripts “0” here denote the $\zeta=0$ values). For anti-symmetric modes, we
fix $\delta j_{R0}=\delta j_{\varphi 0}=0$, $\delta b_{R0}=1$, but $\delta
b_{\varphi 0}$ as the other paramter. When obtaining these type of modes,
parameters $\nu$ and $\delta j_{\varphi 0}$ (symmetric) and $\delta b_{\varphi
0}$ are adjusted, so that $\delta b_{R}|_{\zeta=\zeta_{1}}=\delta
b_{\varphi}|_{\zeta=\zeta_{1}}=0$. The upper limit of integration $\zeta_{1}$
is located at where the steady-state $\mu_{R}$ changes its sign. This approach
of setting up boundary conditions constrains the perturbations that they do
not introduce shears in magnetic fields at $\zeta_{1}$, to avoid electric
current sheets with infinitesimal thickness causing discontinuities. It also
yields zero disturbances to the Maxwell stresses at the wind bases so that the
steady-state wind solutions above the disk surfaces are unaffected.

##  3 Stabilities of Simplified Models

This section applies the methods in §2 to simplified systems, which verify
that our theories comply with existing studies, and also help us to acquire
acquaintances with the behavior of the perturbations in PPDs.

###  3.1 Local ideal MHD limit

The perturbation theories derived in §2.2 quantify the finite conductivities
of plasmas with non-zero dimensionless diffusivity $\alpha$, and can thus be
restricted to the ideal MHD limit by taking $\alpha\rightarrow 0$, yielding,

$\begin{split}&\partial_{\tau}\delta\mu_{R}=\frac{2\partial_{\zeta}\delta
b_{R}}{\beta_{0}\varrho}+2\delta\mu_{\varphi},\
\partial_{\tau}\delta\mu_{\varphi}=\frac{2\partial_{\zeta}\delta
b_{\varphi}}{\beta_{0}\varrho}-\frac{\delta\mu_{R}}{2},\\\
&\partial_{\tau}\delta b_{R}=\partial_{\zeta}\delta\mu_{R},\
\partial_{\tau}\delta
b_{\varphi}=\partial_{\zeta}\delta\mu_{\varphi}-\frac{3}{2}\delta
b_{R}.\end{split}$ (13)

At any specific altitude, $\beta_{0}\varrho$ is fixed. One usually studies the
local behaviors by assuming that the coefficients can be treated as constants
and Fourier transforms are applicable, i.e., $\delta
x\propto\mathrm{e}^{\mathrm{i}(k\zeta-\omega t)}$, and the dispersion relation
can be written as,

$\begin{bmatrix}\mathrm{i}\omega&&\mathrm{i}k&\\\
-3/2&\mathrm{i}\omega&&\mathrm{i}k\\\
2\mathrm{i}k/(\beta_{0}\varrho)&&\mathrm{i}\omega&2\\\
&2\mathrm{i}k/(\beta_{0}\varrho)k&-1/2&\mathrm{i}\omega\end{bmatrix}\begin{bmatrix}\delta
b_{R}\\\ \delta b_{\varphi}\\\ \delta\mu_{R}\\\
\delta\mu_{\varphi}\end{bmatrix}=0.$ (14)

This equation requires a vanishing determinant for non-trivial solutions.
Defining $\tilde{k}^{2}\equiv 2k^{2}/(\beta_{0}\varrho)$, the dispersion
relation reads,

$\omega^{4}-(2\tilde{k}^{2}+1)\omega^{2}-3\tilde{k}^{2}+\tilde{k}^{4}=0.$ (15)

Since $\omega^{2}$ is always real, instability arises from positive imaginary
parts of $\omega$, which requires

$\omega^{2}<0\ \Rightarrow\ 2\tilde{k}^{2}+1<(16\tilde{k}^{2}+1)^{1/2}\ ,$
(16)

or $\tilde{k}<\sqrt{3}$. This result is identical to the MRI dispersion
relation assuming isothermal, radially local (ignoring all radial derivatives
and wavenumbers) conditions (e.g. Balbus & Hawley, 1991).

###  3.2 Constant diffusivities

Figure 2: Local-mode growth rates ($\nu\equiv\mathrm{max}\\{{\rm
Im}[\omega]\\}$) assuming constant diffusivities, showing the distribution in
the space spanned by the modulus of dimensionless wavenumber $|k|$ and the
dimensionless altitude $\zeta$ (§3.2). The colormaps are in symmetric
logarithmic scales: red for positive $\nu$, and blue for negative $\nu$,
separated by solid black curves indicating the $\nu=0$ contours. The
dimensionless diffusivity $\alpha_{\mathrm{O}}$ and the plasma $\beta$ of
models are presented to the upper-left of each panel.

A more meaningful simplification is to conduct calculations with constant,
non-zero diffusivities. Similar to the ideal MHD case, we insert the
perturbations into eq. (11) using the vertically local approximation and the
Fourier transform scheme. Assuming that all dimensionless diffusivities
($\alpha_{\mathrm{O},\mathrm{H},\mathrm{A}}$) are constants of $\zeta$, the
equations reduce to

$\begin{split}&\begin{bmatrix}\mathrm{i}\omega&&a_{\varphi
R}&a_{\varphi\varphi}\\\ -3/2&\mathrm{i}\omega&a_{RR}&a_{R\varphi}\\\
&\mathrm{i}k&1&\\\ \mathrm{i}k&&&-1\end{bmatrix}\begin{bmatrix}\delta b_{R}\\\
\delta b_{\varphi}\\\ \delta j_{R}\\\ \delta j_{\varphi}\end{bmatrix}=0\ ;\\\
&a_{\varphi R}\equiv\mathrm{i}k\alpha_{\varphi R}-2\xi(\zeta+\mathrm{i}k),\\\
&a_{\varphi\varphi}\equiv\mathrm{i}k\alpha_{\varphi\varphi}+\xi\omega(k-\zeta),\\\
&a_{RR}\equiv-\mathrm{i}k\alpha_{RR}-\xi\omega(k-\zeta),\\\
&a_{R\varphi}\equiv-\mathrm{i}k\alpha_{R\varphi}-\xi(\zeta+\mathrm{i}k)/2.\end{split}$
(17)

Zero determinant of the matrix yields the dispersion relations. It is
straightforward to verify that the dispersion relation reduces to eq. (15)
once $\alpha_{\mathrm{O},\mathrm{H},\mathrm{A}}\rightarrow 0$ and
$\zeta\rightarrow 0$. An extra fact allows us to simplify the analyses
further, that one shall have $\alpha_{\mathrm{O}}\gg\alpha_{\mathrm{H}}$ and
$\alpha_{\mathrm{O}}\gg\alpha_{\mathrm{A}}$ inside most regions below the disk
surfaces (e.g. Xu & Bai, 2016), leading to,

$\alpha_{RR}\sim\alpha_{\varphi\varphi}\sim\alpha_{\mathrm{O}},\
\\{\alpha_{RR},\alpha_{\varphi\varphi}\\}\gg\\{\alpha_{\varphi
R},\alpha_{R\varphi}\\}.$ (18)

We use $\alpha_{\mathrm{O}}$ and $\beta_{0}$ as the parameters of models.
Similar to §3.1, Figure 2 plots $\nu\equiv\mathrm{max}\\{{\rm Im}[\omega]\\}$
in the space spanned by the modulus of wavenumber $|k|$ and $\zeta$ for models
with different $\alpha_{\mathrm{O}}$ and $\beta_{0}$. It is understood that
$k\lesssim 10^{0}$ is unphysical in disks since the surfaces of a typical PPD
is usually 2 to 3 scale heights from the equatorial plane, yet we still
analytically prolongate our solutions to $|k|=10^{-1}$ for completeness.

In Figure 2 we observe that, at each altitude $\zeta\lesssim 2$, there is a
critical $k_{\mathrm{crit}}$ that splits $\nu>0$ and $\nu<0$. The location of
this $k_{\mathrm{crit}}$ is mostly independent of $\beta_{0}$, and even
insensitive to $\zeta$, but is related to $\alpha_{\mathrm{O}}$. Numerically
we fit and identify that the critical $k$ roughly follows a simple power-law,
$k_{\mathrm{crit}}\simeq 3\alpha_{\mathrm{O}}^{-1/2}/2$. When this number is
comparable to unity, the unstable modes can grow, as the modes with vertical
wavenumbers smaller than $2\pi/h$ do not exist. This criterion about the
growth of instabilities agrees with the conclusions in e.g., Turner et al.
(2007) and Ilgner & Nelson (2008) semi-quantitatively, although it is
noteworthy that this result relies on the constant-diffusivity assumption.

Figure 3: Non-ideal MHD profiles for thei fiducial model (§4.1) at
$(R/\mathrm{AU})\in\\{0.5,1,2\\}$ presented in different columns. Top row: the
dimensionless diffusivities ($\alpha_{\mathrm{O}}$, $\alpha_{\mathrm{H}}/b$,
and $\alpha_{\mathrm{A}}/b^{2}$) and Elsasser numbers ($\Lambda_{\mathrm{O}}$,
$\Lambda_{\mathrm{H}}$, ${\rm Am}$). Bottom row compares the key MHD profiles
($\varrho$, $b_{R}$, $b_{\varphi}$, $\mu_{R}$) of the semi-analytic solutions
to the simulation results [marked with “(Sim.)”] used in Nemer et al. (2020),
Fang et al. (2023b). Note that in the bottow row, different physical
quantities are distinguished by colors, while dashed and solid line shapes
indicate negative and positive values, respectively.

##  4 Fiducial Model: Magnetic Diffusivites Calculated by Simulations

In realistic magnetized protoplanetary disks (PPDs) models, the diffusivity
profiles are vertically stratified, with a dynamical range of several orders
of magnitude. In contrast to more straightforward local modes discussed in
section §3, these profiles require solving eigenvalue problems (as per eq. 11)
using semi-analytic methods described in section §2.2. This section’s
calculations rely on the non-ideal MHD profiles observed in the fiducial
simulation adopted by Nemer et al. (2020) and Fang et al. (2023b).

###  4.1 Comparisons of steady-state solutions

In order to verify the validity of our semi-analytic approach, we compare the
simulation profiles with the steady-state solutions yielded by the scheme
described in §2.1.2, taking the vertical distributions of magnetic
diffusivities at different radii as the input. It is noticed that the
$|b_{R}|\ll 1$ near the mid-plane, where small variations may not affect the
overall behavior of the solution, yet will still change the apparent profiles
in the comparisons. Here we notice that, if we define
$\chi_{z}\equiv\partial_{\ln R}\ln B_{R}$, then,

$\partial_{\zeta}b_{R}=j_{\varphi}+\frac{\chi_{z}}{\mu_{k}},$ (19)

and the $\chi_{z}/\mu_{k}$ term will be rather important near the equatorial
plane where $|j_{\varphi}|\simeq|\partial_{\zeta}b_{R}|\ll 1$ and $|b_{R}|\ll
1$. If there is no radial transport of magnetic fluxes, $\chi_{z}$ should be a
constant $\chi_{z}=-5/4$ at all radii for the power-law setup of the fiducial
simulation (Wang et al., 2019; Nemer et al., 2020), yet during the evolution
of the simulation run, radial transport of fluxes makes $\chi_{z}$ vary. This
subsection fits $\chi_{z}$ to be $\chi_{z}=0.2$ for $R=0.5~{}\mathrm{AU}$, and
$\chi_{z}=-1$ for $R=1~{}\mathrm{AU}$ and $2~{}\mathrm{AU}$, so that they
yield semi-analytic solutions that roughly match the mid-plane $b_{R}$
profiles of the simulation.

Figure 3 presents the comparisons to the simulation results at
$(R/\mathrm{AU})\in\\{0.5,1,2\\}$, confirming that our semi-analytic approach
indeed yields consistent results. The steady-state solutions should be
sufficiently regular at different radii when used as the foundation of
perturbation theories. At the same time, our numerical tests have found that
varying $\chi_{z}$ will not have visible impacts on the eigenvalues.

Figure 4: Eigensolutions at disk radii $(R/\mathrm{AU})\in\\{0.5,1,2\\}$
presented in different columnes (see also §4.2.1). The top row shows the phase
diagrams for the symmetric perturbation modes in the space spanned by the
dimensionless $\\{\nu_{\rm sym}\\}\times\\{\delta j_{\varphi 0}/\delta
j_{R0}\\}$. The solid and dashed contours indicate the tracks of $\delta
b_{R}|_{\zeta_{1}}=0$ and $\delta b_{\varphi}|_{\zeta_{1}}=0$ ($\zeta_{1}$ is
the altitude of the wind base), respectively, and their intersections indicate
eigensolutions. Note that the negative-$\nu$ and positive-$\nu$ regions,
separated by a vertical solid line, have different scales in $\nu$. The middle
row is similar to the top row but shows the phase diagrams for the anti-
symmetric modes in the $\\{\nu_{\rm asy}\\}\times\\{\delta b_{R0}/\delta
b_{\varphi}0\\}$ phase space. The bottom row presents the eigensolutions with
the largest growth rates for both modes. Note that the symmetric modes always
have negative growth rates.

###  4.2 Eigenmodes of perturbations

Once the steady-state solutions have been obtained, the growth rates of
perturbation modes can be calculated. For the symmetric modes, adjusting two
parameters is necessary to make $\delta_{R,\varphi}$ vanish simultaneously at
$\zeta=\zeta_{1}$, the top of the accretion layer. Because the perturbative
problems are linear, one can always set $\delta
j_{R0}=-[\partial_{\zeta}\delta b_{\varphi}]_{\zeta=0}=1$, and treat $\delta
j_{\varphi 0}=[\partial_{\zeta}\delta b_{R}]_{\zeta=0}$ as one of the
parameters. For the anti-symmetric modes, one sets $\delta b_{R0}=1$, and
$\delta b_{\varphi 0}$ is the parameter. The other parameter in both cases is
always the eigenvalue $\nu$.

####  4.2.1 Eigenmodes at various disk radii

We elaborate the calculations for $R/\mathrm{AU}\in\\{0.5,1,2\\}$ as examples.
Plotting the track of the zero points of $\delta b_{R1}\equiv\delta
b_{R}|_{\zeta_{1}}$ and $\delta b_{\varphi 1}\equiv\delta
b_{\varphi}|_{\zeta_{1}}$ in the phase space spanned by
$\\{\nu\\}\otimes\\{\delta j_{\varphi 0}/\delta j_{R0}\\}$ (symmetric modes)
or $\\{\nu\\}\otimes\\{\delta b_{R0}/\delta b_{\varphi 0}\\}$ (anti-symmetric
modes), the proper eigen modes should locate at the intersections of the
“root” curves for $\delta b_{R1}$ and $\delta b_{\varphi 1}$ (Figure 4).

The first attempts to obtain eigenmodes intentionally set
$\alpha_{\mathrm{H}}=0$ as the fiducial cases. For both types of modes, we
observe that the modes with $|\nu|\lesssim 10^{-1}$ have the least number of
nodes; modes with higher spatial frequencies all lie in the $\nu<0$ half-
space. In the $\nu>1$ domain, the traces of roots seem to converge, but more
detailed examinations confirm that they never touch each other.

At this radius, the maximum $\nu$ for symmetric modes is $\nu\simeq-0.36<0$
for $R=1~{}\mathrm{AU}$. Recovery of dimensions finds that this mode decays at
$\sim 0.4~{}\mathrm{yr}$ per e-fold. This value is tiny compared to the disk
lifetime ($\gtrsim 10^{6}~{}\mathrm{yr}$), let alone other symmetric modes
with higher spatial frequency and more negative $\nu$. Therefore, no
instability will occur through the symmetric modes. In contrast, an anti-
symmetric mode with $\nu\simeq 0.0349>0$ exists. Amplitudes of this growing
mode will increase by $\sim 4.5~{}\mathrm{yr}$ per e-fold and should leave the
linear stage rather quickly.

It is noteworthy that the Hall effect does not efficiently manipulate this
instability. Using the $R=1~{}\mathrm{AU}$ vertical Hall diffusivity profile
($\alpha_{\mathrm{H}}/b$) obtained during the simulation and converting it
into the dimensionless Hall diffusivity, we test the situations that $B_{z}$
is parallel and anti-parallel to the axis of disk rotation, respectively. The
resulting eigenvalues are affected by a tiny fraction ($\nu=0.0355$ for
parallel fields and $\nu=0.0348$ for anti-parallel fields). This consequence
seems to be in odd to some previous studies concluding that the Hall effect
could lead to significant asymmetries across the equatorial plane (e.g. Bai &
Stone, 2017; Bai, 2017). We nevertheless point out that the seeming oddity
comes from the nature of the spontaneous breaking of symmetries: the Hall
effect could be the “first push” that determines the direction, and the
subsequent amplification of asymmetry is dominated by the instabilities
described in this section.

####  4.2.2 Scalings of eigenvalues

Figure 5: Dimensional growth rates of anti-symmetric ($\tilde{\nu}_{\rm asy}$)
and symmetric ($\tilde{\nu}_{\rm sym}$) modes (dashed line for negative
values), for the fiducial model through the range of disk radius
$0.4<(R/\mathrm{AU})<4$. Note that $\tilde{\nu}_{\rm sym}<0$ holds for all
radii. The lines for power-laws $R^{-3/2}$ and $R^{-2}$ are also presented for
comparisons.

When the scheme described in §4.2.1 is applied to all relevant radii in the
disk model in §4.1, stability conditions of the disk can be studied in terms
of the reflection symmetry. To simplify the problem, we take the
representative $\chi_{z}\simeq-1/2$ in what follows to produce the steady-
state solutions throughout the disk.

After multiplying by the local Keplerian angular velocity to recover the
dimensions, we present the dimensional growth rates as functions of disk
radius in Figure 5–$\tilde{\nu}_{\rm sym}$ for the symmetric mode, and
$\tilde{\nu}_{\rm asy}$ for the anti-symmetric mode. It is evident that
throughout all radii in the model $\tilde{\nu}_{\rm sym}<0$ but
$\tilde{\nu}_{\rm asy}>0$, thus symmetric modes always decay, but at least one
of the anti-symmetric modes grows. In other words, the reflective symmetric
morphology of the wind-driven accretion model is unstable, and the instability
always initiates in an anti-symmetric shape.

For both $\tilde{\nu}_{\rm asy}$ and the $|\tilde{\nu}_{\rm sym}|$ at
$R\gtrsim 0.6~{}\mathrm{AU}$, one can roughly read their scaling from Figure
5,

$\dfrac{|\tilde{\nu}_{\rm sym}|}{\mathrm{yr}^{-1}}\sim
2\times\left(\dfrac{R}{\mathrm{AU}}\right)^{-2},\ \ \dfrac{\tilde{\nu}_{\rm
asy}}{\mathrm{yr}^{-1}}\sim
0.2\times\left(\dfrac{R}{\mathrm{AU}}\right)^{-2}.$ (20)

The growth and decay rates decrease with radius faster than $\Omega_{\rm\,K}$.
For larger radii, the growth of anti-symmetric modes could be slow. If the
scaling relation in eq. (20) can be generalized to all concerned radii in a
typical PPD, then $\tilde{\nu}_{\rm asy}^{-1}$ could be as slow as $\sim
1.5\times 10^{5}~{}\mathrm{yr}$ at $R\sim 200~{}\mathrm{AU}$, presumably the
“typical” outer radius of PPDs: the anti-symmetric instability will still take
place and grow through the $\sim 10^{6}~{}\mathrm{yr}$ lifetime of a PPD. The
growth is challenging to identify in numerical simulations available at such a
large distance, and hence, the necessity of the analytic approach is
emphasized.

###  4.3 Evolution of the anti-symmetric modes

For the breaking of reflection symmetry, we have identified that a possible
evolution manner is a quasi-steady state, whose most prominent characteristic
is the asymmetric wind-driven accretion flow. To construct the morphologies
analytically, we slightly modify the methods in §2.1.2 by setting non-
vanishing $b_{R}$ and $b_{\varphi}$ at the mid-plane and integrating to both
$+\zeta$ and $-\zeta$ directions. When obtaining such models, the matching
onto ideal MHD wind solutions is carried out only on the wind-launching side,
and the other side is left free. This scheme will leave the mid-plane $b_{R0}$
and $b_{\varphi 0}$ unconstrained. These two degrees of freedom reserve the
tracks of growth for the anti-symmetric modes, allowing a selection of
$(b_{R0},b_{\varphi 0})$ reflects a snapshot in the path of growth before
saturation. In general, the evolution of instabilities could be a complicated
issue that requires excessive numerical experiments to understand, and we
leave these studies to future works.

Figure 6 illustrates the MHD profiles obtained at $R=1~{}\mathrm{AU}$, in
which we set $b_{\varphi}|_{\zeta=0}=-14$, and the mid-plane value of $b_{R}$
is fixed at $b_{R}|_{\zeta=0}=0.45$. In this model, the accretion mainly
occurs on the $\zeta<0$ side, while the disk wind is launched efficiently only
on the $\zeta>0$ side. Over the disk surface on the $\zeta<0$ side,
suppression of gas density is observed, and $\mu_{R}$ is always negative. Such
a combination of fluid parameters indicates an absence of wind over that
surface.

The one-sided wind launching feature implies that the disk wind will exhibit
an apparent “asymmetric” shape. Such phenomenon has been observed in various
non-ideal MHD simulations of PPDs, in which the asymmetries developed to
different extents (e.g. Bai 2017; Hu et al. 2019; Riols et al. 2020; Hu et al.
2022; X. Hu et al., in prep.). Calculations in this work have confirmed that
these asymmetric morphologies should be physically plausible, and we should
also expect to discover such asymmetries observationally (see discussions in
§6.2).

Figure 6:  Similar to the lower row of Figure 3 but showing a non-symmetric
solution using the diffusivity profiles at $R=1~{}\mathrm{AU}$ as an example.
Colors distinguish different physical quantities, while the line shapes
indicate the signs (solid for positive, dashed for negative).

##  5 Stability in Various MHD Models

The wind-driven accretion problem described in §ef sec:fid-model is a typical
PPD model, yet countless situations still need to be discussed. The diversity
of PPDs, including their mass distributions, host star properties, and
thermochemical conditions, can all affect the diffusivity profiles
significantly. To study those parameters, this section explores different disk
conditions, including gravitation, magnetization, and diffusivity profiles. We
will manipulate one parameter at a time to simplify and clarify the
elaborations.

###  5.1 Functional forms of the vertical diffusivity profiles

For typical PPDs, the vertical stratification of diffusivities can be
qualitatively divided into three regions, (1) the weak coupling mid-plane
($\alpha_{\mathrm{O}}\sim 10^{1}$) in the $\zeta\lesssim 2$ region, (2) the
transition layer near the disk surface
($10^{1}\gtrsim\alpha_{\mathrm{O}}\gtrsim 10^{-4}$, $2\lesssim\zeta/h\lesssim
4$), and (3) the highly ionized wind above. As we have verified that the Hall
effect is more critical as a “first-push” than an amplification mechanism, the
following discussions will focus on the Ohmic and ambipolar diffusivities. We
choose to parameterize $\alpha_{\mathrm{O}}$ and $\alpha_{\mathrm{A}}/b^{2}$,
as they are insensitive to magnetic field strengths. Simple log-linear
functions with caps parameterize the vertical profiles of diffusivities
($\zeta^{\prime}\equiv\mathrm{max}\\{(\zeta-\zeta_{t}),\ 0\\}$),

$\begin{split}&\log_{10}\alpha_{\mathrm{O}}=\log_{10}\alpha_{\mathrm{O}0}+\psi_{\mathrm{O}}\,\zeta^{\prime},\\\
&\log_{10}(\alpha_{\mathrm{A}}/b^{2})=\log_{10}(\alpha_{\mathrm{A}}/b^{2})_{0}+\psi_{\mathrm{A}}\,\zeta^{\prime},\end{split}$
(21)

where $\alpha_{\mathrm{O}0}$ and $(\alpha_{\mathrm{A}}/b^{2})_{0}$ are the
values at the equitorial plane, $\psi_{\mathrm{O}}$ and $\psi_{\mathrm{A}}$
are the slopes in the logarithmic space, and $\zeta_{t}$ marks the transition
altitude. In §4 we already found that $\varrho$ profiles are almost unaffected
by the accretion and remain almost identical to the static disk. In addition,
although the radial derivative of $B_{z}$ may affect the $b_{R}$ near the mid-
plane in steady-state solutions, we have nevertheless confirmed that this
hardly affects the growth rates of perturbation modes. Therefore, we take the
approximations that $\rho=\exp(-\zeta^{2}/2)$ and $\chi_{z}=0$ (§4.1) for all
models in this section to make the results independent of $R$.

A reference model is constructed and illustrated in the left column of Figure
7. This model has $\alpha_{\mathrm{O}0}=10^{1}$ ,
$(\alpha_{\mathrm{A}}/b^{2})_{0}=10^{-3}$, $\psi_{\mathrm{O}}=-2$,
$\psi_{\mathrm{A}}=0.5$, and $\zeta_{t}=2$, which qualitatively resembles the
$R=1~{}\mathrm{AU}$ slice of the fiducial model in §4.1. The consequent
steady-state solution is also quite similar. For the perturbation modes, the
symmetric mode has decay rate $\nu_{\rm sym}\simeq-0.095$, which has the same
sign, but the absolute value is different by a $\sim 3$ factor. The growth
rate of the anti-symmetric mode is quite close to the fiducial model,
$\nu_{\rm asy}=0.039$.

###  5.2 Impacts of MHD parameters

Figure 7: Left column illustrate the model for reference with the setup in
§5.1, showing the diffusivity profiles (upper panel), the steady-state
solution (middle panel), and the perturbation modes (lower pannel). The growth
rates are $\nu_{\rm sym}=-0.10$ and $\nu_{\rm asy}=0.035$. Right two columns
present the dependence of $\nu_{\rm sym}$ (orange lines) and $\nu_{\rm asy}$
(blue lines) on different paramters: $\beta_{0}$, $\zeta_{t}$,
$\alpha_{\mathrm{O}0}$, and $(\alpha_{\mathrm{A}}/b^{2})_{0}$ (see §5.2). Note
that negative values are shown in dashed lines. The vertical dotted lines are
plotted in each panel indicating the value taken by the reference model.

The reference model described in §5.1 is modified to study the influence of
different parameters. We first confirm that $\mu_{k}$, which reflects the
relative importance of orbital motion, only plays a secondary role: taking
$\mu_{k}\in[4,\ 10^{2}]$ only leads to minor variations in $\nu$, such that
$-0.063\gtrsim\nu_{\rm sym}\gtrsim-0.10$ and $0.03\lesssim\nu_{\rm
asy}\lesssim 0.04$. In contrast, the magnetization parameterized by
$\beta_{0}$ is a significant parameter. As shown in Figure 7, the stronger the
disk’s magnetization (or smaller $\beta_{0}$), the faster the perturbation
modes grow or decay. If one sets $\beta_{0}=10^{3}$ based on the reference
model, such a strong magnetization yields $\nu_{\rm sym}\simeq-1.0$ and
$\nu_{\rm asy}\simeq 0.38$, leading to an e-fold increase of anti-symmetric
amplitude every half orbit period. With $\beta_{0}\lesssim 10^{5}$, the growth
or decay rates scale roughly as $\sim\beta_{0}^{-1/2}$, or roughly linearly
with field strength.

Different components of magnetic diffusivities control the stability problem
in different ways. The ambipolar diffusion parameter
$(\alpha_{\mathrm{A}}/b^{2})_{0}$ affects both $\nu_{\rm sym}$ and $\nu_{\rm
asy}$, yet the absolute values of both $\nu$ vary rather reluctantly as the
parameter stay in the range
$(\alpha_{\mathrm{A}}/b^{2})_{0}\in[10^{-4},10^{-2}]$. The situation is
similar about the dependence of $\nu_{\rm asy}$ upon $\alpha_{\mathrm{O}0}$
and $\zeta_{t}$. When we examine the symmetric modes, however, it is found
that $\nu_{\rm sym}$ becomes positive when $\alpha_{\mathrm{O}0}\lesssim 0.6$,
or $\zeta_{t}\lesssim 1.4$. As the conductivity in the adjacency of the mid-
plane becomes higher, by either lowering the $\alpha_{\mathrm{O}0}$ or setting
a lower transition altitude $\zeta_{t}$, the stronger coupling between field
and fluids eventually leads to the instability even for the symmetric modes.
The threshold value for $\alpha_{\mathrm{O}0}$ also matches the vertically
local mode calculations in §3.2 semi-quantitatively. Meanwhile, $\nu_{\rm
asy}$ does not change sign throughout the subset of parameter space we
explored and illustrated in Figure 7, neither in the search over a broader
span of the parameter space that is not elaborated here. This fact indicates
the ubiquity of the instability for the anti-symmetric modes.

##  6 Discussions and Summary

This study explores the instability that leads to reflection symmetry breaking
over the equatorial planes in protoplanetary disks (PPDs) that go through
wind-driven accretion. We develop linear perturbation theories using steady-
state solutions and evaluate eigenmodes and their corresponding eigenvalues
with different symmetries. Symmetric perturbations always decay under typical
MHD profiles, whereas at least one of the anti-symmetric modes has a positive
growth rate. The growth of anti-symmetric modes leads to tilted disk
kinematics, where the accretion flow occurs on one side, and the disk wind is
launched on the other. We also investigate the effects of varying magnetic
diffusivity profiles on the reflection symmetry breaking, finding that the
anti-symmetry is most prominent in regions of the disk with poorly ionized
mid-plane where fields and fluids are weakly coupled and intermediate surface
layers where the ionization and field-fluid coupling are intermediate.

###  6.1 Qualitative behavior of eigenvalues

Qualitative discussions on the behaviors of eigenvalues are necessary to
develop physical insights into the findings of this paper and guide future
explorations. Eq. (12) can be recast in the integrated form,

$\left[(A+\xi V_{1})\delta
b\right]_{\zeta_{0}}^{\zeta_{1}}=V_{0}\int_{\zeta_{0}}^{\zeta_{1}}\delta
b\;\mathrm{d}\zeta.$ (22)

Since we are most interested in the modes with the least nodes, especially the
prospectively anti-symmetric modes with anti-symmetry, the approximation
$\nu\ll 1$ can be safely adopted, hence $\xi\simeq 2/(\beta_{0}\varrho)$. The
bottom two rows in the matrix form of eq. (22) are indeed equivalent to
$\delta j_{R}=-\partial_{\zeta}\delta b_{\varphi}$ and $\delta
j_{\varphi}=\partial_{\zeta}\delta b_{R}$, while the top two rows should be
analysed for $\nu$,

$\begin{split}&\nu\mathcal{B}_{R}=\left[\alpha_{\varphi\varphi}\delta
j_{\varphi}+(\alpha_{\varphi R}-2\xi)\delta
j_{R}\right]_{\zeta_{0}}^{\zeta_{1}},\\\
&\nu\mathcal{B}_{\varphi}=\left[\alpha_{RR}\delta
j_{R}+(\alpha_{R\varphi}-\xi/2)\delta
j_{\varphi}\right]_{\zeta_{0}}^{\zeta_{1}}+\dfrac{3}{2}\int_{\zeta_{0}}^{\zeta_{1}}\delta
b_{R}\ \mathrm{d}\zeta\ ;\\\
&\mathcal{B}_{R}\equiv\int_{\zeta_{0}}^{\zeta_{1}}\delta b_{R}\
\mathrm{d}\zeta-(\xi\delta j_{\varphi})_{\zeta_{0}}^{\zeta_{1}},\\\
&\mathcal{B}_{\varphi}\equiv-\int_{\zeta_{0}}^{\zeta_{1}}\delta b_{R}\
\mathrm{d}\zeta-(\xi\delta j_{R})_{\zeta_{0}}^{\zeta_{1}}.\end{split}$ (23)

####  6.1.1 Anti-symmetric modes

For the anti-symmetric modes, the equation for $\delta\tilde{b}_{R}$ is
typically easier to analyze by locating $\zeta_{0}=0$ (at which $\delta
j_{R}=\delta j_{\varphi}=0$) and $\zeta_{1}$ at the first zero point of
$\delta b_{R}$ (likely the wind base). Without loss of generality, the
linearity of the equations allows one to set $\delta b_{R}>0$ throughout all
altitudes in the concerned modes, which in turn yields $\delta b_{\varphi}<0$,
$\delta j_{R}|_{\zeta_{1}}=-[\partial_{\zeta}\delta
b_{\varphi}]_{\zeta_{1}}<0$, and $\delta j_{\varphi}|_{\zeta_{1}}<0$ due to
the geometries of disks and fields. The coefficient $\mathcal{B}_{R}$ on is
therefore positive definite, and the sign of $\nu$ depends on the right-hand
side.

The cases that we concern usually have $|\delta
b_{\varphi}|_{\zeta_{1}}\gg|\delta b_{R}|_{\zeta_{1}}$ due to the
amplification of toroidal fields by rotation, expecting that $|\delta
j_{R}|_{\zeta_{1}}\gg|\delta j_{\varphi}|_{\zeta_{1}}$. In the meantime, the
wind base usually have $\alpha_{\mathrm{A}}/b^{2}\sim\alpha_{\mathrm{O}}$,
therefore $|b_{\varphi}|_{\zeta_{1}}\gg|b_{R}|_{\zeta_{1}}$ usually leads to
$|\alpha_{\varphi R}|_{\zeta_{1}}\gg|\alpha_{\varphi\varphi}|_{\zeta_{1}}$,
which have been confirmed by all cases involved in §4. These facts focus on
the competition between $\alpha_{\varphi R}$ and $2\xi$ at $\zeta_{1}$. Since
$b_{R}$ and $b_{\varphi}$ have different signs at the wind base due to the
disk geometries and MHD configurations, $\alpha_{\varphi R}|_{\zeta_{1}}>0$
always holds. Our numerical experiments have found $\alpha_{\varphi
R}|_{\zeta_{1}}<2\xi|_{\zeta_{1}}$ throughout the fiducial model, which is
consistent with the finding that $\nu>0$, i.e. there are always unstable anti-
symmetric modes (Figure 5). In order to eliminate the instability, the
inequality $\alpha_{\varphi R}|_{\zeta_{1}}<2\xi|_{\zeta_{1}}$ has to be
reverted.

####  6.1.2 Symmetric modes

The symmetric modes are easier to analyze by still setting $\zeta_{0}=0$ but
locate $\zeta_{1}$ at the first zero point of $\delta
j_{R}=-\partial_{\zeta}\delta b_{\varphi}$. With such selection, it is easy to
verify that $\delta j_{R,\varphi}|_{\zeta_{0}}>0$. Also, taking $\delta
b_{R}>0$ at $\zeta>0$ due to linearity, one can prove that $\mathcal{B}_{R}$
is still positive. On the right-hand side,
$[\alpha_{\varphi\varphi}j_{\varphi}]^{\zeta_{1}}_{\zeta_{0}}$ is definitely
negative. This time, the competition of $\alpha_{\varphi R}$ versus $2\xi$
usually declares the victory of the former for sufficiently weakly ionized
(thus $\alpha_{\mathrm{O}}$ is large), weakly magnetized disks: because now
the values at $\zeta_{0}$ matter, at which $\varrho=1$, and $\xi=2/\beta_{0}$
is tiny for large $\beta_{0}$. When the fluid-field coupling at the mid-plane
becomes stronger, the inequality $\alpha_{\varphi R}<2\xi$ may no longer hold,
leading to the instability of the symmetric modes. The explorations in §5.2
have already witnessed this situation.

####  6.1.3 Diffusivities and stabilities

The discussions elaborated above can relate to the physical picture of
electric currents and diffusivities. The symmetric modes have non-negligible
$|\delta j_{R,\varphi}|=|\partial_{\zeta}b_{\varphi,R}|$ near $\zeta=0$, which
requires sufficient conductivity to develop. However, weakly ionized PPDs do
not have sufficient conductivity there, and the highly resistive plasmas at
the mid-plane damp and inhibit the growth of symmetric modes. Such damping
depends on $\alpha_{\mathrm{O}}$ more than $\alpha_{\mathrm{A}}/b^{2}$, as the
former is generally much greater than the latter at $\zeta\sim 0$, which is
consistent with the emphases on Ohmic resistivity in the stability analyses
stabilities in e.g. Turner et al. (2007) and Ilgner & Nelson (2008).

Anti-symmetric modes, in contrast, bend their field lines at much higher
altitudes (e.g., Figure 4), where the diffusivities are sufficient to support
the electric currents required. Ambipolar diffusion is of greater importance
at this time since $\alpha_{\mathrm{O}}$ decreases drastically near the disk
surfaces. At the same time, the product $|b_{R}b_{\varphi}|$ has much greater
absolute values in the term $\alpha_{\varphi
R}\simeq-(\alpha_{\mathrm{A}}/b^{2})b_{R}b_{\varphi}$ (note again that
$\alpha_{\mathrm{A}}/b^{2}$ is primarily independent of fields and that
$\alpha_{\mathrm{H}}/b$ is negligible due to the lack of the
$|b_{R}b_{\varphi}|$ factor). Qualitatively, this is also consistent with the
emphasis on ambipolar diffusion in the analyses of disk stabilities in e.g.,
Bai & Stone (2011).

The semi-quantitative assumptions and approximations on diffusivities involved
for the perturbative analyses have been examined in the elaborations in §4 and
§5. At lower disk densities and higher ionization fractions (e.g., magnetized
transitional PPDs; Wang & Goodman 2017b), the stability analyses could lead to
different conclusions, converging to the situations discussed in §3.

###  6.2 Predictions of observations

The anti-symmetric modes in PPDs could result in the absence of a magnetized
wind on one side of the disk (§4.3). Some young stellar objects (YSOs) have
been observed to launch asymmetric winds over both disk surfaces, for example,
HH 30 (Burrows et al., 1996; Stapelfeldt et al., 1999). Some other
observational evidence supports asymmetric winds (see a review in Pascucci et
al. 2023). High spectral resolution [O I] $6300~{}\mathrm{\AA}$ emission line
has been used to trace the disk winds (e.g., Rigliaco et al. 2013; Simon et
al. 2016; McGinnis et al. 2018; Fang et al. 2018; Banzatti et al. 2019; Fang
et al. 2023a). The line profile of [O I] typically consists of two types of
components: a high-velocity component (HVC) and a low-velocity component
(LVC). While HVCs are produced in extended jets (Lavalley-Fouquet et al.,
2000; Bacciotti et al., 2000; Woitas et al., 2002), LVCs most likely trace MHD
winds (Fang et al., 2018; Banzatti et al., 2019; Fang et al., 2023a, b).
Though the majority (70%) of LVCs are blueshifted, 21% of LVCs are red-shifted
by more than $\gtrsim 1~{}{\rm\,km}~{}\mathrm{s}^{-1}$ (Fang et al., 2023a).
These redshifted LVCs can be explained if the winds are launched mainly over
the disk surface facing away from the observer. Future observational studies
are also desired to directly reveal the morphologies of disk winds to confirm
the existence of the asymmetries.

The “tilted” accretion may be related to the asymmetric feeding of the central
protostar, which is potentially related to the recycling of the accreted
materials and the formation of asymmetric jets observed on various PPDs (e.g.
Lavalley-Fouquet et al., 2000; Noriega-Crespo et al., 2002; Flores-Rivera et
al., 2023). Nevertheless, the innermost accretion regions and the formation of
jets in PPDs are complicated, and future studies on the jet morphologies of
the jets are required to reveal and confirm the connections between
asymmetries in jets and accretion layers. Meanwhile, accretion lifted to the
surface layers also reduces the damping of differential rotation and can
induce magnetorotational instabilities (MRI). Once MRI develops in the surface
layer, the magnetized wind will be suppressed, and the wind launching area
shall exhibit a “truncation” subsequently. Recent observations have been
conducted with both spatial and spectral resolutions, which figured out the
range of radii for PPD outflows and examined different wind launching
mechanisms (e.g. Fang et al., 2023b). Similar high-resolution observations are
necessary to confirm these theoretical predictions and understand the complex
physical processes that shape magnetized wind structures. For instance,
observations that identify the wind launching regions with spatial resolution
could establish a direct link between tilted accretion, possible surface MRI,
and the suppression of wind launching in PPDs. Both theoretical modeling and
high-resolution observations are vital for unraveling the interplay between
magnetic fields, turbulence, and gravitational instability.

###  6.3 Future works

This study focuses on examining the vertical modes with axisymmetry in
relation to the wind-driven accretion in protoplanetary disks. This approach
does not account for the radial dependence of the problem. A possible
extension of this work would involve studying the radial modes, which play a
crucial role in the radial transport of magnetic fields. Modulated by the non-
ideal MHD features, the radial flux transport and the accretion flow are
intertwined and can result in the formation and movement of disk substructures
such as rings and gaps. Extending the analytic research on PPD controlled by
non-ideal MHD effects, another prospective area of investigation involves
examining the breaking of the disk axisymmetry by azimuthal modes. The close
relationship between the azimuthal modes and angular momentum transport in
protoplanetary disks can lead to complex phenomena like spiral arms, vortices,
and disk warps. Through studying these phenomena, we can better understand how
angular momentum transporting is related to different perturbation modes and
how that impacts the evolution and dispersal of protoplanetary disks.

Existing theories on PPD substructure gaps often attribute their formation to
the influence of existing planets. Nevertheless, a natural mechanism that can
produce these structures without the influence of readily formed planets would
pave its way to a broader and more generic application. Hence, the proposed
study of co-evolved radial, azimuthal, and vertical modes could provide
crucial insights into the natural formation of PPD substructures and
subsequent evolution over time. Additionally, investigating the interactions
of radial modes with the vertical structures and instabilities, possible
surface MRI, and even gravitational instabilities on PPD substructure
formation could also yield valuable insights for future research. Given that
substructures in protoplanetary disks are key to understanding the formation
and evolution of planetary systems, this line of research is vital in
unlocking some of the mysteries surrounding the early stages of planetary
formation.

L. Wang and S. Xu acknowledge the computation resources provided by the Kavli
Institute of Astronomy and Astrophysics at Peking University. We thank our
collegues for helpful discussions (alphabetical order of the last names): Can
Cui, Gregory Herczeg, Xiao Hu, Haifeng Yang, Xinyu Zheng.

## References

  * Bacciotti et al. (2000) Bacciotti, F., Mundt, R., Ray, T. P., et al. 2000, ApJ, 537, L49
  * Bai (2013) Bai, X.-N. 2013, ApJ, 772, 96
  * Bai (2017) —. 2017, ApJ, 845, 75
  * Bai & Stone (2011) Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144
  * Bai & Stone (2013) —. 2013, ApJ, 769, 76
  * Bai & Stone (2017) —. 2017, ApJ, 836, 46
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1
  * Banzatti et al. (2019) Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76
  * Burrows et al. (1996) Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
  * Fang et al. (2023a) Fang, M., Pascucci, I., Edwards, S., et al. 2023a, ApJ, 945, 112
  * Fang et al. (2018) —. 2018, ApJ, 868, 28
  * Fang et al. (2023b) Fang, M., Wang, L., Herczeg, G. J., et al. 2023b, Nature Astronomy, 7, 905
  * Flores-Rivera et al. (2023) Flores-Rivera, L., Flock, M., Kurtovic, N. T., et al. 2023, A&A, 670, A126
  * Harrison et al. (2021) Harrison, R. E., Looney, L. W., Stephens, I. W., et al. 2021, ApJ, 908, 141
  * Hu et al. (2022) Hu, X., Li, Z.-Y., Zhu, Z., & Yang, C.-C. 2022, MNRAS, 516, 2006
  * Hu et al. (2019) Hu, X., Zhu, Z., Okuzumi, S., et al. 2019, ApJ, 885, 36
  * Ilgner & Nelson (2008) Ilgner, M., & Nelson, R. P. 2008, A&A, 483, 815
  * Lavalley-Fouquet et al. (2000) Lavalley-Fouquet, C., Cabrit, S., & Dougados, C. 2000, A&A, 356, L41
  * McGinnis et al. (2018) McGinnis, P., Dougados, C., Alencar, S. H. P., Bouvier, J., & Cabrit, S. 2018, A&A, 620, A87
  * Nemer et al. (2020) Nemer, A., Goodman, J., & Wang, L. 2020, ApJ, 904, L27
  * Noriega-Crespo et al. (2002) Noriega-Crespo, A., Cotera, A., Young, E., & Chen, H. 2002, ApJ, 580, 959
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880
  * Pascucci et al. (2023) Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 567
  * Rigliaco et al. (2013) Rigliaco, E., Pascucci, I., Gorti, U., Edwards, S., & Hollenbach, D. 2013, ApJ, 772, 60
  * Riols et al. (2020) Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95
  * Simon et al. (2013a) Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ, 775, 73
  * Simon et al. (2013b) Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ, 764, 66
  * Simon et al. (2016) Simon, M. N., Pascucci, I., Edwards, S., et al. 2016, ApJ, 831, 169
  * Stapelfeldt et al. (1999) Stapelfeldt, K. R., Watson, A. M., Krist, J. E., et al. 1999, ApJ, 516, L95
  * Turner et al. (2014) Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
  * Turner et al. (2007) Turner, N. J., Sano, T., & Dziourkevitch, N. 2007, ApJ, 659, 729
  * Vlemmings et al. (2019) Vlemmings, W. H. T., Lankhaar, B., Cazzoletti, P., et al. 2019, A&A, 624, L7
  * Wang et al. (2019) Wang, L., Bai, X.-N., & Goodman, J. 2019, ApJ, 874, 90
  * Wang & Goodman (2017a) Wang, L., & Goodman, J. 2017a, ApJ, 847, 11
  * Wang & Goodman (2017b) —. 2017b, ApJ, 835, 59
  * Wardle (2007) Wardle, M. 2007, Ap&SS, 311, 35
  * Wardle & Königl (1993) Wardle, M., & Königl, A. 1993, ApJ, 410, 218
  * Watson & Stapelfeldt (2004) Watson, A. M., & Stapelfeldt, K. R. 2004, ApJ, 602, 860
  * Woitas et al. (2002) Woitas, J., Ray, T. P., Bacciotti, F., Davis, C. J., & Eislöffel, J. 2002, ApJ, 580, 336
  * Xu & Bai (2016) Xu, R., & Bai, X.-N. 2016, ApJ, 819, 68

##  Appendix A Extra terms in the steady-state solution

The general guideline in simplifying the equations for the local wind-driven
accretion model is to ignore the terms with $R$ on their denominators or with
derivatives of $R$. In §2.1.2 the inclusion of
$\partial_{R}(v_{\varphi}B_{R})$ (as part of the $\partial_{R}E_{z}$ term) has
been explained, and this appendix section discusses two other terms with
$\mu_{k}$ on the prospectively important denominators.

One of them is associated with $J_{z}$, which can be expressed as the
derivatives of magnetic fields,

$\dfrac{4\pi
J_{z}}{c}=\dfrac{\partial_{R}(RB_{\varphi})-\partial_{\varphi}B_{R}}{R}=\dfrac{B_{\varphi}}{R}+\partial_{R}B_{\varphi}.$
(A1)

Note that $\partial_{\varphi}=0$ strictly holds due to the axisymmetry. While
it is unable to examine whether $\partial_{R}B_{\varphi}$ plays an essential
role with radially local models, the rest ($B_{\varphi}/R$) could still make
the terms involving $J_{z}B_{\varphi}$ important with sufficiently large
$|B_{\varphi}|$. In the dimensionless form, this product reads
$j_{z}b_{\varphi}\simeq b_{\varphi}^{2}/\mu_{k}$. As one can observe from the
numerical and semi-analytic results in §4.1 and §5, $|b_{\varphi}|$ could
reach $\sim 10^{0}-10^{1}$, making the $j_{z}b_{\varphi}$ term important with
$\mu_{k}\sim 30$ (eq. 4) or even smaller at larger disk radii. Therefore,
although $j_{z}\sim 0$ in eq. (6) essentially holds in most cases, we still
include the $j_{z}$ term when the steady-state solutions are compared to the
simulation results in §4.1.

The other term is related to the gravitational force. At $z=0$, a local model
assuming negligible radial pressure gradient has
$\partial_{R}\Phi=v_{k}^{2}/R$, where $v_{k}$ is the Keplerian velocity. In
the derivation of the radial momentum equation in eq. (5), this equality is
assumed to be true at all altitudes. However, one can prove that
$|\partial_{R}\Phi-v_{k}^{2}/R|$ may be no longer negligible at sufficiently
high $z$,

$\partial_{R}\Phi-\dfrac{v_{k}^{2}}{R}=\dfrac{GM_{*}}{R^{2}}\left[\left(1+\dfrac{z^{2}}{R^{2}}\right)^{-3/2}-1\right]\simeq\dfrac{v_{k}^{2}}{R}\left(\dfrac{3z^{2}}{2R^{2}}\right).$
(A2)

This leads to an extra term $g_{R}$ in the radial momentum equation (eq. 5),
using $z/R=\zeta/\mu_{k}$, and
$\partial_{z}\rightarrow(\mu_{k}/R)\partial_{\zeta}$,

$g_{R}\equiv\dfrac{R(\partial_{R}\Phi-
v_{k}^{2}/R)}{\mu_{k}c_{s}^{2}}\simeq\dfrac{v_{k}^{2}}{\mu_{k}c_{s}^{2}}\left(\dfrac{3z^{2}}{2R^{2}}\right)=\dfrac{3\zeta^{2}}{2\mu_{k}}.$
(A3)

When we consider the disk surface, $\zeta\sim 3-4$, and the $g_{R}$ could also
reach the order of unity. Similar to the terms related to $j_{z}$, the $g_{R}$
is also included to reach better agreements with numerical simulation results.
One may notice that a similar correction can also be applied to the vertical
momentum equation, yet the other terms dwarf the value of the resulting
correction.

##  Appendix B Electric current equations for steady-state solutions

Assuming $\mu_{z}\rightarrow 0$, eqs. (5) yield the dimensionless velocity
components $\mu_{R,\varphi}$,

$\mu_{\varphi}=\dfrac{j_{z}b_{\phi}-j_{\varphi}}{\beta_{0}\varrho}-\dfrac{g_{R}}{2},\quad\mu_{R}=\frac{4(j_{z}b_{R}-j_{R})}{\beta_{0}\varrho}.$
(B1)

Combining these $\mu_{R,\varphi}$ expressions with eqs. (6) and (7), one gets
the equations which shall be solved for $j_{R,\varphi}$,

$\begin{bmatrix}\alpha_{RR}&\alpha_{R\varphi}+1/(\beta_{0}\varrho)\\\
\alpha_{\varphi
R}-4/(\beta_{0}\varrho)&\alpha_{\varphi\varphi}\end{bmatrix}\begin{bmatrix}j_{R}\\\
j_{\varphi}\end{bmatrix}=\begin{bmatrix}\tilde{\varepsilon}_{R}-g_{R}/2+j_{z}b_{\varphi}/(\beta_{0}\varrho)\\\
-4j_{z}b_{R}/(\beta_{0}\varrho)\end{bmatrix},$ (B2)

where $\tilde{\varepsilon}_{R}\equiv\varepsilon_{R}+\mu_{k}$ is the
dimensionless electric field in the Keplerian rotating frame.

##  Appendix C Matching the ideal MHD wind solutions

Bai et al. (2016) introduced a method of constructing magnetized isothermal
wind models above disk surfaces, which is briefly summarized here. This method
is based on the conserved quantities along the magnetic field lines, including
the Bernoulli parameter $H$ (following from the energy conservation), the
mass-field flux ratio $k$ (following from the continuity equation; not to be
confused with dimensionless wavenumbers in §3), and the field line angular
velocity parameter $\omega$, reading,

$H=\dfrac{B_{\mathrm{p}}^{2}}{2k^{2}x^{2}}+\dfrac{\omega^{2}R^{2}}{2}\left[\left(\dfrac{R_{\mathrm{A}}^{2}/R^{2}-1}{x-1}\right)^{2}-1\right]+h-\dfrac{GM_{*}}{(R^{2}+z^{2})^{1/2}}\
;\quad k\equiv\dfrac{4\pi\rho
v_{\mathrm{p}}}{B_{\mathrm{p}}},\quad\omega=\Omega-\dfrac{kB_{\varphi}}{4\pi\rho
R}.$ (C1)

Here $M_{*}$ is the stellar mass, $x\equiv 4\pi\rho/k^{2}$ is the density
parameter, the subscripts “p” in $v_{\mathrm{p}}$ and $B_{\mathrm{p}}$
indicate the poloidal components of velocities and magnetic fields, $R$ and
$z$ are the cylindrical coordinates of the spatial point on the current field
line, and $R_{\mathrm{A}}$ marks the radius of the Alfvénic point on the same
field line. The thermodynamics of gas elements are involved by the specific
enthalpy $h=c_{s}^{2}\ln(\rho/\rho_{0})$, where $\rho_{0}$ is the wind-base
mass density (and the subscript “0” in general marks the quantities at the
wind base in this appendix). The poloidal field $B_{\mathrm{p}}$ in the wind
is prescribed as,

$B_{\mathrm{p}}=B_{\mathrm{p}0}\dfrac{1+q}{(R/R_{0})+q(R/R_{0})^{2}},$ (C2)

in which the parameter $q$ controls the poloidal fields’ transition from
parallel to diverging. In practice, we find that $q=0.5$ can yield accretion
and wind solutions that match the simulation results best. The vertical
coordinate $z$ is related to $R$ along the designated field line by a
straight-line geometry, $z=z_{0}+(R-R_{0})\tan\theta$. It has been found by
inspecting the simulations that the pitch angle of field lines varies in the
wind; $\theta$ should reflect the overall inclination angles of field lines
instead of using the wind-base value ($\theta=\cot^{-1}b_{R}$). In this
manuscript, we take $\theta=\pi/6$ unless specifically noted. Physically
plausible solutions are obtained by,

$H=E,\quad\dfrac{\partial H}{\partial x}=0,\quad\dfrac{\partial H}{\partial
R}=0.$ (C3)

In order to construct consistent wind-driven accretion models, the accretion
solution described in §2.1 has to be connected to the solution described by
eq. (C1) at the wind base, by matching key physical quantities. We first
convert eqs. (C1), (C2) into the dimensionless form,

$\begin{split}&\mathcal{H}\equiv\dfrac{H}{c_{s}^{2}}=\left(\dfrac{b_{\mathrm{p}0}^{2}}{2\tilde{k}^{2}x^{2}}\right)\left(\dfrac{1+q}{\tilde{R}^{2}+q\tilde{R}^{2}}\right)+\dfrac{\tilde{\omega}^{2}\tilde{R}^{2}}{2}\left[\left(\dfrac{\tilde{R}_{\mathrm{A}}^{2}/\tilde{R}^{2}-1}{x-1}\right)^{2}-1\right]+\ln
x-\dfrac{\mu_{k}}{\tilde{r}(R)}\ ;\
b_{\mathrm{p}0}\equiv\dfrac{B_{\mathrm{p}0}}{B_{z0}},\
\tilde{R}\equiv\dfrac{R}{R_{0}},\
\tilde{R}_{\mathrm{A}}\equiv\dfrac{R_{\mathrm{A}}}{R_{0}};\\\
&\tilde{k}\equiv\dfrac{kc_{s}}{B_{z0}}=\left[\left(\dfrac{2}{\beta_{0}\varrho}\right)^{-1}\dfrac{\mu_{R}}{b_{R}}\right]_{\rm
wb},\quad x_{\rm
wb}=\left[\left(\dfrac{2}{\beta_{0}\varrho}\right)\tilde{k}\right]_{\rm
wb}^{-1},\quad\tilde{\omega}\equiv\dfrac{\omega
R_{0}}{c_{s}}=\left[(\mu_{\varphi}+\mu_{k})-\left(\dfrac{2}{\beta_{0}\varrho}\right)\tilde{k}b_{\phi}\right]_{\rm
wb},\\\ \end{split}$ (C4)

where the subscripts “wb” indicate that the variables take the wind-base
values. Note that the assumption $\mu_{z}=0$ in the steady-state accretion.
For the consistency of such matching, we assume that
$\mu_{z}/\mu_{R}=b_{z}/b_{R}$ at the maching point at which the gas has
relatively high conductivities and the velocity field is likely to be largely
parallel to the magnetic field. Note that this approach only affects the value
of $\tilde{k}$ by approximating $\mu_{\mathrm{p}}/b_{\mathrm{p}}$ with
$\mu_{R}/b_{R}$ at the wind base, and our numerical experiments have confirmed
that $\mu_{z}$ only have secondary effects on the disk accretion and wind
solutions as long as the field inclination angle $\theta$ has been prescribed.
The following procedures are taken to obtain a set of solutions with matched
wind and accretion profiles:

  1. 1.

Select a $\varepsilon_{R0}$ (the radial component of the mid-plane
dimensionless electric field; see §2.1.2) and integrate the steady-state
solution to the wind base;

  2. 2.

Calculate the values of $b_{\mathrm{p}0}$, $\tilde{k}$, $\tilde{\omega}$, $x$,
and $\mathcal{H}_{\rm wb}$ at the wind base;

  3. 3.

Solve for a physically plausible wind solution, using $b_{\mathrm{p}0}$,
$\tilde{k}$, $\tilde{\omega}$ calculated in Step 2, by solving the algebraic
equations $\partial\mathcal{H}/\partial
x=\partial\mathcal{H}/\partial\tilde{R}=0$ for the regularity conditions at
both slow and fast magnetosonic points;

  4. 4.

Adjust $\tilde{R}_{\mathrm{A}}$ so that the values of $\mathcal{H}$ at slow
and fast magnetosonic points are equal;

  5. 5.

Adjust $\varepsilon_{R0}$ and conduct Steps 1 through 4 iteratively, so that
the $\mathcal{H}$ value in Step 4 equals to the $\mathcal{H}_{\rm wb}$ in Step
2.