© ESO, 2016

1. Introduction

The evolution of protoplanetary disks (PPDs) is one of the keys to understand planet formation. There are still several unsolved problems, one of which is the dispersal of PPDs (Haisch et al. 2001; Hernández et al. 2008; Takagi et al. 2014, 2015). The evolution and dispersal of PPDs have been extensively studied in the framework of viscously accreting discs that undergo photoevaporation by the irradiation from the central star (e.g., Shu et al. 1993; Hollenbach et al. 2000; Alexander et al. 2006; Kimura et al. 2016).

In addition to the viscous accretion and the photoevaporation, the role of magnetically driven disc winds has recently been received new attention. Suzuki & Inutsuka (2009) and Suzuki et al. (2010) proposed that vertical outflows driven by magnetohydrodynamical (MHD) turbulence might be a viable mechanism that disperses the gas component of PPDs; turbulence is triggered by magnetorotational instability (MRI; Velikhov 1959; Chandrasekhar 1961; Balbus & Hawley 1991), and the Poynting flux associated with the MHD turbulence drives vertical outflows. The idea of MHD turbulence-driven outflow has also been extended by considering various effects, such as a stronger magnetic field (Bai & Stone 2013a), a large-scale magnetic field (Lesur et al. 2013), and the dynamics of dust grains (Miyake et al. 2016), whereas its mass flux is still quantitatively uncertain (Fromang et al. 2013).

Although Suzuki et al. (2010) considered mass loss to be the sole role of the disc wind, the disc wind in reality also carries off the angular momentum (Blandford & Payne 1982; Pelletier & Pudritz 1992; Ferreira et al. 2006; Salmeron et al. 2011). In particular, a dead zone, which is an MRI-inactive region because of the insufficient ionization, is supposed to form in a PPD (Gammie 1996; Sano et al. 2000). In a dead zone the level of the excited turbulence is low, and it is not sufficient to sustain the observed mass accretion onto the central star. In these circumstances, the extraction of the angular momentum by the disc wind possibly plays a primary role in driving mass accretion (Bai & Stone 2013b; Simon et al. 2013). Bai et al. (2016) and Bai (2016) investigated the global evolution of PPDs in such a wind-driven accretion state, by also taking the effect of external heating into account, and reported that a large portion of the mass is removed by the disc wind in comparison to the accreting mass.

A critical open question concerning the disc wind from PPDs is that the mass-loss rate. At the later stage of the evolution, a wind footpoint that is determined by the irradiation from a central star is expected to primarily control the mass-loss rate (Bai et al. 2016; Bai 2016). On the other hand, at the earlier stage when the surface density is high, viscous heating plays an essential role in determining the thermal properties of PPDs (e.g., Ruden & Lin 1986; Nakamoto & Nakagawa 1994; Hirose & Turner 2011; Oka et al. 2011; Bitsch et al. 2015). To investigate the time evolution from the early epoch, we here take the effect of viscous heating in the global evolution of PPDs into account in addition to the loss of mass and angular momentum by the disc wind. We focus in particular on the conditions that create a density structure that is very different from the structure of classic viscously accreting discs, which may help solving long-standing problems such as the radial migration of pebbles, boulders, and protoplanets. For this goal, we evaluate the mass-loss rate from the global energetics of PPDs; the kinetic energy of the vertical outflow is mainly supplied from the gravitational accretion energy. This strategy is different from the method adopted by Bai (2016), in which the mass-loss rate was estimated based on the local profile of magnetically driven wind with external heating. A comparison between the two models is provided in Sect. 4.4.

2. Model

2.1. Basic definitions

We investigated the time evolution of PPDs with magnetically driven disc winds. Suzuki et al. (2010) solved the evolution of PPDs with MRI-triggered disc winds under simplified assumptions: the temperature is locally constant with time, and the disc wind only contributes to the mass loss without removing additional angular momentum. In this paper, we relaxed these assumptions to treat more realistic evolution of PPDs. We considered the heating by viscous accretion (Shakura & Sunyaev 1973; Nakamoto & Nakagawa 1994; Hueso & Guillot 2005) and the effect of disc wind torque on mass accretion (Blandford & Payne 1982; Pelletier & Pudritz 1992; Salmeron et al. 2011; Bai & Stone 2013b)

Throughout this paper, we assume that each annulus at radial distance r from a central star almost rotates with Keplerian frequency, Ω_K$\mathrm{\Omega }\mathrm{\approx }{\mathrm{\Omega }}_{\mathrm{K}}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathit{\star }}}{{\mathit{r}}^{\mathrm{3}}}}\mathit{,}$(1)where G is the gravitational constant and M_⋆ is the mass of the central star. We considered a central star with solar mass, M_⋆ = M_⊙. We defined a vertical scale height, H, of a disc $\mathit{H}\mathrm{=}\frac{\sqrt{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}}{\mathrm{\Omega }}\mathit{,}$(2)where c_s is the sound speed. Temperature T and c_s are related through ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }{\mathit{m}}_{\mathrm{H}}}\mathit{,}$(3)where k_B is the Boltzmann constant, m_H is the proton mass, and we assume mean molecular weight, μ = 2.34 (Hayashi 1981). A different definition for the scale height from ours, c_s/ Ω (without $\sqrt{\mathrm{2}}$), is sometimes used in literatures.

2.2. Evolution of surface density

We treated the time evolution of the radial profile of surface density, Σ = ^∫dzρ, of a disc (1 + 1 D model), while basic formula transformation is done in cylindrical coordinates, (r,φ,z). The time evolution of Σ(r) can be expressed as (see Appendix A for the derivation)

$\begin{array}{ccc}& & \frac{\mathit{\partial }\mathrm{\Sigma }}{\mathit{\partial t}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}[\frac{\mathrm{2}}{\mathit{r}\mathrm{\Omega }}\{\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}\mathrm{\int }\mathrm{d}\mathit{z}\left(\mathit{\rho }{\mathit{v}}_{\mathit{r}}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi }}\right)\right)\\ & & \end{array}$(4)

where δv_φ = v_φ−rΩ is the deviation from the background rotation, and the subscript w stands for disc wind (see below). The [···] parenthesis of the second term represents radial mass flow,

$\begin{array}{ccc}& & \mathrm{-}\mathit{r}\mathrm{\Sigma }{\mathit{v}}_{\mathit{r}}\mathrm{=}\frac{\mathrm{2}}{\mathit{r}\mathrm{\Omega }}\{\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}\mathrm{\int }\mathrm{d}\mathit{z}\left(\mathit{\rho }{\mathit{v}}_{\mathit{r}}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi }}\right)\right)\\ & & \end{array}$(5)

which is derived from the radial balance of angular momentum (Appendix A), and the third term of Eq. (4) denotes the mass loss by the disc wind. The second term consists of the rφ and φz components of Reynolds and Maxwell stresses. The rφ component represents the mass accretion (or decretion) induced by the transport of angular momentum through MHD turbulence. We used the following parametrization based on the α-prescription introduced by Shakura & Sunyaev (1973): $\mathrm{\int }\mathrm{d}\mathit{z}\left(\mathit{\rho }{\mathit{v}}_{\mathit{r}}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi }}\right)\mathrm{\equiv }\mathrm{\int }\mathrm{d}\mathit{z\rho }{\mathit{\alpha }}_{\mathit{r\phi }}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\equiv }\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{,}$(6)where $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ is the mass-weighted vertical average of α_rφ. α_rφ is a nondimensional parameter normalized by gas pressure ($\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$) that describes the transport of angular momentum. We considered α_rφ to originate from the MHD turbulence induced by MRI. $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{(}\lesssim \mathrm{1}\mathrm{)}$ depends on physical conditions of PPDs, such as the ionization and the strength of poloidal magnetic field; see Sect. 2.6 for our adopted values. Although we did not separate the contributions from the Reynolds stress (ρv_rδv_φ) and from the Maxwell stress (− B_rB_φ/ 4π> 0), the latter usually dominates the former by a factor of ~4 in accretion discs with MRI turbulence (Sano et al. 2004; Pessah et al. 2006; Hawley et al. 2011) .

$\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ is an effective turbulent viscosity, and it is mathematically related to viscosity, ν, appeared in a hydrodynamical equation, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{\nu r}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\mathrm{\approx }\frac{\mathrm{3}}{\mathrm{2}}\mathit{\nu }\mathrm{\Omega }\mathit{,}$(7)where the second equality comes from the condition of the Keplerian rotation. The definition of α is not consistent throughout the literature; for example, ν ≈ α_tHc_s, is often used conventionally (e.g., Balbus & Hawley 1998). These two α’s are related by ${\mathit{\alpha }}_{\mathrm{t}}\mathrm{\approx }\frac{\sqrt{\mathrm{2}}}{\mathrm{3}}\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, where note again that the definition of H (Eq. (2)) is also not consistent in the literatures.

The φz component of the second term of Eq. (4) indicates the mass accretion induced by the angular momentum loss with magnetized disc winds, which was not taken into account in Suzuki et al. (2010). The term of (··· )_w represents the sum of the angular momentum flux density carried away by the magnetized outflows from the top and bottom surfaces of a disc. While Reynolds (ρδv_φv_z) and Maxwell (− B_φB_z/ 4π( > 0)) stresses contribute to the φz stress as well, the latter usually dominates in magnetized accretion discs (e.g., Pelletier & Pudritz 1992), similarly to the rφ component. This magnetic braking effect needs to be evaluated in the wind region where it operates; this is the reason why the subscript w is necessary in this term. To incorporate the effect of the wind torque into the 1+1D (t–r) model, α_φz needs to be evaluated by physical quantities at the midplane, and we adopted a similar parametrization to the rφ component, ${\left(\mathit{\rho \delta }{\mathit{v}}_{\mathit{\phi }}{\mathit{v}}_{\mathit{z}}\mathrm{-}\frac{{\mathit{B}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{z}}}{\mathrm{4}\mathit{\pi }}\right)}_{\mathrm{w}}\mathrm{\equiv }\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathit{\phi z}}{\mathrm{)}}_{\mathrm{w}}\mathrm{\equiv }\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathit{,}$(8)where we define the nondimensional stress, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$, normalized by density, ${\mathit{\rho }}_{\mathrm{mid}}\mathrm{(}\mathrm{=}\mathrm{\Sigma }\mathit{/}\mathrm{(}\sqrt{\mathit{\pi }}\mathit{H}\mathrm{\left)}\mathrm{\right)}$, at the midplane,

The third term, (ρv_z)_w, of Eq. (4) represents the sum of the mass loss by the vertical outflows from the upper and lower disc surfaces. Suzuki et al. (2010) introduced the nondimensional mass flux normalized by the density and the sound speed at the midplane: $\mathrm{(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}\mathrm{=}{\mathit{C}}_{\mathrm{w}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}{\mathrm{)}}_{\mathrm{mid}}\mathit{.}$(9)We model C_w in Sect. 2.3.

Substituting Eqs. (6), (8), and (9) into Eq. (4), we finally have $\begin{array}{ccc}\frac{\mathit{\partial }\mathrm{\Sigma }}{\mathit{\partial t}}& \mathrm{-}& \frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\frac{\mathrm{2}}{\mathit{r}\mathrm{\Omega }}\left\{\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\right\}\right]\\ & & \end{array}$(10)We solved this equation for different sets of the three parameters, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$, and C_w. We note that Bai (2016) recently derived essentially the same equation in a different form using mass-loss rate and mass accretion rate instead of the above three-dimensionless parameters.

2.3. Mass-loss rate by disc winds: energetics

We assumed that the energy of the disc wind originates from gravitational accretion. Then, the mass flux of the disc wind, C_w, is constrained by $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ and $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$. A starting point for this energetics constraint is the conservation equation of total MHD energy (e.g., Balbus & Hawley 1998),

$\begin{array}{ccc}& & \frac{\mathit{\partial }}{\mathit{\partial t}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{\mathit{p}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }}\right]\mathrm{+}{\nabla }{·}[{v}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{\mathit{\gamma p}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\right)\\ & & \end{array}$(11)

where p is the gas pressure, γ is a ratio of specific heats, $\mathrm{\Phi }\mathrm{=}\mathrm{-}\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathit{/}\mathit{r}\mathrm{=}\mathrm{-}{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathrm{\approx }\mathrm{-}{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}$ is the gravitational potential by a central star, and F_ot is other contributions to energy flux in addition to the MHD energy, such as thermal conduction and radiative heating or cooling. We considered thin discs with nearly Keplerian rotation (Eq. (1)), and hence, we can assume $\mathit{r}\mathrm{\Omega }\mathrm{\gg }{\mathit{v}}_{\mathit{r}}\mathit{,\delta }{\mathit{v}}_{\mathit{\phi }}\mathit{,}{\mathit{v}}_{\mathit{z}}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathit{,B}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$, and safely neglect the terms concerning gas pressure. Leaving dominant terms in Eq. (11) we finally obtained an approximated energy equation as (Appendix B; Eq. (B.8))

$\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\left(\mathrm{-}\mathrm{\Sigma }\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}[\mathit{r}\mathrm{\Omega }\left\{\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\right\}& & \\ \mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega \Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{]}}\mathrm{+}\mathrm{(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}{\mathit{E}}_{\mathrm{w}}\mathrm{+}{\mathit{F}}_{\mathrm{rad}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(12)

where E_w is the specific total energy of the gas in the disc wind; (ρv_z)_wE_w is the energy carried away by the disc wind. F_rad is radiation loss from the top and bottom surfaces, ${\mathit{F}}_{\mathrm{rad}}\mathrm{=}\mathrm{2}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{surf}}^{\mathrm{4}}\mathit{,}$(13)where σ_SB is the Stefan-Boltzmann constant and T_surf is the temperature at the disc surfaces. We here neglected the energy gain by the irradiation from a central star (Kusaka et al. 1970; Dullemond et al. 2002; Davis 2005) and other external sources. The effect of stellar irradiation was taken into account later when we estimated the temperature.

Equation (12) contains two terms with $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$; the first term in {···} denotes the liberated gravitational energy by mass accretion, and second term outside {···} represents heating by turbulent dissipation, which phenomenologically corresponds to viscous heating. The wind torque, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$, does not contribute to this effective viscous heating because the disc wind does not transport angular momentum within the disc but simply removes it, although $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ contributes to the mass accretion.

Using Eq. (10), we can eliminate the time derivative term of Eq. (12) to derive an energetics constraint on the disc wind (Appendix B):

$\begin{array}{ccc}& & \mathrm{(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}\left({\mathit{E}}_{\mathrm{w}}\mathrm{+}\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{+}{\mathit{F}}_{\mathrm{rad}}\mathrm{=}\\ & & \\ & & \end{array}$

The physical meaning of Eq. (14) is that the energy carried away by disc winds (first term on the left-hand side; l.h.s. hereafter) and radiation (second term on the l.h.s.) is determined by the gravitational energy liberated by accretion (first term on the right-hand side; r.h.s. hereafter) and effective viscous heating (second term on the r.h.s.). We used the Keplerian rotation (Eq. (1)) to transform Eqs. (14) to (15). The term with $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ includes contributions from the gravitational accretion and from the effective viscous heating.

Suzuki et al. (2010) assumed that ${\mathit{E}}_{\mathrm{w}}\mathrm{\ge }\frac{\mathrm{3}}{\mathrm{2}}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}$ is the condition to drive the vertical outflow to a large distance (Eq. (22) of Suzuki et al. 2010). However, we adopt E_w ≥ 0, because this is the sufficient condition for the wind material to reach z ⇒ ∞ (${\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathit{>}\mathrm{0}$ in Eq. (B.5)). Following this consideration, we derived the mass flux of the disc wind that satisfies the energetics constraint with E_w = 0 from Eqs. (14) and (15) in a nondimensional form: $\begin{array}{ccc}{\mathit{C}}_{\mathrm{w}\mathit{,}\mathrm{e}}\mathrm{+}\frac{\mathrm{2}{\mathit{F}}_{\mathrm{rad}}}{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}{\mathrm{)}}_{\mathrm{mid}}}& \mathrm{=}& \frac{\mathrm{2}}{{\mathit{r}}^{\mathrm{3}}\mathrm{\Omega }\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}{\mathrm{)}}_{\mathrm{mid}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}{\mathit{c}}_{\mathrm{s}}}{\mathit{r}\mathrm{\Omega }}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\\ & & \\ & \mathrm{=}& \\ & \mathrm{=}& \mathrm{3}\sqrt{\mathit{\pi }\mathit{/}\mathrm{2}}{\mathit{h}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{+}\sqrt{\mathrm{2}}\mathit{h}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathit{,}\end{array}$where C_w,e stands for the mass flux constrained by the energetics. We here used $\mathrm{\Sigma \Omega }\mathrm{=}\sqrt{\mathrm{2}\mathit{\pi }}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}{\mathrm{)}}_{\mathrm{mid}}$, and for the last equality we introduced an aspect ratio, $\mathit{h}\mathrm{\equiv }\mathit{H}\mathit{/}\mathit{r}\mathrm{=}\sqrt{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}\mathit{/}\mathit{r}\mathrm{\Omega }$.

It is crucial to determine the fractions of the energy transferred to the disc winds (first term on the l.h.s. of Eq. (16)) and to the radiation loss (second term). Following the standard accretion disc model (Shakura & Sunyaev 1973), the available energy from the viscous accretion is transferred to the radiation. In the magnetocentrifugal driven wind model (Blandford & Payne 1982), the angular momentum carried by disc winds is directly related to the wind mass-loss rate. Based on these models, we may infer that the $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ term in Eq. (17) regulates the F_rad term and the $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ term determines C_w,e. However, the situation is not this simple, because disc winds can be launched solely by the $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ term, which was shown by local shearing box simulations with zero-wind torque, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{=}\mathrm{0}$ (Suzuki & Inutsuka 2009). MRI excites MHD turbulence and the associated Poynting flux drives vertical outflows. The original energy source in this mechanism is the energy released by the gravitational accretion.

Despite these complicated problems, we adopted two different strategies to determine C_w,e and F_rad in this paper. The first strategy is that F_rad is equal to the effective viscous heating and all the liberated gravitational energy is transferred to the disc winds. The first corresponds to the first line on the right-hand side of Eq. (16), and the second corresponds to the second line, and then, $\begin{array}{ccc}{\mathit{C}}_{\mathrm{w}\mathit{,}\mathrm{e}}& \mathrm{=}& \\ {\mathit{F}}_{\mathrm{rad}}& \mathrm{=}& \end{array}$where we avoided negative values of C_w,e and F_rad.

