Toward realistic simulations of magneto-thermal winds from weakly-ionized protoplanetary disks
Protoplanetary disks (PPDs) accrete onto their central T Tauri star via magnetic stresses. When the effect of ambipolar diffusion (AD) is included, and in the presence of a vertical magnetic field, the disk remains laminar between 1-5au, and a magnetocentrifugal disk wind forms that provides an important mechanism for removing angular momentum. We present global MHD simulations of PPDs that include Ohmic resistivity and AD, where the time-dependent gas-phase electron and ion fractions are computed under FUV and X-ray ionization with a simplified recombination chemistry. To investigate whether the mass loading of the wind is potentially affected by the limited vertical extent of our existing simulations, we attempt to develop a model of a realistic disk atmosphere. To this end, by accounting for stellar irradiation and diffuse reprocessing of radiation, we aim at improving our models towards more realistic thermodynamic properties.
Interpreting observational properties of T Tauri systems  is closely tied to understanding the complex dynamical evolution of gaseous protoplanetary disks (PPDs), both in terms of their chemistry, and in terms of the microphysics that govern the evolution of the embedded magnetic fields. Moreover, with PPDs being the birth sites for extrasolar planets, a sound physical picture of the dust  and planetesimal  evolution is needed to ultimately provide the building blocks for a comprehensive theory of planet formation.
The fundamental drivers of disk evolution are mass loss processes and redistribution of angular momentum . In sufficiently ionized parts of the disk, the latter can be achieved by turbulent stresses. When accounting for the ionization structure of a typical PPD, large parts of the disk, however, remain laminar owing to the effect of ambipolar diffusion (AD). In this situation, angular momentum is primarily transported via a magnetocentrifugal disk wind [5, 6]. While this picture is further complicated when accounting for the Hall effect [7, 8], it has nevertheless become clear that the thermal structure of the disk plays an important role in setting the mass loading of the disk wind , and hence the timescale on which the system evolves [10, 11].
As in our previous work , we are solving the single-fluid MHD equations including Ohmic resistivity and ambipolar diffusion, that is, the electromotive force resulting from the mutual collision of ions and neutrals. The diffusion coefficients, , and are specified by means of a look-up table, which has been obtained by solving a simple chemical ionization/recombination network. The resulting electromotive forces stemming from the Ohmic and ambipolar diffusion terms are then given by
with being the unit vector along , and where the double vector product results in an additional minus sign, such that positive coefficients and signify diffusion of magnetic fields, with the latter only being sensitive to currents perpendicular to the field direction.
For the typical number densities in the context of PPDs, the frictional coupling between the ions and neutrals happens so quickly (compared to dynamical timescales of interest) that we assume the charged species move at their terminal velocity with respect to the neutrals, that is, where the Lorentz force is balanced by the drag. Because of the low degree of ionization, we formulate our continuity equation and momentum conservation in terms of the neutral component, although it does experience the Lorentz force as mediated by particle collisions.
We use a modified version of the nirvana-iii finite volume Godunov code [12, 13]. The code adopts a total-energy formulation111 Note that we do not include the radiation energy density, , in the total energy (as it is sub-dominant in the problem we consider). We however retain the radiation flux in the momentum equation for consistency. with conserved variables , , and . Together with the conservation of radiation energy, , and defining the total pressure, , as the sum of the gas and magnetic pressures, the system of equations we solve reads
where is the radiation flux, is the radiation pressure tensor, and are the Rosseland and Planck mean opacities, respectively, is the radiation density constant, is the gas temperature, and the gas pressure is , where is chosen as appropriate for an ideal diatomic gas. The gravitational term represents the point-mass potential of the solar-mass star at the center of our spherical-polar coordinate system, and represents external irradiation heating due to the star.
2.1 Flux limited diffusion
The above equations can only be solved once the radiation flux and pressure tensor are specified. An attractive method for obtaining the Eddington tensor that relates with is to solve for the steady-state but angle-dependent radiation intensity and integrate the first (for ) and second (for ) order moments directly . In the interest of maintaining minimum algorithmic complexity and computational expense, we instead adopt the classical ad hoc closure
that specifies in terms of a diffusive flux with diffusion coefficient , where
is a dimensionless number specifying how abruptly the radiation energy density varies compared to the length scale defined by the optical extinction coefficient , and where (for ) is a limiter function  that guarantees in regions of low optical depth, that is, in regions where the diffusion approximation is not valid. In the optically thick limit, (for ), which corresponds to the Eddington approximation.
The described approach has its known shortcomings over characteristics-based methods , but in combination with stellar irradiation heating (determined using a simplified ray tracing algorithm) it has been deemed an acceptable compromise [16, 17, 18] in terms of being able to incorporate radiation thermodynamics in fully dynamical 3D simulations, whereas more accurate Monte Carlo methods are comparatively expensive and have to deal with statistical noise .
We currently treat both irradiation and diffuse redistribution in the gray approximation, but multigroup approaches are straightforward , especially for the irradiation component [17, 18], wherein computational expenditure scales with the number of frequency bins, and where increased realism can be achieved for the thermal structure of the outer disk (). For our PPD model, we precompute look-up tables for mean opacities, and , using D. Semenov’s opacity.f , where a fixed dust-to-gas mass ratio of is assumed. To account for depletion of small dust by grain growth, we enable the reduction of the obtained opacities by a scale factor. Since the opacity tables combine contributions from dust grains and gas molecules, this is only valid for temperatures below the dust sublimation threshold.
2.2 Reduced speed of light approximation
The thermodynamic coupling of the gaseous matter with the radiation field is described by the terms in eqns. (5) and (6), respectively. Subsuming the term in eqn. (6) with the definition of the diffusive flux (8) amounts to a diffusion equation for the redistribution of radiation energy. Both effects can be combined into the subsystem
which we solve by means of operator splitting. Unlike in ref. , we do not include the term in this subsystem but instead treat it as a regular source term when updating the main hyperbolic system of equations. Since, for a large diffusion coefficient , the parabolic system (11) becomes stiff, the most common approach is to use implicit methods to solve it. Especially in view of applications employing adaptive mesh refinement, we have chosen to avoid an implicit update for , as it demands costly non-local communication patterns, and a potentially expensive matrix inversion.
Instead, we integrate the diffusion part of (11) in a time-explicit fashion, and use the reduced speed of light approach  to ameliorate the strict time-step constraints. This method has recently been employed in the context of simulations of the interstellar medium . The approximation is valid as long as the radiative timescales resulting from the adopted artificial value of (with const.) are still short compared to any relevant dynamical timescales. In the context of PPDs, we are mainly interested in the role of radiative effects in setting the consistent temperature structure of the disk, and true radiation hydrodynamic effects are believed to only be of minor importance during the T Tauri phase .
The approximation is introduced by amplifying the left-hand-side of (11) by a factor of . It is crucial that, because all other terms remain unaffected, this implies the modification does not alter the late-time steady-state solution, where , but only changes the timescale on which this solution is achieved. On practical grounds, we multiply (11) by , which implies replacing for in the radiation matter coupling, and attenuating the diffusion term by a factor of . This illustrates how the method works to weaken the stiffness of the diffusion term. To explicitly integrate (11), we employ the second-order accurate Runge-Kutta-Legendre (RKL2) super-time-stepping scheme , that is already used for the updates of the other parabolic terms (such as, viscosity, thermal conduction, resistive and ambipolar diffusion) in the code.
2.3 Radiation matter coupling
Even with the reduced speed, , the radiation-matter-coupling term can itself contribute a severe timestep constraint in regions of high opacity, where the coupling becomes stiff. Ignoring, for the moment, the diffusion term in (11), the coupling amounts to an ordinary differential equation, similar to production/destruction equations that are common in other fields of science. While these are typically solved by explicit Runge-Kutta (RK) methods, higher order predictor/corrector schemes do not typically guarantee positivity or conservation of, in our case, the energy .
There however exists a class of modified RK methods  that employ the so-called Patankar trick – an implicit weighting of the production/destruction terms with the ratio of the evolved quantity before and after the update. It can be shown that for a single-step update, such a weighting precludes that negative values are obtained. Moreover, since the weighting factors are overall symmetric, conservation is guaranteed. Specifically, we use the MPRK scheme given by eqn. (27) in ref. . The resulting update is formally implicit, but can algebraically be manipulated into fully explicit form, that is,
for the predictor step, and
for the corrector step, where is the forward-Euler predictor value of , and is the time-averaged state.222The same conventions, of course, apply for the terms , and . Compared to implicit methods, that demand iterative root-finding, or so-called -schemes (see, e.g., section 3.4 in ref. ) the method presented here offers a relatively inexpensive, parameter-free non-iterative alternative.
2.4 Irradiation heating
Our existing global disk models [27, 28] have either assumed a locally-isothermal temperature , with being the cylindrical radius from the star, or have used an adiabatic equation of state with a Newtonian cooling term in the energy equation that reinstated the profile on a specified timescale (typically a short fraction of the local orbital period).
Compared to this, even relatively basic models  of dust absorption and re-radiation of star light in the disk surface, obtain a much more complex temperature structure within the PPD – with superheated surface layers and cool interiors. To account for such effects, we include a frequency-integrated radially attenuated irradiation flux
from the central star with effective temperature , and , and where the optical depth, is obtained by integrating along radial rays. Following previous work [16, 17, 18], we obtain the irradiation energy source term by computing the negative divergence of (14) over each grid cell. In regions of low optical depth, we use the integral formulation
with instead, as this formulation has been found to produce a more accurate solution on the discretized mesh, when differences across cells are small .
3.1 Radiative transfer test problems
To verify our implementation of the radiation-matter-coupling term, we have performed a simple one-zone model , where the thermal energy density, , is initially out of balance with the radiation energy density . In Fig. 1, we show three cases with , , and , which all converge to the same final equilibrium state. For the purpose of plotting the curves, we have limited the timestep artificially to sample timescales shorter than the radiative equilibration timescale. We have, however, tested that even for numerical time steps somewhat larger than the coupling timescale the scheme remains stable, as expected from the implicit-like integration scheme (see sect. 2.3) that we use.
A standard test case for assessing the interplay of the radiation-matter-coupling with the radiation diffusion are radiative shocks, for which there exist semi-analytic solutions in simplified situations . In Fig. 2, we plot the solution of the case from ref. , using grid cells in the direction, as well as two levels of adaptive mesh refinement (triggered by gradients in the thermal energy, ), which are shown as gray shaded areas in the plot. As seen in the lower panel of Fig. 2, apart from the shock location (which has been shifted by with respect to the semi-analytic solution), all quantities agree to within accuracy. We have also successfully performed the harder test, which we omit here for brevity. An excellent description of the radiative shock test, including the precise values used for initial conditions, can be found in sect. 4.4 of ref. .
3.2 Preliminary global MHD simulations with irradiation
Returning to the original motivation for implementing radiative physics in our modified version of the nirvana-iii code, we conclude by presenting a preliminary snapshot of a global axisymmetric protoplanetary disk simulation including Ohmic resistivity, ambipolar diffusion, radiation diffusion, and stellar irradiation. The ultimate goal of these simulations is to study how the mass-loading of the magnetocentrifugal wind depends on the disk thermal structure.
As a simple proof-of-concept, we present a close-up of the inner disk in a simulation covering seven (initial) pressure scale heights in latitude (see Fig. 3). The basic disk setup is very similar to the one used in ref. , and we have additionally used opacities that correspond to a dust-depletion of a factor of ten compared to the typical interstellar abundance. Similar to our previous simulations, the magnetic field lines (white) bend outward, and in the upper disk layers, where the matter is sufficiently coupled to the magnetic field, a magnetocentrifugal disk wind ensues (black vectors).
Iso-contour lines (gray) of the radiation temperature illustrate the disk’s thermal structure that deviates noticeably from the constant-on-cylinders radial temperature profile , that we have assumed previously. The presented preliminary run used a reduced-speed-of-light factor . Further tests will show whether the chosen time-explicit framework is powerful and efficient enough to be of use in realistic situations.
The author thanks Jon Ramsey for many useful discussions, for carefully reading this manuscript, and for providing his code to compute semi-analytic reference solutions of radiative shocks. Dmitry Semenov is acknowledged for providing helpful clues in connection with his opacity code. The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 638596).
-  Williams J P and Cieza L A 2011 ARA&A 49 67–117 (Preprint 1103.0556)
-  Johansen A, Blum J, Tanaka H, Ormel C, Bizzarro M and Rickman H 2014 PPVI 547
-  Gressel O, Nelson R P and Turner N J 2012 MNRAS 422 1140–1159 (Preprint 1202.0771)
-  Turner N J, Fromang S, Gammie C, Klahr H, Lesur G, Wardle M and Bai X N 2014 PPVI 411–432
-  Bai X N and Stone J M 2013 ApJ 769 76 (Preprint 1301.0318)
-  Gressel O, Turner N J, Nelson R P and McNally C P 2015 ApJ 801 84 (Preprint 1501.05431)
-  Lesur G, Kunz M W and Fromang S 2014 A&A 566 A56 (Preprint 1402.4133)
-  Simon J B, Lesur G, Kunz M W and Armitage P J 2015 MNRAS 454 1117–1131 (Preprint 1508.00904)
-  Bai X N, Ye J, Goodman J and Yuan F 2016 ApJ 818 152 (Preprint 1511.06769)
-  Bai X N 2016 ApJ 821 80 (Preprint 1603.00484)
-  Suzuki T K, Ogihara M, Morbidelli A, Crida A and Guillot T 2016 A&A (Preprint 1609.00437)
-  Ziegler U 2004 JCoPh 196 393–416
-  Ziegler U 2011 JCoPh 230 1035–1063
-  Jiang Y F, Stone J M and Davis S W 2012 ApJS 199 14 (Preprint 1201.2223)
-  Levermore C D and Pomraning G C 1981 ApJ 248 321–334
-  Bitsch B, Crida A, Morbidelli A, Kley W and Dobbs-Dixon I 2013 A&A 549 A124 (Preprint 1211.6345)
-  Kuiper R and Klessen R S 2013 A&A 555 A7 (Preprint 1305.2197)
-  Ramsey J P and Dullemond C P 2015 A&A 574 A81 (Preprint 1409.3011)
-  Noebauer U, Sim S, Kromer M, Röpke F and Hillebrandt W 2012 MNRAS 425 1430 (Preprint 1206.6263)
-  González M, Vaytet N, Commerçon B and Masson J 2015 A&A 578 A12 (Preprint 1504.01894)
-  Semenov D, Henning T, Helling C, Ilgner M and Sedlmayr E 2003 A&A 410 611–621
-  Gnedin N Y and Abel T 2001 NewA 6 437–455 (Preprint astro-ph/0106278)
-  Skinner M A and Ostriker E C 2013 ApJS 206 21 (Preprint 1306.0010)
-  Hartmann L 1998 Accretion Processes in Star Formation
-  Meyer C D, Balsara D S and Aslam T D 2012 MNRAS 422 2102–2115
-  Burchard H, Deleersnijder E and Meister A 2003 AppliedNumMath 47 1 – 30 ISSN 0168-9274
-  Nelson R P, Gressel O and Umurhan O M 2013 MNRAS 435 2610–2632 (Preprint 1209.2753)
-  Gressel O, Nelson R P, Turner N J and Ziegler U 2013 ApJ 779 59 (Preprint 1309.2871)
-  Chiang E I and Goldreich P 1997 ApJ 490 368–376 (Preprint astro-ph/9706042)
-  Bruls J H M J, Vollmöller P and Schüssler M 1999 A&A 348 233–248
-  Turner N J and Stone J M 2001 ApJS 135 95–107 (Preprint astro-ph/0102145)
-  Lowrie R B and Edwards J D 2008 Shock Waves 18 129–143