Free Access
Issue
A&A
Volume 564, April 2014
Article Number A22
Number of page(s) 15
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201321911
Published online 28 March 2014

© ESO, 2014

1. Introduction

Current models of protoplanetary (PP) disks are predicated on the idea that significant regions of the disk are too poorly ionized to sustain magneto-rotational-instability (MRI) turbulence. These PP disks are thought to comprise a turbulent body of plasma (the “active zone”) enveloping a region of cold quiescent gas, in which accretion is actually absent (the “dead zone”) (Gammie 1996; Armitage 2011). These models posit a critical inner radius (~1 au) within which the disk is fully turbulent and beyond which the disk exhibits turbulence only in its surface layers, for a range of radii (110 au) (but see Bai & Stone 2013, for a complication of this picture).

The inner boundary between the MRI-active and dead regions is crucial for several key processes. Because there is a mismatch in accretion across the boundary, a pressure maximum will naturally form at this location, which (a) may halt the inward spiral of centimetre to metre-sized planetesimals (Kretke et al. 2009); and (b) may excite a large-scale vortex instability (“Rossby wave instability”) (Lovelace et al. 1999) that may promote dust accumulation, hence planet formation (Barge & Sommeria 1995; Lyra et al. 2009; Meheut et al. 2012). On the other hand, the interface will influence the radial profiles of the dead zone’s thermodynamic variables – temperature and entropy, most of all. Not only will it affect the global disk structure and key disk features (such as the ice line), but the interface will also control the preconditions for dead-zone instabilities that feed on the disk’s small adverse entropy gradient, such as the subcritical baroclinic instability and double-diffusive instability (Lesur & Papaloizou 2010; Latter et al. 2010).

Most studies of the interface have been limited to isothermality. This is a problematic assumption because of the pervasive interpenetration of dynamics and thermodynamics in this region, especially at the midplane. Temperature depends on the turbulence via the dissipation of its kinetic and magnetic fluctuations, but the MRI turbulence, in turn, depends on the temperature through the ionization fraction, which is determined by thermal ionization (Pneuman & Mitchell 1965; Umebayashi & Nakano 1988). Because of this feedback loop, the temperature is not an additional piece of physics that we add to simply complete the picture; it is instead at the heart of how the interface and its surrounding region work. One immediate consequence of this feedback is that much of the midplane gas inward of 1 au is bistable: if the gas at a certain radius begins as cold and poorly ionized, it will remain so; conversely, if it begins hot and turbulent, it can sustain this state via its own waste heat (Latter & Balbus 2012). Thus two stable states are available at any given radius in the bistable region. This complicates the question of where the actual location of the dead zone boundary lies. It also raises the possibility that the boundary is not static, and may not even be well defined. Similar models have also been explored in the context of FU Ori outbursts Zhu et al. (2010, 2009a,b).

In this paper we simulate these dynamics directly with a set of numerical experiments of MRI turbulence in PP disks. We concentrate on the inner radii of these disks (~0.11 au) so that our simulation domain straddles both the bistable region and the inner dead-zone boundary. Our aim is to understand the interactions between real MRI turbulence and thermodynamics, and thus test and extend previous exploratory work that modelled the former via a crude diffusive process (Latter & Balbus 2012). Our 3D magnetohydrodynamics (MHD) simulations employ the Godunov code RAMSES (Teyssier 2002; Fromang et al. 2006), in which the thermal energy equation has been accounted for and magnetic diffusivity appears as a function of temperature (according to an approximation of Saha’s law). We focus on the optically thick disk midplane, and thus omit non-thermal sources of ionization. This also permits us to model the disk in the cylindrical approximation.

This initial numerical study is the first in a series, and thus lays out our numerical tools, tests, and code-checks. We also verify previous published results for fully turbulent disks and disks with static dead zones (Dzyurkevich et al. 2010; Lyra & Mac Low 2012). We explore the global radial profiles of thermodynamic variables in such disks, and the nature of the turbulent temperature fluctuations, including a first estimate for the magnitude of the turbulent thermal flux. We find that the dead zone can be effectively heated by density waves generated at the dead zone boundary, and we estimate the resulting wave damping and heating. Our simulations indicate that the dead zone may be significantly hotter than most global structure models indicate because of this effect. Another interesting result concerns the dynamics of the dead zone interface itself. We verify that the interface is not static and can migrate from smaller to larger radii. All such simulated MRI “fronts” ultimately stall at a fixed radius set by the thermodynamic and radiative profiles of the disk, in agreement with Latter & Balbus (2012). These fronts move more quickly than predicted because they propagate not via the slower MRI turbulent motions but by the faster density waves. Finally, we discuss instability near the dead zone boundary.

The structure of the paper is as follows. In Sect. 2 we describe the setup we used for the MHD simulations and give our prescriptions for the radiative cooling and thermal ionization. We test our implementation of thermodynamical processes in models of fully turbulent disks in Sect. 3. There we also quantify the turbulent temperature fluctuations and the turbulent heat flux. Section 4 presents results of resistive MHD simulations with dead zones. Here we look at the cases of a static and dynamic dead zone separately. Rossby wave and other instabilities are briefly discussed in Sect. 5; we subsequently summarize our results and draw conclusions in Sect. 6.

2. Setup

We present in this paper a set of numerical simulations performed using a version of the code RAMSES (Teyssier 2002; Fromang et al. 2006) which solves the MHD equations on a 3D cylindrical and uniform grid. Since we are interested in the inner-edge behaviour at the mid-plane, we discount vertical stratification and work under the cylindrical approximation (Armitage 1998; Hawley 2001; Steinacker & Papaloizou 2002). In the rest of this section we present the governing equations, prescriptions, main parameters, and initial and boundary conditions. In order to keep the discussion as general as possible, all variables and equations in this section are dimensionless.

2.1. Equations

Since we are interested in the interplay between dynamical and thermodynamical effects at the dead zone inner edge, we solve the MHD equations, with (molecular) Ohmic diffusion, alongside an energy equation. We neglect the kinematic viscosity, that is much smaller than Ohmic diffusion, because of its minor role in the MRI dynamics. However dissipation of kinetic energy is fully captured by the numerical grid. We adopt a cylindrical coordinate system (R,φ,Z) centred on the central star:

where ρ is the density, v is the velocity, B is the magnetic field, and P is the pressure. Φ is the gravitational potential. In the cylindrical approximation, it is given by Φ =  −GM/R where G is the gravitational constant and M is the stellar mass. In RAMSES, E is the total energy, i.e. the sum of kinetic, magnetic and internal energy eth (but it does not include the gravitational energy). Since we use a perfect gas equation of state to close the former set of equations, the latter is related to the pressure through the relation eth = P/(γ − 1) in which γ = 1.4. The magnetic diffusivity is denoted by η. As we only consider thermal ionization, η will depend on temperature T; its functional form we discuss in the following section. Associated with that resistivity is a resistive flux η that appears in the divergence of Eq. (3). Its expression is given by Eq. (23) of Balbus & Hawley (1998). The symbol denotes radiative losses and it will be described in Sect. 2.3. Dissipative and radiative cooling terms are computed using an explicit scheme. This method is valid if the radiative cooling time scale is much longer than the typical dynamical time (see Sect. 2.4). Finally, we have added a source term in the continuity equation that maintains the initial radial density profile ρ0 (Nelson & Gressel 2010; Baruteau et al. 2011). It is such that (5)The restoring time scale τρ is set to 10 local orbits and prevents the long term depletion of mass caused by the turbulent transport through the inner radius.

2.2. Magnetic diffusion

In this paper, we conduct both ideal (η = 0) and non-ideal (η ≠ 0) simulations. In the latter case, we treat the disk as partially ionized. In such a situation, the magnetic diffusivity is known to be a function of the temperature T and ionization fraction xe (Blaes & Balbus 1994): (6)The ionisation fraction can be evaluated by considering the ionization sources of the gas. In the inner parts of PP disks, the midplane electron fraction largely results from thermal ionization (neglecting radioisotopes that are not sufficient by themselves to instigate MRI). Non thermal contributions due to ionization by cosmic (Sano et al. 2000), UVs (Perez-Becker & Chiang 2011) and X-rays (Igea & Glassgold 1999) are negligible because of the large optical thickness of the disk at these radii. In fact, the ionization is controlled by the low ionization potential (Ei ~ 1 − 5 eV) alkali metals, like sodium and potassium, the abundance of which is of order 10-7 smaller than hydrogen (Pneuman & Mitchell 1965; Umebayashi & Nakano 1988). Then xe can be evaluated using the Saha equation. When incorporated into Eq. (6), one obtains the following relation between the resistivity and the temperature: (7)with T a constant. Because of the exponential, η varies very rapidly with T, and, in the context of MRI activation, can be thought of as a “switch” around a threshold TMRI. In actual protoplanetary disks, it is well known that that temperature is around 103 K (Balbus & Hawley 2000; Balbus 2011). Taking this into account, and for the sake of simplicity, we use in our simulations the approximation: (8)where TMRI is the activation temperature for the onset of MRI. Note that by taking a step function form for η we must neglect its derivative in Eq. (4).

2.3. Radiative cooling and turbulent heating

As we are working within a cylindrical model of a disk, without explicit surface layers, we model radiative losses using the following cooling function: (9)This expression combines a crude description of radiative cooling of the disk (ρσT4) as well as irradiation by the central star (), in which σ should be thought of as a measure of the disk opacity (see also Sect. 2.4). In the absence of any other heating or cooling sources, T should tend toward Tmin, the equilibrium temperature of a passively irradiated disk. In principle, Tmin is a function of position and surface density, but we take it to be uniform for simplicity. In practice, we found it has little influence on our results. Turbulence or waves ensure that the simulated temperature is always significantly larger than Tmin. The parameter σ determines the quantity of thermal energy the gas is able to hold. For computational simplicity and to ease the physical interpretation of our results, we take σ to be uniform. In reality, σ should vary significantly across the dust sublimation threshold (~1500 K) (Bell & Lin 1994), but an MRI front will always be much cooler (~1000 K). As a consequence, the strong variation of σ will not play an important role in the dynamics. Finally, we omit the radiative diffusion of heat in the planar direction, again for simplicity though in real disks it is an important (but complicated) ingredient.

Our prescription for magnetic diffusion, Eq. (8), takes care of the dissipation of magnetic energy when T < TMRI. However, if the disk is sufficiently hot the dissipation of magnetic energy (as with the kinetic energy) is not explicitly calculated. In this case energy is dissipated on the grid; our use of a total energy equation ensures that this energy is not lost but instead fully transferred to thermal energy. It should be noted that though this is perfectly adequate on long length and time scales, the detailed flow structure on the dissipation scale will deviate significantly from reality.

2.4. Initial conditions and main parameters

Our 3D simulations are undertaken on a uniformly spaced grid in cylindrical coordinates (R,φ,Z). The grid ranges over R ∈  [R0,8   R0], φ ∈  [0/4] and Z ∈  [ − 0.3   R0,0.3   R0]. Here the key length scale R0 serves as the inner radius of our disk domain and could be associated with a value between 0.20.4 au in physical units. Each run has a resolution [320,80,80], which has been shown in these conditions to be sufficient for MHD turbulence to be sustained over long time scales (Baruteau et al. 2011). For completeness, we present a rapid resolution study in Appendix D that further strengthens our results.

Throughout this paper, we denote by X0 the value of the quantity X at the inner edge of the domain, i.e. at R = R0. Units are chosen such that: where Ω stands for the gas angular velocity at radius R. Thus all times are measured in inner orbits, so a frequency of unity correspond to a period of one inner orbit. We will refer to the local orbital time at R using the variable τorb.

The initial magnetic field is a purely toroidal field whose profile is built to exhibit no net flux and whose maximum strength corresponds to β = 25: (10)The simulations start with a disk initially in approximate radial force balance (we neglect the small radial component of the Lorentz force in deriving that initial state). Random velocities are added to help trigger the MRI. Density and temperature profiles are initialized with radial power laws: (11)where p and q are free parameters. Pressure and temperature are related by the ideal gas equation of state: (12)where P0 depends on the model and is specified below. We choose p =  −1.5, which sets the initial radial profile for the density. During the simulation, it is expected to evolve on the long secular time scale associated with the large-scale accretion flow (except possibly when smaller scale feature like pressure bump appear, see Sect. 4). Quite differently, the temperature profile evolves on the shorter thermal time scale (which is itself longer than the dynamical time scale τorb). It is set by the relaxation time to thermal equilibrium1.

As shown by Shakura & Sunyaev (1973), turbulent heating in an accretion disk is related to the component of the turbulent stress tensor T: (13)For MHD turbulence T amounts to (Balbus & Papaloizou 1999): (14)The ⟨ . ⟩ notation stands for an azimuthal, vertical and time average, δvφ is the flow’s azimuthal deviation from Keplerian rotation. The last scaling in Eq. (13)assumes such deviations are small, i.e. the disk is near Keplerian rotation. Balancing that heating rate by the cooling function ℒ ~ ρσT4 gives a relation for the disk temperature (here we neglect Tmin in the cooling function): (15)Using the standard α prescription for illustrative purposes, we write T = α ⟨ P ⟩ which defines the Shakura-Sunyaev α parameter, a constant in the classical α disk theory (Shakura & Sunyaev 1973). Since Ω ∝ R-1.5, Eq. (15) suggests that T ∝ R-0.5 in equilibrium for uniform σ values. We thus used q =  −0.5 in our initialization of T.

We can use a similar reasoning to derive an estimate for the thermal time-scale. Using the turbulent heating rate derived by Balbus & Papaloizou (1999) and the definition of α, the mean internal energy evolution equation can be approximated by (16)which yields immediately (17)For α = 0.01, the heating time scale is thus about 25 local orbits. The cooling time-scale τcool is more difficult to estimate, but near equilibrium should be comparable to τheat.

In order to obtain a constraint on the parameter σ we use the relations , where cs is the sound speed and H is the disk scale height. At R = R0 approximate thermal equilibrium, Eq. (15), can be used to express σ as a function of disk parameters: (18)Now setting α ~ 10-2 (and moving to numerical units) ensures that σ only depends on the ratio H0/R0 = c0/(R0Ω0). We investigate two cases: a “hot” disk with H0/R0 = 0.1 which yields σ = 1.1 × 10-4, and a “cold” disk with H0/R0 = 0.05 and thus σ = 2.7 × 10-5. The two parameter choices for σ will be referred to as the σhot and σcold cases in the following. Using the relation between pressure and sound speed respectively yields P0 = 7.1 × 10-3 and P0 = 1.8 × 10-3 for the two models. The surface density profile in the simulations is given by a power law: (19)In resistive simulations, the parameters Tmin, TMRI and η0 (when applicable) need to be specified for the runs to be completely defined. We used Tmin = 0.05   T0, which means that the temperature in the hot turbulent innermost disk radius is 20 times that of a cold passive irradiated disk. As discussed in Sect. 2.3, Tmin is so small that it has little effect on our results. The values of TMRI and η0 vary from model to model and will be discussed in the appropriate sections.

In order to assess how realistic our disk model is, we convert some of its key variables to physical units. This can be done as follows. We consider a protoplanetary disk orbiting around a half solar mass star, with surface density Σ ~ 103 g cm-2 and temperature T ~ 1500 K at 0.1 au. (Given the radial profile we choose for the surface density, this would correspond to a disk mass of about seven times the minimum mass solar nebular within 100 au of the central star.) For this set of parameters, we can use Eq. (18) to calculate σ = 1 × 10-9 in cgs units for model σcold, using α ~ 10-2. In the simulations, the vertically integrated cooling rate is thus given by the cooling rate per unit surface: (20)This value can be compared to the cooling rate of a typical α disk with the same parameters (Chambers 2009): (21)where σb is the Stefan-Boltzman constant and τ = κ0Σ/2 in which κ0 = 1 cm2 g-1 stands for the opacity. Using these figures, we obtain , i.e. a value that is close (given the level of approximation involved) to that used in the simulations2 (see Eq. (20)). We caution that this acceptable agreement should not be mistaken as a proof that we are correctly modelling all aspects of the disk’s radiative physics. The cooling function we use is too simple to give anything more than an idealized thermodynamical model. It only demonstrates the consistency between the thermal time scale we introduce in the simulations and the expected cooling time scale in proto-planetary disks.

2.5. Buffers and boundaries

Boundary conditions are periodic in Z and φ while special care has been paid to the radial boundaries. Here vR, vz, Bφ, Bz are set to zero, BR is computed to enforce magnetic flux conservation, vφ is set to the Keplerian value, and finally temperature and density are fixed to their initial values.

As is common in simulations of this kind, we create two buffer zones adjacent to the inner and outer boundaries in which the velocities are damped toward their boundary values in order to avoid sharp discontinuities. The buffer zones extend from R = 1 to R = 1.5 for the inner buffer and from R = 7.5 to R = 8 for the outer one. For the same reason, a large resistivity is used in those buffer zones in order to prevent the magnetic field from accumulating next to the boundary. As a result of this entire procedure, turbulent activity decreases as one approaches the buffer zones. This would occasionally mean the complete absence of turbulence in the region close to the inner edge because of inadequate resolution. To avoid that problem, the cooling parameter σ gradually increases with radius in the region R0 < R < 2R0. In Appendix A, we give for completeness the functional form of σ we used. These parts of the domain that corresponds to the inner and outer buffers are hatched on all plots of this paper.

3. MHD simulations of fully turbulent PP disks

We start by describing the results we obtained in the “ideal” MHD limit: η = 0 throughout the disk. As a consequence, there is no dead zone and the disk becomes fully turbulent as a result of the MRI. The purpose of this section is (a) to describe the thermal structure of the quasi-steady state that is obtained; (b) to check the predictions of simple alpha models; and (c) to examine the small-scale and short-time thermodynamic fluctuations of the gas, especially with respect to their role in turbulent heat diffusion. Once these issues are understood we can turn with confidence to the more complex models that exhibit dead-zones.

3.1. Long-time temperature profiles

We discuss here the results of the two simulations that correspond to σhot and σcold. We evolve the simulations not only for a long enough time for the turbulence transport properties to reach a quasi-steady state but also for the thermodynamic properties to have also relaxed. Thus the simulations are evolved for a time much longer than the thermal time of the gas, τheat ≈ 25 local orbits. Here, we average over nearly 1000 inner orbits (about 2τheat at the outermost radius R = 7).

In Fig. 1 we present the computed radial temperature profiles and the radial profiles of α for the two simulations. Note that α is clearly not constant in space. In both models, α is of order a few times 10-2 and decreases outward. We caution here that this number and radial profiles are to be taken with care. The α value is well known to be affected by numerical convergence (Beckwith et al. 2011; Sorathia et al. 2012) as well as physical convergence issues (Fromang et al. 2007; Lesur & Longaretti 2007; Simon & Hawley 2009). This may also influence the temperature because it depends on the turbulent activity.

thumbnail Fig. 1

Temperature profiles ⟨ T ⟩ averaged over 900 inner orbits (red curves). Black plain lines show their corresponding theoretical profiles. Black dashed lines show their corresponding initial profiles. Top panel: σcold case. Bottom panel: σhot case. The subframes inserted in the upper right of each panel show the alpha profiles obtained in both cases (red lines). The dashed lines show the analytic profiles used to compute the theoretical temperature profile (see text), given by α = 6.3 × 10-2 − 7.1 × 10-3R/R0 and α = 4.6 × 10-2 − 5.7 × 10-3R/R0 respectively.

Open with DEXTER

The disk temperature T rapidly departs from the initialized power law given by Eq. (11). In both plots the averaged simulated temperature decreases faster than R-0.5. This is partly a result of α decreasing with radius, and a consequent reduced heating with radius. The disk aspect ratios are close to the targeted values: H/R ~ 0.15 and 0.08 in model σhot and σcold respectively, slightly increasing with radius in both cases.

One of our goals is to check the implementation of the thermodynamics, and to ensure that we have reached thermal equilibrium. To accomplish this we compare the simulation temperature profiles with theoretical temperature profiles computed according to the results of Balbus & Papaloizou (1999). The theoretical temperature profile can be deduced from thermal equilibrium, Eq. (15), in which we now include the full expression for the cooling function (i.e. including Tmin): (22)Using the simulated α profiles as inputs, we could then calculate radial profiles for T. In fact, we used linear fits of α (shown as a dashed line) in Eq. (22), for simplicity. The theoretical curves are compared with the simulation results in Fig. 1. The overall good agreement validates our implementation of the source term in the energy equation and also demonstrates that we accurately capture the turbulent heating. It is also a numerical confirmation that turbulent energy is locally dissipated into heat in MRI-driven turbulence and that our simulations have been run sufficiently long to achieve thermal equilibrium. During the steady state phase of the simulations, the mass loss at the radial boundaries is very small: the restoring rate is on average on the domain.

3.2. Turbulent fluctuations of temperature

We now focus on the local and short time evolution of T, by investigating the fluctuations of density and temperature once the mean profile has reached a quasi steady state. We define the temperature fluctuations by , where an overline indicates an azimuthal and vertical average. The magnitude of the fluctuations are plotted in Fig. 2. In the σcold case, the fluctuations range between 4 and 8% while they vary between 3 and 5% in the σhot case. The smaller temperature fluctuations in the latter case possibly reflect the weaker turbulent transport in the latter case, while also its greater heat capacity. The tendency of the relative temperature fluctuations to decrease with radius is due to α decreasing outward.

thumbnail Fig. 2

Temperature standard deviation profiles averaged over 900 inner orbits. The plain line show the result from the σhot case and the dashed curve show the result from the σcold run.

Open with DEXTER

thumbnail Fig. 3

Blue dots map the scaling relation between temperature and density fluctuations δT and δρ in the σhot simulation. 900 events are randomly selected in the whole domain (except buffer zones) among 900 inner orbits. The red line shows the linear function (slope γ − 1)) corresponding to the adiabatic scaling.

Open with DEXTER

These temperature fluctuations can be due to different kinds of events that act, simultaneously or not, to suddenly heat or cool the gas. Such events can, for example, be associated with adiabatic compression or magnetic reconnection. To disentangle these different possibilities, we note that turbulent compressions, being primarily adiabatic, satisfy the following relationship between δT and the density fluctuations δρ: (23)In Fig. 3 we plot the distribution of a series of fluctuation events randomly selected during the σhot model in the plane. The dot opacity is a function of the event’s radius: darker points stand for outer parts of the disk. The events are scattered around the adiabatic slope (represented with a red line), but display a large dispersion that decreases as radius increases. The difference between the disk cooling time and its orbital period is so large that the correlation should be much better if all the heating/cooling events were due to adiabatic compression. This implies that temperature fluctuations are mainly the result of isolated heating events such as magnetic reconnection. This suggestion is supported by the larger dispersion at short distance from the star, where turbulent activity is largest. However, a detailed study of the heating mechanisms induced by MRI turbulence is beyond the scope of this paper. Indeed, the low resolution and the uncontrolled dissipation at small scales – important for magnetic reconnection – both render any detailed analysis of the problem difficult. Future local simulations at high resolution and with explicit dissipation will be performed to investigate that question further.

thumbnail Fig. 4

Thermal diffusivity’s radial profile averaged over 900 inner orbits. Dashed lines show the mean thermal diffusivity for two homogeneous αT (αT = 0.02 and αT = 0.004). The grey area delimit one standard deviation around the mean thermal diffusivity.

Open with DEXTER

As shown in Balbus (2004), MHD turbulence may also induce thermal energy transport through correlations between the temperature fluctuations δT and velocity fluctuations δv. For example, the radial flux of thermal energy is given by Eq. (5) of Balbus (2004): (24)By analogy with molecular thermal diffusivity, we quantify the turbulent efficiency to diffuse heat using a parameter κT: (25)To make connection with standard α disk models, we introduce the parameter αT which is a dimensionless measure of κT: (26)as well as the turbulent Prandtl number PR = α/αT that compares turbulent thermal and angular momentum transports. We measure the turbulent thermal diffusivity in the σhot model and plot in Fig. 4 its radial profile and statistical deviation. For comparison, we also plot two radial profiles of κT that would result from constant αT, chosen such that they bracket the statistical deviations of the simulations results. They respectively correspond to PR ~ 0.8 and PR ~ 3.5 (assuming a constant α = 0.02). The large deviation around the mean value prevents any definitive and accurate measurement of PR, but our data suggest that PR is around unity. This is relatively close to PR = 0.3 that Pierens et al. (2012) used to investigate, with a diffusive model, how the turbulent diffusivity impacts on planet migration.