In the second choice we left the uncertainty to a parameter, ϵ_rad, that determines the fractional energy to the radiation loss: $\begin{array}{ccc}{\mathit{C}}_{\mathrm{w}\mathit{,}\mathrm{e}}& \mathrm{=}& \\ & \mathrm{=}& \mathrm{(}\mathrm{1}\mathrm{-}{\mathit{ϵ}}_{\mathrm{rad}}\mathrm{)}\left[\mathrm{3}\sqrt{\mathit{\pi }\mathit{/}\mathrm{2}}{\mathit{h}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{+}\sqrt{\mathrm{2}}\mathit{h}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\right]\\ {\mathit{F}}_{\mathrm{rad}}& \mathrm{=}& \end{array}$Since the first method is an extreme limit for the maximum disc wind flux, we sought the other extreme limit of great radiation loss in the second method; we adopted ϵ_rad = 0.9. We name the first case (Eqs. (18) and (19)) strong DW and the second case (Eqs. (20) and (21) with ϵ_rad = 0.9) weak DW from here on; DW stands for disc wind.

On the other hand, local MHD shearing box simulations also give the mass flux of disc winds (Suzuki & Inutsuka 2009; Suzuki et al. 2010). We constrained the mass flux of the local simulations, C_w,0, by the energetics of the global accretion to give the C_w that we use in our calculations, ${\mathit{C}}_{\mathrm{w}}\mathrm{=}\begin{array}{c}\min \end{array}\mathrm{\left(}{\mathit{C}}_{\mathrm{w}\mathit{,}\mathrm{0}}\mathit{,}{\mathit{C}}_{\mathrm{w}\mathit{,}\mathrm{e}}\mathrm{\right)}\mathit{,}$(22)where the adopted C_w,0 is presented in Sect. 2.6.

2.4. Temperature: viscous heating and radiative equilibrium

By referring to the terms concerning $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ in Eq. (14), the viscous heating rate can be scaled as ~$\mathrm{\Sigma \Omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$. Since Σ decreases with t and $\mathrm{\Omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ has a negative dependence on r, the viscous heating is anticipated to play a primary role in determining the temperature in the inner region (≲10 au) and at the early stage of the evolution of a PPD. As Σ decreases with the dispersal of the gas component, the disc evolves passively by the illumination from the central star. A number of works have been published that treat this problem with detailed models that include viscous heating and stellar irradiation (e.g., Garaud & Lin 2007; Oka et al. 2011; Bitsch et al. 2015).

If the viscous heating is more effective in a PPD than stellar irradiation, then the temperature at the midplane, T_mid, will be higher than T_surf in Eq. (13). On the other hand, if the viscous heating is ineffective and the stellar irradiation dominates, then T_surf will be higher than T_mid. The radiative transfer needs to be solved to determine the vertical temperature profile. However, since our main focus here is to investigate the roles of magnetically driven disc winds, we adopt the simple prescription for the temperature that was introduced by Nakamoto & Nakagawa (1994). We defined T_vis as the temperature at the midplane determined by viscous heating, $\mathrm{2}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{vis}}^{\mathrm{4}}\mathrm{=}\left(\frac{\mathrm{3}}{\mathrm{8}}{\mathit{\tau }}_{\mathrm{R}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\tau }}_{\mathrm{P}}}\right){\mathit{F}}_{\mathrm{rad}}$(23)where τ_R and τ_P are the Rosseland mean optical depth and the Planck mean optical depth measured at the midplane. τ_R is estimated from the surface density and the Rosseland mean opacity, κ_R, (Hueso & Guillot 2005) as ${\mathit{\tau }}_{\mathrm{R}}\mathrm{=}{\mathit{\kappa }}_{\mathrm{R}}\mathrm{\Sigma }\mathit{/}\mathrm{2}\mathit{,}$(24)where we use $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{R}}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \end{array}& & \end{array}$(25)based on the opacity of dust grains (Nakamoto & Nakagawa 1994, see also Baillié et al. 2015). The Planck mean optical depth can be approximated as ${\mathit{\tau }}_{\mathrm{P}}\mathrm{=}\begin{array}{c}\max \end{array}\mathrm{\left(}\mathrm{2.4}{\mathit{\tau }}_{\mathrm{R}}\mathit{,}\mathrm{0.5}\mathrm{\right)}$(26)(Nakamoto & Nakagawa 1994; Hueso & Guillot 2005), where we give the lower bound on τ_P to obtain the pre-factor of Eq. (23), ${}^{\mathrm{(}}\frac{\mathrm{3}}{\mathrm{8}}{\mathit{\tau }}_{\mathrm{R}}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}{\mathit{\tau }}_{\mathrm{P}}}}^{\mathrm{)}}\mathrm{\Rightarrow }\mathrm{1}$, for the optically thin limit.

We can also define the temperature under the radiation equilibrium, which is determined by the irradiation from the central star, ${\mathit{T}}_{\mathrm{req}}\mathrm{=}{\mathit{T}}_{\mathrm{1}\mathrm{au}}{\left(\frac{\mathit{r}}{\mathrm{1}\mathrm{au}}\right)}^{\mathit{p}}\mathrm{·}$(27)We adopted T_1au = 280 K and p = −1/2 based on the simple radiative equilibrium for the original minimum mass solar nebula (MMSN) (Hayashi 1981; Hayashi et al. 1985). We note that a slightly different scaling is derived, when the geometry of a flared disc is taken into account (Chiang & Goldreich 1997; Chiang & Youdin 2010).

When a PPD becomes optically thin and the viscous heating is ineffective, not only T_surf but also T_mid approaches T_req. To take both viscous heating and stellar irradiation into account, we tool the sum of these two temperatures, ${\mathit{T}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{vis}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{req}}^{\mathrm{4}}\mathit{,}$(28)for the representative z-averaged temperature, T, to estimate c_s in Eq. (3).

2.5. Initial and boundary conditions

We calculated the evolution of Σ of the initial profile, ∝r^{− 3/2} (Hayashi 1981; Hayashi et al. 1985), with a cut-off radius r_cut${\mathrm{\Sigma }}_{\mathrm{int}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{1}\mathrm{au}}{\left(\frac{\mathit{r}}{\mathrm{1}\mathrm{au}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\exp \left(\mathrm{-}\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{cut}}}\right)\mathrm{·}$(29)The original MMSN by Hayashi (1981) considered Σ_1au = 1.7 × 10³ g cm^-2 with a sharp cut-off at 36 au, which gives the initial disc mass, M_disc,int = 0.013 M_⊙. We adopted a ten times larger Σ_1au = 1.7 × 10⁴ g cm^-2 but slightly smaller r_cut = 30 au in this paper, which gives M_disc,int = 0.11 M_⊙. Mass accretion rates are observationally obtained as a function of time (Gullbring et al. 1998; Hartmann et al. 1998; Ricci et al. 2010; Manara et al. 2016), which corresponds to the age of the central stars, while the MMSN corresponds to a late stage of the evolution. Therefore, we chose the massive initial disc to directly compare our results to these observations.

We solved Eq. (10) to track the time evolution of Σ in the region from r_in = 0.01 au to r_out = 10⁴ au with grid spacing, $\mathrm{\Delta }\mathit{r}\mathrm{\propto }\sqrt{\mathit{r}}$. At the inner and outer boundaries, r = r_in and = r_out, we imposed $\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathrm{\Sigma }{\mathit{r}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{0}$, which corresponds to the zero-torque boundary condition (Lynden-Bell & Pringle 1974); the $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ term in Eq. (10), $\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}$, is zero for a constant $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ and ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\propto }{\mathit{r}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ (Eq. (27) with p = −1/2).

2.6. Parameters

The free parameters of our model are turbulent viscosity, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, disc wind mass flux, C_w,0, and disc wind torque, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$. We would like to note that, although we here call $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ turbulent viscosity, large-scale magnetic fields possibly contribute to $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ in realistic situations (Turner & Sano 2008; Johansen et al. 2009).

2.6.1. Turbulent viscosity – $\overline{{\mathit{\alpha }}_{{r}\mathit{\phi }}}$

We compared two cases with spatially uniform $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-3}$, and 8 × 10^-5. $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-3}$ was adopted from the result of local shearing box MHD simulations with sufficient ionization (Suzuki et al. 2010, see also e.g. Sano et al. 2004; Sai et al. 2013) in which MHD turbulence is fully developed by the MRI. When the ionization is not sufficient and non-ideal MHD effects such as resistivity, Hall diffusion, and ambipolar diffusion are important, a magnetically inactive dead zone forms (Gammie 1996) and $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ is smaller (Sano et al. 1998; Lesur & Longaretti 2007; Simon et al. 2011; Okuzumi & Hirose 2011; Flock et al. 2012; Gressel et al. 2015). We adopted $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-5}$ for such MRI-inactive circumstances. Although we assumed constant $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ for simplicity, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ would be spatially dependent on r and evolve with time in realistic situations, because a dead zone generally forms only in the inner region and its size shrinks with time (e.g., Sano et al. 2000; Suzuki et al. 2010; Dzyurkevich et al. 2013). For future elaborate studies, we need to take this spatially and time-dependent $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ into account.

2.6.2. Disc wind mass flux – C_w,0

