Radiation thermo-chemical models of protoplanetary disks. Grain and polycyclic aromatic hydrocarbon charging

Context. Disks around pre-main-sequence stars evolve over time by turbulent viscous spreading. The main contender to explain the strength of the turbulence is the Magneto-Rotational-Instability (MRI) model, whose efficiency depends on the disk ionization fraction. Aims. We aim at computing self-consistently the chemistry including PAH charge chemistry, the grain charging and an estimate of an effective value of the turbulence alpha parameter in order to find observational signatures of disk turbulence. Methods. We introduced PAH and grain charging physics and their interplay with other gas-phase reactions in the physico-chemical code ProDiMo. Non-ideal magnetohydrodynamics effects such as Ohmic and ambipolar diffusion are parametrized to derive an effective value for the turbulent parameter alpha_eff . We explored the effects of turbulence heating and line broadening on CO isotopologue sub-millimeter lines. Results. The spatial distribution of alpha_eff depends on various unconstrained disk parameters such as the magnetic parameter beta_mag or the cosmic ray density distribution inside the protoplanetary disks. The inner disk midplane shows the presence of the so-called dead-zone where the turbulence is quasi-inexistent. The disk is heated mostly by thermal accommodation on dust grains in the dead-zone, by viscous heating outside the dead-zone up to a few hundred astronomical units, and by chemical heating in the outer disk. The CO rotational lines probe the warm molecular disk layers where the turbulence is at its maximum. However, the effect of turbulence on the CO line profiles is minimal and difficult to distinguish from the thermal broadening. Conclusions. Viscous heating of the gas in the disk midplane outside the dead-zone is efficient. The determination of alpha from CO rotational line observations alone is challenging.


Introduction
Pre-main-sequence stars (TTauri and HerbigAe stars) are surrounded by planet-forming disks in a Keplerian rotation (Williams & Cieza 2011;Espaillat et al. 2014). The disks are massive at the early stage reaching gas masses of 10 −2 M or even higher and disappear after a few million years (Alexander et al. 2014). The disk material can evaporate, form giant planets, or fall onto the star. Gas can accrete from the inner rim of disks to the star at a rateṀ of 10 −9 -10 −7 M yr −1 (Calvet et al. 2000;Muzerolle et al. 2004;Johns-Krull et al. 2000;Hartmann et al. 1998), although episodic accretion rate of up to 10 −5 M yr −1 can occur for young disks (Audard et al. 2014). Accretion shocks onto the stellar surfaces result in ultraviolet excess emission. Gas accretion is also traced by optical emission lines (Mendigutía et al. 2012(Mendigutía et al. , 2011Garcia Lopez et al. 2006).
Magneto-rotational instability is an ideal magnetohydrodynamics (MHD) phenomenon that is effective only if the gas is sufficiently ionized to couple dynamically to the magnetic fields (Gammie 1996;Bai 2015). Gas turbulence efficiency is characterized by the factor α so that the large-scale turbulence is ν = αc s h, where c s is the gas sound speed and h is disk pressure scale-height representing the fluid typical scale. For subsonic turbulence, the value of α should be lower than unity, with an usually assumed value of 0.01.
Many non-ideal Magnetohydrodynamics (MHD) dissipation effects can restrict the development of MRI turbulence or even suppress it (Jin 1996). The Ohmic and ambipolar diffusion depend on the abundances of the charge carriers among other factors. Many MHD simulations were carried out to study the extent of "dead zones" defined as regions where dissipation overcomes the MRI turbulence (Wardle & Ng 1999;Sano & Stone 2002;Bai & Stone 2011;Gressel et al. 2015). On the other hand, Hall diffusion revives the MRI under certain conditions (MRI-Hall; Lesur et al. 2014;Wardle & Salmeron 2012;Sano & Stone 2002;Balbus & Terquem 2001).
The charge carriers in disks are either electrons, atomic and molecular ions, polycyclic aromatic hydrocarbons (PAHs) or dust grains (Bai & Stone 2011). The abundances of the carriers vary throughout the disk due to ionization processes (photoionization by UV and X-ray photons, cosmic rays), recombination with an electron, and charge exchanges. The rates of these processes are functions of the gas composition, the gas and dust temperatures, the UV, X-ray, and density fluxes inside the disks (Rab et al. 2017). Therefore, the computations of the abundances require a detailed gas and dust physical and chemical modelling together with the continuum and line radiativetransfer. Bai (2011) and Perez-Becker & Chiang (2011) explored the role of grains in the efficiency of MRI but did not consider the effects of radiation on the grain charging and considered PAHs as small grains and not as macro-molecules.
Ionization fraction in realistic protoplanetary disks was modelled with different levels of sophistication. The aim of those studies was to determine the extent of the dead-zone, the area in disks where the ionization fraction is too low for the magnetic field to couple to the neutral gas, and hence for MRI to sustain turbulence. No observational evidence exists either for or against the presence of a dead-zone yet. The dead-zone encompasses the region of planet formation in disks. Interestingly one of the questions in planet-formation is the influence of low turbulence on the growth of planetesimals. Ilgner & Nelson (2006a) considered viscous heating and radiative cooling and tested the effects of different chemical networks on the ionization fraction in disks. The source of ionization is stellar X-rays. Subsequently, Ilgner & Nelson (2006c) studied the effects caused by X-ray flares. The effect of turbulent mixing on gas chemistry was explored in Ilgner & Nelson (2006b); Heinzeller et al. (2011). Fromang et al. (2002) modelled the ionization fraction and the size of the dead-zone in disks. Perez-Becker & Chiang (2011) focused on the UV ionized layers of disks and argued that this layer cannot sustain the accretion rate generated at larger radii. Dzyurkevich et al. (2013) studied the dependence of disk parameters such as temperature, surface density profile, and gas-to-dust mass ratio on the MRI efficiency. The sources of ionisations in the inner disks are discussed by Desch & Turner (2015). They found that the ambipolar diffusion controls the location of the dead zone. Ivlev et al. (2016) considered the ionization and the effect of dust charging in protoplanetary disks. The role played by dust grains in the ionization fraction in disks has also been discussed in Salmeron & Wardle (2008). Simon et al. (2011) performed non-ideal MHD simulations, predicting turbulence velocities in the disk midplane ranging from 0.01 times the sound speed in the dead-zone to 0.1 outside.
Observational constraints on the turbulence factor α in protoplanetary disks have been obtained for the TW Hya (Teague et al. 2016;Flaherty et al. 2018) and HD 163296 (Flaherty et al. , 2017 disks by determining the contribution of the turbulence broadening to the total line width using ALMA data. Teague et al. (2016) found a turbulence value of (0.2 − 0.4) c s in the TW Hya disk, while Flaherty et al. (2018) constrained α to be between 0.04 c s and 0.13 c s . Flaherty et al. (2015) provided a low upper limit of v turb < 0.03 c s for the turbulence in the HD 163296 disk. Earlier studies with the IRAM Plateau de Bure interferometer provided limits of v turb ≤ (0.3 − 0.5) c s (Dartois et al. 2003;Piétu et al. 2007). Hughes et al. (2011) constrained the value v turb ≤ 0.1 c s for the TW Hya disk and v turb ≤ 0.4 c s for the HD 163296 disk using data obtained by the Smithsonian Millimeter Array. The accuracy of the estimated values is limited by the uncertainties in inverting the protoplanetary disk thermal structure from the observations. Hartmann & Bae (2018) argue that a low viscosity disk model driven by hydrodynamic turbulence is compatible with the observed disk mass accretion rates and the measured low turbulence width.
In this paper, we explore further the effects of a detailed treatment of the physics and chemistry of PAHs and grain charging on the disk ionization. For this purpose, we implemented MRIturbulence heating and cooling in the photo-chemical disk code ProDiMo. We considered far-ultraviolet from the star and the accretion excess, X-ray, and cosmic rays as sources of ionization. We used CO isotopologue rotational lines as potential tracers of turbulence in disks because CO is the most abundant molecule in disks after H 2 , is widespread through the disks, and because its chemistry is well understood. The numerical study of MHD processes in triggering disk turbulence is beyond the scope of this paper. Instead, we parametrize the onset and effects of turbulence on the gas by simple empirical formula. This study focuses on the role of the gas and grain chemistry in the disk's ionization, which in turns can affect the non-ideal MHD coefficients. The scope of the ProDiMo code is to be able to derive disk parameters such as α by matching high signal-to-noise observations carefully considering the relevant disk physico-chemical processes.
The paper is organized as follows: the ProDiMo code is introduced in Sect. 2.1; a brief discussion on the gas-phase and PAH chemistry are given in Sect. 2.2 and 2.3; the dust charging physics implemented in the code is described in Sect. 2.4; the prescription of our MRI-driven turbulence model is provided in Sect. 3; analytical results on the influence of non-ideal MHD are shown in Sect. 4; the protoplanetary disk model and the model results are presented in Sect. 5 and Sect. 6; our findings are discussed in Sect. 7 and we conclude in Sect. 8.

ProDiMo
ProDiMo is a code built to model the gas and dust grain physics and chemistry (Woitke et al. 2009;Kamp et al. 2010;Woitke et al. 2016). It has been used to model disk Spectral Energy Distributions (SEDs, Thi et al. 2011), water deuteration chemistry ) CO rovibrational emissions including UVfluorescence (Thi et al. 2013), and many Herschel observations from the GASPS large programme. X-ray physics are implemented Meijerink et al. 2012;Aresu et al. 2011). ProDIMo has been designed to run within a few CPUhours per model such that automatic fittings of the observed continuum emission and of thousands of gas lines are feasible. Such fitting procedure requires running thousands of models, which would be too time-consuming with a full 3D non-ideal MHD radiation chemical code. As such, ProDiMo does not solve the gas hydrodynamic equations. The disk density structures are parametrized. Only the vertical hydrostatic structure can be selfconsistently re-adjusted with the gas temperature. In this study we chose to model disks with fixed vertical hydrostatic structures. A detailed discussion of the different physics and their implementations are given in the articles listed above. Here we summarize the main features relevant to the modelling of the MRI-turbulence.
In our chemical modelling, we included gas and ice species as well as PAHs. The grain charge is computed for grains of mean radius a . The grains can be up to four times positively or negatively charged. The photoionization and photodissociation rates are computed from the cross-sections and UV field calculated from 2D continuum radiative transfer (van Dishoeck & Visser 2011). Self-shielding is taken into account. The gas temperature at each location in the disk is computed by balancing the heating rate with the cooling rate contrary to the previous study of MRI in disks where the gas temperature follows a power-law distribution in radius T ∝ r −p and the disk is isothermal in the vertical direction or considered viscous heating only (Hughes et al. 2011;Flaherty et al. 2015;Teague et al. 2016).
In ProDiMo heating agents include photoelectrons ejected from PAHs (Bakes & Tielens 1994) and dust grains (photoelectric effects), chemical heating, photoionization heating, and in this study release of viscous energy. Atomic and molecular lines can heat or cool the gas by absorption or emission. At high densities, the gas and the dust grains exchange energy by thermal contact. The thermal accommodation can heat or cool the gas depending on the sign of ∆T =T gas -T dust (Burke & Hollenbach 1983). If ∆T is positive, the gas-dust accommodation will cool the gas and vice-versa. The codes generates detailed line fluxes and profiles from non-LTE radiative transfer, that can be compared directly to observations.
Knowledge of the precise ionization fraction at each location in the disks is paramount for the onset of MRI-driven turbulence. A few additional features have been implemented to improve the physics and chemistry that regulate the charge distribution in protoplanetary disks. Most of the features concern a better treatment of the charging of the PAHs and dust grains.