thumbnail Fig. 5

z = 0 planar maps of radial magnetic fields fluctuation (on the left) and density fluctuations (right plot) from the σhot simulation at t = 900 + τcool(R = 7R0).

Open with DEXTER

Turbulent diffusion must be compared to the radiative diffusion of temperature. In the optically thick approximation, the radiative diffusivity is T3. One can establish, using Eq. (15), that the dimensionless measure of radiative diffusivity αrad, defined similarly to αT (see Eqs. (25) and (26)), is equal to a few α. The turbulent transport of energy is then comparable to radiative transport in PP disks (see also Appendix A in Latter & Balbus 2012). Nevertheless, radiative transport is neglected in our simulations. As indicated by the discussion above, this shortcoming should be rectified in future. Note that radiative MHD simulations are challenging and Flock et al. (2013) performed the first global radiative MHD simulations of protoplanetary disks only very recently. Previous studies had been confined to the shearing box approximation (Hirose et al. 2006; Flaig et al. 2010; Hirose & Turner 2011) because of the inherent numerical difficulties of radiative simulations.

4. MHD simulations of PP disks with a dead zone

We now turn to the non-ideal MHD models in which the disk is composed of an inner turbulent region and a dead zone at larger radii. To make contact with previous work (Dzyurkevich et al. 2010; Lyra & Mac Low 2012), we first consider the case of a resistivity that is only a function of position. Such a simplification is helpful to understand the complex thermodynamics of the dead zone before moving to the more realistic case in which the resistivity is self-consistently calculated as a function of the temperature using Eq. (8).

4.1. The case of a static interface

At t = 900, model σhot has reached thermal equilibrium. We then set η = 0 for R < 3.5, and η = 10-3 for R ≥ 3.5, and restart the simulation. The higher value of η corresponds to a magnetic Reynolds number Rm ~ csH/η ~ 10. It is sufficient to stabilize the MRI, and as expected, the flow becomes laminar outward of R = 3.5 in a few orbits. The structure of the flow after 460 inner orbits (which roughly amounts to one cooling time at R = 7) is shown in Fig. 5. The left panel shows a snapshot of the radial magnetic field. The turbulent active region displays large fluctuations which can easily be identified. The interface between the active and the dead zone, defined as the location where the Maxwell stress (or, equivalently, the magnetic field fluctuations) drops to zero, stands at R = 3.5 until the end of this static non-ideal simulation.

As soon as the turbulence vanishes, the sharp gradient of azimuthal stress at the interface drives a strong radial outflow. As identified by Dzyurkevich et al. (2010) and Lyra & Mac Low (2012), we observe the formation of a pressure and a density maximum near the active/dead interface. In Fig. 6, we present the mean pressure perturbation we obtained in that simulation. The growth time scale of that structure is so short that the source term in the continuity equation (see Sect. 2.1) has no qualitative impact because it acts on the longer time scale τρ. This feature persists even after we freeze the temperature and use a locally isothermal equation of state as shown by the black solid line. We conclude that pressure maxima are robust features that form at the dead zone inner edge independently of the disk thermodynamics.

thumbnail Fig. 6

Pressure profile averaged over 200 orbits after t = 900 + τcool(R = 7) in σHot case. The black solid line show the pressure profile 200 inner orbits after freezing the temperature in the domain. The black dashed line plots the radial temperature profile in the ideal case (i.e. without a dead zone).

Open with DEXTER

thumbnail Fig. 7

The radial profile of temperature (solid red line) averaged over 200 orbits after t = 900 + τcool(R = 7R0). The red dashed line describes the temperature profile from the ideal case and the wave-equilibrium profile is shown by the plain black line. The black dashed lines show the threshold value TMRI used in Sect. 4.3 and (5/4)TMRI. The grey area shows the deviation of temperature at 3 sigmas. The vertical line shows the active/dead zone interface.

Open with DEXTER

We next turn to the temperature’s radial profile. Because we expect the dead zone to be laminar, we might assume the temperature to drop to a value near Tmin = 0.05 (see Sect. 2.3). However, the time averaged temperature radial profile plotted in Fig. 7 shows this is far from being the case. The temperature in the bulk of the dead zone levels off at T ~ 0.3, significantly above Tmin. This increased temperature cannot be attributed to Ohmic heating, since the magnetic energy is extremely small outward of R = 3.5 (at R = 4, the magnetic energy has dropped by a factor of 102 compared to R = 3.5, while at R = 5 it has been reduced by about 103). In contrast, there are significant hydrodynamic fluctuations in the dead zone, as shown on the right panel of Fig. 5. For example, density fluctuations typically amount to δρ/ρ of the order of a few percent at R = 5. The spiral shape of the perturbation indicate that these perturbations could be density waves propagating in the dead zone, likely excited at the dead/active interface. Such a mechanism would be similar to the excitation of sound waves by the active surface layers observed in vertically stratified shearing boxes simulations (Fleming & Stone 2003; Oishi & Mac Low 2009) and also the waves seen in global simulations such as ours (Lyra & Mac Low 2012) but performed using a locally isothermal equation of state. However, as opposed to Lyra & Mac Low (2012), the Reynolds stress associated with this waves is smaller in the dead zone than in the active zone, and amounts to α ~ 10-4. Such a difference is likely due to the smaller azimuthal extent we use here, as it prevents the appearance of a vortex (see Sect. 5.1 and Lyra & Mac Low 2012).

4.2. Wave dissipation in the dead-zone

It is tempting to use Eq. (22) to estimate the temperature that would result if the wave fluctuations were locally dissipated into heat (as is the case in the active zone). However, we obtain T ~ 0.1 which significantly underestimates the actual temperature. This is probably because the density waves we observe in the dead zone do not dissipate locally. We now explore an alternative model that describes their effect on the thermodynamic structure of the dead zone.

We assume that waves propagate adiabatically before dissipating in the form of weak shocks in the dead zone; this assumption is consistent with measured Mach numbers (~0.1) in the bulk of our simulated dead zone. A simple model for the weak shocks, neglecting dispersion, yields a wave heating rate : (27)where f stands for the wave frequency (Ulmschneider 1970; Charignon & Chièze 2013). The derivation of that expression is presented in Appendix B. Balancing that heating by the local cooling rate ⟨ ρ ⟩ σT4 yields an estimate of the temperature’s radial profile. (28)where the relation between the pressure and the disk scale height has been used. This expression provides an estimate of the dead-zone temperature. If one assumes that the waves are excited around the dead-zone interface, the frequency should be of order the inverse of the correlation time τc of the turbulent fluctuations at that location. Such a short period for the waves, much shorter than the disk cooling time, is consistent with the hypothesis that these waves behave adiabatically. Local simulations of Fromang & Papaloizou (2006) suggest τc/Torb ~ 0.15. Given the units of the simulations, this means that we should use f ~ 1. Taking δρ/ρ ~ 0.1 and H/R ~ 0.1, this gives T ~ 0.25, at R = 5, which agrees relatively well with the measured value. The agreement certainly supports the assumption that waves are generated at the interface.

In order to make the comparison more precise, we plot in Fig. 7 the radial profile of the temperature, computed using Eq. (28)3. The agreement with the temperature we measure in the simulation is more than acceptable in the bulk of the dead zone. There are however significant deviations at R ≥ 7, i.e. next to the outer boundary. Associated with these deviations, we also measured large, but non wave-like values of the quantity δρ/ρ at that location (see Fig. 5).

In order to confirm that the outer-buffer zone is responsible for these artefacts, we ran an additional simulation with identical parameters but with a much wider radial extent. The outer-radial boundary is located at R = 22, and thus the outer-buffer zone extends from R = 21.5 to R = 22. The additional simulation shows that indeed the suspicious temperature bump moves to the outer boundary again. The temperature in the bulk of the dead zone, on the other hand, remains unchanged and is in good agreement with the model σhot. Further details of this calculation are given in Appendix C.