The mass flux of disc winds, C_w,0, was also adopted from the local simulations. C_w,0 is controlled by the density at the wind onset region, which is located at the upper regions where the magnetic energy becomes comparable to the thermal energy. For the MRI turbulence, depending on the net vertical magnetic field, the density at the wind footpoint is ≈10^-5−10^-4 times the density at the midplane, which gives C_w,0 ≈ 10^-5−10^-4. Here, we add a note of caution: the local simulations might overestimate the mass-loss rate of the disc winds because the returning mass to the simulation box cannot be properly taken into account. Suzuki et al. (2010) reported that the mass flux is reduced by a factor of 2–3 in simulations with a larger vertical box size. Fromang et al. (2013) also pointed out that the reduction factor could be as large as ~10, but their numerical scheme and other detailed set-up were different from those used in Suzuki et al. (2010). These results show that we must choose C_w,0 carefully from the local simulations.

When we take the face value of the local simulations assuming the ideal MHD condition, C_w,0 ≈ 4 × 10^-5 for the weak vertical magnetic field (Suzuki & Inutsuka 2009). We here set a more conservative value, C_w,0 = 2 × 10^-5, for the MRI-active cases with $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-3}$. If a dead zone is formed, then the mass flux of the disc winds is slightly reduced, but it does not become as low as $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ because the disc winds are driven from the surface regions with sufficient ionization; C_w,0 is only moderately weakened by a factor of a few. We adopted C_w,0 = 1 × 10^-5 for $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-5}$. Moreover, the actual mass flux, C_w, is constrained by the energetics, Eq. (22). We also assumes, in the same way as $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, constant C_w,0 for simplicity. While in realistic situations it would depend on r and vary with time, it does not change as much as $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$.

2.6.3. Disc wind torque – $\overline{{\mathit{\alpha }}_{\mathit{\phi }\mathit{z}}}$

We tested two types of the parametrization for the wind torque: (i) constant $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{=}\mathrm{1}\mathrm{\times }{\mathrm{10}}^{-4}$; and (ii) density dependent with a cap, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{=}\begin{array}{c}\min \end{array}\left({\mathrm{10}}^{-5}{\left(\frac{\mathrm{\Sigma }}{{\mathrm{\Sigma }}_{\mathrm{int}}}\right)}^{-0.66}\mathit{,}\mathrm{1}\right)\mathit{.}$(30)We name (i) constant torque and (ii) Σ-dependent torque from now on. $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ was estimated by local MHD simulations by Bai (2013), who reported $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{~}{\mathrm{10}}^{-5}\mathrm{-}{\mathrm{10}}^{-3}$ with a positive dependence on the strength of the net vertical magnetic field, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{\propto }\mathrm{(}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{8}\mathit{\pi }\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}{\mathrm{)}}^{\mathrm{0.66}}$. ρ_mid is proportional to Σ, while B_z is determined by the inward dragging and outward diffusion of magnetic flux (Lubow et al. 1994; Okuzumi et al. 2014; Guilet & Ogilvie 2014, see also Sect. 4.4). If B_z decreases with the dispersal of gas (decrease of Σ), then $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ will stay approximately constant, which corresponds to (i) constant torque; if B_z does not decrease that much, then $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ has a negative dependence on Σ and will increase with time, which corresponds to (ii) Σ-dependent torque. We tested these two extreme limits for the effect of the wind torque affected by the evolution of the vertical magnetic flux.

3. Results

In this section, we present the properties of the time evolution of PPDs in the MRI-active and MRI-inactive conditions.

3.1. MRI-active cases

Fig. 1 Comparison of time evolutions of four MRI-active PPDs with $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-3}$. The four cases are (i) strong DW + Σ-dependent torque (red); (ii) strong DW + zero-torque (green); (iii) weak DW + zero-torque (purple); and (iv) no DW (black), summarized in Table 1. Top: radial profiles of temperatures, T, at t = 0 (dotted lines), 105 (dashed lines), and 106 (solid lines) years. We note that the initial temperatures of the four cases are almost the same and that the red and green solid lines at t = 106 yr overlap at T = Treq (Eq. (27)). Bottom: radial profiles of surface densities, Σ, at t = 0 (dotted lines), 105 (long dashed lines), 106 (solid lines), and 107 (short dashed lines) years. We note that the radial range of the top panel is more zoomed-in than the radial range of the bottom panel.

In this subsection we show results of four cases of MRI-active PPDs, which are summarized in Table 1. The first three cases take disc winds into account. The magnetic braking by the disc winds is only considered in the first case. The last case (no DW) does not take disc winds into account by substituting ϵ_rad = 1 and C_w,0 = 0 in Eqs. (20)–(22).

Figure 1 compares radial profiles of T and Σ of these four cases. The top panel compares the evolution of the temperatures of these four cases. The initial temperature profiles in 0.1 au ≲ r ≲ 5 au, are kept more or less constant ≲1500−2500 K because dust grains sublimate and the opacity drops above that temperature (Eq. (25); see also Baillié et al. 2015). Furthermore, the initial profiles are almost the same for the four cases, except for different energetics constraints on C_w and wind torques, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$. In particular, the weak DW case (adopting Eqs. (20) and (21); purple dotted line) gives a very similar profile to those of the strong DW cases (adopting Eqs. (18) and (19); red and green dotted line), which needs explanation. In the inner region, ≲10 au, T ≈ T_vis (Eq. (28)) in these cases, and then T is mainly determined from F_rad by Eq. (23). Recalling Σ_int ∝ r^{− 3/2}, we derive $\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma \Omega }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{\approx }\frac{\mathrm{3}}{\mathrm{2}}\mathrm{\Sigma \Omega }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ for ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{~}{\mathit{r}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$. Since the $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$(0 or =10^-5) term in Eq. (21) is negligible in comparison to the $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{(}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-3}\mathrm{)}$ term, both strong DW and weak DW conditions give similar F_rad in Eqs. (19) and (21), and accordingly, the initial temperatures of these cases are similar each other.

In the no DW case (black lines) the viscous heating region (T_vis>T_req) survives until a later time although its size shrinks. In contrast, the temperatures decrease more rapidly in the other cases with disc winds. In the two strong DW cases (red and green lines), the temperatures are mainly determined by T_req in the entire region after t ≳ 10⁶ yr because the surface densities decrease rapidly by the disc winds in the inner region to give T_req ≫ T_vis, while T_vis is no longer negligible in the weak DW case (purple lines) at t = 10⁶ yr.

Fig. 2 Comparison of nondimensional mass flux of disc winds, Cw, of the three MRI-active cases except for the no DW case in Table 1. at t = 0 (dotted lines) and 106 yr (solid lines).

The bottom panel of Fig. 1 compares the evolution of the surface densities. The disc winds reduce Σ particularly in small r regions (Suzuki et al. 2010). A comparison between the two zero-torque cases (green and purple lines) shows the difference between the strong DW and weak DW conditions. As expected, the strong DW case shows faster decrease of Σ because of the higher disc wind mass flux, C_w, which is shown in Fig. 2. At t = 0 the strong DW case (green dotted line) gives quite small C_w ≈ 0 below the displayed range of Fig. 2 because $\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{\approx }\mathrm{0}$ for Σ_int ∝ r^{− 3/2} in Eq. (18). However, as Σ decreases in an inside-out manner and the Σ profile changes, C_w increases and at t = 10⁶ yr this case (green solid line) yields larger C_w than the weak DW case (purple solid line), in which C_w instead decreases with time owing to the decrease in temperature (∝${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$; Eq. (21)). We note that C_w = 0 in the outer region, r> 90 au, of the strong DW + zero-torque case because the gas moves outward ($\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathit{<}\mathrm{0}$ in Eq. (18)) in the outer region and the gravitation energy is not released. In realistic situations, however, a moderate level of external heating by stellar irradiation or other sources would cause disc winds to be launched by relaxing the energetics constraint (see Sect. 4.1), because the gas is only weakly bound by the gravity in the outer region.

The non-zero wind torque also reduces Σ faster (red lines in the bottom panel of Fig. 1) by the enhanced accretion and disc wind mass-loss. A comparison between the red and green lines in Fig. 2 indicates that the removal of angular momentum by the φz stress additionally contributes to the gravitational energy by the accretion to enhance C_w (Eq. (18)). As a result, C_w is not constrained by the energetics, C_w,e, in the almost entire region but is determined by C_w,0( = 2 × 10^-5) at t = 10⁶ yr (red solid lines). The constant C_w = C_w,0 implies faster dispersal of Σ for smaller r because the mass-loss timescale becomes proportional to the Keplerian time (Suzuki et al. 2010; Ogihara et al. 2015a,b), and the slope of Σ is positive in the inner region. The slope of Σ is again negative in the very inner region, r< 0.1 au, at later time, t ≳ 10⁷ yr. This is because $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ is constrained by the cap value = 1 (Eq. (30) there.

Fig. 3 Mass-loss rate by disc wind, Ṁz, (solid lines) and mass accretion rate induced by the rφ stress, Ṁr,rφ (dashed lines) and by the φz stress, Ṁr,rφ (dotted lines) at t = 106 yr of the four MRI-active cases in Table 1. Ṁz = 0 for the no DW case and Ṁr,φz = 0 for the zero-torque cases.

Figure 3 presents the radial profile of the mass-loss rate by disc winds (solid lines), $\mathit{Ṁ}\mathit{z}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathit{r}}^{{\mathit{r}}_{\mathrm{out}}}\mathit{r}\mathrm{d}\mathit{r}\mathrm{(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathit{r}}^{{\mathit{r}}_{\mathrm{out}}}\mathit{r}\mathrm{d}\mathit{r}{\mathit{C}}_{\mathrm{w}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}{\mathrm{)}}_{\mathrm{mid}}\mathit{,}$(31)and the mass accretion rate, $\mathit{Ṁ}\mathit{r}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi r}\mathrm{\Sigma }{\mathit{v}}_{\mathit{r}}\mathit{,}$(32)at t = 10⁶ yr. Here, Ṁ_r can be separated into two parts, Ṁ_r = Ṁ_r,rφ + Ṁ_r,φz, following Eq. (5) with help of Eqs. (6) and (8) (see also Simon et al. 2013): mass accretion induced by the rφ stress (dashed lines), $\mathit{Ṁ}\mathit{r,r\phi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi }}{\mathit{r}\mathrm{\Omega }}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}$(33)and that by the φz stress (dotted line) $\mathit{Ṁ}\mathit{r,\phi z}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{\Omega }}\mathit{r}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\mathit{.}$(34)We note that Ṁ_z(r) in our definition is the total mass loss outside r, while the disc wind mass loss at r is sometimes defined as the mass lost inside r (e.g., Owen et al. 2011; Bai et al. 2016; Bai 2016). We chose our definition to show how $\dot{{\mathit{M}}_{\mathit{r}}}$ is converted into Ṁ_z as mass accretes inward.