Gas phase chemistry
The gas-phase chemistry includes photodissociation, ionneutral, neutral-neutral, as well as a few colliders and three-body reactions. The chemical network is discussed in details in Kamp et al. (2017a). At high gas temperatures, ionization by collisions with hydrogen atoms or molecules and with electrons is possible.
with the rates provided by Hollenbach & McKee (1980). For the thermal ionization by hydrogen, Hollenbach & McKee (1980) suggest using the rate for electrons scaled by a factor 1.7×10 −4 . Collisional ionization of metals can also occur where M is a neutral metal (Fe, Mg, S, Na, Si, ...) and N is either H, H 2 , or He.

PAH charge exchange chemistry
PAHs can become the main charge carriers in disks because of their abundances and electron affinity. Details on the PAH chemistry can be found in Kamp et al. (2017a) and Thi et al., submitted. In protoplanetary disks, their abundances are lower by a factor f PAH =10 −1 -10 −3 compared to their interstellar abundance of 3×10 −7 (Tielens 2008). We chose to use the pericondensed circumcoronene (C 54 H 18 ) as typical PAHs that are large enough to remain unaffected by photodissociation in disks around HerbigAe stars (Visser et al. 2007). The circumcoronene can be once negatively-charged (PAH − ) and three times positively charged by absorbing a UV photon with energy below 13.6 eV or by charge exchange reactions (PAH + , PAH 2+ , PAH 3+ , see Table 1). The effective radius of a PAH is computed by (Weingartner & Draine 2001b) a PAH = 10 −7 N C 468 where N C is the number of carbon atoms in the PAH. The radius for the circumcoronene is a PAH (C 54 H 18 )=4.686 × 10 −8 cm.
The PAH ionization potential can either be taken from the literature when they are measured or estimated (Weingartner & Draine 2001b) where W 0 is the work function assumed to be 4.4 eV (7.05×10 −12 erg), and Z PAH is the charge of the PAH. The ionization potentials (I.P.) for circumcoronene are listed in Table 1.
PAHs are not formed or destroyed in our chemical network and only exchange charges with other positively-charged species (for examples H + , He + , Mg + , Fe + , Si + , S + , HCO + , ...). Chemical reaction rates involving PAHs are highly uncertain. Most of the rates are extrapolations from a few existing laboratory or theoretical rates. PAH freeze-out is presented in Kamp et al. (2017a).
On disk surfaces, PAH are mostly positively charged because of the photoelectric effect. We compute the PAH photoejection rate using detailed PAH cross-sections and disk UV fields computed by radiative transfer (Woitke et al. 2016). A free electron can also attach on PAHs. A detailed discussion on PAH charging is provided in Thi et al., submitted.

Dust grain charging
Dust grains with radius can be major charge carriers in the interstellar medium in general and in protoplanetary disks in particular. Dust grain charging is studied in the field of plasma physics (Mishra & Misra 2015).
We implemented the silicate dust grain charging physics of Draine & Sutin (1987), Weingartner & Draine (2001b), and Umebayashi & Nakano (1980) with a couple of differences. We considered an average grain of radius a ≡ a and a geometric cross-section σ dust = πa 2 . The charge of the grain is Z d with a discrete charge distribution function f (Z d ) with a minimum charge Z min and maximum charge Z max , such that Z max Z min f (Z d ) = 1.

Silicate dust grain ionization potential
Detailed quantum mechanical calculations of the work function of oxide and silicate clusters show variations for different silicate compositions (Rapp et al. 2012). They noticed that silicates have an increase in the work function by 2 eV compared to the oxides (MgO, FeO) and attributed this effect to the presence of the silicon. When a grain is covered by a thick water ice mantle, the work function will be modified. Water ice has a work function of 8.7 eV (Rapp 2009;Baron et al. 1978), comparable to that of the magnesium-rich silicate. Therefore, the work function W 0 for bulk silicate is assumed to be typically 8 eV whether the grain is coated by an icy mantle or not. When a grain is positively charged, the Coulomb attraction between the grain and an electron effectively increases the effective work function by an extra work function that depends on the grain charge Z d . The effective work function W eff , equivalent to the ionization potential for an atom or molecule, becomes (Weingartner & Draine 2001b) where W c is the extra work function defined by where e is an elementary charge.When the grain is positivelycharged (Z d > 1). The value of the constant added to the grain charge is the correction due to the finite size of a perfectly spherical grain and is controversial. A value of 3/8 has also been proposed instead of 1/2 (Wong et al. 2003).