In Appendix C we show that the density wave theory can also account for the density fluctuations as a function of radius. As waves propagate, their amplitudes decrease, not least because their energy content is converted into heat. Overall, the good match we have obtained between the wave properties’ radial profiles (both for the temperature and density fluctuations) and the analytical prediction strongly favour our interpretation that density waves are emitted at the active/dead zone interface and rule out any numerical artefacts that might be associated with the domain outer boundary.

thumbnail Fig. 8

Turbulent thermal flux profile averaged over 200 inner orbits after t = 900 + τcool(R = 7) from the σhot run. The central dotted region locates the percolation region behind active/dead zone interface.

Open with DEXTER

thumbnail Fig. 9

Maps of radial magnetic field at different time t = 1600 + t. t is given in inner orbital time.

Open with DEXTER

thumbnail Fig. 10

Space-time diagrams showing the turbulent activity evolution in the MHD simulation (left panel) and in the mean field model (right panel). In the case of the MHD simulation, α is averaged over Z and φ. In the case of the mean field model, the different white lines mark the front position for the additional models described in Sect. 4.3. The solid line represents the case of remnant heating in the dead zone, the dashed line no remnant heating, and the dotted line the case of constant thermal diffusivity over the whole domain. Note the different vertical scales of the two panels (see discussion in the text for its origin).

Open with DEXTER

The effect of the waves can be viewed as a flux of thermal energy carried outward. Such a thermal flux can be quantified using Eq. (24) and its radial profile in model σhot is shown in Fig. 8. It is positive and decreasing in the dead zone which confirms that thermal energy is transported from the active zone and deposited at larger radii. While it is rapidly decreasing with R, the important point of Fig. 8 is that it remains important in a percolation region outward of the dead zone inner edge (located at R = 3.5). Most of the wave energy that escapes the active region is transmitted to the dead zone over a percolation length, the size of the percolation region. We define the percolation region as the region inside the dead zone where the thermal transport amounts to more that half its value in the active zone. As shown by Fig. 8, the percolation length l ~ 0.5 in model σhot, which translates to about one scale height at that location.

Though the heat flux is greatest within about H of the dead zone interface, the action of the density waves throughout our simulated dead zone keeps temperatures significantly higher than what would be the case in radiative equilibrium. This would indicate that much of the dead zone in real disks would be hotter than predicted by global structure models that omit this heating source. In particular, this should have important implications for the location of the ice line, amongst other important disk features.

4.3. The case of a dynamic interface

We now move to the main motivation of this paper. At t = 1600 we restart model σhot but close the feedback loop between turbulence and temperature. The resistivity η is now given by Eq. (8) with η0 = 10-3. We set the temperature threshold to the value TMRI = 0.4, i.e. slightly above the typical temperature in the dead zone (which is about 0.25 as discussed above). This set-up ensures that at least half of the radial domain is “bistable”, i.e. gas at the same radius can support either one of two quasi-steady stable states, a laminar cold state, and a turbulent hot state.

As seen in Fig. 7, the temperature exceeds TMRI in the region 3.5 ≤ R ≤ 4.5 at restart. Thus, we would naively expect the interface to move to R ~ 4.5. As shown by snapshots of the radial magnetic field fluctuations in the (R,φ) plane at three different times in Fig. 9, indeed the front moves outwards. In doing so, it retains its coherence. However, as shown by the third panel of Fig. 9, the turbulent front does not stop at R ~ 4.5 but moves outward all the way to R ~ 5.56. This is confirmed by the left panel of Fig. 10, in which a space-time diagram of α indicates that the front reaches its stagnation radius in a few hundred orbital periods4. The following questions arise: can we predict the stalling radius? And can we understand this typical time scale? What can it tell us about the physical mechanism of front propagation?

thumbnail Fig. 11

α profiles averaged over 100 orbits after the fronts have reached their final position (α values are listed on the left axis). Left panel: σhot simulation which TMRI = 0.46. Middle and right panels: σcold simulation which TMRI = 0.1 and TMRI = 0.13 respectively. The black plain lines remind the temperature profiles from ideal cases. Black dashed lines show the threshold value TMRI and 5/4TMRI. The temperature values are listed on the right axis.

Open with DEXTER

Stalling radius:

we use the mean field model proposed by Latter & Balbus (2012) to interpret our results. The front-stalling radius can be calculated from the requirement that the integrated cooling and heating across the interface balance each other: (29)in which TDZ and TA stand for the temperature in the dead zone and in the active zone respectively. While the cooling part of the integral is well defined, the heating part Q+ is more uncertain. It is probably a mixture of turbulent heating (with an effective α parameter that varies across the front) and wave heating such as described in Sect. 4.1. For simplicity, we adopt the most naive approach possible and assume that Q+ is a constant within the front (equal to 1.5ΩT) as long as the temperature T is larger than TMRI, and vanishes otherwise. In that situation, we obtain, similarly to Latter & Balbus (2012) that the stalling radius Rc is determined by the implicit relation: (30)In Fig. 7, we draw two horizontal lines that correspond to T = TMRI and T = (5/4)TMRI. The condition given by Eq. (30) is satisfied for R ~ 5.25 which is very close to the critical radius actually reached in simulations (see Figs. 9 and the left panel of Fig. 10).

In order to investigate the robustness of this result, we have carried a further series of numerical experiments. First, we have repeated model σhot using a value of TMRI = 0.46 (instead of TMRI = 0.4). According to Fig. 7, we expect the interface now to propagate over a smaller distance. As shown by the left hand side panel of Fig. 11, this is indeed the case: the front stalls at R ~ 5 which is precisely the position where T = (5/4)TMRI.

We also computed two analogues of those models with σcold. This was done as follows: after the σcold ideal MHD simulations has reached a quasi thermal equilibrium (t = 900), we introduced a static dead zone from R = 3.5 to the outer boundary of the domain and we restarted this simulation. We waited for thermal equilibrium to close the feedback loop at t = 1600. We used TMRI = 0.1. A front propagates until the critical radius R ~ 6.5. It is as fast as the front of the σhot case. We show in the middle panel of Fig. 11 that, in accordance with the previous result, the front stops close to the critical radius predicted by Eq. (30). We repeated the same experiment using TMRI = 0.13 and show on the right panel of Fig. 11 that the front stalling radius is again accounted for by the same equation.

To summarize, the front stops close to the predicted value in every simulations we performed. This shows that the equilibrium position of a dead zone inner edge can be predicted as a function of the disk’s radiative properties and thermal structure. It is robust, and in particular does not depend on the details of the turbulent saturation in our simulations (such as the radial α profile).

Front speed and propagation:

the previous result is based on dynamical systems arguments. A front is static at a given radial location if the attraction of the active turbulent state balances the attraction of the quiescent state. However, the argument does not help us identify the front propagation speed. In order to understand its dynamics, we ran a set of mean field simulations akin to the “slaved model” in Latter & Balbus (2012) but using parameters chosen to be as close as possible to model σhot. We solve the partial differential equation for the temperature (we have dropped the overbar here for clarity): (31)To solve that equation, we have used the heating term measured during the ideal MHD simulation σhot (see Sect. 3). The thermal transport of energy is modeled by a simple diffusive law that is supposed to account for turbulent transport. As discussed in Sect. 4.1, the dead zone is heated by waves excited at the dead/active interface. It is thus highly uncertain (and one of the purpose of this comparison) whether a diffusive model adequately describes heat transport of this fashion. The value for κT is chosen so heat diffusion gives the same mean flux of thermal energy as measured in our simulation (shown in Fig. 8). We have simplified the radial profile in that plot and have assumed that the thermal flux vanishes outward of the percolation region but is uniform in the active zone and in the percolation region: (32)The front location Rf(t) is evaluated at each time step. It is the smallest radius where T < TMRI. As in the MHD simulation, we used TMRI = 0.4 and Tmin = 0.05. When Lp = 0, we have found the dead zone inner edge is static regardless of the value of κ0. As expected, the front displacement requires the active and the dead zone to be thermally connected. We have thus used a percolation length Lp = 0.5 and a value of κ0 = 2.5 × 10-4 in the active part of the disk which matches the thermal flux at the outer edge of the percolation region (see Fig. 8). The simulation is initialized with the same physical configuration as the MHD simulation: a dead zone extends from R = 3.5 to R = 8 and is in thermal equilibrium. As shown on the right panel of Fig. 10, we find that a front propagates outward. The critical radius Rc ≃ 5.5 − 6 is in agreement with the results of the simulation and with the argument based on energy conservation detailed above. However, there is clearly a difference in the typical time scale of the propagation: the front observed in the MHD simulation propagates five time faster than the front obtained in the mean field simulation.