Fig. 4 Time evolution of Ṁz (solid) and Ṁr = Ṁr,rφ + Ṁr,φz (dashed) at r = 0.0225 au of the four MRI-active cases in Table 1.

The no DW case shows a spatially uniform accretion rate, Ṁ_r,rφ = 1.5 × 10^-8M_⊙ yr^-1 in r< 10 au (black dashed line). When disc winds are taken into account, the mass accretion rate decreases with decreasing r as the mass is lost by the disc winds. When we evaluate Ṁ at r = 0.0225 au (≈4.8 R_⊙), which is one grid point outside r_in = 0.01 au and approximately twice the radius of typical T Tauri stars, Ṁ_r,rφ is reduced to 2.5 × 10^-9M_⊙ yr^-1 in the weak DW case (purple dashed line). Instead, the mass is largely lost by the disc winds, Ṁ_z = 1.0 × 10^-8M_⊙ yr^-1 at r = 0.0225 au (purple solid line). This situation is more drastic in the strong DW + zero-torque case, and Ṁ_z ≈ 100Ṁ_r,rφ (green lines) at r = 0.0225 au. We note that Ṁ might have to be evaluated at a slightly outer location when the inner disc is truncated by the magnetosphere of the central star (e.g., Shu et al. 1994; Hirose et al. 1997; Dyda et al. 2015); in this case, Ṁ_r is not as small as the above evaluated values.

The strong DW + Σ-dependent torque case (red lines) gives very small Ṁ_r,rφ at r = 0.0225 au because Σ is small there (Fig. 1). On the other hand, the accretion by the φz stress is non-zero only in this case of the four cases displayed in Fig. 3, and Ṁ_r,φz is still kept = 2.4 × 10^-9M_⊙ yr^-1 at r = 0.0225 au because $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ increases to ≈0.1 in the inner region; the disc is in a wind-driven accretion phase.

Figure 4 compares the time evolutions of Ṁ_z (solid) and Ṁ_r = Ṁ_r,rφ + Ṁ_r,φz (dashed) at r = 0.0225 au of these four cases. The obtained t−Ṁ_r trends can be directly compared to the observed distribution in the t−Ṁ_r plane (Gullbring et al. 1998; Hartmann et al. 1998; Ricci et al. 2010; Manara et al. 2016). Although Ṁ_r of the strong DW + zero-torque case is smaller than the observed lower edge (Ṁ_r ~ 10^-9M_⊙ yr^-1 at t = 10⁶ yr), Ṁ_r of the other three cases are well inside the observed range.

3.2. MRI-inactive cases

Fig. 5 Same as Fig. 1 but for four MRI-inactive cases with $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-5}$. The four cases are (i) strong DW + Σ-dependent torque (red); (ii) weak DW + Σ-dependent torque (blue); (iii) strong DW + constant torque (grey); and (iv) strong DW + zero-torque (green), summarized in Table 2. The initial temperatures of the three strong DW cases (red, grey, and green dotted lines) are the same and the red and grey solid lines at t = 106 yr overlap at T = Treq (Eq. (27)).

We present results of four MRI-inactive cases, which are summarized in Table 2. We focus on effects of the wind torque on the evolution of PPDs in this subsection. Figure 5 compares radial profiles of T and Σ. The temperatures (top panel) of these cases are systematically lower than the temperatures of the MRI-active cases (the top panel of Fig. 1) because smaller $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ gives smaller F_rad (Eqs. (19) and (21)) and accordingly lower T_vis (Eq. (23)).

Fig. 6 Same as Fig. 2 but for the MRI-inactive cases of Table 2.

Smaller $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ also leads to slower evolution; when the MRI-active and MRI-inactive cases are compared, which adopt the same strong DW + zero-torque parameters (green lines in Figs. 1 and 5), the decrease of Σ is much slower in the MRI-inactive case. This is first because the accretion itself is slower owing to the smaller $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ and second because the disc wind mass flux is strongly constrained by the energetics to give smaller C_w (Fig. 6 in comparison to Fig. 2).

Fig. 7 Comparison of $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ at t = 0 (dotted lines), 105 (long dashed lines), 106 (solid lines), and 107 (short dashed lines) years of the non-zero-torque cases of Table 2.

The evolution of Σ is largely affected by non-zero wind torque $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$, because its effect is relatively important for lower turbulent viscosity, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$. The addition of the spatially constant $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{=}{\mathrm{10}}^{-4}$ (constant torque, grey lines) greatly reduces Σ. The two Σ-dependent torque cases (red and blue lines) give positive slopes of Σ in the inner region, which we explain below.

Figure 7 presents the time evolution of $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ for the Σ-dependent torque cases. $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ increases with time from the inside to the outside as Σ decreases in an inside-out manner. As a result, the disc wind mass flux, C_w, is not constrained by the energetics (Eqs. (18) and (20)) but is chosen to be the constant C_w,0( = 10^-5) in Eq. (22) (Fig. 6), which leads to the inside-out dispersal of the gas. In addition, the accretion is faster for smaller r because $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ is larger for smaller r. The positive slopes of Σ can be explained by the combination of these effects.

Although we assumed spatially uniform $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ and C_w,0, they are also anticipated to depend on the strength of net vertical magnetic field, ${\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{8}\mathit{\pi }\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}$. In this case, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ and C_w,0 could inversely correlate with Σ (Suzuki et al. 2010), which additionally enhances the positive slopes of the surface densities.

Fig. 8 Same as Fig. 3 but for the MRI-inactive cases of Table 2.

Fig. 9 Same as Fig. 4 but for the MRI-inactive cases of Table 2.

Figure 8 compares Ṁ_z (Eq. (31); solid), Ṁ_r,rφ (Eq. (33); dashed), and Ṁ_r,φz (Eq. (34); dotted) of the MRI-inactive four cases at t = 10⁶ yr. In the zero-torque case (green lines) the mass is dominantly lost by the disc winds, Ṁ_z ≈ 100Ṁ_r,rφ at r = 0.0225 au. In the constant torque case (grey lines) the mass accretion is mainly driven by the φz stress, the total accretion rate is also largely dominated by the mass loss by the disc winds. Adopting the Σ-dependent torque condition changes the situation; the mass accretion is driven by the φz stress, and the accretion rate is well above 10^-9M_⊙ yr^-1, that is, the weak DW case gives Ṁ_r,φz ≈ Ṁ_z ≈ 5 × 10^-9M_⊙ yr^-1 at r = 0.0225 au.

Figure 9 shows the time evolution of Ṁ_z (solid) and Ṁ_r = Ṁ_r,rφ + Ṁ_r,φz (dashed) at r = 0.0225 au of the same four cases. Ṁ_r(0.0225 au) of the cases of zero or constant torque (green and grey dashed lines) are smaller than the observed range of t−Ṁ_r (Gullbring et al. 1998; Hartmann et al. 1998; Ricci et al. 2010; Manara et al. 2016). On the other hand, Ṁ_r(0.0225au)’s of the Σ-dependent torque cases (red and blue dashed lines) are consistent with the observed t−Ṁ_r. Although the mass accretion rate of the strong DW case is lower than the wind mass-loss rate (red lines), it is not so small; Ṁ_r(0.0225 au) = 6.0 × 10^-9M_⊙ yr^-1 at t = 10⁵ yr, 1.7 × 10^-9M_⊙ yr^-1 at 10⁶ yr, and 5.2 × 10^-10M_⊙ yr^-1 at 10⁷ yr.

4. Discussion

4.1. Uncertainties

Our model has the three free parameters, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, C_w,0, and $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$. Since these parameters are not yet tightly constrained by observations or theoretical calculations, we calculated the evolution of PPDs in the wide ranges of the parameters to test various possibilities (Sect. 3). Uncertainties of the three parameters is largely attributed to the uncertainty of the initial distribution and to the evolution of the poloidal magnetic flux because these three parameters depend on the vertical magnetic field strength (Suzuki et al. 2010; Okuzumi & Hirose 2011; Bai & Stone 2013b).