Photoejection and photodetachment
An electron can be ejected upon absorption of a photon with energy Θ higher than the ionization potential of the grain. The photoejection process concerns positively-charged and neutral grains, while the photodetachment process concerns negativelycharged grains. The photodetachment cross-sections for negatively-charged grains and for neutral and positively-charged grains respectively are calculated following the prescription of Weingartner & Draine (2001b).
The photoejection yield η b for the silicate bulk varies as function of the photon energy Θ following equation 17 of Weingartner & Draine (2001b) where W 0 is the silicate bulk photoejection yield assumed to be 8 eV, see also (Kimura 2016). The bulk yield is enhanced for very small grains by a factor η small according to equation 13 of Weingartner & Draine (2001b). The effective yield becomes η eff = η b η small . For large grains (x = 2πa/λ >5) we used an analytical fit to the experimental data of Abbas et al. (2006) for the effective yield where x 0 = 2.5, λ 0 =120 nm. The formula reproduces the experimental data within factor of a few. The photoelectric yield at three wavelengths for different grain sizes is shown in Fig. 1. The energy threshold for photoejection of a grain of charge Z d reads hν pe = W eff .
When a grain is negatively charged the photodetachment threshold energy is (Weingartner & Draine 2001b)  where the electron affinity is The band gap E bg = 0 for metals and semimetals while it can be several eVs for other materials. Weingartner & Draine (2001b) assumed (W 0 − E bg ) = 3 eV for a silicate grain (a band gap of 5 eV). The value for water ice is much lower at 0.8 eV (do Couto et al. 2006 For a singly negatively-charged grain, Z d =-1 and a = 1 µm, hν pd (Z d , a) 2.9 eV or λ pd 0.4 µm for a pure silicate grains and hν pd (Z d , a) 0.7 eV or λ pd 1.66 µm for pure water ice grains in the near infrared. The threshold energy for photodetachment is much lower than for photoejection. For disks, the dust extinction is lower at longer wavelengths and at the same time the stellar luminosity is higher in the blue. Both effects combine to make photodetachment very efficient. The combined rate for photoejection and photodetachment is where Q abs is the frequency-dependent absorption efficiency, J ν is the specific mean intensity at the frequency ν computed in by the continuum radiative transfer.

Electron attachment
The electron attachment rate coefficient onto neutral grains is k e,n = n e S e 8kT e πm e σ dust [electrons s −1 ], where T e is the electron temperature assumed equal to the gas kinetic temperature, S e is the sticking coefficient for electron attachment, assumed to be 0.5 (50% of the encounters are assumed elastic), σ dust is the geometrical cross-section of a grain of radius a, and m e is the mass of an electron. Umebayashi & Nakano (1980) discussed theoretically the sticking coefficient of electrons on grains. They found that the sticking probability depends on the surface composition and is in most cases larger than 0.3. For positively-charged grains (Z d > 0) the electron recombination rate is enhanced by Coulomb attraction where W c is the extra work function (see eq. 6). On the contrary, the recombination rate is lower if the grain is negatively charged

Thermionic emission
We considered the thermionic emission of electrons by hot dust grains following the Richardson-Dushman theory (Ashcroft & Mermin 1976;Sodha 2014) where T d is the dust temperature, r is the so-called reflection parameter. It corresponds to the fraction of electrons that have enough energy to escape at grain surfaces but do not do so. It is assumed to have a value of 0.5. The thermionic emission work W therm is equal to W eff for positively-charged grains and to EA + E min for negatively-charged grains.

Collisional electron detachment
An atomic hydrogen or molecular hydrogen impinging onto the grain can either transfer momentum or energy to the grain surface such that an electron is ejected.
where the work W cd is equal to W eff for positively-charged grains and to EA(Z d + 1) + E min for negatively-charged grains.

Charge exchange between ions and dust grains
Gas-phase species and dust grains can exchange charges. Weingartner & Draine (2001a) considered the effect of ioncharged grains on the ionization of the interstellar gas. Cations can recombine with an electron from the grains. Likewise, a positively-charged grain can capture an electron from atoms and molecules with ionization potential lower than that of the grain. For negatively-charged grains, the non-dissociative exchange reaction can be written as a suite of reactions (n ≥ 0) gr n− + X g + → gr (n−1)− + X s + ∆E ads (X g ) + ∆E Coulomb (19) gr n− + X s + → gr (n−1)− + X s + ∆E rec + ∆E relax where X + g and X + s are gas-phase and surface ions respectively, or gr n− + X s + → gr (n−1)− + X g + ∆E rec − ∆E ads (X g ).
The gas-phase species X + g acquires the approach energy that is the sum of the adsorption energy (∆E ads (X g )) and the Coulomb energy ∆E Coulomb . ∆E rec is the energy released/required by the recombination reaction and corresponds to the ionization potential minus the electron affinity for negatively-charged grains or the work function of the solid for neutral grains. The excess energy can be radiated away or transfer to the surface (∆E relax < 0). It can also be used to break the weak bond (∆E ads (X s ) ∼ 0.1-0.2 eV) between the species X s and the surface. Therefore, we assumed that the neutral atom immediately leaves the grain surface after the recombination. For molecular ion recombination, extra outcomes are possible. The recombination results first into an excited species When the excess energy of the electronically excited neutral species AH * cannot transfer efficiently to the grain surface or radiated away (because the species has a low dipole moment), the recombination is dissociative, where the molecule dissociates into species that leave immediately the surface carrying the excess energy as kinetic energy; or one of the products remain on the grain surface AH * s → A s + H g + ∆E sdiss,rec .
Here the heavier product A tends to easily transfer the excess energy to the surface or possesses many degrees of freedom and, therefore, can radiate efficiently the excess energy away. On the other hand, the inefficient transfer of energy from the light hydrogen atom to a heavy surface element results in the atom leaving the surface.
Cations can also exchange their positive charge with a neutral or positively-charged grain with charge +n with (n ≥ 0) gr n+ + X g + → gr (n+1)+ + X g .
For large silicate grains, Aikawa et al. (1999) performed classical computation for the recombination of HCO + with negativelycharged grains and concluded that the recombination results in the dissociation of the molecules. Another example is given by Both recombination and charge exchange reactions proceed with the rate where E therm is an energy equal to the endothermicity of the reaction. For an exothermic charge exchange, the energy is null (E therm =0). The ion temperature T ion is equal to the gas thermal temperature T gas and we assumed S ion =1, similar to Weingartner & Draine (2001b). The term in parentheses is positive for negatively-charged grains, enhancing cation recombinations. At the same time, it ensures that cation exchanges with highly-positive grains are prevented due to a repulsive potential. The energetic barrier is defined by for neutral and positively-charged grains. We assume that the adsorption and Coulomb energy are negligible. For negativelycharged grains where IP neu is the ionization energy of the neutral species. Again for Z d =-1 and a = 1 µm, EA(Z d + 1, a) 2.99 eV and E min (Z d < 0, a) 2.9 eV, thus The reality is probably much more complex. The Coulomb interaction acquired by the approaching ion can lead to the tunnelling of a surface electron through the grain work function, resulting in a recombination in the gas-phase (Tielens & Allamandola 1987). Positively-charged grains (n >1) can transfer their charge to a neutral gas-phase species, The positive grains can induce a field that polarizes the neutral species. Since the grain has a large geometrical cross-section, the rate is assumed to be the largest value between a Langevin-type rate and that of a gas impinging on a neutral grain where the reduced mass is basically that of the gas-phase species µ = m n and the polarizability is α pol = 10 −24 where the barrier term becomes in practically all cases, the geometrical rate dominates, Anions can also react with positively-charged grains However, we did not have anions in the chemical network at this current stage. are endothermic, while reaction with HCO + is possibly exothermic.
Charge exchange reactions will result in the positive charges being carried by the species with the lowest ionization potential, whose value is lower than the energy required to detach an electron from a neutral grain. The cations with low ionization potential will preferably recombine with free electrons and not with negatively-charged grains because of the velocity and also on the extra work required to remove an electron from a grain (∼2.7 eV).

Minimum and maximum grain charges
The maximum positive charge that a grain can acquire is determined by the highest energy of the photons if the dominant ionization process is photoejection. Assuming a grain radius a in micron and W 0 =8 eV and a maximum photon energy of hν max =13.6 eV, an elementary charge e in c.g.s. of 4.803206815 × 10 −10 , and the conversion 1 eV = 1.60217657 × 10 −12 erg, the maximum charge, assuming that photo-ejection is the dominant electron ejection process, reads For the minimum grain charge, we adopted again the formalism of Weingartner & Draine (2001b) where U ait V 3.9 + 1200(a/µm) + 0.0002(µm/a) for carbonaceous, 2.5 + 700(a/µm) + 0.0008(µm/a) for silicate (38) For a 1 micron radius silicate grain, U ait 702 V and Z min −487499. This large value corresponds to the maximum negative charges on a silicate grain of radius a without considering any electron ejection process including charge exchange processes. In practice, the maximum amount of negative charges on a grain is limited by the balance between electron recombination 6 and electron emission. Assuming a balance between electron attachment and thermionic emission in the absence of UV photons (k RD (Z d < 1) = k e,− ), we can derive an upper limit to the amount of negative charges on a grain from formula 14 and 17. where Assuming that S e = 1, r=0.5, and T d = T g = T , and For dense regions with n e = 10 12 cm −3 and T = 100 K, the term B is greater than unity. The term kT ln B 2.37 × 10 −5 T ln (5.27 × 10 −8 T −3/2 n e ) eV is always negligible compared to W 0 (=8 eV), thus a strong lower limit to the grain charge can be derived Assuming elementary charge in c.g.s of e = 4.80321 × 10 −10 statC and W 0 =8 eV, E bg =5 eV, and f ∼1, we obtain In this study, we did not include photoejection due to X-ray photons, and we adopted in our numerical models the limits

Dust charge estimates for dense UV-obscured regions
In disk regions where photoemission can be neglected (even photodetachment requires photons in the optical blue domain), one can estimate the grain charging at equilibrium by balancing the electron recombination with the charge exchange between the ions and the grains. Here we also neglect the effects of UV photons created by interaction of cosmic rays with the gas. The grain charge is the solution of the non-linear transcendental equation assuming equal sticking coefficients S e = S i (Evans 1994) where we can set x = Z d e 2 /akT . Using the medium global neutrality n elec = n ion + Z d n d , the equation becomes A representation of the ratio between charges on grain surfaces and free electrons in an obscured region is shown in Fig. 2 for a 1 micron radius grain. Grains are negatively charged because of the difference in velocity between the electrons and the ions. Assuming n ion ∼ n elec an approximate solution to the equation is Ratio between the amount of negative charges on grains and free electrons for a grain of radius 1 micron and different total ionization fraction (from 10 −14 to 10 −8 ) as function of the gas kinetic temperature.
where n nucl is the average number of nucleons of the cations. If the main ion is a molecular ion with 19 nucleons (HCO + ), the approximation becomes The estimated negative charge is much lower than the limit set by the balance between electron recombination and thermionic emission. The approximation is valid when the charge carriers are dominated by the free electrons and by the ions. Interestingly, the dust negative charging increases with the gas temperature.
Using the limit on the number of negative charges on a dust grain (see formula 44), we can derive the maximum fractional abundance of negative charges on grains One can also use the estimate of the dust charge in obscured disk regions (equation 49). The estimated fractional abundance of negative charges on grains becomes This last equation together with Fig. 2 shows that, unless the ionization fraction is below 10 −13 , the negative charges on large silicate grains (with a radius greater than one micron) are negligible in UV-obscured regions. Dust charging is stochastic by nature and dust grains exhibit a distribution of charge states. In addition, the charge on a given grain can fluctuate over time (Piel 2010;Cui & Goree 1994;Matthews et al. 2013). Extensive numerical studies suggest that the dust charge distribution about its equilibrium value Z d when photons are absent has a standard deviation of δZ d = 0.5|Z d | 1/2 . The simultaneous existence of positively and negatively charged grains is possible for small grains of 0.1 µm in radius or smaller at 10 K.
The numerical implementation of the simultaneous computation of the grain charge distribution, the PAH restricted chemistry, and the gas-phase chemistry is explained in more details in Appendix C. When photon-induced electron ejection is present, the number of negative charges on grains will be even lower.

MRI-driven turbulence prescription for hydrostatic protoplanetary disk models.
In this section, we describe our model non-ideal MHD driven turbulence gas heating and cooling and line broadening.

Ideal MHD value of α
The ideal-MHD MRI value for the turbulence parameter α, which we call α ideal , can be evaluated from the results of detailed local 3D-MHD simulations. These simulations have shown that α ideal can be related to β mag , which is the ratio of the thermal P therm (r, z) to magnetic pressures P mag (r, z), by the function where δ is a parameter between 0.5 and 1. The β mag (beta magnetic) parameter is related to other gas parameters where E th and E mag are the thermal and magnetic energy respectively, ρ is the gas mass density, c s = kT/µ n is the sound speed with mean molecular mass µ n = 2.2 amu (atomic mass unit), B z (r) is the vertical component of the magnetic field, and v A is the Alfvèn speed in the disk vertical direction defined by The values for α ideal are bound between 2×10 −3 and 0.1 for β mag between 1 and 10 6 . A high value of β mag means that thermal motions in the plasma is important while β mag < 1 means that the dynamic is dictated by the magnetic field. Alternatively, we can use the prescription of Armitage et al. (2013) log where the value for the constants are: A = −1.9, B = 0.57, C = 4.2 and D = 0.5. Both functionals are plotted in Fig. 3 with δ = 0.5. We assume that β mid does not vary with radius in the midplane. This is equivalent to assume that the magnetic field is accreted with the gas (Guilet & Ogilvie 2014) β mag (r, z) = β mid P therm (r, z) P therm (r, 0) .
where the constant β mid is a free parameter of the disk models. The vertical component of the magnetic field varies with radius but constant in the vertical direction (see Fig. E.1). In both formulations, we set α ideal , the ideal-MHD MRI value of α in the absence of resistivities, to a very low value (10 −10 ), if β mag ≤ 1, i.e. when the flow is directed by the magnetic field. Since the gas thermal pressure decreases with height, β mag will also decrease. The ideal MRI value of α ideal will thus increase together with the disk height until β mag reaches 1, whose location can be considered as the base of a MHD-driven wind. The different plasma resistivities in the case of non-ideal MHD will decrease the value of α ideal to an effective value α eff . In the subsequent sections, we will describe the different non-MHD resistivities before presenting the prescription for the value of α eff that is used in the chemical disk models.

Non-ideal MHD resistivities
In the absence of strong X-ray radiation, the gas in protoplanetary disks is mostly neutral with a maximum ionization fraction reaching a few 10 −4 when all the carbon is ionized in the upper atmospheres (weakly-charged plasma). Hydrogen is predominantly molecular. The ionization fraction at disk surfaces can be much high when X-ray ionization is included, such that a significant fraction of hydrogens is in the form of protons. We considered a fluid that is globally neutral with velocity v(r) and density ρ(r). A charged species (electron, atom, molecule, PAH, or dust grain) j has mass m j , charge Z j , number density n j , and drift velocity relative to neutral v j . The current density J including all the charge-bearing species j (electrons, ions, PAHs, dust grains) with charge Z j is For a weakly ionized gas, the Lorentz force and the drag with the neutral gas force dominate over the inertia, gas pressure and gravitational forces in the equation of motion for the charge species where is the drag coefficient (cm 3 g −1 s −1 ) with σv j being the average collision rate between the charge species and the neutral gas of average mass m. The drag forces measure the rate of momentum exchange via the collisions. Using the gas global neutrality j n j Z j = 0 we obtain The inversion of equation 58 leads to the generalized Ohm's law, which characterizes the non-ideal MHD effects (Wardle & Ng 1999;Norman & Heyvaerts 1985) where E and E ⊥ denote the decomposition of E into vectors parallel and perpendicular to the magnetic field B respectively, and theˆmeans unit vector. σ O is the Ohm conductivity, also referred as the conductivity parallel to the field (σ O = σ ), σ H is the Hall conductivity, and σ P is the Pedersen conductivity. In cgs units, the conductivities have units of s −1 . The knowledge of the conductivities (or equivalently the resistivities) is central to the efficiency of MRI-driven turbulence in disks. The conductivities depend on the charge carriers and hence on the disk chemistry. Conversely, the decay of the turbulence provides energy to heat the gas and turbulence affects the line widths, hence the line cooling.

Ohm resistivity
We consider a neutral gas of number density n <H> (cm −3 ) at gas temperature T gas (K) with a mean mass m n (grams), and mass density ρ n (grams cm −3 ). We define for each charged species j ( j=e, atomic and molecular ions, PAH, gr, where gr stands for dust grains) the units-free plasma β j parameter (which should not be confused with β mag ) as where Z j is the charge of species j of mass m j (in grams), and B is the magnetic field strength, for which we only consider the vertical component B z (the units are the Gauss, denoted Gs, with 1 Gs = 1 cm −1/2 g 1/2 s −1 ) . e is the elementary charge (in units of statC). Finally, γ j = σv j /(m n + m j ) is the drag frequencies of the positively-or negatively-charged species j of number density n j with the neutrals. If a species j has |β j | 1, this means that it is well coupled with the neutral gas. On the other hand, if |β j | 1, the charged species is well coupled with the magnetic field. Using the terms introduced above, the Ohm electrical conductivity is defined as The collision rates σv jn , which appear in the drag frequencies term γ j , are thermal velocity-averaged values of the crosssections σ (without subscript). The electron-neutral collision rate is σv e,n ≈ 8.28 × 10 −10 √ T cm 3 s −1 .
The atomic and molecular ions have average collision rates with neutral species following the Langevin rate σv ion,n ≈ 1.9 × 10 −9 cm 3 s −1 .
Specifically for the HCO + -H 2 system, Flower (2000) provided the formula σv ion,n ≈ 8.5 × 10 −10 T 0.24 cm 3 s −1 , which shows a weak temperature-dependency. At 100 K, the rate becomes 2.5 × 10 −9 . For PAH ions, we use the largest value between the temperature-independent Langevin and the geometrical rate where a PAH is the radius of the PAH. The average collision rate between the negatively-charged grains with the neutral gas species is approximated by σv dust,n ≈ πa 2 8kT gas πm n 1/2 We have assumed that the grains are basically immobile w.r.t the gas. For a fully molecular gas m n ≈2.2 amu = 3.65317 ×10 −24 grams, The Ohm diffusivity is characterized by the dimensionless Elsasser number, which is defined as the ratio between Lorentz and Coriolis force where is the angular speed, which is in the case of a protoplanetary disk assumed to be in Keplerian rotation, where G is the gravitational constant, M * is the mass of the central star (M is the mass of the Sun), and R and r the distance in the disk from the star in cm and in au respectively. One au (astronomical units) corresponds to 1.4959787 × 10 13 cm. The Ohm resistivity is The Ohm Elsasser number can be rewritten as the sum of the contribution of each charged species to the Ohm conductivity: For MRI-driven turbulence to fully develop, the turbulence development timescale τ ideal has to be shorter than the damping by ohmic diffusion timescale τ damp : This criterion translates to In the approximation of vertical isothermal disk (h = c s /Ω) and in the condition of v A = c s (β mag =2, α ideal = 1), the criterion becomes This simple estimate is supported by non-ideal MHD simulations (Sano & Stone 2002). Although the disks are not isothermal in the vertical direction, we will use this criterion in our models. The ohmic Elsasser number criterion can be seen as the rate of dissipation of the magnetic energy (∼ B 2 ) compared to the energy dissipated by ohmic resistivity.

Ambipolar diffusion
The difference between the ion and neutral (atoms and molecules, PAHs, and dust grains) velocities is responsible for ambipolar diffusion (Bittencourt 2004). The Elsasser number for ambipolar diffusion is where the ambipolar resistivity η A is (Wardle & Ng 1999) where σ P the Pedersen conductivity, σ ⊥ is the total conductivity perpendicular to the field, and η O is the Ohm conductivity. The D factor prevents efficient ambipolar diffusion in strongly ionized gas, D = ρ neu /ρ gas where ρ neu and ρ gas are the mass density of the neutral gas and the total gas respectively (Pandey & Wardle 2008). The total perpendicular conductivity is where the Pedersen conductivity is The Hall conductivity σ H is defined by When all the charged species are well coupled to the neutral gas (|β j | 1 for all j), σ O ≈ σ P |σ H | the conductivity is scalar (isotropic) and the regime is resistive. When all the charged species are well coupled to the field, then σ O σ P |σ H | ambipolar diffusion dominates (Wardle & Ng 1999). The ambipolar resistivity becomes In addition to the neutral gas, we included ambipolar conductivity between neutral PAHs, neutral grains and the ions. Based on intensive numerical simulations, Bai & Stone (2011) parametrized the maximum value of α for a given Am The maximum possible value for α is shown in Fig. 4. For Am = 1, max(α) 10 −2 . It should be noted that the formula provides an upper limit on the value of α due to ambipolar diffusion. The ambipolar diffusion prevents the MRI to fully developed at disk surfaces where the neutral-charges coupling is inefficient because of the low density. In the midplane, the gas density is high but the ionization fraction is low at least in the inner disk region. Ambipolar diffusion is most likely the most efficient in the intermediate-height disk layers.

Non-ideal MHD functional form for α and analytical results
In this section, we use the model discussed above to propose an empirical non-ideal MHD MRI-driven formulation for α eff that can be used on physico-chemical protoplanetary disk codes like ProDiMo . Sano & Stone (2002) proposed a limitation to the efficiency of the MRI-driven turbulence when the Elsasser Ohmic number Λ Ohm is larger than one The Elsasser Ohmic criterion tests the coupling between the magnetic field and the charged particles in the disk. For a constant magnetic field strength in the disk, the coupling is more efficient at the disk surfaces, where the gas density is low and the abundances of charged particles are high. Wardle & Salmeron (2012) argued that the Hall diffusion can overcome the resistive damping of the MRI in certain circumstances. Numerical simulations of the Hall diffusivity give a more complex picture (Lesur et al. 2014), however. In light of the uncertainties, we decided to not include the effects of the Hall diffusivity in our estimate of the effective viscosity parameter α eff . The adopted functional form for α eff is when β mag Λ Ohm > 1 and α eff 0, otherwise. It has been shown that, when β mag Λ Ohm < 1, all the perturbation modes are stabilized and thus no MRI-driven turbulence can exist (Jin 1996). Our formula attempts to encompass effects from numerous detailed MHD simulations. As such the effective α eff should only be considered as a practical simplification to be used in non hydrodynamic codes like ProDiMo. With more MHD simulations in the future, the approximation will certainly evolve.
The free parameters of the MRI-model are the value of β mid and the general disk parameters (see Table D.1). All other parameters and variables are derived either from β mid or from the chemistry, gas and dust thermal balance and radiative transfer. While both Ohm and ambipolar diffusion contributions to the effective turbulence depend on the available charges, the two effects differ on the other dependencies. The Ohmic Elsasser num-T gas =100K, n H =10 10 cm -3 min[χ(e)]=10 -3 ζ Effective turbulence coefficient α eff as function of β mag for gas of density of n <H> =10 10 cm −3 and 10 12 cm −3 at T =100 K in a disk midplane. The minimum free electron abundance is 10 −3 with respect to the total number of negative charges (i.e. all the other charges are on dust grains). The left panels show the ideal-MHD, Ohmic, and ambipolar diffusion contributions to α eff , while the two right panels show the resulting α eff including the all mode damping criterion of Jin (1996). Lines at 10 −10 on the right panels mean complete no MRI-driven turbulence for total ionisation fraction of 10 −11 for all values of β mag .
ber is sensitive to the magnetic field strength and the gas temperature, while the ambipolar Elsasser number depends on the absolute gas density. Consequently, in low-mass disks and at disk upper surfaces, MRI efficiency will be limited predominately by ambipolar diffusion. In the disk midplane, the ionization is weak because even the cosmic ray flux can be attenuated (Rab et al. 2017). Ionization is proportional to the gas density (∝ ζn <H> ), whereas the electron recombination depends on the square of the density (∝ k rec n ion n e ). This translates to an inverse density dependency for the number of ions in the midplane. Therefore, Am can also be small in the disk midplane where ionization sources are scarce. Both Elsasser numbers increase with radius with power 1.5. But if the gas density decreases as r −2.5 than the ambipolar diffusion Elsasser number will decrease.
To illustrate those dependencies, we plotted in Fig. 5 the effective turbulence as the function of the parameter β mag for a gas and dust temperature of 100 K and for two gas densities (n <H> =10 10 cm −3 and 10 12 cm −3 ). The panels show the effects of the variation of the Ohmic and ambipolar diffusion resistivity on the value of α eff as well as the value of ideal MHD-MRI α ideal . At low and medium densities and for strong magnetic field strength (i.e. with low values of β mag ), the effective turbulence is dominated by the ambipolar diffusion term. The Ohmic restric-tion on α eff occurs for a large value of β mag . The value of α eff is larger than 10 −3 down to a total ionization fraction of ∼ 10 −8 -10 −10 depending on the gas density. At disk surfaces, where the ionization fraction is high and the value of β mag small, the efficiency of MRI turbulence is probably limited by the ambipolar diffusion resistivity to a maximum value of ∼0.25.

Disk model
We chose to model a typical protoplanetary disk with a total gas+solid mass of 10 −2 M (Woitke et al. 2016). The disk extends from r in =1 au to r out =600 au with a tapered outer disk. The dust is composed of compact spherical grains composed of Olivine, amorphous carbon, and Trolite (FeS). The dust grains size distribution follows a power-law with index -3.5 from a min =0.1 µm to a min =3000 µm. The global gas-to-dust mass ratio has the standard value of 100. Dust grains can settle with a turbulent mixing parameter α settle of 0.01. The PAH abundance depletion factor compared to the interstellar medium value f PAH is set to 0.01 ( f PAH = 1 corresponds to an abundance of 3×10 −7 ). The disk flares with an index of 1.15. The gas scale-height is 1 au at 10 au (10%). The cosmic ray flux is assumed to be constant throughout the disk at 1. varying the X-ray fluxes, lowering the cosmic ray flux and its attenuation in disks (Cleeves et al. 2013(Cleeves et al. , 2014, and modelling the effects of stellar particles (Rab et al. 2017). The protoplanetary disk parameters are summarized in Table D.1.
The key parameter of our parametric MRI model is β mid . We modelled the disks with β mid in the disk midplane from 10 2 to 10 6 . For comparison, we also ran a model with no (MRI) turbulence, i.e. the line broadening is purely thermal (v turb =0 km s −1 ) and other models with a constant turbulence width v turb =0.05, 0.1, 0.15 and 0.2 km s −1 and no turbulence heating throughout the disk (passive disk models). The central object is either a TTauri star (M * = 0.7 M , T eff = 4000 K, L * = 1.0 L ) or a HerbigAe star (M * = 2.3 M , T eff = 8600 K, L * = 32 L ). The star shows excess UV and strong X-ray emission. The disk parameters are summarized in Table D.1.
The disk thermal balance has been modified to include the effects of turbulence heating and turbulence line broadening.
The release of turbulent energy in the gas at each disk location is (Hartmann 1998) The turbulence will effect the gas temperature, which in turn will change the chemistry, the value of β mag , the ohmic and ambipolar Elsasser number. In an analytical analysis in the case of ideal MRI-turbulence, the energy release is Adopting δ=0.5, E acc ∝ P therm Ω.
Assuming that the pressure term can be described by a barotropic law P therm ∝ n γ ,  Fig. 8. The Elsasser Ohm number is shown in the upper panels for the DIANA typical disk and. The white contours corresponds the location of the total charge in the disk. The middle panels show the ambipolar diffusion number in the disk models. The location where the C and C + abundances are equal are overplotted in red. The criterion for all modes to be damped is shown in the lower panels. The lefts panels are models with β mid =10 4 and the right panels are models with 10 6 .  Fig. 9. The effective α eff value and turbulence velocity over sound speed ratio for the DIANA typical disk and β mid =10 4 (left panels) and 10 6 (right) panels. The white contours in the upper panels show the location where the Elsasser Ohm number Λ Ohm is unity. The level where there is no MRI in the inner disk midplane, the so-called "dead-zone", is shown in red (Jin=1). For β 1/2 mid Λ Ohm < 1, MRI is entirely suppressed and α eff tends to zero.The contour where β = 1 is also shown in the upper panels. The contours in the low panels indicate the gas-phase CO abundances. the gas heating rate is The value of γ can vary typically from 1 to 3. It is interesting to note that in our prescription, the turbulence heating efficiency depends on the thermodynamics of the gas. The effects of increasing the gas temperature on the value of α eff are not easy to predict. α ideal is proportional to T −δ , Λ Ohm is proportional to √ T . The dependence of Am on T is not direct, but assuming electronic recombination rates ∝ T −1/2 , an increase of the gas temperature should result in a decrease of the ion recombination rates and thus a higher value for Am. Therefore, as the gas temperature increases due to viscous dissipation, the resistivities decrease but at the same time the ideal-MHD α eff decreases as well.
The turbulent velocity is sub-sonic and is related in our model to α eff and the sound speed c s by v turb = √ α eff c s (Hughes et al. 2011). Another relationship such as v turb =α eff c s has been proposed    Fig. 10. The upper panels show the main heating sources and the lower panels the main cooling sources at the different location in the disk. The left panels are for the β mid =10 4 disk model, while the right panels correspond to the β mid =10 6 disk model. In the upper panels, the contour of the Jin's criterion is shown in white contour. In the panels, the black contours correspond to A V =10. ∆v = (v 2 th + v 2 turb ) 1/2 , where v th is the thermal broadening. The turbulent velocity will affect the line optical depths as 1/∆v. In turn, the cooling term is affected because of the change in the line cooling rates. The MRI-turbulence does not only affect the heating but also the cooling of the gas.
The present model only considers viscous accretion, through the alpha parameter. However, it is known that, in the presence of a poloidal magnetic field, discs can be subject to MHD wind which can also drive accretion (e.g., ). We have purposely neglected this process due to the lack of robust prescriptions for wind-driven accretion.

Disk model results
The disk total charge (a.k.a. ionization fraction) is shown in the upper-left panel of Fig. 6 for a disk model with β mid =10 4 . The figures show the different disk structures with radius r in au in the horizontal-axis and z/r in the vertical axis.
The disk ionization fraction (or total charge) is governed by the balance between ionization due to UV, X-ray, and cosmic rays in the disk and the recombination and charge exchange reactions. Fig. 6 shows that the free electrons and the gas-phase cations are the main negative and positive charge carriers except in a small region in the inner disk midplane.
The dust grain average charge is displayed in the upper-left panel of Fig. 7, while the fractional contribution of the charges (positive and negative) on grains to the total charge is exhibited in the upper-right panel. From the upper disk atmosphere to the midplane, the dust grains have first lost hundreds of electrons due to the photoelectric effect. When the UV field decreases, the free electrons (up to a relative abundance of a few 10 −4 ) from atomic carbon ionization stick to the grains, rendering them negatively-charged. The total charge reaches the value of the elemental abundance of carbon (1 to 3 ×10 −4 ). As carbon becomes neutral at lower altitude, the recombination with the free electrons compensates almost exactly for the photoejection of the electrons from the grains and the grain charge fluctuates around zero. The remaining sources of electrons are the atoms with ionization potential lower than 13.6 eV (Fe, Mg, Si, S). Therefore, gas depletion of the metal species (Fe, Mg, Si), in addition to that of Sulfur will determine the electron fraction χ(e) in the intermediate disk heights. It should be noted that the dust grains remain at temperatures T dust well below the sublimation temperatures (T subl. ∼ 1500 K) even though the gas temperature can reach up to ∼ 10,000 K (Fig. F.1). We have not considered the destruction of dust grains in warm gas due to chemisputtering.
In the midplane inner disk region (r < 1 au with total charge fraction < 10 −12 ), the radiation field is too weak to eject the ex-cess electrons on the grains. The grains are negatively charged after the attachment of the free electrons created by the interaction of cosmic rays or X-ray photons with the neutral gas, consistent with the analytical approximation (see Sect. 2.5). The negatively-charged dust grains become the major dominant negative charge carrier. The cations remain the main positive charge carrier.
In the lower-left panel of Fig.7, we see that the contribution of the positively-charged PAHs to the total charge is small. The contribution of the negatively-charged PAHs is relatively high but not dominant in the region above the freeze-out region of the PAHs. The relative contribution of PAHs as charge carriers depends on their abundances in disks. In the typical model, the PAH depletion factor f PAH is between 0.01 and 1 (Woitke et al. 2016). In the disk midplane, the PAHs are frozen onto the grains.
The Elsasser Ohm number, the ambipolar diffusion number, and the complete mode damping criterion are shown in Fig. 8 for the disk models with β mid =10 4 (left panels) and β mid =10 6 (right panels). The Elsasser Ohm number is the most sensitive to the value of β mag . The ambipolar Elsasser number Am distributions are similar for the β mid =10 4 and β mid =10 6 . In unobscured regions, the electron fractional abundance χ(e) is relatively high in the range 10 −8 -10 −4 . Due to the limit on negative charges on grains (eq. 50), in a gas with total ionization fraction greater than a few 10 −12 , the free electrons will be the main contributor to the Ohm Elsasser conductivity. In the inner disk midplane, the total charge is low, resulting in low values for the Ohm and ambipolar Elsasser numbers.
The approximated Ohm Elsasser distribution for a disk model with β mid = 10 4 is shown in Fig. G.1.
The disk distribution for the effective turbulent parameter α eff is shown in Fig. 9 for β mid of 10 4 (the left panels) and 10 6 (the right panels). The disk turbulence efficiency structure can be divided into zones. The zone limits are defined by the transitions between the MRI-driven region, the Ohm diffusion limited deadzone according to the total damping criterion Jin ≡ β 1/2 mag Λ Ohm =1, the Ohm diffusion restricted MRI region (Λ Ohm =1), and the location of the disk where β mag =1 in the disk atmospheres. The decrease of α eff with the disk radius stems from the decrease in the gas pressure, hence of β mag (see the middle panels of Fig. E.1). The outer disk midplane decrease of α eff stems from the increase of β mag . The ambipolar diffusion does not restrict the value of α eff in the disk. At the top disk surfaces, the hydrogen gas is ionized due to the X-rays. As the protons recombined and the gas becomes more neutral, the main ions are the C + cations. Below the fully ionized carbon layer, the ambipolar diffusion number Am plummets by two orders of magnitude. Despite the drop, Am remains higher than unity for most of the disk. The ambipolar diffusion is also affected by the metal abundances. At large radii, the vertical column density is low such that the gas can be ionized and Am increases again.
The lower panels of Fig. 9 show the turbulence over the sound speed velocity ratios. The location of the gas-phase CO is overplotted as contours. For β mid =10 4 , CO will be emitted significantly in the turbulent warm molecular disk layers with v turb /c s between 0.25 and 0.4 . For the β mid =10 6 model, v turb /c s is between 0.1 and 0.3. For completeness, the sound speed c s for the two models are plotted in Fig. E.1. The CO emission comes mostly from gas with sound speed ∼0.3 -0.5 km s −1 . Fig. 10 shows the main heating and cooling processes for the disk models with β mid =10 4 (left panels) and β mid =10 6 (right panels). In the dead-zone where β 1/2 mag Λ Ohm < 1, the gas heating occurs mainly via the gas thermal accommodation on warm dust grains. Outside the dead-zone, viscous heating is the dominant heating process in both disk models from the midplane to the intermediate layers (z/r ∼ 0.2) till disk radii of a few hundred au, where the density, the Keplerian rotation, and α eff are low (see eq. 85). The outer disk is heated by the energy released by H 2 formation on dust grains and by the chemical reactions. The disk atmosphere region where PAHs are not frozen onto dust grains in the β mid =10 4 model can be heated by the photoelectric effect. The main gas heating processes of passive disks in the atmospheres are hits by cosmic rays and photoelectric effects from dust grains and PAHs (Bergin et al. 2007) and chemical heating (Woitke et al. 2016). Towards the midplane the densities are high enough for the gas and the dust to become thermally coupled. In an active disk, these processes are minor as soon as the magnetic field strength is not too weak (or β mid > 10 6 ). In the regions with β mag < 1, the gas is heated by the X-ray Coulomb heating. The main cooling agents in the dead-zone are the water vapor and molecular hydrogen rovibrational lines. In the viscous heating dominated region, the cooling occurs via thermal accommodation on dust grains. Outside the viscous dominated region, CO isotopologue rotational lines are the main coolants in the warm molecular layers, and the [OI] and [CII] fine-structure lines are the main cooling lines in the atomic disk atmosphere. In the X-ray dominated region, [OIII] and Lyman α are the main cooling lines. This is consistent with the conclusions of Najita & Adámkovics (2017) who modeled the effects of heat generated by turbulence decay and concluded that warm CO and H 2 O lines can trace the mechanically heated gas.
The actual CO line profiles computed in non-LTE for CO, 13 CO and C 18 O J=2-1 and J=3-2 are shown in Fig. 11 for the TTausi and HerbigAe disks seen with a 45 • inclination. The line flux increases with more turbulent gas because of stronger heating and that more extended emitting area per velocity resolution elements.The signature of turbulence broadening on the line profiles is not obvious, although the CO lines should probe the most turbulent regions of protoplanetary disks. The CO profiles for the disk model with β mid =10 6 and that for the model with no turbulence are similar for both the J=2-1 and J=3-2 transitions. For both the CO J=2-1 and CO J=3-2 lines computed in models with a fixed value of v turb , one can constraint v turb down to 0.05 km s −1 (see the red-dotted lines in Fig. 11). Fig. 12 shows the variation of the peak flux, FWHM, and the profile peak separation as function of β mid for the CO lines for the TTauri disk models. The ratio between the peak and line centre flux (peak-totrough) for the TTauri and HerbigAe disk models are shown in Fig. 13 for the full MRI models (upper panels) and for disk models with a fixed value of the turbulence width (lower panels). The choice of β mid influences mostly the ratio between the peak and line centre flux, consistent with the results of the parametric disk model of Flaherty et al. (2015). The trend of decreasing peak-tocenter ratio with increasing disk turbulence velocity shows the same amplitude (0.2 to 0.4) than the trend of a lower ratio for high values of β mid . The peak-to-center flux ratios of the 13 CO transitions (J=3-2 and J=2-1) show the strongest sensitivity to β mid or v turb . A low peak-to-trough ratio indicates a high level of turbulence. The 12 CO lines are emitted in the region where v turb /c s reaches its maximum. The optically thinner 13 CO lines are emitted in the layers closer to the midplane where v turb /c s is lower. The peak flux is mostly emitted in the outer disk, where the C 18 O lines can be optically thin. Therefore, the peak fluxes for the C 18 O lines are low and the peak-to-trough ratios for the C 18 O lines are less sensitive to the turbulence broadening than the optically thick 12 CO and 13 CO lines.
One possible parameter that influences the peak-to-trough ratios for the 13 CO lines is the abundance of 13 CO compared to 12 CO. The lower-right panel in Fig. 13 shows the effect of a global 12 CO/ 13 CO ratio of 150 instead of the standard value of 70 (Henkel et al. 1994;Wilson 1999) in the Interstellar Medium (ISM). A higher value for 12 CO/ 13 CO in Young Stellar Objects (YSOs) environments compared to the ISM value is consistent with the observations of the 12 CO/ 13 CO ratio towards a number of YSOs by Smith et al. (2015), who found a range of 86 to 165. The higher value for 12 CO/ 13 CO may reflect the less efficient 13 CO self-shielding against photodissociation compared to 12 CO (Miotello et al. 2016). Other parameters such as the disk inclination or the central star mass affect the value of the peakto-trough ratios.

Discussion
Gas and dust in passive disks are heated by the conversion of the stellar radiation (from the X-ray to the IR) to thermal energies (Jonkheid et al. 2004;Bergin et al. 2007;Woitke et al. 2009;Bruderer et al. 2012). Passive disks show an increase of the temperatures from the midplane to the disk surfaces because of the extinction by dust grains for the UV and optical photons and by the gas for the X-ray photons. The warm layers reach temperatures of a few hundred to a few thousand Kelvin. In the inner 10-20 au, the turbulence is low due to the low abundance of charge species. Mass accretion has to occur in the low gas densities, hot and relatively ionized upper disk layers above the dead zone to compensate for the decrease of mass accretion in the midplane. This result contrasts with vertical isothermal and relatively cool disk models where even extremely high values of α eff 7 6 5 4 3 2 1 log(β mid ) are not enough to compensate for the low accretion in the midplane (Perez-Becker & Chiang 2011). Our modelling supports the idea put forward by Gammie (1996) that accretion would occur mainly through the active surface layers.
Beyond 10-20 au in the midplane, the ionization is high enough such that the MRI-turbulence can fully develop. Accretion heating becomes important and since the density is high, gas-grain thermal accommodation is the main source of gas cooling. The dust grains are thus indirectly partly viscously heated. The rate of gas-grain accommodation depends on the density-squared and thus will become inefficient at large radii. The dust grains act as a thermostat for the gas through the gasdust thermal accommodation. When the gas is heated by the turbulence viscosity, the dust grains receive extra energy from the gas and radiate it. When there is no turbulence heating, the dust grains transfer some of their energy to the gas.
The value of α eff reaches a maximum of 0.1, which translates to a turbulent width of 0.4 times the sound speed with values between 0.2 and 0.4 c s for a β mid =10 4 in the midplane. The computed turbulent line widths assuming β mid =10 4 are consistent with the limit found by Teague et al. (2016) in the TW Hya disk but not with the low value determined by Flaherty et al. (2018). Flaherty et al. (2015Flaherty et al. ( , 2017 also found a low value in the protoplanetary disk around the HerbigAe star HD 163296 (v turb < 0.03 c s ). Simon et al. (2018) found that a low ionization rate and a low magnetic field is needed to explain a weak turbulence in the outer regions of protoplanetary disks. Teague et al. (2016) found a nearly constant ratio between the turbulence width and the local sound speed of v turb = 0.2 c s . The TW Hya disk is seen at a low inclination so that there is no possibility to use the line peak-to-trough ratio to constrain the turbulence. Flaherty et al. (2015) observed a peak-to-trough value of ∼2.0 for the 12 CO J=3-2 line. In our MRI and fixed v turb models, the peak-to-trough has a maximum value of 1.9 for the HerbigAe disk model with 45 • inclination. In the ALMA observations presented by Flaherty et al. (2015), the 12 CO J=3-2 and J=2-1 profiles show larger peak-to-trough ratios than the 13 CO profiles (see the right panels of Fig. 13) in contrast to our model results. A higher 12 CO/ 13 CO ratio can decrease slightly the peakto-trough value, which is still not sufficient to explain the observed value (see the lower-right panel of Fig. 13).
An alternative to using high-quality CO lines, the ionization structure of protoplanetary disks can be constrained by matching emission lines from molecular ions. Cleeves et al. (2015) argued using emission lines from HCO + and N 2 H + that the TW Hya disk has a dead-zone extending up to 50-65 au. Their disk model requires a very low cosmic ray ionization rate lower than 10 −19 s −1 . Their dead-zone is defined by the criterion R EM ≤ 3300, where R EM is the magnetic Reynolds number (Balbus & Terquem 2001). HCO + and N 2 H + are molecular ions whose theoretical abundances depend sensitively on the treatment of ionization and recombination in the gas (Rab et al. 2017) as well as of the N 2 self-shielding for N 2 H + Visser et al. (2018). In disks models, molecular ions are confined to the upper disk molecular layers. Using molecules other than CO would require a detailed understanding of the chemistry of those species in disks. Kamp et al. (2017b) showed that the N 2 H + abundance and line fluxes in disks vary significantly when different standard gas-phase chemical networks are used. Finally, although the HCO + emissions may be strong, N 2 H + emissions are faint in disks. Guilloteau et al. (2012) used CS as a tracer of turbulence because CS is much heavier than CO such that its thermal width is smaller. The low signal-to-noise data allowed them to constrain turbulence to be between 0.4 c s an 0.5 c s . ALMA observations of CS in the DM Tau disk reveals a non-thermal line width between 0.05 and 0.15 km s −1 at 300 au (Semenov et al. 2018). Sulfur chemistry is complex and CS lines are weak even in massive disks Dutrey et al. (2011).
The magnetic Reynolds number and the Elsasser Ohmic number are then related by assuming c s = hΩ. Flaherty et al. (2017) used DCO + lines in addition to CO and C 18 O to constrain α, and found that the limits given from the three tracers are consistent with each other. When turbulence is low, it becomes difficult to distinguish between the CO line profiles of MRI models with β mid < 10 5 , i.e. that CO line profiles can be used to constrain α down to ∼ 5 × 10 −3 . The CO abundances can also vary over the lifetime of the disks (Yu et al. 2017), which may affect the CO line profile. The CO rotational lines that are seen with large beams in the sub-millimeter domain probe a large disk volume but may not  Fig. 13. Peak flux to zero velocity flux ratio for standard TTauri disk model (left panels) and the HerbigAe disk model (right panels) as function of the β mid parameter (upper panels) and of v turb (lower panels). All the models are viewed with an inclination of 45 • (0 • means that the disk is seen face-on).
be sensitive to the location in the protoplanetary disks with the highest v turb /c s ratio, although high spatial resolution CO maps can be obtained with ALMA.
Alternatively, the CO fundamental rovibrational line (v=1-0, ∆J = ±1) emission area matches the location of the maximum of the v turb /c s ratio (see the lower panels in Fig. 9).
The knowledge of the gas and dust temperature structure is central to the estimate of the line thermal broadening. Fig. F.1 shows that the disk dust and gas temperature structures may be complex. In particular, the gas temperature structure depends on β mag . Because of the dead-zone and the efficiency of the gasdust thermal accommodation at high densities, the inner disk midplane gas temperature does not change much with accretion heating.
The charge in protoplanetary disks depends on many parameters that were fixed in our models. In the UV-dominated regions, the electrons are provided by the ionization of neutral carbon atoms, while the abundance of the gas-phase metals play an important role in controlling the ionization fraction in the UV-shielded molecular region. Therefore, the ionization level in the UV-shielded region depends strongly on the gas-phase metal abundances. In the disk areas where T dust is higher than the sublimation temperature of the main silicate dust grains, the gas phase is enriched by a large quantity of silicon, magnesium, iron, sodium, ..., reaching the stellar photospheric abundances. Those elements can still survive at high gas temperatures in form of metal oxides like SiO, MgO, FeO if UV is not present. But since the main UV opacity, i.e. the dust, is no more present, UV photons permeate the dust-free gas and can dissociate efficiently the oxides unless efficient high-temperature formation of those oxides is possible like for water vapor (Thi & Bik 2005). If the gas chemistry is at equilibrium, one can use the Saha equation to compute the ionization fraction assuming that the collisional ionization controls the ionization fractions (Desch & Turner 2015).
Negative grains are formed by attachment of an electron or charge-exchange with an anion (in our simple chemical network there is only one anion H − ). Electrons are first created by ionization of H 2 into H + 2 by an energetic ionization event (X-ray, Cosmic rays, radioactive decay), and subsequently H + 3 . They can recombine either with an abundant atomic or molecular ion or with a dust grains. Cations and negatively-charged grains can interact to neutralize each other.
In the gas phase, the radiative recombinations of metallic ions with electrons are slow because the excess energy has to be radiated away. On the other hand, the excess energy in the ion recombination on negatively-charged grain surfaces is transferred to the grain, making the recombination rates with grains much higher than with an ion. The grain surfaces act as a heat sink. This may result in neutral gas in the disk midplane and thus a much lower ionization fraction.
Many parameters can influence the ionization structure and hence the efficiency of MRI-driven turbulence in protoplanetary disks. The most important are not well constrained: β mag , the metallicity, the depletion of the metals into refractory solids, cosmic ray flux, and spectrum, stellar particle flux, grain properties such as the size distribution, settling and drift. The effect of those key parameters can compensate each other. For example, a strong magnetic field can compensate for a low cosmic ray flux. Likewise, a strong X-ray flux can compensate for a weak cosmic-ray flux. The direct grain ionization by absorption of Xray photons was not considered in this paper and would the subject of a subsequent study.
Our disk model assumes dust settling, which can affect the UV penetration and hence the ionization structure. In our modelling of the dust gain distribution, the vertical settling is parametrized by a constant dust turbulence mixing parameter α settle of 0.01 throughout the disk (Woitke et al. 2016). How this dust mixing parameter α settle relates to the effective turbulence parameter α eff is not clear. Detailed modelling of the dust vertical structure of HL Tau suggests that the dust vertical mixing is weak with a value of a few 10 −4 for α settle . This corresponds to a disk model with β mid ∼ 10 6 .
In our model, the value of α eff effects both the sound speed (because turbulence is an efficient heating agent and changes the cooling line transfer) and the turbulence width. The change in the line cooling efficiency is small because the CO lines (one major gas coolant for the outer disk gas) are highly optically thick. The contribution of the turbulence to the total line width remains small even for high values of α eff . A fundamental difference between our models and the parametric disk models used to constrain the turbulence speed is that the value of α is assumed uniform in the parametric disk models. Future studies are required to model consistently α settle and α eff , as well as the changes in the disk density structure, in the disk chemistry, and in the grain properties (size, drift, ...) over the disk lifetime.

Conclusions
We have implemented a simple parametrized magnetorotational instability (MRI) driven turbulence model in the physicochemical code ProDiMo to constrain the gas heating by gas turbulence decay and the potential presence of a dead-zone. The strength of the turbulence in protoplanetary disks is computed consistently with the ohmic and ambipolar diffusion, whose resistivities depend on the abundance of the charge carriers (free electrons, gas-phase atomic and molecular cations, PAHs, and dust grains). The gas and dust temperatures, chemistry, PAH and grain charges are computed self-consistently together with the continuum and line radiative transfer. The disk models include at the same time active and passive heating processes. The main conclusions are: -The free electrons and gas-phase ions are the main charge carriers in most parts of the disk. -The negatively-charged PAHs are the major charge carriers in the disk atmosphere above the PAH freeze-out zone. -The dust grains are the dominant negative charge carriers in the inner disk midplane where the ionization fraction is very low. -The Ohm resistivity governs the location of the dead-zone in the inner disk midplane. -The ambipolar diffusion resistivity does not affect much the efficiency of the MRI-driven turbulence because the density of ions is relatively high in the whole disk. The last two conclusions may change when other disks are considered. -In the inner 10-20 au, gas may accrete through the warm layer that sits above the dead-zone. -The gas accretion heating dominates in the disk midplane outside the dead-zone with the dust thermal accommodation being the main gas cooling agent, while in the dead zone thermal accommodation is the main source of gas heating and the gas is cooled by molecular line emissions. In the disk upper layers, the gas is mainly heated by the photoejected electrons from PAHs and dust grains as well as by gas absorption of the dust grain infrared emissions and by H 2 formation. -The signatures of the gas turbulence in the CO emission lines are weak and sophisticated modelling of the high signal-tonoise observations is required to derive meaningful estimates of the turbulence parameter α. -The CO line profiles from the MRI disk models differ from the profiles of disk models with a single constant value of α. Simultaneous fits to multiple transitions of CO lines and to other disk tracers would improve the reliability of the estimate of β mid (or v turb ), although the line profiles are shaped by many disk properties that are currently not well constrained.
Future works include the study of the influence of several model parameters on the location of the dead-zone. The important parameters include the cosmic ray flux attenuation in the disk, the stellar particle flux, the size and mass of the disk, the stellar properties. Another venue is the study of the disk evolution and the changes in the value of α eff with time (Bai 2016). Since our model can generate different types of observables (molecular line maps and continuum emission maps), one can use our model outputs as test cases to improve the reliability of inversion models used to derive α from the observations. Finally we will also explore emission lines from species other than CO as turbulence tracers. where m i γ i 1.9 × 10 −9 (A.16) assuming that m i >> m n . In reality, a fraction χ(e) of the negative charges will be electrons locked in grains.
In UV-dominated regions at disk surfaces, the dust grains will be positively charged with the maximum possible positive charge set by formula 36. Formula A.11 reads σ e,O σ dust,O 5.5 × 10 13 a 4 µm χ(e). (A.17)

Appendix B: Effective turbulence coefficient
Fig. B.1 shows the effective turbulence coefficient α e f f as function of the total ionization fraction for a relatively dense gas at density 10 10 cm −3 and temperature of 100 K. The four panels correspond to different values of β mag (10 4 , 10 3 , 10 2 , 10).  Effective turbulence coefficient α eft as function of the total ionization fraction for a gas of density n <H> =10 10 cm −3 at T =100 K and four values of the β mag parameters (10 4 , 10 3 , 10 2 , 10). The values for α ideal as plotted in red short-dashed-lines have been derived assuming δ=0.5. The limits imposed by the am bipolar diffusion resistivity are shown in dashed-lines.
(integer) and f (Z d ) is the normalized discrete distribution (fractional abundance) of grains with charge Z d . The values for Z min and Z max are chosen such that f (Z min ) = 0 and f (Z max ) = 0.
Instead of solving for each individual charge, which will increase the number of differential equations in the system by the number of charges, we only considered three grain pseudo charge-species Z − , Z, and Z + with the abundance [Z − ], [Z], and [Z + ] respectively. The gas-phase species react with one of these three pseudo-species only.
The grain charge distribution f (Z d ) is The total abundance is We defined the pseudo-charge species (Z − , Z − , Z) such that 3) where we have split the grain population into three categories: the negatively-charged grains the positively-charged grains The neutral grains are modelled by the pseudo-species Z of abundance The average grain charge is The pseudo-neutral grain species Z includes the actual neutral grains and the positively-charged grains exchanging charges with cations of index j to become more positively-charged and absorbing a photon and ejecting an electron: Considering only this reaction, the chemical reaction differential equation is where the positively-grain-charge-averaged charge-exchange rate k 0 ex,i is . (C.12) The absorption of a neutral grain or of a positively-charged grain leads to the ejection an electron: The rate for the pseudo species Z is k 0 pe , whose derivation follows (C.14) We obtain Likewise, the pseudo-neutral grains model the thermal emission of electrons actual neutral and positively-charged grains: The rate of thermionic emission is . (C.18) The pseudo-neutral species Z is also being used to model electron attachment reactions where k 0 e = k e,n f (0) The electron recombinations of positively-charged grains are modelled by the pseudo-reaction Z + + e − → Z.

(C.22)
Consider the case of the electron recombination with a grain with n positive charges and rate k re (n), the sum of all the recombination chemical rate reactions reads The pseudo-negative charge species Z − models the reactions leading to neutralize the ions A of index j: The generic reaction is written as . (C.29) The pseudo-species Z − also participate to the photodetachment reaction Z − + hν → Z + e − , (C.30) with the rate . (C.31) Finally Z − can emit electrons due to the thermionic effect . (C.33) To illustrate that our concept is valid, we can assume first that th the grains can have at most one positive or negative charge, in other words, Z min =-1 and Z max =1. Then the reactions in Table C.1 correspond to the actual reactions as presented in Table C.2. The definition of the rates above depends on the charge distribution f (Z d ), which has to be known or solved simultaneously with the chemical reactions. In a steady-state treatment, gr + A + → gr + + A 0 0 k j ex,0 = k j gr,ion (0) 2 gr + hν → gr + + e − 0 0 k 0 pe = k pe (0) 3 gr → gr + + e − 0 0 k 0 th = k th(0) 4 gr + e − → gr − 0 0 k 0 e = k e,n 5 gr + + e − → gr 1 1 k + e = k e,+ (1) 6 gr − + A + → gr + A -1 -1 k j ex,− = k j gr,ion (−1) 7 gr − + hν → gr -1 -1 k − pe = k pe (−1) 8 gr − → gr + e − -1 -1 k − th = k th (−1) the chemistry and the grain charge distribution are solved alternatively.
Using an estimate of the electron and ion abundances, the grain charge distribution f (Z d ) is calculated by iteration. The average grain charge < Z d > is determined by balancing the photoemission (photoejection and photodetachment) rates and the thermionic rate with the recombination and charge exchange rates. Then all the chemical rates between a gas-phase species and the pseudo-species Z − , Z, and Z + can be derived. The electron abundance is computed from the global gas neutrality considering the charges on gas-phase species, PAHs, and dust grains. The grain charge distribution f (Z d ) is subsequently reevaluated knowing the abundances of the cations and electrons. This iterative process ensures that the grain charges and the gascharge chemistry are computed self-consistently. Especially the global neutrality of the medium is enforced. We assume that a dust grain can have charges from Z min and Z max , which are numerical free parameters. The results do not depend on the exact values for Z min or Z max provided that Z min ≤ Z d − 10 and Z max ≥ Z d + 10. We assumed in this work for the numerical reason that Z min = −Z max . If the chemistry is solved timedependently, the grain charge fractional distribution has to be solved together with the chemical reaction differential equations. Solving a time-dependent grain charge distribution entails many challenges, whose resolution goes beyond the scoop of this paper.