thumbnail Fig. 12

Vortensity ((V − Vinit)/Vinit) map in the disk (R,φ) plane. The left panel shows the midplane vortensity in the σhot simulation with the small azimuthal extend. The right panel shows the midplane vortensity in the σhot simulation which has φ ∈  [0/2].

Open with DEXTER

In order to check the sensitivity of this results to some of the uncertain aspects of the mean field model, we have run three additional mean field simulations similar to that described above. In these additional models, we change one parameter while keeping all others fixed. In the first model, we have changed the minimum value of the temperature Tmin = 0.3 (this is a crude way to model the effect of the dead zone heating at long distances from the interface). In the second model, we have used κ0 = 5 × 10-4. In the last model, we have used Lp =  +∞, thereby extending turbulent heat transport to the entire radial extent of the simulation. The results are shown on the right panel of Fig. 10: in all cases, the front propagates outward slower than observed in the MHD simulation.

A slow front velocity is thus a generic feature of the diffusive approximation to heat transport. But conversely, it reveals that the fast fronts mediated by real MRI turbulence are controlled by non-diffusive heat transport. This in turn strongly suggests that the front moves forward via the action of fast non-local density wave heating, and not via the slow local turbulent motions of the MRI near the interface, as originally proposed by Latter & Balbus (2012). The transport via density waves occurs at a velocity of order cs. Thus, it takes a time Δtwaves ~ H/cs ≃ τorb/2π for thermal energy to be transported through the percolation length Lp ~ H. This is shorter than the typical diffusion time over the same distance Δtdiff = H2/κ > 100/2πτorb at R = 3.5, for the values of κ used in Eq. (31).

5. Instability and structure formation

Other issues that can be explored through our simulations are instability and structure formation at the dead/active zone interface and throughout the dead-zone. The extremum in pressure at the interface is likely to give rise to a vortex (or Rossby wave) instability (Lovelace et al. 1999; Varnière & Tagger 2006; Meheut et al. 2012). On the other hand, the interface will control both the midplane temperature and density structure throughout the dead zone; it hence determines the magnitude of the squared radial Brunt-Väisälä frequency ⟨ NR ⟩. The size and sign of this important disk property is key to the emergence of the subcritical baroclinic instability and resistive double-diffusive instability in the dead zone (Lesur & Papaloizou 2010; Latter et al. 2010; Klahr & Bodenheimer 2003; Petersen et al. 2007b,a). In this final section we briefly discuss these instabilities, leaving their detailed numerical analysis for a future paper.

5.1. Rossby wave instability

The Rossby wave instability has been studied recently by Lyra & Mac Low (2012) with MHD simulations that use a locally isothermal equation of state and a static dead zone. Vortex formation mediated by the Rossby wave instability was reported: a pressure bump forms at the interface which then triggers the formation of a vortex. It is natural to wonder how these results are modified when better account is made of the gas thermodynamics.

In our simulations of a static dead zone (Sect. 4.1), we also find that pressure maxima form. Even with a T dependent η, the bumps survive as the interface travels to its stalling radius. Moreover, the amplitude of the pressure bump is not modified compared to those observed in simulations of static dead zones. On account of the reduced azimuthal extend of our domains, no large-scale vortices formed. (The reduced azimuth was chosen so as to minimize the computational cost of our simulations.) In order to observe the development of the Rossby wave instability we performed one run identical to the σhot simulation except we extended the azimuthal domain: φ =  [0/2]. The η was a given function of position as in Sect. 4.1, thus the dead/active zone interface was fixed. The right panel of Fig. 12 shows a late snapshot of this run, in which a vortex has appeared similar to those observed by Lyra & Mac Low (2012). It survives for many dynamical time scales. Consistent with the results of Lyra & Mac Low (2012), we measure α ≃ 0.01 in the dead zone. This is two orders of magnitude larger than in the σhot model. The dead zone is consequently much hotter in this simulation, which means the vortex plays a crucial role in both accretion and the thermal physics of the dead zone. Simulations are currently underway to investigate the robustness of this results and the vortex survival when we close the feedback loop (setting η = η(T)). This will be the focus of a future publication.

5.2. Subcritical baroclinic and double-diffusive instabilities

Another potentially interesting feature of the interface is the strong entropy gradient that might develop near the interface. It could also impact on the stability of the flow on shorter scales, giving rise potentially to the baroclinic instability (Lesur & Papaloizou 2010; Raettig et al. 2013) or to the double-diffusive instability (Latter et al. 2010). Both instabilities are sensitive to the sign and magnitude of the entropy gradient, which is best quantified by the Brunt-Väisälä frequency NR: (33)In general, we find that ⟨ NR ⟩ 2 takes positive values of order 1020% of the angular frequency squared. Negative values sometimes appear localized next to the interface. Positive values rule out both the baroclinic instability and the double-diffusive instability, and indeed, we see neither in the simulations. However, we caution against any premature conclusions about the prevalence of these instabilities in real disk. First, the source term in the continuity equation might alter the density radial profile and, consequently, the Brunt-Väisälä frequency in the dead zone. Second, it is well known that these instabilities are sensitive to microscopic heat diffusion. We do not include such a process explicitly in our simulations. Instead, there is numerical diffusion of heat by the grid, the nature of which may be unphysical. Both reasons preclude definite conclusions at this stage. Dedicated and controlled simulations are needed to assess the existence and nonlinear development of these instabilities in the dead zone.

6. Conclusion

In this paper, we have performed non-ideal MHD simulations relaxing the locally isothermal equation of state commonly used. We have shown the active zone strongly influences the thermodynamic structure of the dead zone via density waves generated at the interface. These waves transport thermal energy from the interface deep into the dead zone, providing the dominant heating source in its inner 20H. As a consequence, the temperature never reaches the very low level set by irradiation. In the outer regions of the dead-zone, however, temperatures will be set by the starlight reprocessed by the disk’s upper layers (Chiang & Goldreich 1997; D’Alessio et al. 1998). It is because the wave generation and dissipation is located at the midplane that waves should so strongly influence the thermodynamic structure of the dead zone. Note also that density waves generated by the dead/active zone edge are stronger than those excited by the warm turbulent upper layers of the disk as seen in stratified shearing box (Fleming & Stone 2003).

Another result of this paper concerns the dynamical behaviour of the dead-zone inner edge. We find the active/dead interface propagates over several H (i.e. a few tens of an au) in a few hundreds orbits. All the simulated MRI fronts reached a final position that matched the prediction made by a mean field approach (Latter & Balbus 2012), which appeals to dynamical systems arguments. As the gas here is bistable, it can fall into either a dead or active state; the front stalls at the location where the nonlinear attraction of the active and dead states are in balance. In contrast, the mean-field model fails to correctly predict the velocity of the simulated fronts. We find that a diffusive description of the radial energy flux yields front speeds that are too slow. In fact, the simulations show that fronts move rapidly via the efficient transport of energy by density waves across the interface. Fronts do not propagate via the action of the slower MRI-turbulent motions.

In addition, we have used our simulations to probe the thermal properties of turbulent PP disks. We have constrained the turbulent Prandtl number of the flow to be of order unity. We have also quantified the turbulent fluctuations of temperature: they are typically of order a few percent of the local temperature. However, their origin – adiabatic compression vs. reconnection – is difficult to assess using global simulations. In-depth dedicated non-isothermal shearing-box simulations will help to distinguish the dominant cause of the temperature fluctuations. Finally, we have made a first attempt to estimate the radial profile of the radial Brunt-Väisälä frequency NR in the dead zone. This quantity is the key ingredient for the development of both the subcritical baroclinic instability and the resistive double-diffusive instability. We find that NR can take both positive and negative values at different radii; but we caution that these preliminary results require more testing with dedicated simulations.