The evolution of poloidal magnetic flux in accretion discs has been studied by a number of groups (Lubow et al. 1994; Rothstein & Lovelace 2008; Guilet & Ogilvie 2012; Suzuki & Inutsuka 2014) and has recently been specifically applied to PPDs (Okuzumi et al. 2014; Guilet & Ogilvie 2014; Takeuchi & Okuzumi 2014). Accreting gas drags the vertical magnetic field inward, while the vertical field also possibly diffuses outward by magnetic diffusivity, which consists of both effective turbulent resistivity and non-ideal MHD effects (Sect. 2.6). The radial motion of the vertical magnetic flux is determined by the balance between these inward dragging and outward diffusion. The direction of the magnetic flux itself is still uncertain, which depends on the initial configuration of the poloidal magnetic field, in addition to the combination of accretion and magnetic diffusion.

One future possibility is that we finally obtain a universal tendency for the time evolution of vertical magnetic fields. In this case, we can constrain our free parameters, and evolutions of surface densities will not show a variety but converge to a unified trend. On the other hand, if the evolution of the poloidal magnetic flux is different in different PPDs, depending on physical circumstances, such as initial magnetic flux and disc mass, and stellar irradiation, which controls the non-ideal MHD effects through the ionization, then the evolutions of surface densities are also different in different PPDs as shown so far, which should lead to a wide variety of the subsequent planet formation processes and final exoplanet systems.

At present, the unified picture of the evolution of the poloidal magnetic field is not well understood at all, and therefore it is worth pursuing various possibilities. Our calculations took the effect of the evolution of the vertical magnetic field in the wind torque into account; the two cases of constant $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ and Σ-dependent $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ correspond to the case in which the magnetic energy decreases in the same manner as the decrease of the surface density and the case with the preserved magnetic flux, respectively. The Σ-dependent torque cases show a runaway behavior of the gas dispersal in an inside-out manner; once the gas is dispersed, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ increases, which further accelerates the dispersal of the gas. This is the main reason why the positive slope of Σ is produced. Although we did not consider this effect, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$ and C_w,0 depend similarly on Σ, which causes an additional runaway dispersal of the gas (Suzuki et al. 2010, see also Sect. 3.2). The case with constant $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ even gives the moderately positive slope (Fig. 5). Within the two cases we tested, the positive slope of Σ on r is not peculiar, but a common feature. However, we should note that our calculations do not cover all the possible distributions and evolutions of the vertical magnetic field. Therefore, it would be premature to conclude that the positive slope of Σ is a natural outcome of the accretion induced by the magnetically driven disc wind. For example, when the outward diffusion of vertical magnetic field is effective and the magnetic flux is dispersed more rapidly than the gas, the effect of the wind torque is suppressed with time. In this case, the Σ profile would maintain a normal negative slope.

We now discuss other ambiguities of the mass flux of the disc winds, in addition to the uncertainty of the vertical magnetic field. At the moment, the mass flux, C_w,0, is available only from local MHD simulations (e.g. Suzuki & Inutsuka 2009; Fromang et al. 2013; Bai & Stone 2013a). As discussed in Sect. 2.6, these local simulations may overestimate the mass flux. Although we adopted the conservative C_w,0 by reducing the simulation results by half (see Sect. 2.6), it might be even lower (Fromang et al. 2013). We here briefly discuss how the results are affected and particularly focus on the slope of the surface density when C_w,0 is smaller.

As shown in Figs. 2 and 6,C_w is already constrained by the energetics. In most cases except for the MRI-inactive cases with Σ-dependent torque, the energetics constraint already suppresses C_w in the inner region. Therefore, adopting a smaller C_w,0 does not affect C_w in the inner region but reduces C_w in the outer region, which suppresses the gas dispersal there. Hence, the slope of Σ would be more positive in these cases. On the other hand, in the MRI-inactive cases with Σ-dependent torque, the energetics constraint suppresses C_w at the relatively outer location, r ~ 10 au. In these cases, a smaller C_w,0 reduces C_w in the inner region. As a result, the obtained large positive Σ slopes in these cases (Fig. 6) would be reduced to moderately positive ones.

When we determined the mass flux of the disc winds, we applied the energetics constraint from the gravitational accretion without external heating or momentum inputs (Sect. 2.3; Eq. (22)). This treatment is expected to give a reasonable constraint at the early phase when viscous heating dominates the radiative heating or other effects from the central star. However, at the later time this is not the case because the surface density decreases and the viscous heating becomes relatively unimportant. Effects of external heating or momentum inputs need to be considered. They weaken the energetics constraint to give a larger C_w in the region with C_w,e<C_w,0 (see Sect. 4.5).

4.2. Radial drift of pebbles and boulders

Although calculations still include uncertainties that mainly stem from the ambiguity of the evolution of poloidal magnetic fields, the positive slopes of the surface densities obtained in Sect. 3 are a possible consequence of the evolution of PPDs with disc winds, as discussed in Sect. 4.1. These positive slopes raise various interesting implications for planet formation. In this and the next subsections, we demonstrate how the obtained Σ profiles affect the solid component of PPDs by studying cases that show large positive slopes of Σ.

The first example is the radial drift of solid bodies through gas drag. In general the rotation velocity of the gas in PPDs is slightly slower than the local Keplerian velocity because of the radial pressure gradient force. On the other hand, solid particles rotate with Keplerian velocity without the support from the gas pressure. As a result, the solid particles feel a head wind from the gas, which causes them to drift inward. Considering the momentum balance, solid particles with nondimensional stopping time ≈1, which corresponds to a meter-sized spherical boulder at 1 au of the MMSN, experience the radial drift most severely (Weidenschilling 1977; Nakagawa et al. 1986), and their drift timescale in the midplane is given by ${\mathit{\tau }}_{\mathrm{dr}\mathit{,}\max }\mathrm{\approx }\frac{\mathrm{1}}{\mathit{\eta }{\mathrm{\Omega }}_{\mathrm{K}}}\mathit{,}$(35)where η is pressure gradient force normalized by the twice of centrifugal force, $\mathit{\eta }\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{mid}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathrm{mid}}}{\mathit{\partial r}}\frac{\mathrm{1}}{\mathrm{2}\mathit{r}{\mathrm{\Omega }}_{\mathrm{K}}^{\mathrm{2}}}\mathrm{·}$(36)In the usual condition, η ~ 10^-3−10^-2> 0, which causes solid particles to drift inward. Smaller η leads to slower inward drift; if η< 0, the direction of the drift is opposite and solid particles move outward.

Figure 10 shows η of the two MRI-inactive ($\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-5}$) cases with Σ-dependent torque of Table 2 (red and blue lines; the same as in Figs. 5–9) in comparison to the no disc wind (no DW) case with the same $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-5}$ (black lines). We here derive p_mid from Σ by ${\mathit{p}}_{\mathrm{mid}}\mathrm{=}{\mathit{\rho }}_{\mathrm{mid}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{\Sigma \Omega }{\mathit{c}}_{\mathrm{s}}}{\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{·}$(37)The no DW case shows η remains within 10^-3−10^-2, which implies fast inward drift. In contrast, η’s are considerably reduced in the Σ-dependent torque cases. In particular, the red lines (strong DW case) show negative η in part (red lines are truncated between 0.04–0.4 au at t = 10⁵ yr and 1–2 au at t = 10⁶ yr), which indicates that solid particles move outward in this region. As a result, the solid component will accumulate around the outer edge of the negative η region, which offers suitable conditions for planet formation (Kobayashi et al. 2012). Furthermore, this location moves outward with time; the suitable site for the planet formation also moves outward.

Fig. 10 Comparison of normalized pressure gradient force, $\mathrm{-}\left(\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{mid}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathrm{mid}}}{\mathit{\partial r}}\right)\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{r}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\right)}$, of MRI-inactive PPDs at t = 0 (dotted), 105 (solid), and 106 yr (dashed). The MRI-inactive cases with Σ-dependent torque in Table 2, blue lines for weak DW and red lines for strong DW, which corresponds to the red and blue lines in Fig. 5, are compared to the MRI-inactive no DW case with Cw,0 = 0 and $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{-5}$ (black lines).

4.3. Type I migration

Fig. 11 Migration efficiency for Earth-mass planets for MRI-inactive cases with Σ-dependent torque (red for strong DW and blue for weak DW of Table 2) at t = 106 yr (corresponding to the solid red and blue lines in Fig. 5) in comparison to the MRI-inactive no DW case (black). CI > 0 means outward migration.

Another interesting implication of the positive Σ slopes is that an inward migration of low-mass planets (type I migration) can be slowed down or even reversed. The torque for type I migration can be expressed by the sum of Lindblad and corotation torques. The corotation torque is more sensitive to the slope of the gas surface density and can be positive for positive slopes.

Here we estimate the migration rate of Earth-mass planets embedded in MRI-inactive PPDs with the surface densities shown in Fig. 5. We used the formulae of Paardekooper et al. (2011) to calculate the migration timescale, t_a (see Eqs. (8)–(16) in Ogihara et al. (2015a) for details of the formulae). We introduced a parameter of the efficiency of inward type I migration, C_I ≡ −t_a,TTW/t_a, where t_a,TTW is the migration time in a locally isothermal disc derived by a linear analysis by Tanaka et al. (2002). The migration timescale is defined as t_a ≡ a/ (−ȧ); positive migration time means inward migration.

Figure 11 shows the migration efficiency for the Σ-dependent torque case (the red and blue curves in Fig. 5) at t = 10⁶ yr in comparison to the no DW case (black line). The migration rate depends on the planetary mass and the orbital eccentricity; Earth-mass planets with zero eccentricity were considered here. The blue curve shows that the type I migration is slowed down inside a few au by several factors from t_a,TTW. The migration is even reversed (outward migration) between 0.1–0.5 au in the red curve (strong DW case). Thus the disc wind would also play important roles in the late stage of planet formation.

4.4. Comparison to previous work

Recently, Bai (2016) also presented a global evolution model for PPDs with magnetically driven disc winds. However, none of the cases in his model calculations resulted in a surface density with a drastic positive slope relative to r as some of our cases have shown. The two main differences between his setup and ours is the mass-loss rate by the disc wind and the evolution of the vertical magnetic field.

Our calculations, which started from a relatively massive initial disc (M_disc,int = 0.11 M_⊙) to study the evolution from the early stage, neglected the heating by the irradiation from a central star but considered viscous heating, and the mass-loss rate was constrained by the global energetics of the viscous accretion. In contrast, the initial disc mass adopted by Bai (2016) is lower, =0.035 M_⊙, to focus on the later stage of the evolution, and the location of the wind base in the inner region r ≲ 10−30 au is determined from heating by far-ultraviolet (FUV) irradiation from a central star. Here, the penetration depth of the FUV was assumed to be spatially constant. Since the surface density decreases with r, the penetration depth normalized by the scale height is deeper for larger r. Therefore, the mass loss by the disc wind affects the depletion of the gas at outer locations more severely than in our model setting, and consequently a positive slope of Σ was not obtained in the results of Bai (2016).

As for the evolution of the vertical magnetic field, Bai (2016) considered two cases: in the first case the total magnetic flux is preserved with time, and in the second case it decreases in the same manner as the total mass. In both cases, the plasma $\mathit{\beta }\mathrm{=}\mathrm{(}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{8}\mathit{\pi }\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}{\mathrm{)}}^{-1}$ at the midplane was assumed to be spatially uniform. Even in the first case, the vertical magnetic field was redistributed to follow the density profile (Armitage et al. 2013). This spatially uniform β was also adopted in our constant torque setting. In contrast, our Σ-dependent torque assumed the preserved vertical magnetic field at each location, which led to a runaway inside-out dispersal and produced a large positive slope of Σ (Sect. 3), compared to the above-mentioned cases with the spatially uniform β.

4.5. Stellar wind and photoevaporation

We did not take the effects of a central star into account except to determine the radiative equilibrium temperature, T_req (Eq. (27)). However, the stellar wind and irradiation affect the evolution of PPDs.

In our calculations, the mass flux of the disc wind is C_w,e constrained by the energetics of accretion, and it can be smaller than C_w,0 determined by the mass loading expected from the local MHD simulations. When this is the case, gaseous clouds are lifted up by vertical upflows but cannot stream out to large z; they float in the disc atmosphere or return to the disc because they are bound by the gravity of the central star. The stellar wind from the central star would change this situation.

The mass flux of the stellar wind from pre-main sequence stars is much higher, by an order of 4–6, than that of the current solar wind partly because of the energy supply from accretion (Hirose et al. 1997; Matt & Pudritz 2005; Cranmer 2009). Even after the accretion terminates, the mass flux of the stellar wind is expected to be still high because of the high magnetic activity (Wood et al. 2005; Cranmer & Saar 2011; Suzuki et al. 2013). The strong stellar wind would blow away the clouds that are lifted up by the disc winds (see Suzuki et al. 2010, for the energetics). In the framework of our model, the contribution from the stellar wind would increase C_w,e in Eqs. (16) and (17), in the small r region. The increase of C_w in the inner region reduces Σ there, which also produces a larger positive slope of Σ.

In this discussion, we neglected the roles of global magnetic fields that are rooted in the central star and in the PPD. When the field strength is strong enough, the stellar wind region and the disc wind region are separated by a boundary layer formed by magnetospheric ejections (Zanni & Ferreira 2013). In this case the stellar winds will not contribute to driving the disc winds. It depends on the relative strength of the magnetic energy to the sum of the dynamic pressure and the gas pressure whether the interaction between the stellar winds and the disc winds is efficient. When the magnetic energy is weaker, the interaction is stronger, and vice versa.

Photoevaporation by irradiation from the central star or neighbouring stars has been extensively studied as a viable source for dispersing PPDs (e.g., Shu et al. 1993; Hollenbach et al. 2000; Adams et al. 2004). The mass-loss rate by the photoevaporation, which depends on the flux in different spectral ranges, FUV, extreme UV, and X-rays, yields a wide variety of ~10^-10−10^-8M_⊙ yr^-1 (Alexander et al. 2006; Ercolano et al. 2008; Gorti & Hollenbach 2009; Owen et al. 2010; Tanaka et al. 2013). After the mass accretion rate or the mass-loss rate by the disc wind decreases below this level, the photoevaporation would quickly disperse PPDs (e.g. Armitage 2011); our results would be affected at the late stage of the evolution.

However, we expect that the evolution of the Σ profile of a photoevaporating PPD is qualitatively different from our results with the magnetically driven disc wind because the photoevaporation mostly affects the disc dispersal in the outer region where the sound speed of the heated gas exceeds the local escape velocity from the central star. Although the photoevaporation could create an inner hole by the combination with the viscous accretion, the local slope of Σ remains negative except at the inner edge of the hole (e.g. Alexander et al. 2006; Owen et al. 2011). This is in clear contrast to the evolution with the magnetically driven disc wind.

5. Summary

We have studied the global evolution of PPDs by considering viscous heating and magnetically driven disc winds. We constructed a global model from fundamental MHD equations for the time-evolution of PPDs. One of the key features of our model is that the mass-loss rate by the disc wind is derived from both the local MHD shearing box simulations and the global energetics of the gravitational accretion. Our model has three dimensionless parameters, which are turbulence viscosity, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, disc wind mass flux, C_w, and disc wind torque, $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$, and these three parameters are constrained by the above-mentioned global energetics. We performed model calculations in a wide parameter range to cover both MRI-active PPDs and MRI-inactive PPDs with dead zones.

We started our calculations from the relatively massive disc, M_disc.int = 0.11 M_⊙. Initially, the viscous heating dominantly determines the temperature in the inner region <10 au; for instance, T ≃ 1500 K at 1 au, which is much higher than the temperature estimated from the radiative equilibrium. As the surface density decreases with time, the temperature approaches the radiative equilibrium temperature. In the cases that consider the disc wind mass loss, the gas in the inner region is rapidly dispersed before 10⁶ yr, and the viscous heating is negligible in determining the temperature after t ≳ 10⁶ yr, whereas in the no disc wind cases the viscous heating is not negligible even up to several 10⁶ yr.

The mass accretion rates decrease with time as the surface densities decrease, regardless of whether the accretion is induced by turbulent viscosity or wind torque. The obtained accretion rates are consistent with observed accretion rates for a wide range of the adopted parameters.

The three free parameters, $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}$, C_w,0, and $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ still contain ambiguities, arising mainly from the uncertainty of the evolution of vertical magnetic fields. We have pursued various possibilities by testing different combinations of these parameters. The detailed profiles of the temperatures and the surface densities show a wide variety. Since physical properties of a PPD affect planet formation processes that take place in the disc (e.g., Kobayashi et al. 2016), the obtained variety of our PPD calculations would be a source of the observed variety of exoplanet systems (e.g. Howard et al. 2012).

The wind-driven accretion can promote an increase in disc surface density with r in the inner region; this is the case in our calculations for MRI-inactive PPDs when the distribution of the vertical magnetic flux is preserved with time evolution (Sect. 3.2). This large positive slope of the surface density suppresses or reverses the inward drift of pebble- or boulder-sized solids through gas drag (Sect. 4.2) and the inward migration of protoplanets (Sect. 4.3), which is a favourable condition for the planet formation.

Acknowledgments

T.K.S. is supported by the Astrobiology Center Project of National Institute of Natural Sciences (NINS) (Grant Number AB271020, AB281018). A.M., A.C. and T.G. acknowledge support by the French ANR, project number ANR-13–13-BS05-0003-01 project MOJO (Modeling the Origin of JOvian planets). T.K.S. thanks Hiroshi Kobayashi and Shinsuke Takasao for fruitful discussions. The authors thank the referee for many valuable comments.

References

Appendix A: Derivation of the equation for the surface density

In this appendix, we derive Eq. (10) from the conservation equations for angular momentum and mass. Under the axisymmetric approximation, a general MHD expression of the conservation of angular momentum (e.g., Balbus & Hawley 1998) is

$\begin{array}{ccc}& & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}\mathit{\rho r}{\mathit{v}}_{\mathit{\phi }}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[{\mathit{r}}^{\mathrm{2}}\left(\mathit{\rho }{\mathit{v}}_{\mathit{r}}{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi }}\right)\right]\\ & & \end{array}$(A.1)

