{'fraction_non_alphanumeric': 0.03403467180323011, 'fraction_numerical': 0.02181063861312787, 'mean_word_length': 4.429686242960579, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Protostellar Disk Evolution Over Million-Year Timescales with a Prescription for Magnetized Turbulence

Magnetorotational instability (MRI) is the most promising mechanism behind accretion in low-mass protostellar disks. Here we present the first analysis of the global structure and evolution of non-ideal MRI-driven T-Tauri disks on million-year timescales. We accomplish this in a 1+1D simulation by calculating magnetic diffusivities and utilizing turbulence activity criteria to determine thermal structure and accretion rate without resorting to a 3-D magnetohydrodynamical (MHD) simulation. Our major findings are as follows. First, even for modest surface densities of just a few times the minimum-mass solar nebula, the dead zone encompasses the giant planet-forming region, preserving any compositional gradients. Second, the surface density of the active layer is nearly constant in time at roughly 10 g/cm2, which we use to derive a simple prescription for viscous heating in MRI-active disks for those who wish to avoid detailed MHD computations. Furthermore, unlike a standard disk with constant-alpha viscosity, the disk midplane does not cool off over time, though the surface cools as the star evolves along the Hayashi track. The ice line is firmly in the terrestrial planet-forming region throughout disk evolution and can move either inward or outward with time, depending on whether pileups form near the star. Finally, steady-state mass transport is a poor description of flow through an MRI-active disk. We caution that MRI activity is sensitive to many parameters, including stellar X-ray flux, grain size, gas/small grain mass ratio and magnetic field strength, and we have not performed an exhaustive parameter study here.


INTRODUCTION
A typical, low-mass T-Tauri star accretes mass at a rate of 10 −8 M yr −1 -one Jupiter mass every 100,000 years (e.g. Hartmann et al. 1998;Sicilia-Aguilar et al. 2004). Balbus & Hawley (1991) put forth the magnetorotational instability (MRI) as the most likely driver of accretion in T-Tauri disks, followed by Brandenburg et al. (1995), Hawley et al. (1995), Balbus et al. (1996), and Balbus & Hawley (1998). Before a first-principles physical description of angular momentum transport was available, accretion was often modeled using the αprescription (Shakura & Syunyaev 1973), The α-prescription relates turbulent viscosity to length and velocity scales in the disk based on dimensional analysis. In Equation 1, ν is the turbulent viscosity, c s is the russell.landry@gmail.com sound speed and H is the pressure scale height. α is a dimensionless efficiency that is often assumed, without physical motivation, to be uniform throughout the disk. The first numerical investigations of MRI-driven turbulence were local shearing box simulations (Hawley et al. 1995), which treat a box of approximate size 2πH × H × H centered at a given distance from the star. The shearing box is still a useful technique for investigating detailed properties of turbulence. However, more recent investigations computed turbulent viscosity from first principles for global disk models, which were previously the domain of the α-prescription. Such global simulations confirmed that the MRI can lead to turbulence-driven accretion matching observed rates of 10 −8 M yr −1 in either unstratified disks (Hawley 2001;Steinacker & Papaloizou 2002;Lyra et al. 2008) or thin, stratified disks (Sorathia et al. 2010). MRI turbulence is self-sustaining for the simulation timeframe of around 1000 orbits near the inner boundary (Fromang & Nelson 2006;Flock et al. 2011).
In parallel with the development of global accretion disk models came investigations of the behavior of the MRI in non-ideal, partially ionized fluids. Gammie (1996) was the first to point out that protostellar disks have surface layers ionized well enough to couple to magnetic fields, and an interior dead zone where extremely low ionization levels prevent magnetically driven turbulence. Subsequent investigations incorporating Ohmic resistivity into the MHD equations confirmed Gammie's prediction that a high enough resistivity would dampen the growth of the MRI (Jin 1996;Sano & Miyama 1999). Numerical simulations showed a "dead zone" near the midplane and near the star, where the surface density is high enough to shield the disk interior from ionizing radiation (Fleming et al. 2000;Sano et al. 2000). Recent research has revealed that the dead zone is not entirely dead: the shielded interior still experiences a shear stress only about an order of magnitude less than the active layer due to propagating acoustic waves (Fleming & Stone 2003) and smooth, large-scale magnetic fields (Turner et al. 2007;Turner & Sano 2008). Other nonideal effects that affect the MRI include the Hall current (e.g. Wardle 1999;Balbus & Terquem 2001;Sano & Stone 2002) and ambipolar diffusion (e.g. Blaes & Balbus 1994;Mac Low et al. 1995;Kunz & Balbus 2004;Bai & Stone 2011) (See Section 2.1 for a more complete description of each effect). Improved understanding of both non-ideal MRI and ionization in disks (Igea & Glassgold 1999;Ilgner & Nelson 2006) has provided the tools to describe the MRI in gas of any ionization fraction or density that can be found in a protostellar disk.
Detailed short-timescale snapshots have now been constructed of angular momentum transport in a protostellar disk, including the dead zone, turbulent layers and corona (Dzyurkevich et al. 2010;Kretke & Lin 2010;Bai 2011;Flaig et al. 2012). Yet it is possible that MRIactive T-Tauri disks may not be in steady state due to (a) the dead zone and its resulting, radially varying accretion rate and (b) the lack of a protostellar envelope to provide material to maintain steady, inward mass transport. Our protostellar disk snapshot, then, must change over the million-year timescales on which the star/disk system evolves. The goal of this work is to illustrate how disk interior structure and angular momentum transport change over the entire, multi-million-year lifetime of the T-Tauri phase. We extend the work of Armitage et al. (2001) and Zhu et al. (2010), who also performed millionyear simulations, by calculating a radially and vertically varying α based on the ionization state rather than assuming a constant α-value in the active zone. Martin et al. (2012b) also modeled FU Orionis outbursts using time-dependant global simulations of MRI-active disks, including Ohmic resistivity, using 1-D layered models which used the α-prescription for the active layer, and an analytical approximation for the active layer surface density (Martin et al. 2012a). Our models build on their work by adding a vertical dimension to the disk structure.
Our evolving model of a magnetically turbulent T-Tauri disk answers the following questions: 1. How do the relative sizes of the dead zone and active layers change over time?
2. How doesṀ vary with radius and time?
3. How can disk modelers parameterize heating in the active layers and dead zone without resorting to a 3-D MHD simulation?
4. Does the disk midplane heat up or cool off over time?
5. Where is the ice line in an MRI-active disk and how does its location change over time?
Questions 1, 2 and 3 elucidate basic properties of an MRI-turbulent accretion disk. Questions 4 and 5 highlight fundamental ways in which disk models based on the standard viscosity prescription with constant α lead us astray. Note that our disk model does not include photoevaporation, which is another process that operates on million-year timescales that can produce radially varying accretion accretion rates. Our ability to simulate an entire T-Tauri disk lifetime is due to a new MRI activity prescription that allows us to compute the thermal and viscous effects of MRI turbulence without resorting to 3-D magnetohydrodynamic simulations of the turbulence itself. We can thus reduce our computational domain from three spatial dimensions to 1+1 spatial dimensions-1-D vertical structures representing axisymmetric disk annuli that are connected only by a 1-D radial mass transport equation (Dodson-Robinson et al. 2009). Sacrificing information about small-scale turbulent fluctuations, we retain our ability to accurately describe large-scale structures such as the dead zone and active layers while dramatically improving our ability to simulate long timescales.
Our paper is organized as follows. In Section 2 we discuss the basic equations governing the MRI under nonideal MHD conditions and give our prescription for determining turbulent viscosity. We present our method of computing vertical strtucture and mass transport in Section 3. We outline a basic picture of MRI-turbulent disk evolution and answer Questions 1, 2 and 3 in Section 4. In Section 5, we discuss the differences in thermal structure between our model and constant-α disk models, which leads us to answer Questions 4 and 5. We discuss the limitations of our model in Section 6 and present our conclusions in Section 7. Readers who wish to skip over the details of the computations may wish to proceed directly to Section 4.

INSTABILITY-DRIVEN TURBULENCE IN PARTIALLY IONIZED GASES
When an accretion disk is fully ionized and the magnetic field is weak, the entire disk is MRI turbulent. Yet when a disk is only partially ionized, as is the case for a protoplanetary disk, there is an incomplete coupling between the disk gas and the magnetic field, and non-ideal effects become important to the growth of MRI-driven turbulence. In this section, we describe how we treat angular momentum transport from MRI turbulence in non-ideal, partially ionized gases. We begin by describing our turbulence criterion in Section 2.1, then list our method for determining the diffusion regime and resulting turbulent stress in Section 2.2.
2.1. MRI activity criteria in the three non-ideal regimes The non-ideal magnetic induction equation has three extra terms corresponding to the three non-ideal effects (Wardle 2007), Ohmic resistivity, the Hall effect and ambipolar diffusion: In Equation 2, v is the gas velocity, B andB are the magnetic field and magnetic field unit vector, and ⊥ refers to the component perpendicular to B. η O , η H , and η A are the Ohmic, Hall, and ambipolar diffusivities respectively. Ohmic resistivity dominates other nonideal effects when collisions with neutrals cause both electrons and ions to decouple from field lines. When collisional drag is sufficient to decouple ions and grains from the magnetic field, but not electrons, the relative velocity between the ions and electrons is non-negligible and the Hall effect is dominant. In the ambipolar diffusion regime, electrons and ions decouple from the neutral gas and the magnetic field lines are frozen to the charged species and drift through the neutral gas. When Ohmic resistivity is the largest non-ideal effect, the MRI will only occur if the Elsasser number, is at least of order unity (Sano & Stone 2002;Turner et al. 2007). In Equation (3), v Az is the Alfvén speed in the vertical direction and Ω is the Keplerian angular velocity. Physically, Λ is the ratio of the wavelength of maximum growth to the diffusive scale length. The tangled magnetic fields in MRI turbulence usually have a toroidal component with pressure 10 to 30 times greater than the pressure in the vertical component (Turner & Sano 2008), so the Alfvén speed in the vertical direction used to calculate Λ is where v A = B/ √ 4πρ is the total Alfvén speed. The Λ 1 criterion ensures that the most unstable mode can grow more quickly than the charged particles can diffuse across magnetic field lines. Simulations by Sano & Stone (2002) suggest the Hall effect does not change the conditions for turbulence if Ohmic diffusion is also present. Further work is required to determine the growth of the MRI in regimes where the Hall term is much stronger than other terms (Wardle 2012). Ambipolar diffusion arises from the relative motion of ions and neutral particles in the disk gas. In the "strong coupling" limit, in which ion density is negligible and electron recombination time is much smaller than the orbital time 1/Ω, as is generally the case for protoplanetary disks (Bai 2011), ion density cannot be assumed to follow the continuity equation. Instead, the ion density is determined by the ionization-recombination equilibrium, and characterized by the parameter Am (Chiang & Murray-Clay 2007): where ρ i is the ion density and γ is the neutral-ion drag coefficient, In Equation 6, σ ni is the effective cross section for neutral-ion collision and w ni is the relative velocity between neturals and ions. Physically, Am is the ratio of the orbital period to the collisional timescale between ions and neutrals. Since which is equivalent to the Elsasser number Λ in the Ohmic regime. Similarly then, Am is the ratio of the wavelength of maximum growth to the ambipolar diffusive scale length.
In their three-dimensional shearing-box simulations exploring the effect of ambipolar diffusion on MRI turbulence, Bai & Stone (2011) determine that heavily ionized yet tenuous disks can only sustain turbulence when threaded by weak magnetic fields. The magnetic field strength is characterized by the plasma β, the ratio of the gas pressure to the magnetic pressure: The requirement that the magnetic field energy be small in comparison to the gas thermal energy (a "weak" field) restricts MRI turbulence to values of the plasma β that are greater than a minimum (Bai & Stone 2011): The maximum field strength, beyond which the field is too strong to be destabilized for any given field geometry, decreases the more important ambipolar diffusion becomes (smaller Am). However, for a sufficiently weak field, MRI can be sustained even for Am 1. In a protoplanetary disk, ambipolar diffusion dominates in the atmosphere, which is diffuse and highly ionized by stellar X-rays and cosmic rays. Ambipolar diffusion is also important in the outer disk where the surface density is very low. Ohmic resistivity dominates in the dense inner region shielded from ionizing radiation. Ignoring the effects of the Hall diffusivity, which are unlikely to alter either the conditions required for MRI or the strength of the turbulence where Ohmic dissipation is also present (Sano & Stone 2002), a non-ideal protoplanetary disk differs from an ideal accretion disk through the possibility of a dead zone in the inner disk. The dead zone would remain cold and would not efficiently transport angular momentum. The accretion efficiency in the upper, ionized layers of non-ideal disks also lags behind ideal disks due to ambipolar diffusion effects.
In order to compute the disk viscosity and angular momentum transport properties, we need to know (a) whether the MRI is operating, and (b), if so, how strong the turbulence is. Since the MRI growth timescale is roughly 1/Ω, hundreds of thousands of times shorter than our ∼ 1 Myr simulation timescale, we assume the MRI is either fully saturated or completely damped. Our MRI turbulence criterion is therefore equivalent to that of Bai (2011): • If Λ ≥ 1 and β > β min , neutral gas couples to the magnetic field and MRI is saturated; • If Λ < 1 or β ≤ β min , neutral gas decouples from the magnetic field and MRI is damped.
In the next section, we discuss the computation of the magnetic diffusivities that determine Λ and Am.
2.2. The diffusion regime and turbulent stress We cannot apply our turbulence criterion without knowing the values of η O , η H and η A . For a given charged species j, the ratio of Lorentz force to the neutral drag force is where Z j e is the charge of j (negative or positive), B is the magnitude of the magnetic field, m j is the mass of j, c is the speed of light, and γ j is defined according to Equation 6. (Note that β j is not the same as the plasma β of Equation 8.) For each diffusion regime, one can define a conductivity by summing over all charged species: (Wardle 2007). In Equations 11-13, n j is the number density of species j. Finally, one can write the diffusivities according to where σ ⊥ ≡ σ 2 H + σ 2 P . To determine the equilibrium abundances of charged species n j , we solve a simplified set of chemical reactions from Model 4 of Ilgner & Nelson (2006), which we briefly motivate here. The set is derived from the following re-actions: where HCO + is a representative molecular ion, Mg + is a representative metal ion and g is a grain. Here every species (except the energetic particle X-a cosmic ray, Xray or radionuclide decay product) is created in at least one reaction, and destroyed in at least one other. Over the whole set, no species is produced or consumed on balance. The subset producing the ions and electrons reduces to That is, each energetic particle striking a hydrogen molecule yields one ion and one electron. In constructing the conductivity lookup tables we therefore approximate Equations 17-21 by neglecting the fact that the molecular ion contains just one hydrogen atom. Since HCO + is orders of magnitude less abundant than H 2 , forming ions leaves the H 2 density unchanged. Similarly, we don't model CO destruction and reformation because the ion is so much less abundant than the molecule. Equation 22 then becomes The simplified network consists of Eqs. 25, 26, 27 and 23, together with the grain surface reactions described by Ilgner & Nelson (2006). The metal atoms' thermal adsorption and desorption on the grains is included. Ilgner & Nelson (2006) found that this reduced network yields similar results to a detailed version including hundreds of species and thousands of reactions, in the most common situation where the recombination occurs mostly on the grains. The internal grain density, gas/small grain ratio, and grain size used in the chemical reaction network are listed in Table 1. Here we deviate from the standard interstellar gas/small grain mass ratio of 100 and assume some grain growth has occurred, so that 90% of the grain mass is in grains larger than 1µm. Using the standard gas/dust ratio of 100 resulted in no MRI turbulence (see section 4.1). To avoid having to run the chemical reaction network at every timestep of our million-year simulations, we followed the approach of Flaig et al. (2012) and created a look-up table of magnetic diffusivities as a function of temperature T , gas density ρ, ionization rate ζ and plasma β. To generate the look-up table, we ran the chemical reaction network until it reached equilibrium abundances of all species for each combination of T , ρ, ζ and β. We then computed the conductivities and tabulated diffusivities according to Equations 10-16.
The rate coefficient for the reaction H 2 + X → H + 2 + e − is, of course, the ionization rate ζ. To determine ζ, we consider cosmic rays, stellar X-rays, and short-lived radionuclides (τ 1/2 10 8 yr), including 26 Al. Following Umebayashi & Nakano (2009), we take the ionization rate from short-lived radionuclides as ζ R = 7.6 × 10 −19 s −1 and calculate the attentuated cosmic ray ionization as: where ζ surf CR is the unattenuated cosmic ray ionization rate, λ CR is the cosmic ray penetration depth (Umebayashi & Nakano 1981), and Σ 1,2 (z) are the mass columns above and below the vertical height z. Values of ζ surf CR , λ CR and all other numerical inputs to our model are listed in Table 1. Finally, following Bai & Goodman (2009), we calculate the stellar X-ray ionization rate ζ X as: where L X,29 ≡ L X /(10 29 erg s −1 ), and L X is the the stellar X-ray luminosity. We take L X,29 = 20 to match the young solar-mass stars observed in the Orion Nebula (Garmire et al. 2000). Here we keep the stellar X-ray flux constant in time, though it could certainly vary in either a smooth, systematic way with age or stochastically with accretion bursts. All parameters in Equation 29 are listed in Table 1. The first exponential represents attenuation of X-rays by absorption, and the second represents the contribution from scattered X-rays. We show vertical profiles of the ionization rate at two different disk radii in Figure 1. After the ionization rate is determined at each (R, z) zone in our disk model, an interpolation through the look-up table can return the magnetic diffusivities. We can then compute the Elsasser number Λ and Am and determine whether MRI turbulence in the zone is active or not, according to our turbulence criterion (Section 2.1). The last remaining ingredient in our MRI prescription is a rule for determining the strength of the turbulence, where it is present. Where MRI is saturated, we use the scaling relations found between turbulent stress and magnetic field strength in a variety of shearing-box simulations (Hawley et al. 1995 Since there is an accretion flow caused by large-scale magnetic fields even in the dead zone (Turner & Sano 2008;Nelson & Gressel 2010), we set a minimum value of α where the MRI is damped. In the active layer, α is close to its maximum value of ∼ 0.5, set by the cessation of MRI in the strong-field limit of β ≥ 1 (Hawley et al. 1995;Balbus & Hawley 1998). Bai & Stone (2011) found a similar result in the ambipolar regime, α ≈ 0.4 for Am → ∞. The shear stress in the dead zone is an order of magnitude less than in the active layer (Fleming & Stone 2003;Turner et al. 2007;Turner & Sano 2008), and the plasma β at the midplane is typically two to three orders of magnitude higher than at the top of our grid (see Figure 1, which shows profiles of the vertical component of β, β z = 8πP/B 2 z , for two different disk radii). α min must be therefore be roughly four orders of magnitude less than α max to produce an appropriate level of dead zone shear stress. We take α min = 10 −5 (Table 1).

DISK STRUCTURE AND MASS TRANSPORT
With our method of computing turbulent viscosity at any location in a protostellar disk, we may now simulate how the disk re-distributes its mass throughout its multi-million-year existence. Since our viscosity prescription depends on the four inputs T , ρ, β and ζ, which are functions of radius and vertical height (R, z), we must compute the detailed vertical and radial structure of the disk. Our computational setup, similar to the disk models of Dodson-Robinson et al. (2009), is based on the following simplifying assumptions: 1. The disk is axisymmetric and symmetric about the midplane; 2. The disk is geometrically thin, H/R 1; 3. Heat escapes in the vertical direction much faster than it is carried with the gas flow in the radial direction.
Assumption 1 reduces the physical three-dimensional disk to a two-dimensional quadrant with zero flux at the midplane (required for symmetry). Neglecting the azimuthal dimension is a critical step in speeding up the code to allow long-timescale simulations. Assumption 2 allows the vertical and radial dimensions to be decoupled in a 1+1-D framework, so that energy transport proceeds only in the vertical direction. Each radial gridpoint contains an independent vertical structure model. Assumption 2 is valid as long as the vertical sound-crossing time (of order the orbital timescale) is much less than the accretion timescale, which is generally true in T-Tauri disks.
In section 3.1 we describe our vertical structure model, while in Section 3.2 we discuss our mass transport parametrization. Section 3.3 contains a description of our computational methods.
3.1. Vertical structure Hydrostatic equilibrium and the thermal balance between stellar irradiation, radiative cooling and viscous heating govern the vertical structure of our disk. We use the flux-limited diffusion approximation for the transport of viscously generated energy. The accretion flux gradient is determined by the viscous energy generation rate per unit volume (Pringle 1981): In Equation 32, Ω is the height-dependent Keplerian frequency Temperature and pressure are related by the ideal gas equation of state, where µ = 2.33 g mol −1 is the mean molar weight of the disk gas (Table 1). The temperature and pressure gradients are In Equation 35, T acc is the temperature contribution from viscous heating only-stellar irradiation and the ambient molecular cloud also contribute some of the thermal energy. T is the true temperature and includes all heat sources.
To calculate ∇ ≡ d ln T /d ln P , the thermodynamic gradient, we use the Schwarzschild criterion for stability against convection: ∇ ad = 2/7 is the adiabatic thermodynamic gradient for diatomic gas and is the radiative thermodynamic gradient, where a is the radiation density constant and κ is the local Rosseland mean opacity. For full details on how to compute ∇ conv in the disk's convective zone, see Kippenhahn & Weigert (1994). Computing ∇ requires the local Rosseland mean opacity, κ. At low temperatures, T < 700 K, we use the opacities of Semenov et al. (2003) calculated for a 5layered sphere topology. At higher temperatures where molecular gas dominates opacity (T > 1000 K), we use the tables of Ferguson et al. (2005). For 700 K ≤ T ≤ 1000 K, we interpolate between the two tables using a weighted average in log(T ) space. For more details, and plots of the resulting opacities, see Dodson-Robinson et al. (2009).
Integrating the coupled ODEs in Equations 32, 35 and 36 requires two different temperatures: T , the true temperature resulting from all sources of thermal energy, and T acc , the component from viscous heating only. The other heat sources are the central star, which sets an equilibrium temperature component T eq , and the ambient star-forming region, from which long-wavelength radiation that penetrates the disk sets a minimum temperature T amb . To compute T eq , we begin by assuming the disk surface is flared as a result of hydrostatic equilibrium and radiative equilibrium with the star. Following the models developed by Chiang & Goldreich (1997), we calculate the grazing angle θ at which stellar energy enters the disk: where R * and T * are the star's radius and effective tem-perature and T c is a measure of the gravitational potential at the surface of the star: In Equation 40, M * is the star's mass and σ is the Stefan-Boltzmann constant. As the disk evolves, we determine the star's temperature and radius as a function of age from the pre main-sequence evolutionary tracks of D' Antona & Mazzitelli (1994). Since the star's luminosity decreases as it moves down the Hayashi track, the disk flaring becomes less pronounced over time and the disk surface, where T ≈ T eq , cools.
In the direction parallel to the disk midplane, the stellar radiation penetrates to an optical depth τ * ∼ 1. The asterisk subscript denotes that this optical depth is measured at the peak wavelength of the starlight, near 1 µm. Measured perpendicular to the disk midplane, stellar radiation is mostly attenuated by optical depth τ ⊥ * ∼ τ * θ ∼ θ. The equilibrium temperature between the stellar heating and radiative cooling at the disk surface-neglecting viscous heating-is (D'Alessio et al. 2006;Natta et al. 2000): where τ surf is the Rosseland mean optical depth at the disk surface for blackbody radiation at T eq,surf . The top surface of our vertical grid is defined by τ surf = 0.2. Since half of the stellar radiation absorbed by grains at the disk surface is re-radiated into space, the equilibrium temperature with the star as a function of height in the disk-again, neglecting viscous heating-is where τ z is the Rosseland mean optical depth to height z: Finally, the true temperature T simply the flux sum of the individual temperature components: 3.2. Radial diffusion Since the radial and vertical dimensions of our disk model are decoupled, we must treat mass transport as a one-dimensional problem. Yet the MRI-active, partially ionized disk is vertically layered, with the most active accretion occurring at the surface. The key to a successful 1-D description of layered accretion is the fact that vertical re-distribution of mass within an annulus occurs more quickly than accretion in the active layers: 1/Ω R 2 /ν active . By computing a mass-weighted, vertically averaged value of turbulent viscosity in each annulus,ν where Σ is the surface density in the annulus, we can describe mass transport using the radial diffusion equation: In each (R, z) zone, we compute viscosity according to Equation 1. We use a height-dependent modified scale height H, softened into a non-singular form (Milsom et al. 1994).
In MRI-active zones, α is given by Equation 31, while in inactive zones we set α to our chosen value of α min ( Table 1). The sound speed used to compute ν(R, z) is

Computational methods and initial conditions
To initialize the disk evolution model, we compute the vertical structure of a disk with the following features: 1. A surface density profile Σ ∝ R −3/2 , predicted by Zhu et al. (2010) for layered accretion disks in the T-Tauri phase.
2. A pre-main-sequence star with a mass of 0.95 M and an initial age of 0.1 Myr, which roughly coincides with the beginning of the T-Tauri phase (Dunham & Vorobyov 2012). The star will continue to accrete a small amount of mass from the disk during the ∼ 3 Myr T-Tauri phase.
3. An ambient temperature of 20 K (Table 1) to match the typical background temperatures of infrared dark clouds (Peretto et al. 2010).
4. An outer radius of 70 AU, set by the 80 AU solar nebula size limit of Kretke et al. (2012).
(Note that the disk expands from its initial radius.) Kretke et al. show that a solar nebula with R out > 80 AU would excite Kozai oscillations in some of the planetesimals scattered by Jupiter and Saturn, stranding them in stable, high-inclination, low-eccentricity orbits that surveys have not detected.
5. An inner radius of R in = 0.5 AU. The requirement that the disk stay below the dissociation temperature of H 2 , so that the ideal gas equation holds, dictates our choice of R in . Ruden & Lin (1986) find that the exact value of R in does not affect the overall disk structure as long as R out R in . 6. A disk mass of either 0.015M (Model 1) or 0.03M (Model 2). Our disk masses are designed to be comparable to previous global-disk MRI simulations (e.g. Lyra et al. 2008;Dzyurkevich et al. 2010;Bai 2011), almost all of which use minimummass solar nebulae. Note that such low-mass disks are probably not viable giant planet-forming environments (e.g. Thommes et al. 2008). In a forthcoming study, we will examine the evolution of MRI-active, high-mass, planet-forming disks.
We then evolve the disk forward in time using Equation 46. We use fully implicit finite differencing, adjusting the timestep ∆t so that surface density varies by a maximum of 1.0% during a single timstep. The inner boundary R in experiences zero stress, such that matter falls directly from R in onto the star. The disk is allowed to expand freely from the outer boundary R out , with four new zones added to the disk each diffusion time of R 2 out /ν(R out ). At each timestep, we independently calculate the vertical structure for each zone in the radial grid. We begin the vertical structure solver with initial guesses of T acc , ρ and vertical component of the magnetic pressure P B,z = B 2 z /8π at the top of our grid, defined by τ surf = 0.2. We find the height of the grid surface by The accretion flux at the grid surface is F acc (z surf ) = σT 4 acc . We use a fourth-order Runge-Kutta integrator with adaptive stepsize control to integrate the coupled ODEs in Equations 32, 35 and 36 from the surface to the midplane. The vertical magnetic pressure stays constant in height, though it varies with radius. A solved vertical structure model has the properties required to keep the disk symmetric about the midplane, so that the volume densities add up to the surface density in the annulus, and β(z = 0) = 1000.
Equation 52 requires that the plasma β be constant throughout the disk midplane. After turbulence is established, Fromang & Nelson (2006) find midplane values of 100 β 1000 for a wide range of vertical box sizes, resolutions and boundary conditions in global ideal MHD calculations. We use the Newton-Raphson algorithm (Press et al. 1992) to adjust the initial guesses of T acc , ρ and P B,z until a solution is found that satisfies Equations 50, 51, and 52.
The on-off nature of the MRI creates discontinuities in and ν(z) and F acc (z) that can cause the Newton-Raphson algorithm to oscillate between two sets of input parameters that bracket the correct solution of Equations 50, 51 and 52. To avoid such oscillations in the Ohmic regime, we decrease the value of α gradually in the range 1.6 ≥ Λ ≥ 0.4, using a sigmoid function: where k = 1.5 ln(1/(α min β) − 1). We use a similar sigmoid function to smooth α in the ambipolar regime for 0.4β min ≤ β ≤ 1.6β min .

MASS TRANSPORT IN MRI-ACTIVE DISKS
Here we present our simulations of the evolution of magnetically turbulent disks over million-year timescales. In §4.1, we discuss the relative sizes of the dead zone and active layers (Question 1 in Introduction). In §4.2, we demonstrate how the shrinking of the active layer over regions of high density enhances mass pileup in the dead zone. We also analyze the radial mass floẇ M (R) through the disk and show that the disk never reaches a steady state, even on million-year timescales (Question 2 in Introduction). In §4.3, we give a simple prescription for accretional heating due to MRI for use with semi-analytical, non-evolving disk models to predict observables (Question 3 in Introduction). Figure 2 shows viscosity as a function of (R, z) for the inner 30 AU of Model 1 (0.015M , top) and Model 2 (0.03M , bottom). Outside of 30 AU, viscosity is almost independent of height z. The left-hand panels of Figure  2 show the disks after 10 4 years of simulation timea star age of 0.11 Myr, since we began the simulations with a star at the beginning of the T-Tauri phase, age 0.1 Myr (Dunham & Vorobyov 2012 2. Although the vertical heights of both the dead zone and the overall disk shrink with time as the disk loses mass and cools (see Section 5), the radial extent of the dead zone stays approximately constant in time.

Turbulent Structure
The radial size of the dead zone is larger than the ∼ 5 AU typically quoted for disks similar to the minimummass solar nebula (MMSN) (Matsumura & Pudritz 2003;Salmeron & Wardle 2008;Flaig et al. 2012). Here we define the dead zone as the region of the disk where Λ < 1. The primary reason our dead zone is so extensive is the surface density of our disks: the MMSN contains 0.01M within 100 AU of the sun, while our least massive disk contains 0.015M within 70 AU. Our dead zone is also extensive because we treat ambipolar diffusion, as well as Ohmic diffusion. Sano et al. (2000) use the same tenfold dust depletion we do, but treat only Ohmic diffusion and find a dead zone that extends 5-10 AU.
Dead zones with α 10 −5 are required for maintaining compositional gradients such as the ice gradient in the asteroid belt, which would radially diffuse on million-year timescales for fully turbulent disks (Nelson & Gressel 2010). The fact that Uranus and Neptune have atmospheres with higher carbon enrichment than Jupiter and Saturn provides some evidence that the solar nebula had compositional gradients covering the entire giant planetformation region (Encrenaz 2005). Given that the dead zone in a disk of just 0.015M can reach ∼ 15 AU, the outer boundary of the giant planet-formation region in the Nice model (Tsiganis et al. 2005), the solar nebula could have supported such a compositional gradient. The lack of change in the dead zone radius over time suggests that compositional gradients are stable over million-year timescales. Figure 1 gives further insight into the vertical structures of Model 1 (top) and Model 2 (bottom). In each panel, the active layer is shaded with red dots. At 10 AU, both disks show the classical layered structure of an MRIactive zone on tope of a dead zone. At 1.5 AU, however, we see additional complexity in the disk's vertical structure. At 45 mG (Model 1) and 70 mG (Model 2), the magnetic tension is strong enough to prevent MRI from bending the field lines in the low-density surface gas despite a high ionization rate of 10 −10 s −1 , forming a stable corona (e.g. Miller & Stone 2000). Our stabilizing magnetic field values are roughly consistent with those of Salmeron & Wardle (2008), who calculated that unstable MRI modes can only grow for B 80 mG in the presence of 1µm grains. (See §2.2 and Table 1 for more on the grain properties used in our models.) Figure 3 is a zoom-in on the stress coefficient α(R, z), viscosity ν(R, z) and ionization rate ζ(R, z) in the inner 4 AU of Model 2 after 1 Myr. The inactive corona is apparent for radii R 3 AU, but moves above our computed disk surface (the location where τ * ∼ 1) as the magnetic field stength decreases with radius.
Also noticeable in Figures 1 and 3 is a split active layer in the inner part of Model 2. MRI-active regions of high α and ν sandwich a "dead slice" that has reduced stress and viscosity by an order of magnitude. (The fact that α does not immediately plunge to its minimum value in the dead slice is a result of the sigmoid smoothing described in §3.3.) While not present at the start of our simulations, the split active layer appears after only 5000 years of disk evolution and persists until the end of the simulation at 3 Myr. The split active layer extends from roughly 1.2 AU < R < 1.7 AU (Figure 3), though its radial extent shrinks slightly as the disk evolves.
The dead slice, in the part of the disk where ambipolar diffusion is the strongest non-ideal MHD term, is the result of two competing effects. First, MRI in the ambipolar regime requires high ionization: decreasing ζ toward the disk midplane reduces Am and shuts down turbulence. Yet ambipolar diffusion is quenched at high densities: increasing ρ toward the midplane increases Am, favoring turbulence. In the dead slice, the dropping ζ temporarily dominates over the increasing ρ and shuts down the MRI. Model 1 does not have a dead slice because its surface density is about 1/2 that of Model 2: ζ can stay high enough for the MRI to operate until very near the midplane, where Ohmic diffusion begins to dominate. An open question is whether or not a dead slice would be present in a fully 3-D simulation with identical vertical profiles of ρ, ζ and β to our disk: the thickness of the dead slice is smaller than the MRI wavelength by a factor of 2-10, so turbulence would very likely erase it.
We have seen that MRI-active disks may have complex vertical structure, with multiple layers of turbulent and non-turbulent zones. In the next section we explore the overall mass flow through the disk and discuss its evolution on million-year timescales.

Mass Flow
As might be expected for a layered disk with a dead zone, radial mass transport is not in steady state for either Model 1 or Model 2. The left-hand panel of Figure  Fig for both models after 1 Myr of evolution, where v R is a density-weighted, vertically averaged gas radial velocity. The convention in Equation 54 is that v R is positive when gas flows toward the star and negative when gas flows away from the star. In the inner part of the disk, where gas flows inward, the highest accretion rates occur (a) at the outer edge of the dead zone (at 21 AU for Model 1 and 16 AU for Model 2), and (b) at the inner disk boundary. There is a clear drop inṀ (R) associated with the dead zone. (The noisyṀ profile where the disk has two or more layers is a reflection of the root-finding tolerance in our Newton-Raphson algorithm. Fluctuations in the location of the boundary between dead and active zones are random and average out over time.) The potential for the dead zone to become gravitationally unstable, given enough time to accumulate mass, is clear (Gammie 1996;Armitage et al. 2001;Zhu et al. 2010)-though we see only a slow, steady density growth at ∼ 1 AU-4 AU over the course of 3 Myr in Model 2 ( Figure 5). We will examine the potential for gravitational instability in high-mass, MRI-active disks in a forthcoming paper. One important caveat is that the accumulation of mass in the dead zone is unstable to the Rossby wave instability (Hawley 1987;Li et al. 2001;Meheut et al. 2012), which triggers spiral density waves. Without azimuthal information in our model, we cannot model the effect of the Rossby wave instability on dead-zone overdensity. Ultimately, the overdensity may not survive and may break up into large-scale vortices (Papaloizou & Pringle 1985;Lovelace et al. 1999).
At t = 1 Myr, there is turn-around in the accretion flow near 40 AU in each model. The outward mass flow outside 40 AU is necessary for overall conservation of angular momentum as material in the inner disk moves toward the star. The fact that there must be a change in the accretion flow direction also follows directly from the diffusive nature of T-Tauri disks (e.g. Lynden-Bell & Pringle 1974), which lack a circumstellar envelope to feed steadystate inward accretion. In both models, the disk has expanded beyond its initial 70 AU size after 1 Myr. Note that the 40 AU region of the disk does not drain, as the "turn-around radius" where the accretion flow changes direction moves inward with time (right-hand panel of Figure 4). We heartily discourage the use of a steadystate accretion rate for any reasonable description of angular momentum transport in T-Tauri disks. However, for predicting observables to order-of-magnitude using a static, non-time-varying disk model, we suggest approximating the accretional heating in MRI-active disks with a constantṀ ∼ 10 −9 M yr −1 . Figure 4 suggests that at a given time,Ṁ has a modest dependence on disk mass.
How consistent are our modeled accretion rates with observations? An oft-quoted value ofṀ * , the T-Tauri star accretion rate, is 10 −8 M yr −1 (Sicilia- Aguilar et al. 2004). Accretion rates in transitional disks are roughly 3×10 −9 M yr −1 (Espaillat et al. 2012). Our simulations achieve an accretion rate of roughlyṀ * ∼ 10 −9 M yr −1 , more consistent with the median transitional disk accretion rate. Reproducing high accretion rates in the inner regions of MRI-active disks is a delicate balance between magnetic field strength, X-ray ionization, grain size, and gas/small grain mass ratio. For example, we performed simulations using the same disk masses and X-ray luminosity (0.015M and 0.03 M ; 2 × 10 30 erg s −1 , respectively) but a grain size of 0.1µm instead of 1µm and a standard gas/small grain mass ratio of 100, and quenched all MRI activity in the disk entirely. However, as models by Zsom et al. (2011) show grain growth and settling within 1000 years, it is reasonable to assume the gas/small grain mass ratio has evolved from the interstellar value by the T-Tauri phase. In our simulations, a stronger magnetic field (lower value of β at the midplane) also suppressed MRI-driven accretion, which occurs in regions of the disk dominated by ambipolar diffusion. Most simulations of non-ideal MRI-driven accretion have difficulty reaching the 10 −8 M yr −1 benchmark (e.g. Zhu et al. 2010;Bai 2011), though there is a substantial amount of scatter in both T-Tauri and transitional disk accretion rates (Romero et al. 2012).

Simple prescription for accretional heating in an
MRI-active disk Previous studies have often relied on static, nonevolving disk models to connect measured line fluxes, velocities or spectral indices to physical properties of disks (e.g. D' Alessio et al. 2006;Pinte et al. 2010). Unfortunately the constant-α prescription for turbulent angular momentum transport leads to unphysical thermodynamic descriptions of disk midplanes. Predicted midplane temperatures in a constant-α model, where most turbulence is concentrated at the midplane, are too high. The poor match of the constant-α model to realistic protostellar disk accretion is definitely a problem for observations that trace disk midplanes, such as sub-millimeter measurements of continuum emission from large grains or surveys of rare molecules like HCO + or HD. However, the problem also affects observations that trace surface layers, such as Spitzer emission lines, infrared SEDs and sub-millimeter maps of abundant gases like CO. If a significant subset of T-Tauri disks rely on MRI to drive accretion, overestimating their midplane temperatures leads to overpredicting the photosphere and pressure scale heights, underpredicting the optical depths of low-excitation lines such as CO (J = 2 → 1), and possibly underpredicting the rates of grain settling and growth.
Here we use our simulation results to present a simple prescription for α in an MRI-active disk. First we determine the depth of the MRI-active layer. Figure 6 shows the surface density of the active layer as a function of radius for four time snapshots of Model 1 (left) and Model 2 (right). Clearly, the functional form of Σ active (R) is the same for both models and does not vary significantly with time. Near the star, where X-ray ionization is important, Σ active is high but falls off quickly as X-ray irradiation declines (see Equation 29). Over the dead zone, in the region where only cosmic-ray ionization is important, Σ active (R) increases as the magnetic field weakens, shrinking the inactive corona. Σ active (R) reaches its maximum at the outer edge of the dead zone. Here, where the entire vertical column is active, the decreasing depth of the active layer simply reflects the fact that Σ(R) is a decreasing function.
Although Σ active does vary with R in the parts of the A plot ofṀ (R) for both disk models after 1 Myr of disk evolution reveals that low-mass MRI-driven accretion disks have not reached steady state after 1 Myr. As required by conservation of angular momentum, there is a turn-around in the mass flow between inward accretion and outward decretion. The highest accretion rates in the inner disk, where mass moves toward the star, occur on either side of the dead zone. Right: The location of the turn-around radius as a function of time for both models. The turn-around radius moves inward over time. disk with a dead zone, the variation is less than a factor of two for both models presented here. Similarly, Σ active (R) is slightly higher for Model 2 than for Model 1, suggesting that the active column may increase modestly with disk mass. Σ active ≈ 10 g cm −2 is a good approximation for low-mass disks with surface density less than about four times the MMSN. The semi-analytical viscosity prescription is then simple: 1. Set the depth of the active column to Σ active ≈ 10 g cm −2 . Set α ≈ 0.01 in the active column and α ≈ 10 −5 in the dead zone. While there may be a corona on the inner disk surface, it contains very little mass and may safely be ignored as long as the disk model includes stellar heating.
2. Smooth the transition between the active layer and the dead zone if desired.
3. The entire vertical column will be active where Σ ≈ 20 g cm −2 , such that the two active layers meet at the midplane. Outside the outer radius of the dead zone, where Σ ≈ 20 g cm −2 , use α ≈ 0.01.
4. Calculate viscosity ν at each (R, z) in the disk using Equation 1.

THERMAL EVOLUTION OF MRI-ACTIVE DISKS
Now that we understand how mass moves through a disk where accretion is driven by MRI, we turn our attention to the thermal structure and evolution of the disk.
Here we see some important differences from disk models that use a constant-α viscosity prescription. In Section 5.1 we discuss how different parts of MRI-active disks heat up and cool off over time (Question 4 in Introduction). In Section 5.2 we investigate the location of the ice line and its motion through an MRI-active disk (Question 5 in Introduction).

Disk Heating and Cooling
The top panels of Figure 7 show the surface and midplane temperatures of Models 1 and 2 after 10 thousand years of evolution (left) and 1 Myr of evolution (right). First, note the obvious feature that the surface is far hotter than the midplane throughout most of the disk. Disk models with constant-α viscosity prescriptions usually have warm midplanes, T 100 K, in the inner 5 AU. In the inner 1 AU of constant-α models, the midplane temperature can even exceed the surface temperature, approaching 1000 K (e.g. Hersant et al. 2001;D'Alessio et al. 2006;Dodson-Robinson et al. 2009). In our models, there is so little turbulent energy generated in the dead zone that the disk midplane falls to 20 K-the ambient temperature of the remnant molecular cloud surrounding the disk. Here we assume the disk is optically thin to long-wavelength radiation from the ambient cloud and cannot cool below the ambient temperature. Only in the inner ∼ 2 AU does residual mass transfer by large-scale magnetic fields lift the midplane temperature above the minimum value.
Moving outward through the disk, the midplane temperature rises modestly until it equalizes with the surface temperature where the disk becomes optically thin to stellar radiation. Both of our models, despite their modest masses (0.015M and 0.03M ), are optically thick to stellar radiation inside ∼ 70 AU, decoupling the surface and midplane temperatures. Passively heated disk models (e.g. Woitke et al. 2009) are therefore poor approximations to the midplane temperatures of our model disks, as are constant-α models in which T decreases with R at the midplane. Instead, coupling the simple viscosity prescription in Section 4.3 with a radiative transfer scheme is preferable for modeling observables that trace the inner ∼ 70 AU. The decrease in surface viscosity with time is due to the increasing density in the surface layers as the disk cools. Despite the fact that the disk surface cools with time, the tendency of the MRI to pile up mass unevenly leads the temperature to increase with time in certain parts of the disk. The bottom panels of Figure 7 show the temperature increase in the pileup at 1 AU in Model 2. Though modest, the temperature increase does affect the location of the ice line, which we discuss in the next section.

The Ice Line in MRI-Active Disks
Ice forms in protoplanetary disks when the temperature falls below 145-170 K, depending on the water vapor's partial pressure (Podolak & Zucker 2004;Lecar et al. 2006). Observations of the outer asteroid belt place the ice line in today's Solar System at 2.7 AU. Theoretical estimates of the ice line location in the solar nebula place it a minimum of 0.6 AU from the sun (Davis 2005;Garaud & Lin 2007) and a maximum of 6 AU at the beginning of the T-Tauri phase, moving inward as the disk evolves (Dodson-Robinson et al. 2009). In constant-α disks, the ice line moves inward with time as the optically thick disk radiates away its accretion energy. At late times, however, the inner disk loses enough mass to become optically thin to stellar irradiation, causing the ice line to move outward with time (Garaud & Lin 2007;Oka et al. 2011). The more massive a constant-α disk, the higher its midplane temperature will be. Taking the Davis (2005), Garaud & Lin (2007) and Oka et al. (2011) disk models up to a reasonable planet-forming mass of 0.04M or higher (Thommes et al. 2008) would push their ice lines outside the terrestrial planet-forming region, more in agreement with the results of Dodson-Robinson et al. (2009).
In an MRI-active disk, the overall thermal evolution is determined by the dimming of the parent star, which causes viscosity in the active surface layers to decrease with time (see Section 5.1). Figure 8 shows the 2-D structure of the ice line in Model 2 at timepoints 100 yr, 1000 yr, 10,000 yr, 100,000 yr, 1 Myr and 3 Myr. At early times, we recover the "two-branch" structure of the ice line seen by Davis (2005), in which a nearly horizontal ice line divides the hot surface from the cool interior and a midplane ice line separates the warm midplane near the star from the cool midplane far away from the star. After 10,000 yr of evolution, the inner edge of the disk at 0.5 AU cools enough for ice to freeze at the midplane, causing a "pinch-off" in the midplane branch of the ice line that leaves an H 2 O gas bubble at 0.6 AU. This pinchoff is probably a boundary effect: since mass flows freely from our inner boundary at 0.5 AU onto the star, the disk near our inner boundary loses mass and cools quickly.
By 1 Myr, mass loss from the inner edge of our grid combined with less activity in the surface layers cools the disk enough to push the midplane ice line inside 0.5 AU, the inner boundary of our computation. Here our results are still consistent with the low-mass models of Davis (2005), Garaud & Lin (2007) and Oka et al. (2011), who all show the ice line moving inside 1 AU for accretion ratesṀ ∼ 10 −10 M yr −1 . An important difference between our models and those of Davis (2005), Garaud & Lin (2007) and Oka et al. (2011) is that the accretion rate does not have to drop extremely low for the terrestrial planet-forming region to cool enough to freeze ice: both Model 1 and Model 2 have active layers that drivė M ∼ 10 −9 M yr −1 through the disk surface. Our model disks predict colder midplanes in the inner 10 AU than those of Terquem (2008) due to the lower minimum α assumed here (α min = 10 −5 for our work vs. α D = 10 −3 or α D = 10 −4 in Figure 3 of Terquem (2008), where α D the value of α in the dead zone). Unlike Terquem (2008), our disk has a colder midplane than surface because we include the effects of stellar heating.
The MRI can also create local regions that heat up with time not because the disk becomes optically thin, but because mass piles up in the dead zone. Between 1 Myr and 3 Myr of evolution, the pileup at 1 AU of Model 2 heats up enough for the midplane ice line to reappear-a real physical effect. Model 1, which does not show any such pileup, follows a similar evolutionary path as Model 2 up to 1 Myr. In an MRI-active disk, water ice in the terrestrial planet-forming region may be transient. Understanding which parts of the inner disk have water ice available for planet formation requires a careful comparison of the planetesimal growth timescale, the star cooling timescale and the growth timescale of any pileups deposited in the dead zone. In Model 2, the pileup grows and the star dims on a timescale similar to the disk lifetime. We will examine high-mass MRI-active disks, in which dead-zone pileups can grow much more quickly, in a forthcoming paper.

MODEL LIMITATIONS
Although our model incorporates much of the physics of disk evolution during the T-Tauri phases of stars, our long-timescale computation requires a number of simplifications. Our disk accretion model suffers from the following limitations, which may affect our conclusions: 1. We model our disks as 1+1d instead of 3d. We assume radial symmetry, and that the vertical structure is not coupled to the radial mass transport. This ignores the possibility of accretion driven by non-axisymmetric instabilities not represented by our parameterization of stress from MHD turbulence, such as the Rossby wave instability (Hawley 1987;Li et al. 2001;Meheut et al. 2012). We may then underestimate accretion rates in some regions of the disk.
2. We do not model the effects of gravitational instability on momentum transport, as is done by Martin et al. (2012b), but instead limit our simulations to those disks which are gravitationally stable for their entire lifetime. To test whether a disk is gravitationally stable to axisymmetric perturbations, we calculate the Toomre Q parameter at every radial grid point. Stability requires that Note that for disks more massive than those we have simulated, the disk does eventually become gravitationally unstable to axisymmetric perturbations within the dead zone (see Section 4.2).
3. We do not model the effects of Hall diffusivity, which is the least well understood magnetic diffusion regime. However, it is unlikely to alter either the conditions required for MRI or the strength of the turbulence where Ohmic dissipation is also present (Sano & Stone 2002).
4. We make several assumptions about the magnetic field strength and presence of grains. The activity of MRI driven turbulence is very sensitive to these parameters. First, we deviate from the standard interstellar gas/small grain mass ratio of 100 and instead assume a ratio of 1000. This is equivalent to assuming that 90% of the grain mass is either in grains larger than 1µm (Oliveira et al. 2010;Scaife et al. 2012), which do not significantly affect electron density, or has settled below the dead zone. Mohanty et al. (2013) find that grain depletion through growth or settling is required to account for the observed accretion rates of low mass protostars. Without the high gas/small grain mass ratio, the inner Ohmic dead zone and the upper Ambipolar dead zone overlap, producing a passive thermal structure for much of the radial extent of the disk (see Section 4.1).
We also assume a vertically constant magnetic pressure with a midplane plasma β of 1000. This is at the upper end of the range for midplane plasma β found by Fromang & Nelson (2006) in global MHD simulations with saturated turbulence. Without a sufficiently weak magnetic field, the ambipolar diffusivity can be large enough that there is once again an overlap in the dead zones and an at least partially passive disk.
5. We assume a zero stress boundary condition at the inner radius of our disks (see Section 3.3). This creates a non-realistic boundary effect in which mass is rapidly depleted from the inner annuli. One result of this can be seen in Figure 8, as the midplane within 0.6 AU becomes cold due to the loss of surface density.
Due to computational convergence difficulties in our model, at any particular timestep there are a few annuli with vertical structures which do not meet one or more of our midplane boundary conditions (Equations 50, 51, and 52), due to the discontinuous nature of the MRI and opacity. We do not consider this to be a limitation on our simulations, as the unsolved annuli's contribution to the viscosity profile is smoothed with a Gaussian filter before being used to update the surface density profile (Equation 46). Although there are always a few "bad" annuli present, at any particular radius the lack of convergence for the vertical structure persists for only a few timesteps. The unconverged annuli can be seen as occasional incongruous vertical bars in our contour plots (Figures 2, 3, 5, 7, and 8).

CONCLUSIONS
In the Introduction, we asked five questions about the structure and evolution of MRI-active disks. Here we summarize our findings and answer each question: 1. How do the relative sizes of the dead zone and active layers change over time?
The radial size of the dead zone is almost constant in time, while the vertical height of the dead zone shrinks over time. What was surprising about our results was not the evolution of the dead zone, but the complexity of the disk structure. Between 1.2 AU and 1.7 AU, Model 2 has a five-layer structure throughout most of its evolution: inactive corona, active layer, dead slice, active layer, dead midplane (see Section 4.1 and Figures 1 and 3). Note, however, that detailed 3-D simulations would likely show no dead slice since the MRI wavelength in the active layers bracketing the dead slice is of order the dead slice thickness. The lower-mass Model 1 has at most three layers: inactive corona, active layer, dead midplane. Throughout this work, we have seen that increasing disk mass leads to increasing complexity in the disk structure and accretion flow. Finally, the radial size of the dead zone was somewhat higher than predicted in previous work: 16 AU for Model 1 and 21 AU for Model 2. Previous papers reporting a ∼ 5 AU dead zone used the MMSN (e.g. Matsumura & Pudritz 2003;Salmeron & Wardle 2008;Flaig et al. 2012), but the modestly higher surface densities of our disks expanded the dead zone. There is some evidence that the solar nebula had a large dead zone consistent with our findings-the giant planets have an atmospheric composition gradient that, if primordial, would have diffused on million-year timescales if not protected by a dead zone (Nelson & Gressel 2010). One caveat, though, is that we have assumed that the stellar X-ray flux is constant in time, as is the cosmic ray flux. A decreasing stellar X-ray flux, which ionizes mainly the inner ∼ 3 AU of the disk surface, might erase the dead slice over time, while changing the cosmic ray flux as the ambient molecular cloud disperses would certainly change the radial extent of the dead zone.

How doesṀ vary with radius and time?
Throughout the disk evolution, |Ṁ | is highest at the inner and outer boundaries of the dead zone. The lowest |Ṁ | in the inner disk, where gas flows toward the star, is in the middle of the dead zone. |Ṁ | is about 50% higher for Model 2 than Model 1, suggesting that higher-mass disks support higher MRI-driven accretion rates-though the increase inṀ with M disk is modest. |Ṁ (t)| decreases extremely slowly: though the depth of the active layer does not change with time ( Figure 6), the viscosity in the active layer drops modestly as the star cools. As required by conservation of angular momentum, both model disks expand as they evolve, creating a turn-around in the accretion flow. The turn-around radius is almost entirely determined by R out at t = 0 and moves inward as the disk evolves. Here, with R out = 70 AU at t = 0, the turn-around radius eventually reaches 40 AU after 1 Myr of evolution. Note that the turn-around radius moves steadily inward in both models and does not converge toward a particular location (Figure 4). One expects the turn-around radius to move steadily inward because mass must continually join the outward "decretion" flow in order to keep transporting angular momentum outward. Our models predictṀ ∼ 10 −9 M yr −1 in the planet-forming region of the disk, but the exact value ofṀ depends on many free parameters such as grain size, gas/small grain mass ratio and magnetic field strength. We have not attempted an exhaustive parameter study ofṀ as a function of all variables. We merely note that for a gas/small grain mass ratio of 100 and a grain size of 0.1µm, all MRI activity in the disk was suppressed. Likewise, decreasing plasma β at the midplane from 1000 to 100 suppressed MRI turbulence, though not as severely as small grains.
3. How can disk modelers parameterize heating in the active layers and dead zone without resorting to a 3-D MHD simulation?
The value of non-evolving, "snapshot" disk models is inarguable, particularly for modeling observables. Section 4.3 presents a simple modification of the standard, constant-α irradiated disk that approximates the thermal structure of the disk where the dead zone is present. Simply set the surface density of the active layer to 10 g cm −2 and use α ≈ 0.01 (see Figures 1, 3 and 6). For the dead zone, use α ≈ 10 −5 . For numerical models, we recommend smoothing the transition between the dead zone and the active layer. While the depth of the active layer does vary across the dead zone, its variation is at most a factor of two in a given disk. The fact that the active layer depth is almost constant in time for both Model 1 and Model 2 makes our simple prescription applicable to any stage of T-Tauri disk evolution, provided that stellar irradiation is included.
4. Does the disk midplane heat up or cool off with time?
In both models, the midplane temperature varies little with time, but the disk surface cools as the star evolves down the Hayashi track. Despite the modest masses of our model disks, the optical depth of both to stellar irradiation is enough to thermally decouple the surface from the midplane. Most of the dead zone is so lacking in energy generation that it falls to the assumed ambient temperature of surrounding molecular cloud material, 20 K in these models. Since the radial extent of the dead zone changes little with time, the disk midplane temperature remains static except at R 60 AU, where the disk thins enough over time to become optically thin to stellar irradiation (Figure 7).
In the midplane, there are two possible locations where the temperature is not static but increases with time ( Figure 7). The first is at the outer edge of the dead zone. A slight decrease in surface density with time pushes the dead zone boundary modestly inward, allowing material at the edge of the dead zone to become turbulent and heat up. The other location of increasing temperature with time is the pileup at 1 AU of Model 2. Lowermass Model 1 does not develop any such pileups on million-year timescales. To the extent that MRI-driven accretion deposits piles of material in the dead zone, the disk midplane may heat up. Note, however, that the thermal properties of the disk depend on the ionizing radiation it receives. Near the star, the disk structure would evolve if the X-ray flux were to change with time. In future work, it would be interesting to let L X scale with bolometric luminosity on the Hayashi track. A spike in cosmic ray ionization from a nearby supernova would affect the global disk structure (Fatuzzo et al. 2006).

5.
Where is the ice line in an MRI-active disk and how does its location change over time?
Due to the paucity of energy generation in the dead zone, the midplane ice line falls somewhere inside the terrestrial planet-forming region. In our models, the midplane ice line actually moves off the inner edge of our grid at 0.5 AU after 10 5 years of disk evolution, reappearing in Model 2 after the pileup at 1 AU reheats the midplane. Despite the different physics used in computing the ice line location, our results roughly agree with those of Davis (2005), Garaud & Lin (2007) and Oka et al. (2011) in that the midplane ice line is inside 1 AU for most of the disk's evolution. Determining whether and when ice is available for terrestrial planet formation requires carefully comparing the disk cooling timescale, planetesimal growth timescale and the timescale on which MRI-deposited pileups grow. An important difference between our model and standard, constant-α disk models is that the accretion rate does not have to drop extremely low to move the ice line inside 1 AU: throughout their evolution, our disks haveṀ 10 −9 M yr −1 moving through the active layers that sandwich the icy inner disk.
Here we have presented the first analysis of the structure and evolution of an entire MRI-active disk on million-year timescales. While we have chosen to focus on low-mass disks in this work in order to compare with previous studies, we will expand our analysis to include high-mass, planet-forming disks in a forthcoming paper. Already, at M disk = 0.03M -just short of the minimum 0.04M required for giant planet formation (Thommes et al. 2008)-Model 2 exhibits some new features such as the split active layer and the re-heating midplane.
Funding for this work was provided by NASA through grant NNX10AH28G to S.D.R. and N.J.T., and by University of Texas through a startup grant to S.D.R. Computing and visualization support were provided by the Texas Advanced Computing Center, which is funded by the National Science Foundation. N.J.T. carried out his work at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA.