{'fraction_non_alphanumeric': 0.030556449018977163, 'fraction_numerical': 0.01998805311767679, 'mean_word_length': 4.82860203535083, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heating and Cooling Protostellar Disks

We examine heating and cooling in protostellar disks using 3-D radiation-MHD calculations of a patch of the Solar nebula at 1 AU, employing the shearing-box and flux-limited radiation diffusion approximations. The disk atmosphere is ionized by stellar X-rays, well-coupled to magnetic fields, and sustains a turbulent accretion flow driven by magneto-rotational instability, while the interior is resistive and magnetically dead. The turbulent layers heat by absorbing the light from the central star and by dissipating the magnetic fields. They are optically-thin to their own radiation and cool inefficiently. The optically-thick interior in contrast is heated only weakly, by re-emission from the atmosphere. The interior is colder than a classical viscous model, and isothermal. The magnetic fields support an extended atmosphere that absorbs the starlight 1.5 times higher than the hydrostatic viscous model. The disk thickness thus measures not the internal temperature, but the magnetic field strength. Fluctuations in the fields move the starlight-absorbing surface up and down. The height ranges between 13% and 24% of the radius over timescales of several orbits, with implications for infrared variability. The fields are buoyant, so the accretion heating occurs higher in the atmosphere than the stresses. The heating is localized around current sheets, caused by magneto-rotational instability at lower elevations and by Parker instability at higher elevations. Gas in the sheets is heated above the stellar irradiation temperature, even though accretion is much less than irradiation power when volume-averaged. The hot optically-thin current sheets might be detectable through their line emission.


INTRODUCTION
Measurements of the temperatures, densities and flow velocities inside protostellar disks are vitally important for understanding planet formation. New probes of the conditions are made possible by the detections of infrared line emission from molecules including H 2 O, OH, CO, CO 2 , C 2 H 2 and HCN, in the central 3 AU of the disks around T Tauri stars (Najita et al. 2003;Carr et al. 2004;Carr & Najita 2008;Salyk et al. 2008;Pontoppidan et al. 2010). However interpreting the line emission requires some knowledge of the power source. The line excitation temperatures range from a few hundred to a few thousand Kelvin, and can be reached in gas near 1 AU through accretion heating, if the stresses in the disk atmosphere approach the local gas pressure (Glassgold et al. 2004).
The accretion stresses in T Tauri disks likely come from magneto-rotational turbulence (Balbus & Hawley 1998). Magnetic fields couple only to gas that is sufficiently ionized, and while much of the disk is too cold for thermal ionization, the top and bottom layers are ionized by the X-rays from the young star (Igea & Glassgold 1999;Ilgner & Nelson 2006). The turbulent, magnetically-active surface layers are optically-thin, and thus produce line emission directly (Bai & Goodman 2009).
Another important heat source is the light from the central star. The starlight enters the disk at a grazing angle, and so is absorbed high in the atmosphere, well above the surface of optical depth unity for the disk's own longer-wavelength continuum emission (Chiang & Goldreich 1997). The stellar illumination dominates the atmospheric heating at 1 AU in traditional viscous models, while the accretion dominates in the interior, at typical accretion rates around 10 −8 M ⊙ yr −1 (Dullemond et al. 2007).
Here we report first-principles radiation-MHD calculations of the heating and cooling, treating the conversion of gravitational to magnetic energy through magneto-rotational instability, the dissipation of the magnetic energy as heat in the gas, and the emission and escape of the disk radiation. The calculations build on recent work by Flaig et al. (2010) in including the stellar irradiation heating and allowing distinct gas and radiation temperatures. We also choose a much lower surface density of 1 000 g cm −3 , similar to that at 1 AU in the minimum-mass Solar nebula (Hayashi 1981), so that temperatures are too low for collisional ionization. The ionization is dominated by the stellar energetic photons, which are absorbed in the atmosphere, leaving a magnetically-decoupled interior dead zone.

METHODS
We solve the equations of radiation-MHD using the ZEUS code , conserving total energy following Hirose et al. (2006). The direct starlight is treated by explicitly solving the transfer equation, while the disk's own emission is treated using the fluxlimited diffusion module of Turner & Stone (2001). The gas and radiation energy equations are advanced through time in fully-implicit fashion using a generalized Newton-Raphson method ) with multigrid solver. Radiation emission, absorption and diffusion are computed simultaneously, as needed for accuracy when two of them almost cancel. We use opacities for homogeneous spherical silicate grains of normal iron content, from Semenov et al. (2003). Opacities for the direct starlight are averaged over the stellar spectrum, while for the disk's radiation we use the Planck and Rosseland means at the local temperature.
The star has one-half the mass and twice the radius of our Sun, and the spectrum of a 4 000 K blackbody. The starlight enters the domain top and bottom at a grazing angle, 0.05 radians ≈ 3 • from the horizontal. The heating is computed each timestep by solving the transfer equation in the horizontally-averaged density profile. We neglect the small amount of absorption that would occur outside the MHD domain.
The Ohmic resistivity varies in space and time, and is found by interpolating in a precomputed lookup table whose three axes are the density, temperature and local ionization rate. The table holds equilibrium solutions of the simplified Ilgner & Nelson (2006) recombination network including reactions on grains (their Model 4). The grain size of 0.1 µm is chosen to give almost the same geometric cross-section per unit mass as the particle size distribution used for the opacities. The two differ by less than 15%. The ionization rate due to the stellar X-rays and long-lived radionuclides is computed from the horizontally-averaged overlying mass column using the prescriptions in Turner & Sano (2008), taking an X-ray luminosity 2 × 10 30 erg s −1 and a long-lived radionuclide ionization rate 1.4 × 10 22 s −1 at nominal dust abundance. The magnetic diffusivity is inversely proportional to the resulting electron fraction (Blaes & Balbus 1994) and is very high in the weakly-ionized disk interior, reaching 6 × 10 24 cm 2 s −1 (or if short-lived 26 Al were present, 2 × 10 21 ). To avoid short timesteps, we cap the diffusivity, limiting it to no more than 5 × 10 16 cm 2 s −1 . This is below the threshold of 5 × 10 20 for diffusion to overcome the shearing of radial fields under differential rotation (eq. 1 of Turner & Sano 2008), so weak magnetic activity is expected in the dead zone. However this is probably unimportant for the active layers, as we obtain similar MHD results with a cap ten times higher.
T Tauri disk atmospheres appear dust-depleted by factors of 10 1−4 compared with the interstellar medium, based on their mid-infrared colors (Furlan et al. 2006(Furlan et al. , 2009. Atmo-spheric depletion is consistent with modeling of some protostellar disks over wider wavelength ranges, incorporating spatially-resolved data (Duchêne et al. 2003;Pinte et al. 2008). Furthermore, the concentration of solid material in the interior is helpful for planet formation. We therefore assume a dust-to-gas mass ratio 10 −4 and reduce all dust effects 100 times below nominal, including the opacities, the cross-sections for grain surface reactions, and the radionuclide decay ionization. The dust is well-mixed, and thermally coupled to the gas so they share a single temperature, which is allowed to differ from the local radiation temperature. The gas has mean molecular weight 2.3 and an ideal equation of state with γ = 1.4. A density floor of 10 −5 times the initial midplane density is applied to prevent high Alfven speeds from halting the calculations. We checked the results are basically unchanged with a density floor ten times lower. Flow velocities also are limited to no more than 10 times the shear speed, 3 2 Ω times the domain width. The domain is a box of size 0.07 × 0.28 × 0.7 AU along the radial, orbital and vertical directions, centered in the midplane at 1 AU and divided into 32 × 64 × 320 grid cells. The boundaries are shearing-periodic, periodic, and outflow, respectively. Like Flaig et al. (2010), we assume the magnetic fields are vertical on the top and bottom boundaries. Results are similar using the standard ZEUS outflow boundary condition, where the fields on the boundaries are computed from the extrapolated EMFs. We also added a small artificial resistivity near the boundaries for numerical stability (Hirose et al. 2009). The disk radiation on the vertical boundaries is placed in the streaming limit, with flux equal to the product of the radiation energy density and the speed of light.
Initially the gas is isothermal at 280 K and in hydrostatic balance, while the magnetic field is the sum of an 0.02-Gauss uniform vertical component with pressure 3 × 10 5 times less than the midplane gas pressure, and a vertical component varying sinusoidally in the radial direction with maximum pressure 1/270th the midplane gas pressure. Results were similar in a further calculation with the pressure in the sinusoidal component 100 times less. We also compared an ideal-MHD version of our calculation against the results of Flaig et al. (2010), obtaining comparable profiles of the plasma beta parameter, turbulent Mach number and density fluctuations, despite our lower surface density.
Below we relate the MHD results to the classical viscous, irradiated disk model with the same surface density and accretion rate. We follow D' Alessio et al. (1998) in solving the equations of vertical hydrostatic equilibrium, radiative equilibrium in the Eddington approximation, and thermal balance, together with a transfer equation for the illuminating starlight. Scattered radiation is neglected in both the viscous and MHD calculations.

RESULTS
The turbulent layers on the top and bottom are optically-thin to their own radiation, while the magnetically-dead interior is optically-thick. The base of the turbulent layer, defined by unit Elsasser number v 2 Az /ηΩ (Sano & Inutsuka 2001;Turner et al. 2007), has height 0.08 AU and Rosseland mean optical depth 0.39, so the disk photosphere falls within the dead zone. The starlight on the other hand is absorbed in the turbulent layers, dominating the mean heating above 0.11 AU (Figure 1). The turbulent layers are more opaque in the optical than the infrared, and absorb starlight better than they re-radiate, becoming several times hotter than the disk radiation passing through them. The horizontal averaging in the irradiation prescription yields mean temperatures higher than the true irradiation temperature of 410 K found in the viscous model. Below the disk photosphere the gas and radiation temperatures match. The interior, in contrast, receives no direct starlight and is heated by the re-radiation from the hot surface layers. The cooling is almost wholly radiative throughout. Turbulent advection is negligible, and unlike the viscous models in D' Alessio et al. (1998) figure 10, carries no significant heat to the interior. The heating and cooling together yield the temperature profile in Figure 2.
The turbulent layers are supported by magnetic forces above 0.12 AU, while the interior is gas-pressure-supported. The magnetic dissipation is the second-largest contribution to the volume-integrated heating. The magnetic fields are buoyant and rise while dissipating, so the  (red) and corresponding viscous model (blue). Shadings with the same colors indicate the penetration of the starlight to unit optical depth, while green shading marks the part of the disk that is opaque to its own radiation. The gas and dust (solid) is hotter than the disk radiation (dotted) at low optical depths. The profiles are averaged horizontally, and in time from 50 to 200 orbits. accretion heating occurs higher on average than the accretion stress ( Figure 3). By contrast, in the classical model the stress and heating occur at the same place and time, and both result from an effective viscosity, nature unspecified. The dissipation rate in our calculation is comparable to the local pressure p(z) times the shear rate 3 2 Ω, the ratio increasing from 0.4 near the turbulent layer base to 2 at the boundary. Ratio unity is sufficient to produce CO molecular line emission in the models of Glassgold et al. (2004) at lower columns around 10 21 cm −2 . The dissipation is also high relative to the stress in two narrow layers near ±0.07 AU where turbulent magnetic fields are destroyed on entering the dead zone. Overall, the stresses correspond to a mean mass flow rate 1.3×10 −8 M ⊙ yr −1 , within the range among classical T Tauri stars (Muzerolle et al. 2003).
The magnetic field structure varies with height. Near the base of the turbulent layers, current sheets are typically close to horizontal, while high in the atmosphere, current sheets are most often vertical. Near the base, the current sheets separate magnetic fields moving radially inward and outward under the magneto-rotational instability, while higher in the atmosphere, the current sheets separate magnetic fields lying along adjacent orbits and having out-of-phase vertical undulations. The upward slope of one field line annihilates against the downward slope of the next (Figure 4). The undulations come from the Parker instability, whose fastest-growing modes have large radial wavenumber under Keplerian differential rotation. By solving in each grid cell the dispersion relation of Terquem & Papaloizou (1996) eq. 33, describing localized non-axisymmetric disturbances on toroidal magnetic fields, we have verified that MRI modes (k z ≫ k r ) grow fastest near the base, where the field is weak, while Parker modes (k r ≫ k z ) grow fastest high in the atmosphere, where the plasma beta is low and magnetic buoyancy strong. A similar pattern appears in ideal-MHD models of strongly-magnetized disks (Johansen & Levin 2008). The Parker modes grow faster than the orbital frequency in locations where the fields decline steeply with height, including just below reversals in the sign of the dominant toroidal field, and along the tops of magnetic flux tubes. Among the solutions of the dispersion relation are still faster-growing Parker modes with azimuthal wavelengths 2-3 times our domain size. Future experiments with larger boxes are indicated. Since the dissipation is concentrated in thin current sheets (figure 4), it dominates the heating at some places and times. Dissipation is stronger than irradiation in a volume fraction 3.2 × 10 −4 , time-averaging from 50 to 200 orbits and considering the layers above the disk photosphere (|z| > 0.18 AU). Temperatures in the current sheets are up to 50% greater than in the gas heated only by starlight. Finally, mass leaves the top and bottom boundaries in an outflow comparable to those in isothermal calculations by Suzuki et al. (2010).
Most of the energy input between 50 and 200 orbits (97%) comes from stellar irradiation. Of the 3% coming from work done by magnetic stresses, 52% is converted directly to heat via magnetic dissipation, while 26% is first converted to kinetic energy and then dissipated via shocks and viscosity. The remainder crosses the boundaries as magnetic and turbulent kinetic energy and the gravitational potential energy of the small amount of escaping gas. Overall, energy ultimately leaves the box almost entirely by radiation diffusion. The field strength in Gauss is indicated by colors: stronger fields are red, weaker fields blue. Inside the yellow isosurfaces the magnetic dissipation rate exceeds 10 −9 erg cm −3 s −1 . Field lines drawn in green are those passing through a zone of strong magnetic dissipation at top center. Right: close-up of the green field lines, viewed from the star (top) and along the orbit (bottom). The dissipation is concentrated in a current sheet of narrow radial extent, separating two counter-undulating bundles of field lines.
The surface where the starlight is absorbed moves up and down irregularly over intervals of several orbits ( Figure 5) as fluctuations in the magnetic fields lead to changes in the density profile. The fluctuations are associated with episodic reversals in the toroidal field, which lift away from the midplane at a rate that increases with height, producing the butterflywing pattern visible in figure 5. A majority of the volume-integrated accretion heating occurs near these reversals. As in previous stratified MRI turbulence calculations (e.g. Miller & Stone 2000;Hirose et al. 2006;Blaes et al. 2007;Flaig et al. 2010), the density in the magnetically-dominated atmosphere fluctuates horizontally by more than a decade, and the velocity dispersion exceeds the sound speed, approaching the Alfven speed.

DISCUSSION
Diagnosing conditions in planet-forming disks requires an understanding of the power sources for the observed line and continuum emission. As a step in this direction, we carried out energy-conserving radiation-MHD calculations of heating and cooling in a local shearing box. The starlight that dominates the overall heating is absorbed in an extended, magnetically-supported atmosphere. The surface of optical depth unity lies on average 1.5 times higher than in the corresponding hydrostatic model. The fraction of the stellar luminosity intercepted will be larger in proportion, with consequences for interpreting the spectral energy distributions of protostellar disks. The puffy, magnetically-supported atmospheres may account for the puzzling fact that many young stars have near-infrared excesses too large to explain with hydrostatic models, as reviewed by Dullemond & Monnier (2010). Greater dust mass fractions than our assumed value of 10 −4 would raise the opacities, pushing the absorbing surface yet higher, while also increasing the recombination cross-section, making the weakly-ionized dead zone thicker.
The accretion flow occurs in the disk atmosphere, not the interior as assumed in classical viscous models, and the dissipation occurs still higher (figure 1). The viscous model's midplane temperature bump vanishes in the MHD case: relocating the dissipation to the atmosphere leaves the interior cooler and isothermal (figure 2). The magnetic fields are generated in an optically-thin turbulent layer and rise buoyantly as they dissipate. While the total mass column is 1 000 g cm −2 , over half (56%) of the accretion heat is deposited at columns within 1 g cm −2 of the boundaries, and more than one-fifth (23%) within 0.1 g cm −2 . The turbulent layer is optically-thin to its own continuum emission independent of the dust depletion, so long as the grains control both the opacity and the recombination (Bai & Goodman 2009). Because the accretion power heats optically-thin gas, the resulting spectral lines can potentially be used to measure the mass flow rate. The dissipation is spatially inhomogeneous, heating gas locally above the irradiation temperature despite the starlight dominating the volume-integrated power. Even a small disk patch thus produces emission at a range of temperatures.
Since the disk atmosphere is magnetically-supported, its structure varies in response to changes in the strength of the fields. In particular, the height where the starlight is absorbed moves up and down by almost a factor of two over timescales of several orbits, with an amplitude around 10% of the distance to the star. The resulting changes in the stellar illumination of the disk surface could contribute to the widespread low-amplitude infrared variability among young stars with disks (Morales-Calderón et al. 2009;Luhman et al. 2010).
To make the calculations feasible, we simplified the situation in several respects. Among these, we capped the magnetic diffusivity well above the threshold for shutting off the MRI, but well below the value computed from the midplane ionization fraction. The weak midplane magnetic oscillations visible in figure 5 are an artifact of the reduced midplane diffusivity, and would disappear if the cap were lifted above the threshold where diffusion cancels the winding-up of radial fields. The dead zone stresses are unimportant however, as they contribute only 2.6% to the total mean mass accretion rate. Also, we treated the stellar irradiation in a horizontally-averaged fashion and kept the starlight incidence angle fixed despite the changes in the disk thickness, since the shearing-box approximation gives no information on the tilt of the disk surface. Future work could explore the effects of relaxing these simplifications. In addition, our current sheets are artificially broadened to the grid scale by the numerical resistivity. The sheets should be narrower, with still greater accretion heating per unit volume, and thus greater chances of exceeding the stellar irradiation. The magnetic dissipation furthermore heats the gas, which may become hotter than the dust where densities are low enough so heat transfer by gas-dust collisions is slow. Since we assumed good gas-dust thermal coupling, our peak temperatures are lower limits. Future modeling should follow Glassgold et al. (2004), Kamp & Dullemond (2004) and Nomura & Millar (2005) in allowing distinct dust and gas temperatures.