Free Access

This article has an erratum: [https://doi.org/10.1051/0004-6361/201628955e]


Issue
A&A
Volume 596, December 2016
Article Number A74
Number of page(s) 15
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201628955
Published online 02 December 2016

© 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ΩΩK=GMr3,\begin{equation} \Omega \approx \Omega_{\rm K} = \sqrt{\frac{G M_{\star}}{r^{3}}}, \label{eq:Keprot} \end{equation}(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 H=2csΩ,\begin{equation} H = \frac{\sqrt{2}c_{\rm s}}{\Omega}, \label{eq:sclh} \end{equation}(2)where cs is the sound speed. Temperature T and cs are related through cs2=kBTμmH,\begin{equation} c_{\rm s}^2 = \frac{k_{\rm B} T} {\mu m_{\rm H}}, \label{eq:cs} \end{equation}(3)where kB is the Boltzmann constant, mH is the proton mass, and we assume mean molecular weight, μ = 2.34 (Hayashi 1981). A different definition for the scale height from ours, cs/ Ω (without 2\hbox{$\sqrt{2}$}), is sometimes used in literatures.

2.2. Evolution of surface density

We treated the time evolution of the radial profile of surface density, Σ = d, 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)

Σ ∂t 1 r ∂r [ 2 r Ω { ∂r ( r 2 d z ( ρ v r δ v φ B r B φ 4 π ) ) \begin{eqnarray} &&\frac{\partial \Sigma}{\partial t} - \frac{1}{r}\frac{\partial}{\partial r} \left[\frac{2}{r\Omega}\left\{\frac{\partial}{\partial r}\left(r^2 \int {\rm d}z \left(\rho v_r \delta v_{\phi} - \frac{B_rB_{\phi}}{4\pi}\right)\right) \right.\right. \nonumber\\ && \left.\left. \qquad\qquad + r^2 \left(\rho \delta v_{\phi} v_z - \frac{B_{\phi}B_z}{4\pi}\right)_{\rm w} \right\}\right] + (\rho v_z)_{\rm w} = 0, \label{eq:sgmevlori} \end{eqnarray}(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,

r Σ v r = 2 r Ω { ∂r ( r 2 d z ( ρ v r δ v φ B r B φ 4 π ) ) \begin{eqnarray} &&-r \Sigma v_r = \frac{2}{r\Omega}\left\{\frac{\partial}{\partial r} \left(r^2 \int {\rm d}z\left(\rho v_r \delta v_{\phi} - \frac{B_rB_{\phi}}{4\pi}\right)\right) \right. \nonumber\\ &&\left. \qquad\qquad+ r^2 \left(\rho \delta v_{\phi} v_z - \frac{B_{\phi}B_z}{4\pi}\right)_{\rm w} \right\}, \label{eq:rdmsfx} \end{eqnarray}(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 and φz components of Reynolds and Maxwell stresses. The 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): dz(ρvrδvφBrBφ4π)dαcs2Σαcs2,\begin{equation} \int {\rm d}z \left(\rho v_r \delta v_{\phi} - \frac{B_rB_{\phi}}{4\pi}\right) \equiv \int {\rm d}z \rho \alpha_{r\phi} c_{\rm s}^2 \equiv \Sigma \overline{\alpha_{r\phi}} c_{\rm s}^2, \label{eq:rphistress} \end{equation}(6)where α\hbox{$\overline{\alpha_{r\phi}}$} is the mass-weighted vertical average of α. α is a nondimensional parameter normalized by gas pressure (ρcs2\hbox{$\rho c_{\rm s}^2$}) that describes the transport of angular momentum. We considered α to originate from the MHD turbulence induced by MRI. α(1)\hbox{$\overline{\alpha_{r\phi}}({\lesssim} 1)$} 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 (ρvrδvφ) and from the Maxwell stress (BrBφ/ 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) .

α\hbox{$\overline{\alpha_{r\phi}}$} is an effective turbulent viscosity, and it is mathematically related to viscosity, ν, appeared in a hydrodynamical equation, αcs2=νrΩ∂r32νΩ,\begin{equation} \overline{\alpha_{r\phi}} c_{\rm s}^2 = -\nu r \frac{\partial \Omega}{\partial r} \approx \frac{3}{2}\nu\Omega, \label{eq:viscosity} \end{equation}(7)where the second equality comes from the condition of the Keplerian rotation. The definition of α is not consistent throughout the literature; for example, ναtHcs, is often used conventionally (e.g., Balbus & Hawley 1998). These two α’s are related by αt23α\hbox{$\alpha_{\rm t} \approx \frac{\sqrt{2}}{3}\overline{\alpha_{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φvz) and Maxwell (BφBz/ 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 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 (tr) model, αφz needs to be evaluated by physical quantities at the midplane, and we adopted a similar parametrization to the component, (ρδvφvzBφBz4π)w(ρcs2αφz)w(ρcs2)midαφz,\begin{equation} \left(\rho \delta v_{\phi} v_z - \frac{B_{\phi} B_z}{4\pi}\right)_{\rm w} \equiv (\rho c_{\rm s}^2\alpha_{\phi z})_{\rm w} \equiv (\rho c_{\rm s}^2)_{\rm mid} \overline{\alpha_{\phi z}}, \label{eq:phizstress} \end{equation}(8)where we define the nondimensional stress, αφz\hbox{$\overline{\alpha_{\phi z}}$}, normalized by density, ρmid(=Σ/(πH))\hbox{$\rho_{\rm mid} (=\Sigma /(\!\!\sqrt{\pi} H))$}, at the midplane,

The third term, (ρvz)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: (ρvz)w=Cw(ρcs)mid.\begin{equation} (\rho v_z)_{\rm w} = C_{\rm w}(\rho c_{\rm s})_{\rm mid}. \label{eq:Cw} \end{equation}(9)We model Cw in Sect. 2.3.

Substituting Eqs. (6), (8), and (9) into Eq. (4), we finally have Σ∂t1r∂r[2rΩ{∂r(r2Σαcs2)+r2αφz(ρcs2)mid}]\begin{eqnarray} \frac{\partial \Sigma}{\partial t}& -& \frac{1}{r}\frac{\partial}{\partial r} \left[\frac{2}{r\Omega}\left\{\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}} (\rho c_{\rm s}^2)_{\rm mid} \right\}\right] \nonumber\\ &&\qquad\qquad \qquad \qquad\, \qquad+ C_{\rm w}(\rho c_{\rm s})_{\rm mid} = 0. \label{eq:sgmevl} \end{eqnarray}(10)We solved this equation for different sets of the three parameters, α\hbox{$\overline{\alpha_{r\phi}}$}, αφz\hbox{$\overline{\alpha_{\phi z}}$}, and Cw. 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, Cw, is constrained by α\hbox{$\overline{\alpha_{r\phi}}$} and αφz\hbox{$\overline{\alpha_{\phi z}}$}. A starting point for this energetics constraint is the conservation equation of total MHD energy (e.g., Balbus & Hawley 1998),

∂t [ 1 2 ρ v 2 + ρ Φ + p γ 1 + B 2 8 π ] + · [ v ( 1 2 ρ v 2 + ρ Φ + γp γ 1 ) \begin{eqnarray} &&\frac{\partial}{\partial t}\left[\frac{1}{2}\rho v^2 + \rho \Phi +\frac{p}{\gamma -1} + \frac{B^2}{8\pi} \right] +\mbf{\nabla\cdot}\left[\mbf{v}\left(\frac{1}{2}\rho v^2 + \rho \Phi +\frac{\gamma p}{\gamma -1} \right) \right. \nonumber\\ &&\left. \qquad +\frac{\mbf{B}}{4\pi}\times (\mbf{v}\times\mbf{B}) + \mbf{F}_{\rm ot}\right] = 0, \label{eq:totengMHD} \end{eqnarray}(11)

where p is the gas pressure, γ is a ratio of specific heats, Φ=GM/r=r2ΩK2r2Ω2\hbox{$\Phi=-GM_{\star}/r = -r^2\Omega_{\rm K}^2\approx -r^2\Omega^2$} is the gravitational potential by a central star, and Fot 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 rΩvrvφ,vz,cs,B/4πρ\hbox{$r\Omega \gg v_r, \delta v_{\phi}, v_z, c_{\rm s}, B/\!\sqrt{4\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))

∂t ( Σ r 2 Ω 2 2 ) + 1 r ∂r [ r Ω { ∂r ( r 2 Σ α c s 2 ) + r 2 α φz ( ρ c s 2 ) mid } + r 2 ΩΣ α c s 2 ] + ( ρ v z ) w E w + F rad = 0 , \begin{eqnarray} \frac{\partial}{\partial t}\left(-\Sigma \frac{r^2\Omega^2}{2}\right) +\frac{1}{r}\frac{\partial}{\partial r}\left[r\Omega \left\{\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}} (\rho c_{\rm s}^2)_{\rm mid} \right\} \right. \nonumber\\ \left. + r^2\Omega\Sigma\overline{\alpha_{r \phi}} c_{\rm s}^2\right] + (\rho v_z)_{\rm w} E_{\rm w}+ F_{\rm rad} =0, \label{eq:totengring} \end{eqnarray}(12)

where Ew is the specific total energy of the gas in the disc wind; (ρvz)wEw is the energy carried away by the disc wind. Frad is radiation loss from the top and bottom surfaces, Frad=2σSBTsurf4,\begin{equation} F_{\rm rad} = 2\sigma_{\rm SB}T_{\rm surf}^4, \label{eq:SB} \end{equation}(13)where σSB is the Stefan-Boltzmann constant and Tsurf 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 α\hbox{$\overline{\alpha_{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, αφz\hbox{$\overline{\alpha_{\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 αφz\hbox{$\overline{\alpha_{\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):

( ρ v z ) w ( E w + r 2 Ω 2 2 ) + F rad = \begin{eqnarray} &&(\rho v_z)_{\rm w}\left(E_{\rm w}+\frac{r^2\Omega^2}{2}\right) + F_{\rm rad}=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \nonumber\\ &&\,\,\,\,\,\, \frac{\Omega}{r}\left[\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid}\right] -\frac{1}{r}\frac{\partial}{\partial r}\left(r^2 \Sigma\Omega \overline{\alpha_{r\phi}}c_{\rm s}^2\right) \label{eq:engcnt1} \\ &&\,\,\,\,\,\,=\frac{3}{2}\Omega\Sigma\overline{\alpha_{r\phi}}c_{\rm s}^2 + r\Omega \overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid} . \label{eq:engcnt2} \end{eqnarray}

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 α\hbox{$\overline{\alpha_{r\phi}}$} includes contributions from the gravitational accretion and from the effective viscous heating.

Suzuki et al. (2010) assumed that Ew32vz2\hbox{$E_{\rm w}\ge \frac{3}{2} v_z^2$} is the condition to drive the vertical outflow to a large distance (Eq. (22) of Suzuki et al. 2010). However, we adopt Ew ≥ 0, because this is the sufficient condition for the wind material to reach z ⇒ ∞ (vz2>0\hbox{$v_z^2>0$} in Eq. (B.5)). Following this consideration, we derived the mass flux of the disc wind that satisfies the energetics constraint with Ew = 0 from Eqs. (14) and (15) in a nondimensional form: Cw,e+2Fradr2Ω2(ρcs)mid=2r3Ω(ρcs)mid∂r(r2Σαcs2)+2csrΩαφz==3π/2h2α+2hαφz,\begin{eqnarray} C_{\rm w,e} + \frac{2 F_{\rm rad}}{r^2\Omega^2 (\rho c_{\rm s})_{\rm mid}} &=& \frac{2}{r^3\Omega(\rho c_{\rm s})_{\rm mid}}\frac{\partial}{\partial r} (r^2 \Sigma\overline{\alpha_{r\phi}}c_{\rm s}^2) + \frac{2c_{\rm s}}{r\Omega}\overline{\alpha_{\phi z}} \nonumber \\ & &- \frac{2}{r^3\Omega^2(\rho c_{\rm s})_{\rm mid}}\frac{\partial}{\partial r} (r^2 \Sigma\Omega\overline{\alpha_{r\phi}}c_{\rm s}^2) \label{eq:Cwcn1} \\ &=&\frac{3\sqrt{2\pi}c_{\rm s}^2}{r^2\Omega^2}\overline{\alpha_{r\phi}} + \frac{2c_{\rm s}}{r\Omega}\overline{\alpha_{\phi z}} \label{eq:Cwcn2} \\ &=& 3\sqrt{\pi/2}h^2 \overline{\alpha_{r\phi}} + \sqrt{2}h\overline{\alpha_{\phi z}}, \nonumber \end{eqnarray}where Cw,e stands for the mass flux constrained by the energetics. We here used ΣΩ=2π(ρcs)mid\hbox{$\Sigma \Omega = \sqrt{2\pi}(\rho c_{\rm s})_{\rm mid}$}, and for the last equality we introduced an aspect ratio, hH/r=2cs/rΩ\hbox{$h\equiv H/r=\sqrt{2}c_{\rm s}/r\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 α\hbox{$\overline{\alpha_{r\phi}}$} term in Eq. (17) regulates the Frad term and the αφz\hbox{$\overline{\alpha_{\phi z}}$} term determines Cw,e. However, the situation is not this simple, because disc winds can be launched solely by the α\hbox{$\overline{\alpha_{r\phi}}$} term, which was shown by local shearing box simulations with zero-wind torque, αφz=0\hbox{$\overline{\alpha_{\phi z}}=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 Cw,e and Frad in this paper. The first strategy is that Frad 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, Cw,e=Frad=\begin{eqnarray} C_{\rm w,e} &=& \max\left(\frac{2}{r^3\Omega(\rho c_{\rm s})_{\rm mid}} \frac{\partial}{\partial r}(r^2 \Sigma\overline{\alpha_{r\phi}}c_{\rm s}^2) + \frac{2c_{\rm s}}{r\Omega}\overline{\alpha_{\phi z}},0\right) \label{eq:DWlossCw} \\ F_{\rm rad} &=& \max\left(-\frac{1}{r}\frac{\partial}{\partial r} (r^2 \Sigma\Omega\overline{\alpha_{r\phi}}c_{\rm s}^2),0\right), \label{eq:DWlossFrad} \end{eqnarray}where we avoided negative values of Cw,e and Frad.

In the second choice we left the uncertainty to a parameter, ϵrad, that determines the fractional energy to the radiation loss: Cw,e==(1ϵrad)[3π/2h2α+2hαφz]Frad=\begin{eqnarray} C_{\rm w,e} &=& (1-\epsilon_{\rm rad})\left[\frac{ 3\sqrt{2\pi}c_{\rm s}^2}{r^2\Omega^2}\overline{\alpha_{r\phi}} + \frac{2c_{\rm s}}{r\Omega}\overline{\alpha_{\phi z}}\right] \label{eq:RadlossCw} \\ &=& (1-\epsilon_{\rm rad})\left[3\sqrt{\pi/2}h^2 \overline{\alpha_{r\phi}} + \sqrt{2}h\overline{\alpha_{\phi z}}\right] \nonumber \nonumber\\ F_{\rm rad} &=& \epsilon_{\rm rad} \left[\frac{3}{2}\Omega\Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2 + r\Omega \overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid}\right]. \label{eq:RadlossFrad} \end{eqnarray}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, Cw,0, by the energetics of the global accretion to give the Cw that we use in our calculations, Cw=min(Cw,0,Cw,e),\begin{equation} C_{\rm w} = \min(C_{\rm w,0},C_{\rm w,e}), \label{eq:Cwe,0} \end{equation}(22)where the adopted Cw,0 is presented in Sect. 2.6.

2.4. Temperature: viscous heating and radiative equilibrium

By referring to the terms concerning α\hbox{$\overline{\alpha_{r\phi}}$} in Eq. (14), the viscous heating rate can be scaled as ~ΣΩcs2\hbox{$\Sigma \Omega c_{\rm s}^2$}. Since Σ decreases with t and Ωcs2\hbox{$\Omega c_{\rm s}^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, Tmid, will be higher than Tsurf in Eq. (13). On the other hand, if the viscous heating is ineffective and the stellar irradiation dominates, then Tsurf will be higher than Tmid. 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 Tvis as the temperature at the midplane determined by viscous heating, 2σSBTvis4=(38τR+12τP)Frad\begin{equation} 2\sigma_{\rm SB} T_{\rm vis}^4 = \left(\frac{3}{8}\tau_{\rm R} +\frac{1}{2\tau_{\rm P}}\right) F_{\rm rad} \label{eq:visT} \end{equation}(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 τR=κRΣ/2,\begin{equation} \tau_{\rm R} = \kappa_{\rm R}\Sigma/2, \end{equation}(24)where we use κR={\begin{eqnarray} \kappa_{\rm R} = \left\{ \begin{array}{ll} 4.5\left(\frac{T}{150~{\rm K}}\right)^2 {\rm cm^{2}\,g^{-1}} & \quad T<150\; {\rm K}\\ 4.5\;{\rm cm^{2}g^{-1}} & \quad 150\;{\rm K}\le T\le 1500\; {\rm K},\\ 0\;{\rm cm^{2}g^{-1}} & \quad T>1500\; {\rm K} \end{array} \right. \label{eq:kappaR} \end{eqnarray}(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 τP=max(2.4τR,0.5)\begin{equation} \tau_{\rm P}=\max(2.4\tau_{\rm R},0.5) \end{equation}(26)(Nakamoto & Nakagawa 1994; Hueso & Guillot 2005), where we give the lower bound on τP to obtain the pre-factor of Eq. (23), (38τR+12τP)1\hbox{$\left(\frac{3}{8}\tau_{\rm R}+\frac{1}{2\tau_{\rm P}}\right)\Rightarrow 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, Treq=T1au(r1au)p·\begin{equation} T_{\rm req} = T_{\rm 1au} \left(\frac{r}{1{~\rm au}}\right)^{p}\cdot \label{eq:treq} \end{equation}(27)We adopted T1au = 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 Tsurf but also Tmid approaches Treq. To take both viscous heating and stellar irradiation into account, we tool the sum of these two temperatures, T4=Tvis4+Treq4,\begin{equation} T^4 = T_{\rm vis}^4 + T_{\rm req}^4, \label{eq:Tselect} \end{equation}(28)for the representative z-averaged temperature, T, to estimate cs 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 rcutΣint=Σ1au(r1au)3/2exp(rrcut)·\begin{equation} \Sigma_{\rm int} = \Sigma_{1{\rm au}} \left(\frac{r}{\rm 1\; au}\right)^{-3/2} \exp\left(-\frac{r}{r_{\rm cut}}\right)\cdot \label{eq:Sgminit} \end{equation}(29)The original MMSN by Hayashi (1981) considered Σ1au = 1.7 × 103 g cm-2 with a sharp cut-off at 36 au, which gives the initial disc mass, Mdisc,int = 0.013 M. We adopted a ten times larger Σ1au = 1.7 × 104 g cm-2 but slightly smaller rcut = 30 au in this paper, which gives Mdisc,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 rin = 0.01 au to rout = 104 au with grid spacing, Δrr\hbox{$\Delta r \propto \sqrt{r}$}. At the inner and outer boundaries, r = rin and = rout, we imposed ∂r(Σr3/2)=0\hbox{$\frac{\partial}{\partial r}(\Sigma r^{3/2})=0$}, which corresponds to the zero-torque boundary condition (Lynden-Bell & Pringle 1974); the α\hbox{$\overline{\alpha_{r\phi}}$} term in Eq. (10), ∂r(r2Σαcs2)\hbox{$\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2)$}, is zero for a constant α\hbox{$\overline{\alpha_{r\phi}}$} and cs2r1/2\hbox{$c_{\rm s}^2\propto r^{-1/2}$} (Eq. (27) with p = −1/2).

2.6. Parameters

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

2.6.1. Turbulent viscosity – αrφ\hbox{$\overline{\alpha_{\mathsfsl{r}\phi}}$}

We compared two cases with spatially uniform α=8×10-3\hbox{$\overline{\alpha_{r\phi}}=8\times 10^{-3}$}, and 8 × 10-5. α=8×10-3\hbox{$\overline{\alpha_{r\phi}}=8\times 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 α\hbox{$\overline{\alpha_{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 α=8×10-5\hbox{$\overline{\alpha_{r\phi}} =8\times 10^{-5}$} for such MRI-inactive circumstances. Although we assumed constant α\hbox{$\overline{\alpha_{r\phi}}$} for simplicity, α\hbox{$\overline{\alpha_{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 α\hbox{$\overline{\alpha_{r\phi}}$} into account.

2.6.2. Disc wind mass flux – Cw,0

The mass flux of disc winds, Cw,0, was also adopted from the local simulations. Cw,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 Cw,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 Cw,0 carefully from the local simulations.

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

2.6.3. Disc wind torque – αφz\hbox{$\overline{\alpha_{\phi \textit{z}}}$}

We tested two types of the parametrization for the wind torque: (i) constant αφz=1×10-4\hbox{$\overline{\alpha_{\phi z}}=1\times 10^{-4}$}; and (ii) density dependent with a cap, αφz=min(10-5(ΣΣint)-0.66,1).\begin{equation} \overline{\alpha_{\phi z}}=\min\left(10^{-5}\left(\frac{\Sigma}{\Sigma_{\rm int}} \right)^{-0.66},1\right). \label{eq:apzddp} \end{equation}(30)We name (i) constant torque and (ii) Σ-dependent torque from now on. αφz\hbox{$\overline{\alpha_{\phi z}}$} was estimated by local MHD simulations by Bai (2013), who reported αφz~10-510-3\hbox{$\overline{\alpha_{\phi z}}\sim 10^{-5} {-} 10^{-3}$} with a positive dependence on the strength of the net vertical magnetic field, αφz(Bz2/8π(ρcs2)mid)0.66\hbox{$\overline{\alpha_{\phi z}}\propto (B_z^2 /8\pi (\rho c_{\rm s}^2)_{\rm mid})^{0.66}$}. ρmid is proportional to Σ, while Bz 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 Bz decreases with the dispersal of gas (decrease of Σ), then αφz\hbox{$\overline{\alpha_{\phi z}}$} will stay approximately constant, which corresponds to (i) constant torque; if Bz does not decrease that much, then αφz\hbox{$\overline{\alpha_{\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

thumbnail Fig. 1

Comparison of time evolutions of four MRI-active PPDs with α=8×10-3\hbox{$\overline{\alpha_{r\phi}}=8\times 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 Cw,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 Cw and wind torques, αφz\hbox{$\overline{\alpha_{\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, TTvis (Eq. (28)) in these cases, and then T is mainly determined from Frad by Eq. (23). Recalling Σintr− 3/2, we derive 1r∂r(r2ΣΩαcs2)32ΣΩαcs2\hbox{$-\frac{1}{r}\frac{\partial}{\partial r}(r^2 \Sigma \Omega \overline{\alpha_{r\phi}} c_{\rm s}^2) \approx \frac{3}{2}\Sigma \Omega \overline{\alpha_{r\phi}} c_{\rm s}^2$} for cs2~r1/2\hbox{$c_{\rm s}^2 \sim r^{-1/2}$}. Since the αφz\hbox{$\overline{\alpha_{\phi z}}$}(0 or =10-5) term in Eq. (21) is negligible in comparison to the α(=8×10-3)\hbox{$\overline{\alpha_{r\phi}} ({=}8\times 10^{-3})$} term, both strong DW and weak DW conditions give similar Frad 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 (Tvis>Treq) 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 Treq in the entire region after t ≳ 106 yr because the surface densities decrease rapidly by the disc winds in the inner region to give TreqTvis, while Tvis is no longer negligible in the weak DW case (purple lines) at t = 106 yr.

thumbnail 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, Cw, which is shown in Fig. 2. At t = 0 the strong DW case (green dotted line) gives quite small Cw ≈ 0 below the displayed range of Fig. 2 because ∂r(r2Σαcs2)0\hbox{$\frac{\partial}{\partial r} (r^2\Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2)\approx 0$} for Σintr− 3/2 in Eq. (18). However, as Σ decreases in an inside-out manner and the Σ profile changes, Cw increases and at t = 106 yr this case (green solid line) yields larger Cw than the weak DW case (purple solid line), in which Cw instead decreases with time owing to the decrease in temperature (cs2\hbox{$c_{\rm s}^2$}; Eq. (21)). We note that Cw = 0 in the outer region, r> 90 au, of the strong DW + zero-torque case because the gas moves outward (∂r(r2Σαcs2)<0\hbox{$\frac{\partial}{\partial r}(r^2\Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2)< 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.

Table 1

Parameters for MRI-active cases.

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 Cw (Eq. (18)). As a result, Cw is not constrained by the energetics, Cw,e, in the almost entire region but is determined by Cw,0( = 2 × 10-5) at t = 106 yr (red solid lines). The constant Cw = Cw,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 ≳ 107 yr. This is because αφz\hbox{$\overline{\alpha_{\phi z}}$} is constrained by the cap value = 1 (Eq. (30) there.

thumbnail Fig. 3

Mass-loss rate by disc wind, z, (solid lines) and mass accretion rate induced by the 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), z(r)=2πrroutrdr(ρvz)w=2πrroutrdrCw(ρcs)mid,\begin{equation} \dot{M}_z(r) = 2\pi\int_{r}^{r_{\rm out}} r {\rm d}r (\rho v_{z})_{\rm w} =2\pi\int_{r}^{r_{\rm out}} r {\rm d}r C_{\rm w}(\rho c_{\rm s})_{\rm mid}, \label{eq:Mzdef} \end{equation}(31)and the mass accretion rate, r(r)=2πrΣvr,\begin{equation} \dot{M}_r(r) = -2\pi r \Sigma v_r, \label{eq:Mrdef} \end{equation}(32)at t = 106 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 stress (dashed lines), r,rφ(r)=4πrΩ∂r(r2Σαcs2),\begin{equation} \dot{M}_{r,r\phi}(r) = - \frac{4\pi}{r\Omega}\frac{\partial}{\partial r} \left(r^2\Sigma\overline{\alpha_{r\phi}}c_{\rm s}^2\right), \label{eq:Mrrpdef} \end{equation}(33)and that by the φz stress (dotted line) r,φz(r)=4πΩrαφz(ρcs2)mid.\begin{equation} \dot{M}_{r,\phi z}(r) = - \frac{4\pi}{\Omega}r\overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid}. \label{eq:Mrpzdef} \end{equation}(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 Mr˙\hbox{$\dot{M_r}$} is converted into z as mass accretes inward.

thumbnail 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 rin = 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 ≈ 100r,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 αφz\hbox{$\overline{\alpha_{\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 tr trends can be directly compared to the observed distribution in the tr 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 = 106 yr), r of the other three cases are well inside the observed range.

3.2. MRI-inactive cases

Table 2

Parameters for MRI-inactive cases.

thumbnail Fig. 5

Same as Fig. 1 but for four MRI-inactive cases with α=8×10-5\hbox{$\overline{\alpha_{r\phi}}=8\times 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 α\hbox{$\overline{\alpha_{r\phi}}$} gives smaller Frad (Eqs. (19) and (21)) and accordingly lower Tvis (Eq. (23)).

thumbnail Fig. 6

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

Smaller α\hbox{$\overline{\alpha_{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 α\hbox{$\overline{\alpha_{r\phi}}$} and second because the disc wind mass flux is strongly constrained by the energetics to give smaller Cw (Fig. 6 in comparison to Fig. 2).

thumbnail Fig. 7

Comparison of αφz\hbox{$\overline{\alpha_{\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 αφz\hbox{$\overline{\alpha_{\phi z}}$}, because its effect is relatively important for lower turbulent viscosity, α\hbox{$\overline{\alpha_{r\phi}}$}. The addition of the spatially constant αφz=10-4\hbox{$\overline{\alpha_{\phi z}}=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 αφz\hbox{$\overline{\alpha_{\phi z}}$} for the Σ-dependent torque cases. αφz\hbox{$\overline{\alpha_{\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, Cw, is not constrained by the energetics (Eqs. (18) and (20)) but is chosen to be the constant Cw,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 αφz\hbox{$\overline{\alpha_{\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 α\hbox{$\overline{\alpha_{r\phi}}$} and Cw,0, they are also anticipated to depend on the strength of net vertical magnetic field, Bz2/8π(ρcs2)mid\hbox{$B_z^2 /8\pi(\rho c_{\rm s}^2)_{\rm mid}$}. In this case, α\hbox{$\overline{\alpha_{r\phi}}$} and Cw,0 could inversely correlate with Σ (Suzuki et al. 2010), which additionally enhances the positive slopes of the surface densities.

thumbnail Fig. 8

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

thumbnail 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 = 106 yr. In the zero-torque case (green lines) the mass is dominantly lost by the disc winds, z ≈ 100r,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,φzz ≈ 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 tr (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 tr. 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 = 105 yr, 1.7 × 10-9M yr-1 at 106 yr, and 5.2 × 10-10M yr-1 at 107 yr.

4. Discussion

4.1. Uncertainties

Our model has the three free parameters, α\hbox{$\overline{\alpha_{r\phi}}$}, Cw,0, and αφz\hbox{$\overline{\alpha_{\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 αφz\hbox{$\overline{\alpha_{\phi z}}$} and Σ-dependent αφz\hbox{$\overline{\alpha_{\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, αφz\hbox{$\overline{\alpha_{\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, α\hbox{$\overline{\alpha_{r\phi}}$} and Cw,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 αφz\hbox{$\overline{\alpha_{\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, Cw,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 Cw,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 Cw,0 is smaller.

As shown in Figs. 2 and 6,Cw is already constrained by the energetics. In most cases except for the MRI-inactive cases with Σ-dependent torque, the energetics constraint already suppresses Cw in the inner region. Therefore, adopting a smaller Cw,0 does not affect Cw in the inner region but reduces Cw 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 Cw at the relatively outer location, r ~ 10 au. In these cases, a smaller Cw,0 reduces Cw 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 Cw in the region with Cw,e<Cw,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 τdr,max1ηΩK,\begin{equation} \tau_{\rm dr,max} \approx \frac{1}{\eta \Omega_{\rm K}}, \end{equation}(35)where η is pressure gradient force normalized by the twice of centrifugal force, η=1ρmidpmid∂r12rΩK2·\begin{equation} \eta = -\frac{1}{\rho_{\rm mid}}\frac{\partial p_{\rm mid}}{\partial r} \frac{1}{2r\Omega_{\rm K}^2}\cdot \end{equation}(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 (α=8×10-5\hbox{$\overline{\alpha_{r\phi}}=8\times 10^{-5}$}) cases with Σ-dependent torque of Table 2 (red and blue lines; the same as in Figs. 59) in comparison to the no disc wind (no DW) case with the same α=8×10-5\hbox{$\overline{\alpha_{r\phi}}=8\times 10^{-5}$} (black lines). We here derive pmid from Σ by pmid=ρmidcs2=ΣΩcs2π·\begin{equation} p_{\rm mid} = \rho_{\rm mid} c_{\rm s}^2 = \frac{\Sigma\Omega c_{\rm s}}{\sqrt{2\pi}}\cdot \end{equation}(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 = 105 yr and 1–2 au at t = 106 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.

thumbnail Fig. 10

Comparison of normalized pressure gradient force, (1ρmidpmid∂r)/(2rΩ2)\hbox{$-\left(\frac{1}{\rho_{\rm mid}}\frac{\partial p_{\rm mid}}{\partial r}\right) /(2r\Omega^2)$}, 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 α=8×10-5\hbox{$\overline{\alpha_{r\phi}} =8\times 10^{-5}$} (black lines).

4.3. Type I migration

thumbnail 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, ta (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, CI ≡ −ta,TTW/ta, where ta,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 taa/ (−ȧ); 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 = 106 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 ta,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 (Mdisc,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 β=(Bz2/8π(ρcs2)mid)-1\hbox{$\beta = (B_z^2/8\pi(\rho c_{\rm s}^2)_{\rm mid})^{-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, Treq (Eq. (27)). However, the stellar wind and irradiation affect the evolution of PPDs.

In our calculations, the mass flux of the disc wind is Cw,e constrained by the energetics of accretion, and it can be smaller than Cw,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 Cw,e in Eqs. (16) and (17), in the small r region. The increase of Cw 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, α\hbox{$\overline{\alpha_{r\phi}}$}, disc wind mass flux, Cw, and disc wind torque, αφz\hbox{$\overline{\alpha_{\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, Mdisc.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 106 yr, and the viscous heating is negligible in determining the temperature after t ≳ 106 yr, whereas in the no disc wind cases the viscous heating is not negligible even up to several 106 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, α\hbox{$\overline{\alpha_{r\phi}}$}, Cw,0, and αφz\hbox{$\overline{\alpha_{\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

  1. Adams, F. C., Hollenbach, D., Laughlin, G., & Gorti, U. 2004, ApJ, 611, 360 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229 [NASA ADS] [CrossRef] [Google Scholar]
  3. Armitage, P. J. 2011, ARA&A, 49, 195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bai, X.-N. 2013, ApJ, 772, 96 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bai, X.-N. 2016, ApJ, 821, 80 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bai, X.-N., & Stone, J. M. 2013a, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bai, X.-N., & Stone, J. M. 2013b, ApJ, 769, 76 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [Google Scholar]
  10. Baillié, K., Charnoz, S., & Pantin, E. 2015, A&A, 577, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
  12. Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
  13. Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon) [Google Scholar]
  16. Chiang, E., & Youdin, A. N. 2010, Ann. Rev. Earth Planet. Sci., 38, 493 [Google Scholar]
  17. Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cranmer, S. R. 2009, ApJ, 706, 824 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cranmer, S. R., & Saar, S. H. 2011, ApJ, 741, 54 [NASA ADS] [CrossRef] [Google Scholar]
  20. Davis, S. S. 2005, ApJ, 620, 994 [NASA ADS] [CrossRef] [Google Scholar]
  21. Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Dyda, S., Lovelace, R. V. E., Ustyugova, G. V., et al. 2015, MNRAS, 450, 481 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114 [NASA ADS] [CrossRef] [Google Scholar]
  24. Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398 [NASA ADS] [CrossRef] [Google Scholar]
  25. Ferreira, J., Dougados, C., & Cabrit, S. 2006, A&A, 453, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Flock, M., Henning, T., & Klahr, H. 2012, ApJ, 761, 95 [NASA ADS] [CrossRef] [Google Scholar]
  27. Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
  29. Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [Google Scholar]
  30. Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539 [Google Scholar]
  31. Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
  32. Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097 [NASA ADS] [CrossRef] [Google Scholar]
  33. Guilet, J., & Ogilvie, G. I. 2014, MNRAS, 441, 852 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323 [NASA ADS] [CrossRef] [Google Scholar]
  35. Haisch, Jr., K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
  36. Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84 [NASA ADS] [CrossRef] [Google Scholar]
  38. Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  39. Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, eds. D. C. Black, & M. S. Matthews (Tucson: The University of Arizona Press), 1100 [Google Scholar]
  40. Hernández, J., Hartmann, L., Calvet, N., et al. 2008, ApJ, 686, 1195 [NASA ADS] [CrossRef] [Google Scholar]
  41. Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30 [NASA ADS] [CrossRef] [Google Scholar]
  42. Hirose, S., Uchida, Y., Shibata, K., & Matsumoto, R. 1997, PASJ, 49, 193 [Google Scholar]
  43. Hollenbach, D. J., Yorke, H. W., & Johnstone, D. 2000, Protostars and Planets IV (Tucson: The University of Arizona Press), 401 [Google Scholar]
  44. Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15 [NASA ADS] [CrossRef] [Google Scholar]
  45. Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
  47. Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016, MNRAS, 461, 2257 [NASA ADS] [CrossRef] [Google Scholar]
  48. Kobayashi, H., Ormel, C. W., & Ida, S. 2012, ApJ, 756, 70 [NASA ADS] [CrossRef] [Google Scholar]
  49. Kobayashi, H., Tanaka, H., & Okuzumi, S. 2016, ApJ, 817, 105 [NASA ADS] [CrossRef] [Google Scholar]
  50. Kusaka, T., Nakano, T., & Hayashi, C. 1970, Prog. Theor. Phys., 44, 1580 [NASA ADS] [CrossRef] [Google Scholar]
  51. Lesur, G., & Longaretti, P.-Y. 2007, MNRAS, 378, 1471 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235 [NASA ADS] [Google Scholar]
  54. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
  55. Manara, C. F., Fedele, D., Herczeg, G. J., & Teixeira, P. S. 2016, A&A, 585, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Matt, S., & Pudritz, R. E. 2005, ApJ, 632, L135 [NASA ADS] [CrossRef] [Google Scholar]
  57. Miyake, T., Suzuki, T. K., & Inutsuka, S.-I. 2016, ApJ, 821, 3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, ICARUS, 67, 375 [Google Scholar]
  59. Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640 [NASA ADS] [CrossRef] [Google Scholar]
  60. Ogihara, M., Kobayashi, H., Inutsuka, S.-I., & Suzuki, T. K. 2015a, A&A, 579, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Ogihara, M., Morbidelli, A., & Guillot, T. 2015b, A&A, 584, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 [NASA ADS] [CrossRef] [Google Scholar]
  63. Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 [NASA ADS] [CrossRef] [Google Scholar]
  64. Okuzumi, S., Takeuchi, T., & Muto, T. 2014, ApJ, 785, 127 [NASA ADS] [CrossRef] [Google Scholar]
  65. Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 [NASA ADS] [CrossRef] [Google Scholar]
  66. Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [NASA ADS] [CrossRef] [Google Scholar]
  67. Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293 [NASA ADS] [CrossRef] [Google Scholar]
  68. Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117 [NASA ADS] [CrossRef] [Google Scholar]
  69. Pessah, M. E., Chan, C.-K., & Psaltis, D. 2006, MNRAS, 372, 183 [NASA ADS] [CrossRef] [Google Scholar]
  70. Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Rothstein, D. M., & Lovelace, R. V. E. 2008, ApJ, 677, 1221 [NASA ADS] [CrossRef] [Google Scholar]
  72. Ruden, S. P., & Lin, D. N. C. 1986, ApJ, 308, 883 [NASA ADS] [CrossRef] [Google Scholar]
  73. Sai, K., Katoh, Y., Terada, N., & Ono, T. 2013, ApJ, 767, 165 [NASA ADS] [CrossRef] [Google Scholar]
  74. Salmeron, R., Königl, A., & Wardle, M. 2011, MNRAS, 412, 1162 [NASA ADS] [Google Scholar]
  75. Sano, T., Inutsuka, S.-I., & Miyama, S. M. 1998, ApJ, 506, L57 [NASA ADS] [CrossRef] [Google Scholar]
  76. Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef] [Google Scholar]
  77. Sano, T., Inutsuka, S.-I., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321 [NASA ADS] [CrossRef] [Google Scholar]
  78. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  79. Shu, F. H., Johnstone, D., & Hollenbach, D. 1993, Icarus, 106, 92 [NASA ADS] [CrossRef] [Google Scholar]
  80. Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 [NASA ADS] [CrossRef] [Google Scholar]
  81. Simon, J. B., Hawley, J. F., & Beckwith, K. 2011, ApJ, 730, 94 [NASA ADS] [CrossRef] [Google Scholar]
  82. Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013, ApJ, 775, 73 [NASA ADS] [CrossRef] [Google Scholar]
  83. Suzuki, T. K., & Inutsuka, S.-I. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
  84. Suzuki, T. K., & Inutsuka, S.-I. 2014, ApJ, 784, 121 [Google Scholar]
  85. Suzuki, T. K., Muto, T., & Inutsuka, S.-I. 2010, ApJ, 718, 1289 [NASA ADS] [CrossRef] [Google Scholar]
  86. Suzuki, T. K., Imada, S., Kataoka, R., et al. 2013, PASJ, 65, 98 [NASA ADS] [Google Scholar]
  87. Takagi, Y., Itoh, Y., & Oasa, Y. 2014, PASJ, 66, 88 [NASA ADS] [Google Scholar]
  88. Takagi, Y., Itoh, Y., Arai, A., Sai, S., & Oasa, Y. 2015, PASJ, 67, 87 [NASA ADS] [Google Scholar]
  89. Takeuchi, T., & Okuzumi, S. 2014, ApJ, 797, 132 [NASA ADS] [CrossRef] [Google Scholar]
  90. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
  91. Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773, 155 [NASA ADS] [CrossRef] [Google Scholar]
  92. Turner, N. J., & Sano, T. 2008, ApJ, 679, L131 [NASA ADS] [CrossRef] [Google Scholar]
  93. Velikhov, E. P. 1959, Zh. Eksp. Teor. Fiz., 36, 1398 [Google Scholar]
  94. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
  95. Wood, B. E., Müller, H.-R., Zank, G. P., Linsky, J. L., & Redfield, S. 2005, ApJ, 628, L143 [Google Scholar]
  96. Zanni, C., & Ferreira, J. 2013, A&A, 550, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

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

∂t ( ρr v φ ) + 1 r ∂r [ r 2 ( ρ v r v φ B r B φ 4 π ) ] \appendix \setcounter{section}{1} \begin{eqnarray} &&\frac{\partial}{\partial t}(\rho r v_{\phi}) + \frac{1}{r} \frac{\partial}{\partial r}\left[r^2\left(\rho v_r v_{\phi} -\frac{B_r B_{\phi}}{4\pi}\right)\right] \nonumber\\ && \qquad\qquad + \frac{\partial}{\partial z}\left[r\left(\rho v_{\rm \phi} v_z -\frac{B_{\phi}B_z}{4\pi}\right)\right] = 0. \label{eq:angmom} \end{eqnarray}(A.1)

The azimuthal velocity, vφ, is decomposed into the mean Keplerian flow and perturbation, vφ=rΩ+δvφ.\appendix \setcounter{section}{1} \begin{equation} v_{\phi} = r\Omega + \delta v_{\phi}. \label{eq:vphdcmp} \end{equation}(A.2)We use the α prescription (Shakura & Sunyaev 1973) for the second and third terms of Eq. (A.1): ρvrvφBrBφ4π=ρvrrΩ+ρ(vrδvφBrBφ4πρ)\appendix \setcounter{section}{1} \begin{eqnarray} \rho v_r v_{\phi} - \frac{B_r B_{\phi}}{4\pi} &=& \rho v_r r\Omega + \rho \left(v_r \delta v_{\phi} - \frac{B_r B_{\phi}}{4\pi \rho}\right) \nonumber\\ &\equiv& \rho v_r r\Omega + \rho \alpha_{r\phi}c_{\rm s}^2, \label{eq:Mxrp} \end{eqnarray}(A.3)and ρvφvzBφBz4π=ρrΩvz+ρ(δvφvzBφBz4πρ)\appendix \setcounter{section}{1} \begin{eqnarray} \rho v_{\phi} v_z - \frac{B_{\phi}B_z}{4\pi}& =& \rho r\Omega v_z + \rho \left(\delta v_{\phi}v_z - \frac{B_{\phi}B_z}{4\pi \rho}\right) \nonumber\\ &\equiv& \rho r\Omega v_z + \rho \alpha_{\phi z}c_{\rm s}^2. \label{eq:Mxpz} \end{eqnarray}(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 ∂t(Σr3Ω)+∂r[r2Σ(vrrΩ+αcs2)]\appendix \setcounter{section}{1} \begin{eqnarray*} \frac{\partial}{\partial t}(\Sigma r^3\Omega) + \frac{\partial}{\partial r} \left[r^2\Sigma \left(v_r r\Omega + \overline{\alpha_{r\phi}}c_{\rm s}^2\right)\right] \end{eqnarray*}+r2[(ρvz)wrΩ+αφz(ρcs2)mid]=0,\appendix \setcounter{section}{1} \begin{equation} \qquad\qquad +r^2 \left[(\rho v_z)_{\rm w}r\Omega + \overline{\alpha_{\phi z}} (\rho c_{\rm s}^2)_{\rm mid}\right]= 0, \label{eq:amintgz} \end{equation}(A.5)where α=ραdz/Σ\hbox{$\overline{\alpha_{r\phi}}=\int \rho \alpha_{r\phi}{\rm d}z/\Sigma$} is the mass-weighted vertical average. The third term, which represents the angular momentum loss from both surfaces, is derived from [ρrΩvz+ραφzcs2]w=[ρvz]wrΩ+(ρcs2)midαφz,\appendix \setcounter{section}{1} \begin{equation} \left[\rho r \Omega v_z + \rho \alpha_{\phi z}c_{\rm s}^2\right]_{\rm w} =[\rho v_z]_{\rm w} r\Omega + (\rho c_{\rm s}^2)_{\rm mid} \overline{\alpha_{\phi z}}, \end{equation}(A.6)where the subscript w stands for disc wind. αφz\hbox{$\overline{\alpha_{\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 Σ∂t+1r∂r(rΣvr)+(ρvz)w=0.\appendix \setcounter{section}{1} \begin{equation} \frac{\partial \Sigma}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r} (r\Sigma v_r) + (\rho v_z)_{\rm w} =0. \label{eq:masscont} \end{equation}(A.7)By combining Eq. (A.7) multiplied by r3Ω and Eq. (A.5), we have rΣvr∂r(r2Ω)+∂r(r2Σαcs2)+r2αφz(ρcs2)mid=0,\appendix \setcounter{section}{1} \begin{equation} r\Sigma v_r \frac{\partial}{\partial r}(r^2 \Omega) + \frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2\overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid} = 0, \end{equation}(A.8)which determines the accretion rate, rΣvr=2rΩ[∂r(r2Σαcs2)+r2αφz(ρcs2)mid],\appendix \setcounter{section}{1} \begin{equation} r\Sigma v_r = - \frac{2}{r\Omega} \left[\frac{\partial}{\partial r} (r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2\overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid}\right], \label{eq:accvel} \end{equation}(A.9)where we here assumed the Keplerian rotation to derive ∂r(r2Ω)=rΩ2\hbox{$\frac{\partial}{\partial r} (r^2\Omega) = \frac{r\Omega}{2}$}.

By substituting Eq. (A.9) into Eq. (A.7), we finally have the equation for the time evolution of Σ (Eq. (10)): Σ∂t1r∂r[2rΩ{∂r(r2Σαcs2)+r2αφz(ρcs2)mid}]\appendix \setcounter{section}{1} \begin{eqnarray*} \frac{\partial \Sigma}{\partial t} - \frac{1}{r}\frac{\partial}{\partial r} \left[\frac{2}{r\Omega}\left\{\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}} (\rho c_{\rm s}^2)_{\rm mid} \right\}\right] \end{eqnarray*}+(ρvz)w=0,\appendix \setcounter{section}{1} \begin{eqnarray*} \qquad \qquad + (\rho v_z)_{\rm w} = 0, \end{eqnarray*}

Appendix B: Energetics of accretion discs

A general MHD expression of the total energy conservation under the axisymmetric approximation is ∂t[12ρv2+ρΦ+pγ1+B28π]\appendix \setcounter{section}{2} \begin{eqnarray*} \frac{\partial}{\partial t}\left[\frac{1}{2}\rho v^2 +\rho \Phi + \frac{p}{\gamma-1} + \frac{B^2}{8\pi}\right] \end{eqnarray*}+1r∂r[r{vr(12ρv2+ρΦ+γγ1p+Bφ2+Bz24π)\appendix \setcounter{section}{2} \begin{eqnarray*} \qquad + \frac{1}{r}\frac{\partial}{\partial r}\left[r\left\{v_r\left(\frac{1}{2}\rho v^2 +\rho\Phi + \frac{\gamma}{\gamma -1}p + \frac{B_{\phi}^2+B_z^2}{4\pi} \right) \right.\right. \end{eqnarray*}Br4π(vφBφ+vzBz)+Fot,r}]\appendix \setcounter{section}{2} \begin{eqnarray*} \qquad \left.\left.- \frac{B_r}{4\pi}(v_{\phi}B_{\phi}+v_z B_z)+F_{{\rm ot},r}\right\}\right] \end{eqnarray*}+∂z[vz(12ρv2+ρΦ+γγ1p+Br2+Bφ24π)\appendix \setcounter{section}{2} \begin{eqnarray*} \qquad +\frac{\partial}{\partial z}\left[v_z\left(\frac{1}{2}\rho v^2 +\rho\Phi + \frac{\gamma}{\gamma -1}p + \frac{B_r^2+B_{\phi}^2}{4\pi} \right) \right. \end{eqnarray*}Bz4π(vrBr+vφBφ)+Fot,z]=0,\appendix \setcounter{section}{2} \begin{equation} \qquad \left.- \frac{B_z}{4\pi}(v_r B_r + v_{\phi}B_{\phi}) + F_{{\rm ot},z}\right]= 0, \label{eq:toteng} \end{equation}(B.1)where we refer to Eq. (11) for the definition of each variable. Decomposing vφ by Eq. (A.2) and assuming rΩvrvφ,vz,cs,B/4πρ\hbox{$r\Omega \gg v_r, \delta v_{\phi}, v_z, c_s, B/\sqrt{4\pi\rho}$} in a disc, we rewrite Eq. (B.1) with leaving dominant terms. The time-derivative term becomes ∂t[12ρv2+pγ1+ρΦ+B28π]∂t[12ρv2+ρΦ]\appendix \setcounter{section}{2} \begin{eqnarray*} \frac{\partial}{\partial t}\left[\frac{1}{2}\rho v^2 + \frac{p}{\gamma-1} +\rho \Phi + \frac{B^2}{8\pi}\right] \approx \frac{\partial}{\partial t}\left[\frac{1}{2}\rho v^2 +\rho \Phi\right] \end{eqnarray*}∂t[12ρ(rΩ+δvφ)2ρr2Ω2]∂t(12ρr2Ω2),\appendix \setcounter{section}{2} \begin{equation} \qquad \approx \frac{\partial}{\partial t}\left[\frac{1}{2}\rho(r\Omega + \delta v_{\phi})^2 - \rho r^2\Omega^2 \right] \approx \frac{\partial}{\partial t}\left(-\frac{1}{2}\rho r^2\Omega^2\right), \label{eq:tderiv} \end{equation}(B.2)where we set rΩδvφ = 0 after the azimuthal average. The r-derivative term, except for Fot,r, can be approximated as ∂r[r{vr(12ρv2+ρΦ+γγ1p+Bφ2+Bz24π)\appendix \setcounter{section}{2} \begin{eqnarray*} \frac{\partial}{\partial r}\left[r\left\{v_r\left(\frac{1}{2}\rho v^2 +\rho\Phi + \frac{\gamma}{\gamma -1}p + \frac{B_{\phi}^2+B_z^2}{4\pi} \right) \right.\right. \end{eqnarray*}Br4π(vφBφ+vzBz)}]\appendix \setcounter{section}{2} \begin{eqnarray*} \left.\left.\quad - \frac{B_r}{4\pi}(v_{\phi}B_{\phi}+v_z B_z)\right\}\right] \end{eqnarray*}∂r[r{vr(12ρv2+ρΦ)Br4πvφBφ}]\appendix \setcounter{section}{2} \begin{eqnarray*} \qquad \approx \frac{\partial}{\partial r}\left[r\left\{v_r\left(\frac{1}{2}\rho v^2 +\rho\Phi\right) - \frac{B_r}{4\pi}v_{\phi}B_{\phi}\right\}\right] \end{eqnarray*}∂r[r{ρvrr2Ω22+ρrΩ(vrδvφBrBφ4πρ)}]\appendix \setcounter{section}{2} \begin{eqnarray*} \qquad\approx \frac{\partial}{\partial r}\left[r\left\{-\rho v_r \frac{r^2\Omega^2}{2}+\rho r\Omega\left(v_r\delta v_{\phi} -\frac{B_r B_{\phi}}{4\pi\rho}\right)\right\}\right] \end{eqnarray*}=∂r[r{ρvrr2Ω22+ρrΩαcs2}],\appendix \setcounter{section}{2} \begin{equation} \qquad=\frac{\partial}{\partial r}\left[r\left\{-\rho v_r \frac{r^2\Omega^2}{2}+\rho r\Omega\alpha_{r\phi}c_{\rm s}^2\right\}\right], \label{eq:rderiv} \end{equation}(B.3)where the second is derived from vr(v22+Φ)vr[(rΩ+δvφ)22r2Ω2]vrr2Ω22+ρrΩvrδvφ\hbox{$v_r\left(\frac{v^2}{2} + \Phi\right)\approx v_r\left[\frac{(r\Omega+\delta v_{\phi})^2}{2} - r^2\Omega^2\right]\approx - v_r \frac{r^2\Omega^2}{2} + \rho r\Omega v_r \delta v_{\phi}$}, and for the last equality Eq. (A.3) is used. We set the z-derivative term, except for Fot.z, to be

∂z [ v z ( 1 2 ρ v 2 + ρ Φ + γ γ 1 p + B r 2 + B φ 2 4 π ) \appendix \setcounter{section}{2} \begin{eqnarray} &&\frac{\partial}{\partial z}\left[v_z\left(\frac{1}{2}\rho v^2 +\rho\Phi + \frac{\gamma}{\gamma -1}p + \frac{B_r^2+B_{\phi}^2}{4\pi} \right) \right. \nonumber\\ &&\left.\qquad - \frac{B_z}{4\pi}(v_r B_r + v_{\phi}B_{\phi})\right] \equiv \frac{\partial}{\partial z}(\rho v_z E_{\rm w}). \label{eq:zderiv} \end{eqnarray}(B.4)

In the wind region, the kinetic energy will eventually dominate (Pelletier & Pudritz 1992), Ewvz22(z),\appendix \setcounter{section}{2} \begin{equation} E_{\rm w} \approx \frac{v_z^2}{2} \;\; ({z\Rightarrow \infty}), \label{eq:Edwinfty} \end{equation}(B.5)provided that the disc wind is accelerated with increasing z.

By substituting Eqs. (B.2)–(B.4) into Eq. (B.1), we obtain

∂t ( ρ r 2 Ω 2 2 ) + 1 r ∂r [ r { ρ v r r 2 Ω 2 2 + ρr Ω α c s 2 + F ot ,r } ] \appendix \setcounter{section}{2} \begin{eqnarray} && \frac{\partial}{\partial t}\left(-\rho \frac{r^2\Omega^2}{2}\right) +\frac{1}{r}\frac{\partial}{\partial r}\left[r\left\{-\rho v_r \frac{r^2\Omega^2}{2} + \rho r\Omega\alpha_{r\phi} c_{\rm s}^2 + F_{{\rm ot},r} \right\}\right] \\\nonumber &&\qquad+\frac{\partial}{\partial z}(\rho v_z E_{\rm w} + F_{{\rm ot},z})= 0. \label{eq:toteaprx} \end{eqnarray}(B.6)

We integrate Eq. (B.7) from the bottom surface to the top surface along z:

∂t ( Σ r 2 Ω 2 2 ) + 1 r ∂r [ r { Σ v r r 2 Ω 2 2 + Σ r Ω α c s 2 } ] \appendix \setcounter{section}{2} \begin{eqnarray} &&\frac{\partial}{\partial t}\left(-\Sigma \frac{r^2\Omega^2}{2}\right) +\frac{1}{r}\frac{\partial}{\partial r}\left[r\left\{-\Sigma v_r \frac{r^2\Omega^2}{2} + \Sigma r\Omega\overline{\alpha_{r\phi}} c_{\rm s}^2\right\}\right] \nonumber\\ &&\qquad + (\rho v_z)_{\rm w} E_{\rm w}+ F_{\rm rad} =0, \label{eq:totezint1} \end{eqnarray}(B.7)

where (ρvz)wEw and Frad are the energy loss by disc winds and radiation from the top and bottom surfaces. Here Frad is from Fot. By substituting Eq. (A.9) into Eq. (B.7), we have

∂t(Σr2Ω22)+1r∂r[rΩ{∂r(r2Σαcs2)+r2αφz(ρcs2)mid}\appendix \setcounter{section}{2} \begin{eqnarray} &&\frac{\partial}{\partial t}\left(-\Sigma \frac{r^2\Omega^2}{2}\right) +\frac{1}{r}\frac{\partial}{\partial r}\left[r\Omega\left\{ \frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}} (\rho c_{\rm s}^2)_{\rm mid} \right\}\right. \nonumber\\ && \left. \qquad+ r^2\Omega\Sigma\overline{\alpha_{r\phi}} c_{\rm s}^2\right] + (\rho v_z)_{\rm w} E_{\rm w}+ F_{\rm rad} =0, \label{eq:totezint2} \end{eqnarray}(B.8)By multiplying Eq. (10) by r2Ω2/ 2, we have ∂t(Σr2Ω22)r2Ω2∂r[1rΩ{∂r(r2Σαcs2)+r2αφz(ρcs2)mid}]\appendix \setcounter{section}{2} \begin{eqnarray} &&\frac{\partial}{\partial t}\left(\Sigma \frac{r^2\Omega^2}{2}\right) -r^2\Omega^2\frac{\partial}{\partial r} \left[\frac{1}{r\Omega}\left\{\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}} (\rho c_{\rm s}^2)_{\rm mid} \right\}\right] \nonumber\\ &&\qquad + (\rho v_z)_{\rm w}\frac{r^2\Omega^2}{2} =0. \label{eq:sgmro2} \end{eqnarray}(B.9)

By combining Eqs. (B.8) and (B.9), we finally obtain a simple relation for the energetics of disc wind, Eqs. (14) and (15)

( ρ v z ) w ( E w + r 2 Ω 2 2 ) + F rad = Ω r [ ∂r ( r 2 Σ α c s 2 ) + r 2 α φz ( ρ c s 2 ) mid ] 1 r ∂r ( r 2 ΣΩ α c s 2 ) = 3 2 ΩΣ α c s 2 + r Ω α φz ( ρ c s 2 ) mid . \appendix \setcounter{section}{2} \begin{eqnarray} & & (\rho v_z)_{\rm w}\left(E_{\rm w}+\frac{r^2\Omega^2}{2}\right) + F_{\rm rad} \nonumber\\ &&\qquad =\frac{\Omega}{r}\left[\frac{\partial}{\partial r}(r^2 \Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2) + r^2 \overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid}\right] -\frac{1}{r}\frac{\partial}{\partial r}(r^2 \Sigma\Omega \overline{\alpha_{r\phi}}c_{\rm s}^2) \nonumber \\ &&\qquad =\frac{3}{2}\Omega\Sigma\overline{\alpha_{r\phi}}c_{\rm s}^2 + r\Omega \overline{\alpha_{\phi z}}(\rho c_{\rm s}^2)_{\rm mid} . \nonumber \end{eqnarray}

When the disc wind is neglected, (ρvz)w = 0, αφz=0\hbox{$\overline{\alpha_{\phi z}}=0$}, Eq. (15) is simplified to σSBT4=34ΩΣαcs2,\appendix \setcounter{section}{2} \begin{equation} \sigma_{\rm SB} T^4 = \frac{3}{4}\Omega\Sigma \overline{\alpha_{r\phi}}c_{\rm s}^2, \label{eq:visdsk} \end{equation}(B.10)where we use Eq. (13). Since the mass accretion rate is approximated as r=2πΣrvr2πΣr(αcs2/rΩ)\hbox{$\dot{M}_r=-2\pi \Sigma r v_r \approx 2\pi \Sigma r (\overline{\alpha_{r\phi}}c_{\rm s}^2/r\Omega)$}, Eq. (B.10) is rewritten as σSBT4=38πrΩ2=38πGMrr3,\appendix \setcounter{section}{2} \begin{equation} \sigma_{\rm SB} T^4 = \frac{3}{8\pi}\dot{M}_r \Omega^2 = \frac{3}{8\pi} \frac{GM_{\star} \dot{M}_r}{r^3}, \end{equation}(B.11)which is consistent with the expression for the standard accretion disc in the outer region (Shakura & Sunyaev 1973).

All Tables

Table 1

Parameters for MRI-active cases.

Table 2

Parameters for MRI-inactive cases.

All Figures

thumbnail Fig. 1

Comparison of time evolutions of four MRI-active PPDs with α=8×10-3\hbox{$\overline{\alpha_{r\phi}}=8\times 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 the text
thumbnail 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).

In the text
thumbnail Fig. 3

Mass-loss rate by disc wind, z, (solid lines) and mass accretion rate induced by the 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 5

Same as Fig. 1 but for four MRI-inactive cases with α=8×10-5\hbox{$\overline{\alpha_{r\phi}}=8\times 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)).

In the text
thumbnail Fig. 6

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

In the text
thumbnail Fig. 7

Comparison of αφz\hbox{$\overline{\alpha_{\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.

In the text
thumbnail Fig. 8

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

In the text
thumbnail Fig. 9

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

In the text
thumbnail Fig. 10

Comparison of normalized pressure gradient force, (1ρmidpmid∂r)/(2rΩ2)\hbox{$-\left(\frac{1}{\rho_{\rm mid}}\frac{\partial p_{\rm mid}}{\partial r}\right) /(2r\Omega^2)$}, 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 α=8×10-5\hbox{$\overline{\alpha_{r\phi}} =8\times 10^{-5}$} (black lines).

In the text
thumbnail 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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.