Several improvements are possible and are the basis of future work. As discussed in Sect. 5, one obvious extension is to investigate the fate and properties of emergent vortices at the dead-zone inner edge. This can be undertaken with computational domains of a wider azimuthal extent. Such simulations may investigate the role of large-scale vortices, and the waves they generate, on the thermodynamic structure of the dead-zone. They can also observe any feedback of the thermodynamics on vortex production and evolution. We also plan to investigate other magnetic field configurations. For example, vertical magnetic fields might disturb the picture presented here because of the vigorous channel modes that might develop in the marginal gas at the dead-zone edge (Latter et al. 2010). Such an environment may militate against the development of pressure bumps and/or vortices.

Our results, employing the cylindrical approximation, represents a thin region around the PP disk midplane. These results must be extended so that the vertical structure of the disk is incorporated. An urgent question to be addressed is the location of density wave dissipation in such global models. Waves can refract in thermally stratified disks and deposit their energy at the midplane or in the upper layers depending on the type of wave and the stratification (Bate et al. 2002, and references therein). In particular, Bate et al. demonstrate that large-scale axisymetric (and low m) density waves (fe modes) propagate upwards, as well as radially, until they reach the upper layers of the disk, at which point they transform into surface gravity waves and propagate along the disk surface. But this is only shown for vertically polytropic disks in which the temperature decreases with z and for waves with relatively low initial amplitudes. In PP-disk dead zones we expect the opposite to be the case, and it is uncertain how density waves behave in this different environment.

Finally, in this work we have increased the realism of one aspect of the physical problem, the dynamics of the turbulence, (via direct MHD simulations of the MRI), but have greatly simplified the physics of radiative cooling. As a result, the simulations presented here are still highly idealized and several improvements should be the focus of future investigations. For example, our approach completely neglects the fact that dust sublimates when the temperature exceeds ~1500 K. As a result, the opacity drops by up to four orders of magnitude with potentially dramatic consequences for the disk energy budget (leading to a increase of the cooling rate Q of the same order, as opposed to our assumption of a constant σ). The radial temperature profile of our simulations indicates that the sublimation radius should be located 35 disk scaleheights away from the turbulent front. Given that the dominant dynamical process we describe here is mediated by density waves characterized by a fast time scale (compared to the turbulent and radiative time scales), we do not expect the front dynamics to be completely altered. Nevertheless, the relative proximity between the sublimation radius and the turbulent front is still likely to introduce quantitative changes. Clearly, the thermodynamics of that region is more intricate than the simple idealized treatment we use in this work. This further highlights the need, in future work, for a more realistic treatment of radiative cooling (for example using the flux limited diffusion approximation with appropriate opacities, in combination with vertical stratification. Such simulations will supersede our heuristic cooling law with its constant σ parameter. They are an enormous challenge, but will be essential to test the robustness of the basic results we present here.


1

One consequence is that the ratio H/R which measures the relative importance of thermal and rotational kinetic energy, is no longer an input parameter, as in locally isothermal simulations, but rather the outcome of a simulation.

2

With our choice of parameters, we have Q ~ R-2.5 and . This means that the agreement between Q and improves with radius. At small radial distances, dust sublimates and our model breaks down (see Sect. 6).

3

We used the azimuthally and temporally averaged simulation data for H/R and δρ/ρ when computing the temperature using Eq. (28).

4

To check that Eq. (8) is an acceptable simplification of the actual exponential dependence of η with temperature, we ran an additional simulation, with the same parameters as for the σhot but using Eq. (7) for η. We find no difference with our fiducial model: the front propagates equally fast and stops at R ≃ 5.3, i.e. almost the same radius as for the σhot case.

Acknowledgments

J.F. and S.F. acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/20072013)/ERC Grant agreement No. 258729. H.L. acknowledge support via STFC grant ST/G002584/1. The simulations presented in this paper were granted access to the HPC resources of Cines under the allocation x2012042231 and x2013042231 made by GENCI (Grand Equipement National de Calcul Intensif).

References

Online material

Appendix A: Opacity law in the buffer zones

Here, we give the functional form of σ we used to prevent the temperature to drop at the inner edge of the simulation. (A.1)where σ0 stands for σhot or σcold depending on the model. The value of σ is kept constant in the outer buffer.

Appendix B: Wave heating

The large-scale density waves witnessed in our simulations develop weak-shock profiles, which are controlled by a competition between nonlinear steepening and wave dispersion. Keplerian shear may also play a role as it “winds up” the spiral and decreases the radial wavelength; though by the time this is important most of the wave energy has already dissipated.

A crude model that omits the strong dispersion inherent in our large-scale density waves nevertheless can successfully account for the energy dissipation in the simulations. In such a model the density wave profiles are dominated by steepening and can thus be approximated by a sawtooth shape propagating at the sound speed velocity cs (see Fig. B.1). The evolution of the amplitude of such isentropic waves is given by Landau & Lifshitz (1959). In the wave frame of reference, the gas velocity at the shock crest evolves over time as the shock wave dissipates: (B.1)where v0 is the excitation amplitude of the wave and λ0 its wavelength, assumed to be conserved over the wave propagation. The mean mechanical energy embodied in one wave period at time t is given by an integral over radius: (B.2)where λ is the wavelength at time t. We next compute the mechanical energy radial flux through a unit surface: (B.3)where δρ is the difference between the shocked and pre-shocked density. To compute the last equality, we have used the fact that, under the weak shock approximation, the wave evolution is isentropic and thus δρ/ρ = v/cs (Landau & Lifshitz 1959). The wave energy and its flux are related through the following conservation law: (B.4)while the dissipation rate of the wave is expressed using the mechanical energy conversion into heat per unit time: (B.5)where f = cs/λ0 is the wave frequency. We used the density fluctuations in our simulations as the difference between the shocked and pre-shocked density to estimate the local wave heating in Eq. (27).

thumbnail Fig. B.1

Velocity fluctuations profile: series of “teeth” modelling the wave shocks.

Open with DEXTER

In addition to that estimate, the thermal energy flux divergence can be used, through Eq. (B.4), as a way to estimate the radial variation of the wave amplitude: (B.6)This must equal the energy rate released as thermal heat given by D. Combining the last two expression thus provides an expression for the radial decay of the wave amplitude: (B.7)The first term is the geometrical term that described the wave dilution as it propagates cylindrically. The second and third terms are specific to waves propagating in stratified media where mean pressure and sound speed are not uniform. Finally, the last term of the right hand side of this equation reflects the wave damping by shocks. In our simulations, all four terms are of comparable importance.

Appendix C: Model with an extended radial extent

In this section of the appendix we describe the results obtained in the radially extended σhot run. We initiated MHD turbulence in that model in the absence of any dissipative term until the region between R = 1 and R = 3.5 has reached thermal equilibrium (t = 900). At that point, we set η0 = 10-3 for R ≥ 3.5. As expected, a dead zone quickly appeared at those radii. Because of the prohibitic computational coast of that simulation (there are 960 cells in the radial direction!), we were not able to run that model until thermal equilibrium is established at all radii. Indeed, the cooling time becomes very long at large radius. Instead, we followed the time evolution of the temperature at nine locations in the dead zone. Due to the absence of turbulent heating, we found that the temperature slowly decreases with time. Assuming this decrease is due to a combination of cooling (resulting from the cooling function) and wave heating, it can be modelled using the following differential equation (C.1)where the term accounts for the (unknown) wave heating. We fitted the time evolution of the temperature during the duration of the simulation (~1000 orbits) in order to obtain a numerical estimate of at those nine locations. An example of that fit is provided in the insert of Fig. C.1. Using these value, we can obtain an estimate for the equilibrium temperature in the dead zone, as shown in Fig. C.1. The comparison with the estimate of Eq. (28) provided by the black line is excellent everywhere in the dead zone.