The azimuthal velocity, v_φ, is decomposed into the mean Keplerian flow and perturbation, ${\mathit{v}}_{\mathit{\phi }}\mathrm{=}\mathit{r}\mathrm{\Omega }\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}\mathit{.}$(A.2)We use the α prescription (Shakura & Sunyaev 1973) for the second and third terms of Eq. (A.1): $\begin{array}{ccc}\mathit{\rho }{\mathit{v}}_{\mathit{r}}{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi }}& \mathrm{=}& \mathit{\rho }{\mathit{v}}_{\mathit{r}}\mathit{r}\mathrm{\Omega }\mathrm{+}\mathit{\rho }\left({\mathit{v}}_{\mathit{r}}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi \rho }}\right)\\ & \mathrm{\equiv }& \end{array}$(A.3)and $\begin{array}{ccc}\mathit{\rho }{\mathit{v}}_{\mathit{\phi }}{\mathit{v}}_{\mathit{z}}\mathrm{-}\frac{{\mathit{B}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{z}}}{\mathrm{4}\mathit{\pi }}& \mathrm{=}& \mathit{\rho r}\mathrm{\Omega }{\mathit{v}}_{\mathit{z}}\mathrm{+}\mathit{\rho }\left(\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}{\mathit{v}}_{\mathit{z}}\mathrm{-}\frac{{\mathit{B}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{z}}}{\mathrm{4}\mathit{\pi \rho }}\right)\\ & \mathrm{\equiv }& \end{array}$(A.4)We integrate Eq. (A.1) along the vertical direction, z, with Eqs. (A.3) and (A.4) from the bottom surface to the top surface of a disc, and we have $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}\mathrm{\Sigma }{\mathit{r}}^{\mathrm{3}}\mathrm{\Omega }\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{[}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }{\mathrm{(}{\mathit{v}}_{\mathit{r}}\mathit{r}\mathrm{\Omega }\mathrm{+}\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{]}}& & \end{array}$$\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathrm{\left[}\mathrm{\right(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}\mathit{r}\mathrm{\Omega }\mathrm{+}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{\left(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\mathrm{\right]}\mathrm{=}\mathrm{0}\mathit{,}$(A.5)where $\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}\mathrm{=}{}^{\mathrm{\int }}\mathit{\rho }{\mathit{\alpha }}_{\mathit{r\phi }}\mathrm{d}\mathit{z}\mathit{/}\mathrm{\Sigma }$ is the mass-weighted vertical average. The third term, which represents the angular momentum loss from both surfaces, is derived from ${\mathrm{[}\mathit{\rho r}\mathrm{\Omega }{\mathit{v}}_{\mathit{z}}\mathrm{+}\mathit{\rho }{\mathit{\alpha }}_{\mathit{\phi z}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{]}}}_{\mathrm{w}}\mathrm{=}\mathrm{[}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{]}}_{\mathrm{w}}\mathit{r}\mathrm{\Omega }\mathrm{+}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathit{,}$(A.6)where the subscript w stands for disc wind. $\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}$ is the angular momentum loss by the φz component of the stress normalized by the density and the sound speed at the midplane, Eq. (8).

The equation of mass conservation is $\frac{\mathit{\partial }\mathrm{\Sigma }}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathit{r}\mathrm{\Sigma }{\mathit{v}}_{\mathit{r}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}\mathrm{=}\mathrm{0.}$(A.7)By combining Eq. (A.7) multiplied by r³Ω and Eq. (A.5), we have $\mathit{r}\mathrm{\Sigma }{\mathit{v}}_{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\mathrm{=}\mathrm{0}\mathit{,}$(A.8)which determines the accretion rate, $\mathit{r}\mathrm{\Sigma }{\mathit{v}}_{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{2}}{\mathit{r}\mathrm{\Omega }}\left[\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\right]\mathit{,}$(A.9)where we here assumed the Keplerian rotation to derive $\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}\frac{\mathit{r}\mathrm{\Omega }}{\mathrm{2}}$.

By substituting Eq. (A.9) into Eq. (A.7), we finally have the equation for the time evolution of Σ (Eq. (10)): $\begin{array}{ccc}\frac{\mathit{\partial }\mathrm{\Sigma }}{\mathit{\partial t}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\frac{\mathrm{2}}{\mathit{r}\mathrm{\Omega }}\left\{\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Sigma }\overline{{\mathit{\alpha }}_{\mathit{r\phi }}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\overline{{\mathit{\alpha }}_{\mathit{\phi z}}}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}_{\mathrm{mid}}\right\}\right]& & \end{array}$$\begin{array}{ccc} \mathrm{+}\mathrm{(}\mathit{\rho }{\mathit{v}}_{\mathit{z}}{\mathrm{)}}_{\mathrm{w}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$

Appendix B: Energetics of accretion discs

A general MHD expression of the total energy conservation under the axisymmetric approximation is $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{\mathit{p}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }}\right]& & \end{array}$$\begin{array}{ccc} \mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}[\mathit{r}\{{\mathit{v}}_{\mathit{r}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{p}\mathrm{+}\frac{{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\right)& & \end{array}$$\begin{array}{ccc} \mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}}{\mathrm{4}\mathit{\pi }}\mathrm{(}{\mathit{v}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{\phi }}\mathrm{+}{\mathit{v}}_{\mathit{z}}{\mathit{B}}_{\mathit{z}}\mathrm{)}\mathrm{+}{\mathit{F}}_{\mathrm{ot}\mathit{,r}}\}]& & \end{array}$$\begin{array}{ccc} \mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial z}}[{\mathit{v}}_{\mathit{z}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{p}\mathrm{+}\frac{{\mathit{B}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\right)& & \end{array}$$\mathrm{-}\frac{{\mathit{B}}_{\mathit{z}}}{\mathrm{4}\mathit{\pi }}\mathrm{(}{\mathit{v}}_{\mathit{r}}{\mathit{B}}_{\mathit{r}}\mathrm{+}{\mathit{v}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{\phi }}\mathrm{)}\mathrm{+}{\mathit{F}}_{\mathrm{ot}\mathit{,z}}]\mathrm{=}\mathrm{0}\mathit{,}$(B.1)where we refer to Eq. (11) for the definition of each variable. Decomposing v_φ by Eq. (A.2) and assuming $\mathit{r}\mathrm{\Omega }\mathrm{\gg }{\mathit{v}}_{\mathit{r}}\mathit{,\delta }{\mathit{v}}_{\mathit{\phi }}\mathit{,}{\mathit{v}}_{\mathit{z}}\mathit{,}{\mathit{c}}_{\mathit{s}}\mathit{,B}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$ in a disc, we rewrite Eq. (B.1) with leaving dominant terms. The time-derivative term becomes $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\frac{\mathit{p}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }}\right]\mathrm{\approx }\frac{\mathit{\partial }}{\mathit{\partial t}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\right]& & \end{array}$$\mathrm{\approx }\frac{\mathit{\partial }}{\mathit{\partial t}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }\mathrm{(}\mathit{r}\mathrm{\Omega }\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathit{\rho }{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\right]\mathrm{\approx }\frac{\mathit{\partial }}{\mathit{\partial t}}\left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\right)\mathit{,}$(B.2)where we set rΩδv_φ = 0 after the azimuthal average. The r-derivative term, except for F_ot,r, can be approximated as $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial r}}[\mathit{r}\{{\mathit{v}}_{\mathit{r}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\mathrm{+}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{p}\mathrm{+}\frac{{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\right)& & \end{array}$$\begin{array}{ccc} \mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}}{\mathrm{4}\mathit{\pi }}\mathrm{(}{\mathit{v}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{\phi }}\mathrm{+}{\mathit{v}}_{\mathit{z}}{\mathit{B}}_{\mathit{z}}\mathrm{)}\}]& & \end{array}$$\begin{array}{ccc} \mathrm{\approx }\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\mathit{r}\left\{{\mathit{v}}_{\mathit{r}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{\rho }\mathrm{\Phi }\right)\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}}{\mathrm{4}\mathit{\pi }}{\mathit{v}}_{\mathit{\phi }}{\mathit{B}}_{\mathit{\phi }}\right\}\right]& & \end{array}$$\begin{array}{ccc} \mathrm{\approx }\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\mathit{r}\left\{\mathrm{-}\mathit{\rho }{\mathit{v}}_{\mathit{r}}\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{\rho r}\mathrm{\Omega }\left({\mathit{v}}_{\mathit{r}}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}\mathrm{-}\frac{{\mathit{B}}_{\mathit{r}}{\mathit{B}}_{\mathit{\phi }}}{\mathrm{4}\mathit{\pi \rho }}\right)\right\}\right]& & \end{array}$$\mathrm{=}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\mathit{r}\left\{\mathrm{-}\mathit{\rho }{\mathit{v}}_{\mathit{r}}\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{\rho r}\mathrm{\Omega }{\mathit{\alpha }}_{\mathit{r\phi }}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\right\}\right]\mathit{,}$(B.3)where the second ≈ is derived from ${\mathit{v}}_{\mathit{r}}{}^{\mathrm{(}}\frac{{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}{\mathrm{\Phi }}^{\mathrm{)}}\mathrm{\approx }{\mathit{v}}_{\mathit{r}}\left[\frac{\mathrm{(}\mathit{r}\mathrm{\Omega }\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\right]\mathrm{\approx }\mathrm{-}{\mathit{v}}_{\mathit{r}}\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{\rho r}\mathrm{\Omega }{\mathit{v}}_{\mathit{r}}\mathit{\delta }{\mathit{v}}_{\mathit{\phi }}$, and for the last equality Eq. (A.3) is used. We set the z-derivative term, except for F_ot.z, to be