As a sanity check, we test here if the wave amplitude decreases as their energy content is converted into heat. We show in Fig. C.2 the mean density fluctuation profile obtained in the extended σhot simulation. It exhibits two maxima close to the dead zone inner edge located at R(g,1) and R(g,2) and the amplitude decay with radius. The reason why we see two such maxima is not clear but might be due to waves originated at the dead/active interface as well as waves excited as the location of the pressure maximum. In any case, we found that modelling the amplitude of fluctuations as the signature of a combination of two waves generated at R(g,i) gives acceptable results. We use an explicit scheme to numerically integrate Eq. (B.7) from the wave generation locations R(g,i). The two waves are excited with the amplitude measured at R = R(g,i) with the frequency 1/f = 0.15τcool(R(g,i)). We plot the solution thus obtained in Fig. C.2. The good agreement between the analytical solution and the profile gives a final confirmation that waves control the dead zone thermodynamics.

thumbnail Fig. C.1

Temperature profile ⟨ T ⟩ averaged over 200 orbits after t = 900 + τcool(R = 7) obtained in the extended σhot case. The black line show wave-equilibrium temperature profile. Red dots show the extrapolated temperature for 7 radii. The inserted frame show in red the temperature evolution after t = 900 at R = 9.8R0 in the extended σhot simulation. On this subplot, the black plain line shows the accurate fit and the black dashed lines show the two extreme fits used to determine the measurement error.

Open with DEXTER

thumbnail Fig. C.2

Profile of density standard deviation ⟨ δρ ⟩ averaged over 200 inner orbits after t = 900 + τcool(R = 7) in the extended σhot case. The black plain line shows the density fluctuation amplitude deducted from the two waves model. The excitation locations used in the two waves model are shown by blue dashed lines.

Open with DEXTER

Appendix D: Simulation with a higher resolution

Here we present a brief test of the impact of spatial resolution on our results. We have restarted the σhot run from the thermal equilibrium of the ideal MHD case (t = 600) and the static dead zone case (t = 900), with twice as many cells in each direction. The resolution is (640,160,160). We show the averaged temperature profile of both cases in Fig. D.1. Because of the large computational cost, each models are integrated for 200 orbits and time averages are only performed over the last 100 orbits. The temperature profiles of each case are very close to those obtained with the fiducial runs. We conclude that resolution as little impact on our results.

thumbnail Fig. D.1

Profiles of temperature ⟨ T ⟩ averaged over 100 inner orbits after t = 700 and t = 1000 in the highly resolved σhot case (plain lines). The dashed lines remind the temperature profiles obtained with the low resolution run. For both resolutions, the temperature in the ideal case is shown by a red line and the temperature in the static dead zone case is shown by a black line.

Open with DEXTER

All Figures

thumbnail Fig. 1

Temperature profiles ⟨ T ⟩ averaged over 900 inner orbits (red curves). Black plain lines show their corresponding theoretical profiles. Black dashed lines show their corresponding initial profiles. Top panel: σcold case. Bottom panel: σhot case. The subframes inserted in the upper right of each panel show the alpha profiles obtained in both cases (red lines). The dashed lines show the analytic profiles used to compute the theoretical temperature profile (see text), given by α = 6.3 × 10-2 − 7.1 × 10-3R/R0 and α = 4.6 × 10-2 − 5.7 × 10-3R/R0 respectively.

Open with DEXTER
In the text
thumbnail Fig. 2

Temperature standard deviation profiles averaged over 900 inner orbits. The plain line show the result from the σhot case and the dashed curve show the result from the σcold run.

Open with DEXTER
In the text
thumbnail Fig. 3

Blue dots map the scaling relation between temperature and density fluctuations δT and δρ in the σhot simulation. 900 events are randomly selected in the whole domain (except buffer zones) among 900 inner orbits. The red line shows the linear function (slope γ − 1)) corresponding to the adiabatic scaling.

Open with DEXTER
In the text
thumbnail Fig. 4

Thermal diffusivity’s radial profile averaged over 900 inner orbits. Dashed lines show the mean thermal diffusivity for two homogeneous αT (αT = 0.02 and αT = 0.004). The grey area delimit one standard deviation around the mean thermal diffusivity.

Open with DEXTER
In the text
thumbnail Fig. 5

z = 0 planar maps of radial magnetic fields fluctuation (on the left) and density fluctuations (right plot) from the σhot simulation at t = 900 + τcool(R = 7R0).

Open with DEXTER
In the text
thumbnail Fig. 6

Pressure profile averaged over 200 orbits after t = 900 + τcool(R = 7) in σHot case. The black solid line show the pressure profile 200 inner orbits after freezing the temperature in the domain. The black dashed line plots the radial temperature profile in the ideal case (i.e. without a dead zone).

Open with DEXTER
In the text
thumbnail Fig. 7

The radial profile of temperature (solid red line) averaged over 200 orbits after t = 900 + τcool(R = 7R0). The red dashed line describes the temperature profile from the ideal case and the wave-equilibrium profile is shown by the plain black line. The black dashed lines show the threshold value TMRI used in Sect. 4.3 and (5/4)TMRI. The grey area shows the deviation of temperature at 3 sigmas. The vertical line shows the active/dead zone interface.

Open with DEXTER
In the text
thumbnail Fig. 8

Turbulent thermal flux profile averaged over 200 inner orbits after t = 900 + τcool(R = 7) from the σhot run. The central dotted region locates the percolation region behind active/dead zone interface.

Open with DEXTER
In the text
thumbnail Fig. 9

Maps of radial magnetic field at different time t = 1600 + t. t is given in inner orbital time.

Open with DEXTER
In the text
thumbnail Fig. 10

Space-time diagrams showing the turbulent activity evolution in the MHD simulation (left panel) and in the mean field model (right panel). In the case of the MHD simulation, α is averaged over Z and φ. In the case of the mean field model, the different white lines mark the front position for the additional models described in Sect. 4.3. The solid line represents the case of remnant heating in the dead zone, the dashed line no remnant heating, and the dotted line the case of constant thermal diffusivity over the whole domain. Note the different vertical scales of the two panels (see discussion in the text for its origin).

Open with DEXTER
In the text
thumbnail Fig. 11

α profiles averaged over 100 orbits after the fronts have reached their final position (α values are listed on the left axis). Left panel: σhot simulation which TMRI = 0.46. Middle and right panels: σcold simulation which TMRI = 0.1 and TMRI = 0.13 respectively. The black plain lines remind the temperature profiles from ideal cases. Black dashed lines show the threshold value TMRI and 5/4TMRI. The temperature values are listed on the right axis.

Open with DEXTER
In the text
thumbnail Fig. 12

Vortensity ((V − Vinit)/Vinit) map in the disk (R,φ) plane. The left panel shows the midplane vortensity in the σhot simulation with the small azimuthal extend. The right panel shows the midplane vortensity in the σhot simulation which has φ ∈  [0/2].

Open with DEXTER
In the text
thumbnail Fig. B.1

Velocity fluctuations profile: series of “teeth” modelling the wave shocks.

Open with DEXTER
In the text
thumbnail Fig. C.1

Temperature profile ⟨ T ⟩ averaged over 200 orbits after t = 900 + τcool(R = 7) obtained in the extended σhot case. The black line show wave-equilibrium temperature profile. Red dots show the extrapolated temperature for 7 radii. The inserted frame show in red the temperature evolution after t = 900 at R = 9.8R0 in the extended σhot simulation. On this subplot, the black plain line shows the accurate fit and the black dashed lines show the two extreme fits used to determine the measurement error.

Open with DEXTER
In the text
thumbnail Fig. C.2

Profile of density standard deviation ⟨ δρ ⟩ averaged over 200 inner orbits after t = 900 + τcool(R = 7) in the extended σhot case. The black plain line shows the density fluctuation amplitude deducted from the two waves model. The excitation locations used in the two waves model are shown by blue dashed lines.

Open with DEXTER
In the text
thumbnail Fig. D.1

Profiles of temperature ⟨ T ⟩ averaged over 100 inner orbits after t = 700 and t = 1000 in the highly resolved σhot case (plain lines). The dashed lines remind the temperature profiles obtained with the low resolution run. For both resolutions, the temperature in the ideal case is shown by a red line and the temperature in the static dead zone case is shown by a black line.

Open with DEXTER
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.