© W. F. Thi et al. 2019

1 Introduction

Pre-main-sequence stars (e.g., T Tauri and Herbig Ae 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⁻² 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⁻⁹–10⁻⁷ M_⊙ yr⁻¹ (Calvet et al. 2000; Muzerolle et al. 2004; Johns-Krull et al. 2000; Hartmann et al. 1998), although episodic accretion rates of up to 10⁻⁵ M_⊙ yr⁻¹ 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, 2011; Garcia Lopez et al. 2006).

The most promising mechanism that can explain mass accretion in disks is turbulence driven by the magnetorotational instability (MRI; see Balbus & Hawley 1991; Balbus 2011; Fromang & Nelson 2006; Bai 2015; Simon et al. 2013, 2015; Bhat et al. 2017; Béthune et al. 2016; O’Keeffe & Downes 2014; Fleming & Stone 2003; Davis et al. 2010) because molecular viscosity is too weak. Alternatively, a weak turbulence may be generated by purely hydrodynamic instabilities, such as the vertical shear instability (Nelson & Papaloizou 2004; Lin & Youdin 2015), gravitational instability for self-gravitating disks (Gammie 2001; Forgan et al. 2012; Hirose & Shi 2017), vertical shear instability (Flock et al. 2017), zombie vortex instability (Marcus et al. 2015; Lesur & Latter 2016), and baroclinic instabilities (Klahr & Bodenheimer 2003; Lyra & Klahr 2011).

Magnetorotational 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_sh, 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 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 physical and chemical modeling of the gas and dust together with the continuum and line radiative-transfer. 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 to be small grains and not macro-molecules.

Ionization fractions in realistic protoplanetary disk s were modeled at different levels of sophistication. The aim of these 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 oragainst 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 by Ilgner & Nelson (2006b) and Heinzeller et al. (2011). Fromang et al. (2002) modeled 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 ionizations 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 disk s. 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 disk s have been obtained for the TW Hya (Teague et al. 2016; Flaherty et al. 2018) and HD 163296 (Flaherty et al. 2015, 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 lie between 0.04 c_s and 0.13 c_s. Flaherty et al. (2015) provided a low upper limit of v_turb < 0.03c_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.1c_s for the TW Hya disk and v_turb ≤ 0.4c_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 MRI-turbulence heating and cooling in the photochemical disk code PRODIMO. We considered far-ultraviolet emission from the star and the accretion excess, X-rays, 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₂, because it iswidespread through the disks, and because its chemistry is well understood. The numerical study of MHD processes in triggeringdisk turbulence is beyond the scope of this paper. Instead, we parametrize the onset and effects of turbulence on the gas by a simple empirical equations. This study focuses on the role of the gas and grain chemistry in the disk’s ionization, which in turn 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 physico-chemical processes of the disk.

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 Sects. 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 Sects. 5 and 6; our findings are discussed in Sect. 7; and we conclude in Sect. 8.

2 Gas and dust charge exchange reactions

2.1 ProDiMo

PRODIMO is a code built to model the gas and dust grain physics and chemistry (Woitke et al. 2009, 2016; Kamp et al. 2010). It has been used to model disk spectral energy distri- butions (SEDs, Thi et al. 2011), water deuteration chemistry (Thi et al. 2010), CO rovibrational emissions including UV- fluorescence (Thi et al. 2013), and many Herschel observations from the GASPS large program. X-ray physics are implemented (Aresu et al. 2012, 2011; Meijerink et al. 2012). PRODIMO has been designed to run within a few CPU-hours per model such that automatic fittings of the observed continuum emission and of thousands of gas lines are feasible. This 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 self-consistently 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 modeling of the MRI-turbulence.

In our chemical modeling, 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 fields calculated from 2D continuum radiative transfer (van Dishoeck & Visser 2014). 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, which 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 disk s. Most of the features concern a better treatment of the charging of the PAHs and dust grains.

2.2 Gas phase chemistry

The gas-phase chemistry includes photodissociation, ion-neutral, neutral-neutral, as well as a few colliders and three-body reactions. The chemical network is discussed in detail in Kamp et al. (2017).

At high gas temperatures, ionization by collisions with hydrogen atoms or molecules and with electrons is possible, $$\text{H}+{\text{H}}_2\to {\text{H}}^{-}+{\text{H}}_2^{+},$$(1)

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⁻⁴. Collisional ionization of metals can also occur $$\text{M}+\text{N}\to {\text{M}}^{-}+\text{N}+{\text{e}}^{-},$$(2)

where M is a neutral metal (e.g., Fe, Mg, S, Na, Si) and N is H, H₂, or He.

2.3 PAH charge exchange chemistry

Polycyclic aromatic hydrocarbonss 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. (2017) and Thi et al. (2019).

In protoplanetary disks, their abundances are lower by a factor f_PAH = 10⁻¹–10⁻³ compared to their interstellar abundance of 3 × 10⁻⁷ (Tielens 2008). We chose to use the peri-condensed circumcoronene (C₅₄ H₁₈) as typical PAHs, which are large enough to remain unaffected by photodissociation in disks around Herbig Ae 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²⁺, PAH³⁺; see Table 1). The effective radius of a PAH is computed as (Weingartner & Draine 2001b) $$a_{\text{PAH}}={10}^{-7}{\left(\frac{N_{\text{C}}}{468}\right)}^{1/3}\text{\hspace{0.17em}cm},$$(3)

where N_C is the number of carbon atoms in the PAH. The radius for the circumcoronene is a_PAH (C₅₄H₁₈) = 4.686 × 10⁻⁸ cm. The PAH ionization potential can either be taken from the literature when they are measured or estimated using (Weingartner & Draine 2001b) $${\text{IP}}_{\text{PAH}}=W_0+(Z_{\text{PAH}}+0.5)\frac{e^2}{a_{\text{PAH}}}+(Z_{\text{PAH}}+2)\frac{e^2}{a_{\text{PAH}}}\frac{0.3\times {10}^{-8}}{a_{\text{PAH}}}\text{\hspace{0.17em}erg},$$(4)

where W₀ is the work function assumed to be 4.4 eV (7.05 × 10⁻¹² erg) and Z_PAH is the charge of the PAH. The ionization potentials (IP) 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 (e.g., 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. (2017).

On disk surfaces, PAHs 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 a PAH. A detailed discussion on PAH charging is provided in Thi et al. (2019).

2.4 Dust grain charging

Dust grains with radius can be major charge carriers in the interstellar medium in general, and in protoplanetary disk s 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 few differences. We considered an average grain of radius ⟨a⟩ ≡ a and a geometric cross section σ_dust = πa². 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 ${\sum }_{Z_{\text{min}}}^{Z_{\text{max}}}f(Z_{\text{d}})=1$.

2.4.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₀ for bulk silicate is typically assumed to be 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) $$\text{IP}=W_{\text{eff}}=W_0+W_{\text{c}},$$(5)

where W_c is the extra work function defined by $$W_{\text{c}}=\left(Z_{\text{d}}+\frac{1}{2}\right)\frac{e^2}{a},$$(6)

where e is an elementary charge.When the grain is positively charged (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).

2.4.2 Photoejection and photodetachment

An electroncan 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 negatively charged grains. The photodetachment and photoejection 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 Eq. (17) of Weingartner & Draine (2001b) $${\eta }_{\text{b}}(\Theta )=\frac{0.5(\Theta /W_0)}{1+5(\Theta /W_0)},$$(7)

where W₀ 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 Eq. (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 $${\eta }_{\text{eff}}(x,\lambda )=(1-{\text{e}}^{-0.025(x-x_0)})\times {10}^{-0.05(\lambda -{\lambda }_0)-1.3},$$(8)

where x₀ = 2.5 and λ₀ = 120 nm. The equation closely reproduces the experimental data. 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{\nu }_{\text{pe}}=W_{\text{eff}}.$$(9)

When a grain is negatively charged the photodetachment threshold energy is (Weingartner & Draine 2001b) $$h{\nu }_{\text{pd}}(Z_{\text{d}}<0)=EA(Z_{\text{d}}+1,a)+E_{\text{min}}(Z_{\text{d}},a),$$(10)

where the electron affinity is $$EA(Z_{\text{d}},a)=W_0-E_{\text{bg}}+\left(Z_{\text{d}}+\frac{1}{2}\right)\frac{e^2}{a}.$$(11)

The band gap E_bg = 0 for metals and semimetals, while it can be several eVs for other materials. Weingartner & Draine (2001b) assumed (W₀ − 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).

The value for E_min follows the definition of Weingartner & Draine (2001b) $$E_{\text{min}}(Z_{\text{d}}<0,a)=-(Z_{\text{d}}+1)\frac{e^2}{a}{\left[1+{\left(\frac{27\AA }{a}\right)}^{0.75}\right]}^{-1}.$$(12)

For a singly negatively charged grain, assuming Z_d = −1 and a = 1 μm, one gets 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 $$k_{\text{pe}}=\pi a^2{\displaystyle {\int }_{{\nu }_{\text{pe}}}^{{\nu }_{\text{max}}}{\eta }_{\text{eff}}}{Q}_{\text{abs}}{J}_{\nu }\text{d}\nu +\pi a^2{\displaystyle {\int }_{{\nu }_{\text{pd}}}^{{\nu }_{\text{max}}}{\eta }_{\text{pd}}}{Q}_{\text{abs}}{J}_{\nu }\text{d}\nu ,$$(13)

where Q_abs is the frequency-dependent absorption efficiency and J_ν is the specific mean intensity at the frequency ν computed by the continuum radiative transfer.

Fig. 1 Analytical fits to the Abbas et al. (2006) experimental data shown as stars, triangles, and squares at wavelengths 120, 140, and 160 nm, respectively.

2.4.3 Electron attachment

The electron attachment rate coefficient onto neutral grains is $$k_{\text{e},\text{n}}=n_{\text{e}}{S}_{\text{e}}\sqrt{\frac{8k{T}_{\text{e}}}{\pi m_{\text{e}}}}{\sigma }_{\text{dust}}[{\text{electrons s}}^{-1}],$$(14)

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 higher than 0.3.

For positively charged grains (Z_d > 0) the electron recombination rate is enhanced by Coulomb attraction $$k_{\text{e},\text{+}}=k_{\text{e},\text{n}}\left(1+\frac{W_{\text{c}}}{kT}\right),$$(15)

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 (Z_d < 0): $$k_{\text{e},\text{-}}=k_{\text{e},\text{n}}\exp \left(\frac{W_{\text{c}}}{kT}\right).$$(16)

2.4.4 Thermionic emission

We considered the thermionic emission of electrons by hot dust grains following the Richardson–Dushman theory (Ashcroft & Mermin 1976; Sodha 2014) $$k_{\text{RD}}=\frac{4\pi m_{\text{e}}{k}^2}{h^3}{T}_{\text{d}}^2(1-r)\exp \left(-\frac{W_{\text{therm}}}{k{T}_{\text{d}}}\right){\sigma }_{\text{dust}}[{\text{electronss}}^{\text{-1}}],$$(17)

where T_d is the dust temperatureand 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.

2.4.5 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. $$k_{\text{cd}}^{\text{d}}=n_{\text{n}}{\sigma }_{\text{dust}}\sqrt{\frac{8kT}{\pi m_{\text{n}}}}\exp \left(-\frac{W_{\text{cd}}}{kT}\right),$$(18)

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.

2.4.6 Charge exchange between ions and dust grains

Gas-phase species and dust grains can exchange charges. Weingartner & Draine (2001a) considered the effect of ion-charged 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): $${\text{gr}}^{n-}+X_{\text{g}}{}^{+}\to {\text{gr}}^{(n-1)-}+X_{\text{s}}+\Delta E_{\text{ads}}(X_{\text{g}})+\Delta E_{\text{Coulomb}},$$(19) $${\text{gr}}^{n-}+X_{\text{s}}{}^{+}\to {\text{gr}}^{(n-1)-}+X_{\text{s}}+\Delta E_{\text{rec}}+\Delta E_{\text{relax}},$$(20)

where $X_{\text{g}}^{+}$ and $X_{\text{s}}^{+}$ are gas-phase and surface ions, respectively, or $${\text{gr}}^{n-}+X_{\text{s}}{}^{+}\to {\text{gr}}^{(n-1)-}+X_{\text{g}}+\Delta E_{\text{rec}}-\Delta E_{\text{ads}}(X_{\text{g}}).$$(21)

The gas-phase species $X_{\text{g}}^{+}$ acquires the approach energy, which is the sum of the adsorption energy (ΔE_ads(X_g)) and the Coulomb energy ΔE_Coulomb. The parameter Δ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 transferred 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 in an excited species: $${\text{gr}}^{n-}+{\text{AH}}_{\text{g}}^{+}\to {\text{gr}}^{(n-1)-}+{\text{AH}}_{\text{s}}^{*}+\Delta E_{\text{rec}}.$$(22)

When the excess energy of the electronically excited neutral species AH^* cannot transfer efficiently to the grain surface or be radiated away (because the species has a low dipole moment), the recombination is dissociative, $${\text{AH}}_{\text{s}}^{*}\to {\text{A}}_{\text{g}}+{\text{H}}_{\text{g}}+\Delta E_{\text{gdiss},\text{rec}},$$(23)

where the molecule dissociates into species that immediately leave the surface carrying the excess energy as kinetic energy or where one of the products remains on the grain surface $${\text{AH}}_{\text{s}}^{*}\to {\text{A}}_{\text{s}}+{\text{H}}_{\text{g}}+\Delta E_{\text{sdiss},\text{rec}}.$$(24)

Here the heavier product A tends to easily transfer the excess energy to the surface or possesses many degrees of freedom and, therefore, can efficiently radiate 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) $${\text{gr}}^{n+}+X_{\text{g}}{}^{+}\to {\text{gr}}^{(n+1)+}+X_{\text{g}}.$$(25)

For large silicate grains, Aikawa et al. (1999) performed classical computation for the recombination of HCO⁺ with negatively charged grains and concluded that the recombination results in the dissociation of the molecules. Another example is given by $$\text{gr}+{\text{NH}}^{+}\to {\text{gr}}^{+}+\text{N}+\text{H}.$$(26)

Both recombination and charge exchange reactions proceed with the rate $$k_{\text{gr},\text{ion}}=n_{\text{ion}}{S}_{\text{ion}}\sqrt{\frac{8k{T}_{\text{ion}}}{\pi m_{\text{ion}}}}{\sigma }_{\text{dust}}\max \left(0,1-\frac{Z_{\text{d}}{e}^2}{ak{T}_{\text{ion}}}\right){\text{e}}^{-E_{\text{therm}}/k{T}_{\text{ion}}},$$(27)

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 $$\begin{array}{ccc}{E}_{\text{therm}}& =& \Delta E_{\text{rec}}-\Delta E_{\text{ads}}(X_{\text{g}})-\Delta E_{\text{Coulomb}}\\ & \simeq & \max \left[I{P}_{\text{gr}}-I{P}_{\text{neu}},0\right]\end{array}$$(28)

for neutral and positively charged grains. We assume that the adsorption and Coulomb energy are negligible. For negatively charged grains $$E_{\text{therm}}=\max \left[\left(\text{EA}(Z_{\text{d}}+1,a)+E_{\text{min}}(Z_{\text{d}},a)\right)-I{P}_{\text{neu}},0\right]\text{},$$(29)

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 $$E_{\text{therm}}=\max \left[5.89-{\text{IP}}_{\text{neu}},0\right]\text{\hspace{0.17em}in}\text{\hspace{0.17em}eV}.$$(30)

The reality is probably much more complex. The Coulomb interaction acquired by the approaching ion can lead to the tunneling 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, $${\text{gr}}^{n+}+{\text{X}}_{\text{g}}\to {\text{gr}}^{(n-1)+}+{\text{X}}_{\text{g}}^{+}.$$(31)

The positivegrains 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 $$\begin{array}{ll}{k}_{{\text{gr}}^{+},\text{n}}=\hfill & n_{\text{n}}\times \max \left[2\pi \left|Z_{\text{d}}\right|e\sqrt{\frac{{\alpha }_{\text{pol}}}{\mu }},S_{\text{n}}\sqrt{\frac{8k{T}_{\text{n}}}{\pi m_{\text{n}}}}{\sigma }_{\text{dust}}\right]\hfill \\ \hfill & \text{}\text{}\times {\text{e}}^{-E_{\text{therm}}/k{T}_{\text{n}}},\hfill \end{array}$$(32)

where the reduced mass is basically that of the gas-phase species μ =m_n and the polarizability is α_pol = 10⁻²⁴, where the barrier term becomes $$E_{\text{therm}}=\text{max}\left[I{P}_{\text{neu}}-I{P}_{\text{gr}},0\right]$$(33)

in practically all cases, the geometrical rate dominates, $$\begin{array}{l}{k}_{{\text{gr}}^{+},\text{n}}=n_{\text{n}}{S}_{\text{n}}\sqrt{\frac{8k{T}_{\text{n}}}{\pi m_{\text{n}}}}{\sigma }_{\text{dust}}{\text{e}}^{-E_{\text{therm}}/k{T}_{\text{n}}}.\hfill \end{array}$$(34)

Anions can also react with positively charged grains $${\text{gr}}^{n+}+X_{\text{g}}^{-}\to {\text{gr}}^{(n-1)+}+X_{\text{g}}.$$(35)

However, we did not have anions in the chemical network at this current stage. Table 2 gives a few examples of species ionization potentials. Interestingly, charge exchange between neutral dust grains and the chemically important molecular ions H₃ O⁺ and NH${}_4^{+}$ 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).

2.4.7 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, W₀ = 8 eV, a maximum photon energy of hν_max = 13.6 eV, an elementary charge e in cgs of 4.803206815 × 10⁻¹⁰, and the conversion 1 eV = 1.60217657 × 10⁻¹² erg, the maximum charge (assuming that photoejection is the dominant electron ejection process) reads $$Z_{\text{max}}=\left(h{\nu }_{\text{max}}-W_0\right)a/e^2-0.5\simeq 3884\left(\frac{a}{\mu \text{m}}\right).$$(36)

For the minimum grain charge, we adopted again the formalism of Weingartner & Draine (2001b) $$Z_{\text{min}}=\text{int}\left(\frac{U_{\text{ait}}}{14.4}\frac{a}{\AA }\right)+1,$$(37)

where $$\frac{U_{\text{ait}}}{\text{V}}\simeq \{\text{}\text{}\text{}\begin{array}{ll}3.9+1200(a/\mu \text{m})+0.0002(\mu \text{m}/a)\hfill & \text{for carbonaceous}\hfill \\ 2.5+700(a/\mu \text{m})+0.0008(\mu \text{m}/a)\hfill & \text{for silicate}.\hfill \end{array}$$(38)

For a silicategrain 1 micron in radius, U_ait ≃ 702 V and Z_min ≃−487 499. 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 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 Eqs. (14) and (17), $$Z_{\text{d}}\simeq \frac{1}{2}\frac{a}{e^2}(f-2)(W_0-E_{\text{bg}}+kT\ln B)-\frac{1}{2},$$(39)

where $$B=\frac{n_{\text{e}}{S}_{\text{e}}\sqrt{8k{T}_{\text{g}}/\pi m_{\text{e}}}}{(4\pi m_{\text{e}}{k}^2/h^3)T_{\text{d}}^2(1-r)}.$$(40)

Assuming that S_e = 1, r = 0.5, and T_d = T_g = T, $$B\simeq 5.27\times {10}^{-8}{T}^{-3/2}{n}_{\text{e}},$$(41)

and $$f={\left(1+{\left(\frac{2.7\times {10}^{-7}}{a}\right)}^{0.75}\right)}^{-1}.$$(42)

For dense regions with n_e = 10¹² cm⁻³ and T = 100 K, the term B is greater than unity. The term kTlnB ≃ 2.37 × 10⁻⁵Tln(5.27 × 10⁻⁸T^−3∕2n_e) eV is always negligible compared to W₀ (=8 eV), thus a strong lower limit to the grain charge can be derived: $$Z_{\text{d}}>\frac{1}{2}(f-2)\frac{a}{e^2}(W_0-E_{\text{bg}})-\frac{1}{2}.$$(43)

Assuming elementary charge in cgs of e = 4.80321 × 10⁻¹⁰ statC and W₀ = 8 eV, E_bg = 5 eV, and f ~1, we obtain $$Z_{\text{d}}>-2107\left(\frac{a}{\mu \text{m}}\right).$$(44)

In this study, we did not include photoejection due to X-ray photons, and we adopted in our numerical models the limits $$\left|Z_{\text{d}}\right|<4000\left(\frac{a}{\mu \text{m}}\right).$$(45)

2.5 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), we 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) $$\left(1-\frac{Z_{\text{d}}{e}^2}{akT}\right)\exp -\left(\frac{(Z_{\text{d}}+0.5)e^2}{akT}\right)\frac{n_{\text{ion}}}{n_{\text{elec}}}=\sqrt{\frac{m_{\text{ion}}}{m_{\text{elec}}}},$$(46)

where we can set x = Z_de²∕akT. Using the medium global neutrality n_elec = n_ion + Z_dn_d, the equation becomes $$(1-x){\text{e}}^{-x}\simeq \sqrt{\frac{m_{\text{ion}}}{m_{\text{elec}}}}\left(1+x\frac{n_{\text{d}}}{n_{\text{ion}}}\frac{akT}{e^2}\right).$$(47)

A representation of the ratio between charges on grain surfaces and free electrons in an obscured region is shown in Fig. 2 for a grain 1 micron in radius. 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 $$Z_{\text{d}}~-6(2.7+0.4\ln (n_{\text{nucl}}))\left(\frac{a}{\mu \text{m}}\right)\left(\frac{T}{100\text{K}}\right),$$(48)

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 $$Z_{\text{d}}~-23\left(\frac{a}{\mu \text{m}}\right)\left(\frac{T}{100\text{\hspace{0.17em}K}}\right).$$(49)

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’s negative charge increases with the gas temperature.

Using the limit on the number of negative charges on a dust grain (see Eq. (44)), we can derive the maximum fractional abundance of negative charges on grains $$\chi (\text{d})=\frac{\left|Z_{\text{d}}\right|n_{\text{d}}}{n_{\langle \text{H}\rangle }}\le 6\times {10}^{-12}\left(\frac{\mu {\text{m}}^2}{a^2}\right)\left(\frac{100}{gd}\right).$$(50)

We can also use the estimate of the dust charge in obscured disk regions (Eq. (49)). The estimated fractional abundance of negative charges on grains becomes $$\chi (\text{d})=\frac{\left|Z_{\text{d}}\right|n_{\text{d}}}{n_{\langle \text{H}\rangle }}~6.7\times {10}^{-14}\left(\frac{\mu {\text{m}}^2}{a^2}\right)\left(\frac{100}{gd}\right)\left(\frac{T}{100\text{\hspace{0.17em}K}}\right).$$(51)

This last equation together with Fig. 2 shows that, unless the ionization fraction is below 10⁻¹³, the negativecharges on large silicate grains (with a radius greater than 1 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, 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 greater detail in Appendix C. When photon-induced electron ejection is present, the number of negative charges on grains will be even lower.

Fig. 2 Ratio of the number of negative charges on grains to the free electrons for agrain 1 micron in radius, and different total ionization fraction (from 10−14 to 10−8) as a function of the gas kinetic temperature.

Fig. 3 Variation in the value of αideal a function of the value of βmag for two prescriptions. In this paper we adopted δ = 0.5.

3 MRI-driven turbulence prescription for hydrostatic protoplanetary disk models

In this section we describe our model of non-ideal MHD driven turbulence gas heating and cooling and line broadening.

3.1 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 $${\alpha }_{\text{ideal}}={\left(\frac{2}{{\beta }_{\text{mag}}}\right)}^{\delta },$$(52)

where δ is a parameter between 0.5 and 1. The β_mag (beta magnetic) parameter is related to other gas parameters $${\beta }_{\text{mag}}(r,z)\equiv \frac{P_{\text{therm}}(r,z)}{P_{\text{mag}}(r,z)}=8\pi \frac{\rho c_{\text{s}}^2}{B_{\text{z}}{(r)}^2}=2{\left(\frac{c_{\text{s}}}{{\text{v}}_{\text{A}}}\right)}^2=\frac{2}{3}\frac{E_{\text{th}}}{E_{\text{mag}}},$$(53)

where E_th and E_mag are the thermal and magnetic energy respectively, ρ is the gas mass density, $c_{\text{s}}=\sqrt{kT/{\mu }_{\text{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 $$v_{\text{A}}=\frac{B_{\text{z}}}{\sqrt{4\pi \rho }}.$$(54)

The values for α_ideal are bound between 2 × 10⁻³ and 0.1 for β_mag between 1 and 10⁶. A high value of β_mag means that thermal motion 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 {\alpha }_{\text{ideal}}=A+B{\tan }^{-1}\left[\frac{C-\log \beta }{D}\right],$$(55)

where the value for the constants are A = −1.9, B = 0.57, C = 4.2, and D = 0.5. Both prescriptions are plotted in Fig. 3 with δ = 0.5. We assume that β_mid does not vary with radius in the midplane. This is equivalent to assuming that the magnetic field is accreted with the gas (Guilet & Ogilvie 2014) $${\beta }_{\text{mag}}(r,z)={\beta }_{\text{mid}}\left(\frac{P_{\text{therm}}(r,z)}{P_{\text{therm}}(r,0)}\right).$$(56)

Here the constant β_mid is a free parameter of the disk models. The vertical component of the magnetic field varies with radius, but is constant in the vertical direction (see Fig. E.1). In both equations we set α_ideal, the ideal-MHD MRI value of α in the absence of resistivities, to a very low value (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 to be 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 following sections, we describe the different non-MHD resistivities before presenting the prescription for the value of α_eff that is used in the chemical disk models.

3.2 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⁻⁴ 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 higher 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 $$J=e{\sum }_j{n}_j{Z}_j{v}_j.$$(57)

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 $$n_j{Z}_{j}e\left(E+\frac{v_j}{c}\times B\right)=n_j{\gamma }_j\rho m_j{v}_j,$$(58)

where $${\gamma }_j\equiv \frac{{\langle \sigma v\rangle }_j}{m_j+m}$$(59)

is the drag coefficient (cm³ g⁻¹ s⁻¹) with ${\langle \sigma v\rangle }_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 ∑_jn_jZ_j = 0 we obtain $$\frac{J\times B}{c}={\sum }_j{n}_j{\gamma }_j\rho m_j{v}_j.$$(60)

The inversion of Eq. (58) leads to the generalized Ohm’s law, which characterizes the non-ideal MHD effects (Wardle & Ng 1999; Norman & Heyvaerts 1985) $$J={\sigma }_{\text{O}}{E}^{\prime }{}_{\parallel }+{\sigma }_{\text{H}}(\widehat{B}\times E^{\prime }{}_{\perp })+{\sigma }_{\text{P}}{E}^{\prime }{}_{\perp },$$(61)

where $E^{\prime }{}_{\parallel }$ and $E^{\prime }{}_{\perp }$ denote the decomposition of E′ into vectors parallel and perpendicular to the magnetic field B, respectively, and the ˆ means unit vector. The parameter σ_O is the ohmic conductivity, also referred to as the conductivity parallel to the field (σ_O= σ_∥), σ_H is the Hall conductivity, and σ_P is the Pedersen conductivity. In cgs units, the conductivities are in units of s⁻¹. 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 theturbulence provides energy to heat the gas, and turbulence affects the line widths and thus the line cooling.

3.2.1 Ohmic resistivity

We consider a neutral gas of number density n_⟨H⟩ (cm⁻³) at gas temperature T_gas (K) with a mean mass m_n (in grams) and mass density ρ_n (in grams cm⁻³). We define for each charged species j (j = e, atomic and molecular ions, PAH, gr, where gr stands for dust grains) the unit-free plasma β_j parameter (which should not be confused with β_mag) as $${\beta }_j\equiv \frac{\left|Z_j\right|eB}{m_{j}c}\frac{1}{{\gamma }_j{\rho }_{\text{n}}},$$(62)

where Z_j is the charge of species j of mass m_j (in grams); B is the magnetic field strength, for which we only consider the vertical component B_z (in units of gauss, Gs, with 1 Gs = 1 cm^−1∕2 g^1∕2 s⁻¹); and e is the elementary charge (in units of statC). Finally, ${\gamma }_j={\langle \sigma v\rangle }_j/(m_{\text{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 closely coupled with the neutral gas. On the other hand, if |β_j | ≫ 1 the charged species is closely coupled with the magnetic field. Using the terms introduced above, the ohmic electrical conductivity is defined as $${\sigma }_{\text{O}}\equiv \frac{ec}{B}{\sum }_j{n}_j\left|Z_j\right|{\beta }_j={\sum }_j{\sigma }_{\text{O},j}=\frac{e^2}{{\rho }_{\text{n}}}{\sum }_j\frac{Z_j^2{n}_j}{m_j{\gamma }_j}.$$(63)

The collisionrates ${\langle \sigma v\rangle }_{j\text{n}}$, which appear in the drag frequencies term γ_j, are thermal velocity-averaged values of the cross sections σ (without subscript). The electron-neutral collision rate is $${\langle \sigma v\rangle }_{\text{e},\text{n}}\approx 8.28\times {10}^{-10}\sqrt{T}{\text{cm}}^3{\text{s}}^{-1}.$$(64)

The atomic and molecular ions have average collision rates with neutral species following the Langevin rate $${\langle \sigma v\rangle }_{\text{ion},\text{n}}\approx 1.9\times {10}^{-9}{\text{cm}}^3{\text{s}}^{-1}.$$(65)

Specifically for the HCO⁺-H₂ system, Flower (2000) provided the equation $${\langle \sigma v\rangle }_{\text{ion},\text{n}}\approx 8.5\times {10}^{-10}{T}^{0.24}{\text{cm}}^3{\text{s}}^{-1},$$(66)

which shows a weak temperature dependency. At 100 K, the rate becomes 2.5 × 10⁻⁹. For PAH ions, we use the largest value between the temperature-independent Langevin and the geometrical rate $${\langle \sigma v\rangle }_{\text{PAH}\text{\hspace{0.17em}ion},\text{n}}\approx \text{max}\left[1.9\times {10}^{-9},\pi a_{\text{PAH}}^2{\left(\frac{8k{T}_{\text{gas}}}{\pi m_{\text{n}}}\right)}^{1/2}\right]{\text{cm}}^3{\text{s}}^{-1},$$(67)

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 $${\langle \sigma v\rangle }_{\text{dust},\text{n}}\approx \pi a^2{\left(\frac{8k{T}_{\text{gas}}}{\pi m_{\text{n}}}\right)}^{1/2}{\text{cm}}^3{\text{s}}^{-1}.$$(68)

We assume that the grains are basically immobile with respect to the gas. For a fully molecular gas m_n ≈ 2.2 amu = 3.65317 × 10⁻²⁴ grams, $${\langle \sigma v\rangle }_{\text{dust},\text{n}}\approx 3\times {10}^{-3}\left(\frac{a^2}{\mu {\text{m}}^2}\right){\left(\frac{T_{\text{gas}}}{100\text{\hspace{0.17em}K}}\right)}^{1/2}{\text{cm}}^3{\text{s}}^{-1}.$$(69)

The ohmic diffusivity is characterized by the dimensionless Elsasser number, which is defined as the ratio of the Lorentz to the Coriolis force, $${\Lambda }_{\text{Ohm}}\equiv \frac{B_{\text{z}}^2}{4\pi \rho {\eta }_{\text{O}}\Omega }\equiv \frac{v_{\text{A}}^2}{{\eta }_{\text{O}}\Omega }\equiv \left(\frac{4\pi {\sigma }_{\text{O}}}{\Omega }\right){\left(\frac{v_{\text{A}}}{c}\right)}^2,$$(70)

where $$\Omega =\sqrt{\frac{G{M}_{*}}{r^3}}\simeq 2\times {10}^{-7}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{1/2}{\left(\frac{r}{\text{au}}\right)}^{-3/2}\text{rad}{\text{s}}^{-1}$$(71)

is the angular speed, which in the case of a protoplanetary disk is 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 are respectively the distance in the disk from the star in cm and in au. One astronomical unit (1 au) corresponds to 1.4959787 × 10¹³ cm. The ohmic resistivity is η_O = (c²∕4π)∕σ_O. The ohmic Elsasser number can be rewritten as the sum of the contribution of each charged species to the ohmic conductivity: $${\Lambda }_{\text{Ohm}}\equiv \frac{4\pi }{\Omega }{\left(\frac{v_{\text{A}}}{c}\right)}^2({\sigma }_{\text{e},\text{O}}+{\sigma }_{\text{ions},\text{O}}+{\sigma }_{\text{PAH},\text{O}}+{\sigma }_{\text{dust},\text{O}}).$$(72)

For MRI-driven turbulence to fully develop, the turbulence development timescale τ_ideal has to be shorter than the damping by ohmic diffusion timescale τ_damp: $${\tau }_{\text{ideal}}~\frac{h}{V_{\text{A}}}<{\tau }_{\text{damp}}~\frac{h^2}{{\eta }_{\text{O}}}.$$(73)

This criterion translates to $$\frac{h{\text{v}}_{\text{A}}}{{\eta }_{\text{O}}}>1.$$(74)

In the approximation of vertical isothermal disk (h = c_s∕Ω) and in the condition of v_A = c_s (β_mag = 2, α_ideal = 1), the criterionbecomes $${\Lambda }_{\text{Ohm}}\equiv \frac{v_{\text{A}}^2}{{\eta }_{\text{O}}\Omega }>1.$$(75)

This simple estimate is supported by non-ideal MHD simulations (Sano & Stone 2002). Although the disks are not isothermal in the vertical direction, we use this criterion in our models. The ohmic Elsasser number criterion can be seen as the rate of dissipation of the magnetic energy (~B²) compared to the energy dissipated by ohmic resistivity.

Fig. 4 max(α) as a function of the ambipolar diffusion term Am.

3.2.2 Ambipolar diffusion

The difference between the ion and neutral (atoms and mole- cules, PAHs, and dust grains) velocities is responsible for ambi- polar diffusion (Bittencourt 2004). The Elsasser number for ambipolar diffusion is $$\text{Am}\equiv \frac{v_{\text{A}}^2}{{\eta }_{\text{A}}\Omega },$$(76)

where the ambipolar resistivity η_A is (Wardle & Ng 1999) $${\eta }_{\text{A}}=D^2\left(\frac{c^2}{4\pi }\frac{{\sigma }_{\text{P}}}{{\sigma }_{\perp }^2}-{\eta }_{\text{O}}\right),$$(77)

where σ_P the Pedersen conductivity, σ_⊥ is the total conductivity perpendicular to the field, and η_O is the ohmic 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 $${\sigma }_{\perp }=\sqrt{{\sigma }_{\text{H}}^2+{\sigma }_{\text{P}}^2},$$(78)

where the Pedersen conductivity is $${\sigma }_{\text{P}}\equiv \frac{ec}{B}{\sum }_j{n}_j{Z}_j\left(\frac{{\beta }_j}{1+{\beta }_j^2}\right).$$(79)

The Hall conductivity σ_H is defined by $${\sigma }_{\text{H}}\equiv \frac{ec}{B}{\sum }_j\frac{n_j{Z}_j}{1+{\beta }_j^2}.$$(80)

When all the charged species are closely coupled to the neutral gas (|β_j |≪ 1 for all j), then σ_O ≈ σ_P ≫|σ_H|, the conductivityis scalar (isotropic), and the regime is resistive. When all the charged species are closely coupled to the field, then σ_O ≫ σ_P ≫|σ_H| and ambipolar diffusion dominates (Wardle & Ng 1999). The ambipolar resistivity becomes $${\eta }_{\text{A}}=D^2\left({\eta }_{\text{P}}-{\eta }_{\text{O}}\right)\approx D^2{\eta }_{\text{P}}.$$(81)

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: $$\text{max}(\alpha )=\frac{1}{2}{\left[{\left(\frac{50}{{\text{Am}}^{1.2}}\right)}^2+{\left(\frac{8}{{\text{Am}}^{0.3}}+1\right)}^2\right]}^{-1/2}.$$(82)

The maximum possible value for α is shown in Fig. 4. For Am = 1, max(α) ≃ 10⁻². It should be noted that the equation provides an upper limit on the value of α due to ambipolar diffusion.

The ambipolar diffusion prevents the MRI from fully developing 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.

4 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 equation 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 ohmic Elsasser number Λ_Ohm is larger than one: $${\alpha }_{\text{Ohm}}={\alpha }_{\text{ideal}}\times \text{min}\left(1,\frac{v_{\text{A}}^2}{{\eta }_{\text{O}}\Omega }\right)={\alpha }_{\text{ideal}}\times \text{min}\left(1,{\sigma }_{\text{O}}{v}_{\text{A}}^2\Omega \right).$$(83)

The ohmic Elsasser 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; however, numerical simulations of the Hall diffusivity give a more complex picture (Lesur et al. 2014). 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 $$\begin{array}{lll}{\alpha }_{\text{eff}}\hfill & =\hfill & {\alpha }_{\text{ideal}}({\beta }_{\text{mag}})\times \text{min}\left(1,\frac{v_{\text{A}}^2}{{\eta }_{\text{O}}\Omega }\right)\hfill \\ \hfill & \hfill & \times {\left[{\left(\frac{50}{{\text{Am}}^{1.2}}\right)}^2+{\left(\frac{8}{{\text{Am}}^{0.3}}+1\right)}^2\right]}^{-1/2},\hfill \end{array}$$(84)

when $\sqrt{{\beta }_{\text{mag}}}{\Lambda }_{\text{Ohm}}>1$, and α_eff ≃ 0 otherwise. It has been shown that, when $\sqrt{{\beta }_{\text{mag}}}{\Lambda }_{\text{Ohm}}<1$, all the perturbation modes are stabilized and thus no MRI-driven turbulence can exist (Jin 1996). Our equation 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, the thermal balance of gas and dust, and radiative transfer. While both ohmic and ambipolar diffusion contributions to the effective turbulence depend on the available charges, the two effects differ on the other dependencies. The ohmic Elsasser number 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 theelectron recombination depends on the square of the density (∝ k_recn_ionn_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 then the ambipolar diffusion Elsasser number will decrease.

To illustrate these 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¹⁰ cm⁻³ and 10¹² cm⁻³). The panels show the effects of the variation in the ohmic and ambipolar diffusivities on the value of α_eff and α_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 restriction on α_eff occurs for a high value of β_mag. The value of α_eff is higher than 10⁻³ down to a total ionization fraction of ~10⁻⁸–10⁻¹⁰ depending on the gas density. At disk surfaces, where the ionization fraction is high and the value of β_mag low, the efficiency of MRI turbulence is probably limited by the ambipolar diffusion resistivity to a maximum value of ~0.25.

Fig. 5 Effective turbulence coefficient αeff as a functionof βmag for gas of density of n⟨H⟩ = 1010 cm−3 and 1012 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). Left panels: ideal-MHD, ohmic, and ambipolar diffusion contributions to αeff, while two right panels: resulting αeff including the all-mode damping criterion of Jin (1996). Lines at 10−10 in the right panels mean no MRI-driven turbulence for a total ionization fraction of 10−11 for all values of βmag.

5 Disk model

We chose to model a typical protoplanetary disk with a total gas+solid mass of 10⁻² 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 size distribution of the dust grains 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⁻⁷). 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.7×10⁻¹⁷ s⁻¹. Future studies will include varying the X-ray fluxes, lowering the cosmic ray flux and its attenuation in disks (Cleeves et al. 2013, 2014), and modeling 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 modeled the disks with β_mid in the disk midplane from 10² to 10⁶. For comparison, we also ran a model with no (MRI) turbulence, i.e., the line broadening is purely thermal (v_turb = 0 km s⁻¹) and other models with a constant turbulence width v_turb = 0.05, 0.1, 0.15, and 0.2 km s⁻¹ and no turbulence heating throughout the disk (passive disk models). The central object is either a T Tauri star (M_* = 0.7 M_⊙, T_eff = 4000 K, L_* = 1.0 L_⊙) or a Herbig Ae 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) $$E_{\text{acc}}=\frac{9}{4}{\alpha }_{\text{eff}}{P}_{\text{gas}}\Omega$$(85)

and is included as a gas heating agent in the heating-cooling balance to determine the gas temperature. McNally et al. (2014) modeled the energy transfer from the turbulence decay to the gas.

The turbulence will effect the gas temperature, which in turn will change the chemistry, the value of β_mag, and the ohmic and ambipolar Elsasser numbers. In an analytical analysis in the case of ideal MRI-turbulence, the energy release is $$E_{\text{acc}}\propto P_{\text{therm}}^{1-\delta }\Omega .$$(86)

Adopting δ = 0.5, $$E_{\text{acc}}\propto \sqrt{P_{\text{therm}}}\Omega .$$(87)

Assuming that the pressure term can be described by a barotropic law $$P_{\text{therm}}\propto n^{\gamma },$$(88)

the gas heating rate is $$E_{\text{acc}}\propto n_{\langle \text{H}\rangle }^{\gamma /2}\Omega .$$(89)

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. The value of α_ideal is proportional to T^−δ, Λ_Ohm is proportional to $\sqrt{T}$. The dependence of Am on T is not direct, but assuming electronic recombination rates ∝ T^−1∕2, an increase in the gas temperature should result in a decrease in 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 subsonic and is related in our model to α_eff and the sound speed c_s by v_turb = $\sqrt{{\alpha }_{\text{eff}}}$ c_s (Hughes et al. 2011). Another relationship such as v_turb = α_eff c_s has been proposed (Simon et al. 2013). The local line width becomes $\Delta v={(v_{\text{th}}^2+v_{\text{turb}}^2)}^{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 affects not only the heating, butalso the cooling of the gas.

The present model only considers viscous accretion through the alpha parameter. However, in the presence of a poloidal magnetic field, it is known that discs can be subject to MHD wind which can also drive accretion (e.g., Bai & Stone 2013). We have purposely neglected this process due to the lack of robust prescriptions for wind-driven accretion.

6 Disk model results

The disk total charge (also known as the ionization fraction) is shown in the upper left panel of Fig. 6 for a disk model with β_mid = 10⁴. 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. Figure 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⁻⁴) 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⁻⁴). As carbon becomes neutral at lower altitudes, 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⁻¹²), the radiation field is too weak to eject the excess 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 ohmic Elsasser number, the ambipolar diffusion number, and the complete mode damping criterion are shown in Fig. 8 for the disk models with β_mid = 10⁴ (left panels) and β_mid = 10⁶ (right panels). The ohmic Elsasser number is the most sensitive to the value of β_mag. The ambipolar Elsasser number Am distributions are similar for β_mid = 10⁴ and β_mid = 10⁶. In unobscured regions, the electron fractional abundance χ(e) is relatively high in the range 10⁻⁸–10⁻⁴. Due to the limit on negative charges on grains (Eq. (50)), in a gas with total ionization fraction greater than a few 10⁻¹², the free electrons will be the main contributor to the ohmic Elsasser conductivity. In the inner disk midplane, the total charge is low, resulting in low values for the ohmic and ambipolar Elsasser numbers.

We can derive a good approximation to the ohmic Elsasser number in disk regions where σ_O ≃ σ_e,O and χ(e) > 10⁻¹³: $${\Lambda }_{\text{Ohm}}\simeq 1{\left(\frac{T}{100}\right)}^{1/2}\left(\frac{{10}^4}{{\beta }_{\text{mag}}}\right){\left(\frac{r}{\text{au}}\right)}^{3/2}\left(\frac{\chi (\text{e})}{{10}^{-9}}\right).$$(90)

The approximated ohmic Elsasser distribution for a disk model with β_mid = 10⁴ is shown in Fig. G.1.

The disk distribution for the effective turbulent parameter α_eff is shown in Fig. 9 for β_mid of 10⁴ (left panels) and 10⁶ (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 ohmic diffusion limited dead zone according to the total damping criterion $\text{Jin}\equiv {\beta }_{\text{mag}}^{1/2}{\Lambda }_{\text{Ohm}}$ = 1, the ohmic diffusion restricted MRI region (Λ_Ohm = 1), and the location of the disk where β_mag = 1 in the disk atmospheres. The decrease in α_eff with the disk radius stems from the decrease in the gas pressure, hence of β_mag (see middle panels of Fig. E.1). The outer disk midplane decrease in α_eff stems from the increase in β_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⁴, 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⁶ 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⁻¹. Figure 10 shows the main heating and cooling processes for the disk models with β_mid = 10⁴ (left panels) and β_mid = 10⁶ (right panels). In the dead zone where ${\beta }_{\text{mag}}^{1/2}{\Lambda }_{\text{Ohm}}<1$, gas heating occurs mainly via 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) up to disk radii of a few hundred au, where the density, the Keplerian rotation, and α_eff are low (seeEq. (85)). The outer disk is heated by the energy released by H₂ 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⁴ model can be heated by the photoelectric effect. The main gas heating processes of passive disks in the atmospheres are hit 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 is strong enough (or β_mid > 10⁶). 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 water vapor and molecular hydrogen rovibrational lines. In the viscous heating dominated region, 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 & Ádámkovics (2017) who modeled the effects of heat generated by turbulence decay and concluded that warm CO and H₂ O lines can trace the mechanically heated gas.

The actual CO line profiles computed in non-LTE for CO, ¹³CO, and C¹⁸O J = 2–1 and J = 3–2 are shown in Fig. 11 for the T Tauri and Herbig Ae disks seen with a 45° inclination. The line flux increases with more turbulent gas because of stronger heating and the 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 disk s. The CO profiles for the disk model with β_mid = 10⁶ and that for the model with no turbulence are similar for both transitions, J = 2–1 and J = 3–2. For the CO J = 2–1 and the CO J = 3–2 lines computed in models with a fixed value of v_turb, we can constrain v_turb down to 0.05 km s⁻¹ (see the red dotted lines in Fig. 11). Figure 12 shows the variation in the peak flux, FWHM, and the profile peak separation as function of β_mid for the CO lines for the T Tauri disk models. The ratio of the peak to the line center flux (peak-to-trough) for the T Tauri and Herbig Ae 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 of the peak to the line center flux, consistent with the results of the parametric disk model of Flaherty et al. (2015). The trend of decreasing peak-to-center ratio with increasing disk turbulence velocity shows the same amplitude (0.2–0.4) as the trend of a lower ratio for high values of β_mid. The peak-to-center flux ratios of the ¹³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 ¹²CO lines are emitted in the region where v_turb/c_s reaches its maximum. The optically thinner ¹³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¹⁸O lines can be optically thin. Therefore, the peak fluxes for the C¹⁸O lines are low and the peak-to-trough ratios for the C¹⁸O lines are less sensitive to the turbulence broadening than the optically thick ¹²CO and ¹³CO lines.

One possible parameter that influences the peak-to-trough ratios for the ¹³CO lines is the abundance of ¹³CO compared to¹²CO. The lower right panel in Fig. 13 shows the effect of a global ¹²CO/¹³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 ¹²CO/¹³CO in young stellar object (YSO) environments compared to the ISM value is consistent with the observations of the ¹²CO/¹³CO ratio towards a number of YSOs by Smith et al. (2015), who found a range of 86–165. The higher value for ¹²CO/¹³CO may reflect the less efficient ¹³CO self-shielding against photodissociation compared to ¹²CO (Miotello et al. 2016). Other parameters such as the disk inclination or the central star mass affect the value of the peak-to-trough ratios.

Fig. 6 Upperleft: total charge (positive or negative) abundance. The red dashed line corresponds to the location in the disk where the abundance of C+ and C are equal. Upper right: contribution of free electrons to the total charge. Lower left: contribution of the negative ions (H−). Lower right: contribution of the positive ions (H+). The disk model is the DIANA typical disk with βmid = 104.

Fig. 7 Average charge per grain (upper left) and its fractional contribution to the total charge (upper right). The location of the transition from ionized to neutral carbon is overplotted in the upper panels. The contribution of the positively charged PAHs (PAH+, PAH++, and PAH+++) is shown in the lower left panel, while the contribution of the negatively charged PAHs (PAH−) is shown in the lower right panel. The contours show the extinction and the location in the disk where the abundance of gas phase and frozen PAHs are equal. The model is the DIANA typical disk with βmid = 104. Silicate dust grains remains at Tdust < 1500 K at all disk heights (see Fig. F.1).

Fig. 8 Ohmic Elsasser number for the DIANA typical disk (first row). White contours correspond to the location of the total charge in the disk. Ambipolar diffusion Elsasser number is shown in the second row. The location where the C and C+ abundances are equal are overplotted in red. The stability criterion Eq. (74) for all modes to be damped, corresponding to a dead zone, is shown in the lower panels (Jin = 1). Left column:models with βmid = 104. Right column: models with βmid = 106.

Fig. 9 Effective αeff value and turbulence velocity over sound speed ratio for the DIANA typical disk and βmid = 104 (left panels) and βmid = 106 (right panels). The white contours in the upper panels show the location where the ohmic Elsasser 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 ${\beta }_{\text{mid}}^{1/2}{\Lambda }_{\text{Ohm}}<1$, MRI is entirely suppressed and αeff tends to zero. Also shown in the upper panels is the contour where β = 1. The contours in the lower panels indicate the gas-phase CO abundances.

Fig. 10 Main heating sources (upper panels) and main cooling sources (lower panels) at the different locations in the disk. Left panels: βmid = 104 disk model, while right panels: βmid = 106 disk model. Upper panels: the contour of the Jin’s criterion is shown in white contour. The black contours correspond to AV = 10 in all the panels.

Fig. 11 Line profiles of CO and isotopologues from the typical disk with different values of βmid (dashed lines), from the disk model without turbulence (continuous black lines), and from the model with a constant turbulent width of 0.05, 0.1, and 0.2 km s−1 (dotted lines). The disk is seen with an inclination of 45° (0° means that the disk is seen face-on).

Fig. 12 CO line profile statistics as a function of the value of βmid. Left panel: peak fluxes in Jy. Middle panel: FWHM in km s−1. Right panel: double-peaked profile peak separation in km s−1. The value for βmid = 107 corresponds to a model with no MRI turbulence in the entire disk.

7 Discussion

Gas and dust in passive disks are heated by the conversion of 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 in temperature from the midplane to the disk surfaces because of the extinction by dust grains for the UV and optical photons, and by 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 of the hot and relatively ionized upper disk layers above the dead zone to compensate for the decrease in mass accretion in the midplane. This result contrasts with vertical isothermal and relatively cool disk models where even extremely high values of α_eff are not enough to compensate for the low accretion in the midplane (Perez-Becker & Chiang 2011). Our modeling 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 fraction is high enough 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 gas-dust 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⁴ in the midplane. The computed turbulent line widths assuming β_mid = 10⁴ 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. (2015, 2017) also found a low value in the protoplanetary disk around the Herbig Ae 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 ¹²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 Herbig Ae disk model with 45° inclination. In the ALMA observations presented by Flaherty et al. (2015), the ¹²CO J = 3–2 and J = 2–1 profiles show higher peak-to-trough ratios than the ¹³CO profiles (see right panels of Fig. 13) in contrast to our model results. A higher ¹²CO/¹³CO ratio can decrease the peak-to-trough value slightly, which is still not sufficient to explain the observed value (see lower right panel of Fig. 13).

As an alternative to using high-quality CO lines, the ionization structure of protoplanetary disk s can be constrained by matching emission lines from molecular ions. Cleeves et al. (2015) argued using emission lines from HCO⁺ and N₂ H⁺ that the TW Hya disk has a dead zone extending to 50–65 au. Their disk model requires a cosmic ray ionization rate lower than 10⁻¹⁹ s⁻¹. Their dead zone is defined by the criterion R_EM ≤ 3300, where R_EM is the magnetic Reynolds number (Balbus & Terquem 2001). The molecular ions HCO⁺ and N₂ H⁺ have theoretical abundances that depend sensitively on the treatment of ionization and recombination in the gas (Rab et al. 2017) and on the N₂ self-shielding for N₂ H⁺ (Visser et al. 2018). In disk 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. (2017) showed that the N₂ 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₂ 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 ratio allowed them to constrain turbulence to between 0.4 c_s and 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⁻¹ 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 ohmic Elsasser number are then related by $$R_{\text{EM}}={\left(\frac{c_{\text{s}}}{v_{\text{A}}}\right)}^2{\Lambda }_{\text{Ohm}}$$(91)

assuming c_s = hΩ. Flaherty et al. (2017) used DCO⁺ lines in addition to CO and C¹⁸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⁵, i.e., that CO line profiles can be used to constrain α down to ~5 × 10⁻³. 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 submillimeter domain probe a large disk volume, but may not be sensitive to the location in the protoplanetary disk s 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).

Knowing the gas and dust temperature structure is central to estimating the thermal line broadening. Figure 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 gas-dust thermal accommodation at high densities, the inner disk midplane gas temperaturedoes not change much with accretion heating.

The charge in protoplanetary disk s 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, and sodium, among others, reaching the stellar photospheric abundances. These elements can still survive at high gas temperatures in the form of metal oxides like SiO, MgO, and FeO if UV is not present. However, since the main UV opacity, i.e., the dust, is no longer present, UV photons permeate the dust-free gas and can efficiently dissociate the oxides unless the efficient high-temperature formation of these oxides is possible, as it is for water vapor (Thi & Bik 2005). If the gas chemistry is at equilibrium, the Saha equation can be used 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₂ into H${}_2^{+}$ by an energetic ionization event (X-ray, cosmic rays, radioactive decay), and subsequently H${}_3^{+}$. They can recombine either with abundant atomic or molecular ions or with 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 disk s. 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, and grain properties such as the size distribution, settling, and drift. The effect of these key parameters can balance 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 X-ray 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 modeling 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 modeling of the dust vertical structure of HL Tau suggests that the dust vertical mixing is weak with a value of a few 10⁻⁴ for α_settle (Pinte et al. 2016). This corresponds to a disk model with β_mid ~ 10⁶.

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 α_settle and α_eff consistently, and to model the changes in the disk density structure, in the disk chemistry, and in the grain properties (e.g., size, drift) over the disk’s lifetime.

Fig. 13 Ratio of Peak flux to zero velocity flux for a standard T Tauri disk model (left panels) and the Herbig Ae disk model (right panels) as a function of the βmid parameter (upper panels) and of vturb (lower panels). All the models are viewed with an inclination of 45° (0° means that the disk is seen face-on).

8 Conclusions

We implemented a simple parametrized magnetorotational instability (MRI) driven turbulence model in the physico-chemical 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 anddust temperatures, chemistry, and PAH and grain charges are computed self-consistently together with the continuum and line radiative transfer. The disk models include both active and passive heating processes. The main conclusions are as follows:

• 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 ohmic resistivity governs the location of the dead zone in the inner disk midplane.

• The ambipolar diffusion resistivity does not greatly affect 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, by gas absorption of the dust grain infrared emissions, and by H₂ formation.

• The signatures of the gas turbulence in the CO emission lines are weak, and sophisticated modeling of the high signal-to-noise ratio (S/N) 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, and 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), our model outputs can be used as test cases to improve the reliability of inversion models used to derive α from the observations. Finally, in the future we will also explore emission lines from species other than CO as turbulence tracers.

Acknowledgements

We thank ANR (contract ANR-07-BLAN-0221) and PNPS of CNRS/INSU, France, for support. I.K., W.F.T., and P.W. acknowledge funding from the EU FP7-2011 under Grant Agreement no. 284405. We acknowledge discussions with Ch. Pinte and F. Ménard.

Appendix A Comparing the contributions of electrons, ions, and dust grains to the ohmic conductivity

The contribution of the electrons to the ohmic conductivity is $${\sigma }_{\text{e},\text{O}}\simeq \frac{e^2}{m_{\text{e}}{\langle \sigma v\rangle }_{\text{e},\text{n}}}\chi (\text{e}),$$(A.1)

where χ(e) is the free electron fractional abundance. The equation can be rewritten as $${\sigma }_{\text{e},\text{O}}\simeq 3\times {10}^{16}{\left(\frac{T}{100}\right)}^{-1/2}\chi (\text{e})[{\text{s}}^{-1}].$$(A.2)

One can make a comparison between the contribution of each charged species to the ohmic Elsasser number. For instance, the ratio between the contribution from the free electrons to the charged dust grains is of interest: $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}={\left(\frac{Z_{\text{e}}}{Z_{\text{d}}}\right)}^2\left(\frac{n_{\text{e}}}{n_{\text{d}}}\right)\left(\frac{m_{\text{d}}{\gamma }_{\text{d}}}{m_{\text{e}}{\gamma }_{\text{e}}}\right).$$(A.3)

For the electrons, $$m_{\text{e}}{\gamma }_{\text{e}}\simeq 8.28\times {10}^{-9}{T}_{100}^{1/2}\frac{m_{\text{e}}}{m_{\text{n}}},$$(A.4)

(m_e≪ m_n) and for the dust grains $$m_{\text{d}}{\gamma }_{\text{d}}\simeq 3\times {10}^{-3}{a}_{\mu \text{m}}^2{T}_{100}^{1/2},$$(A.5)

where the average mass of a grain is $$m_{\text{d}}=(4/3)\pi {\rho }_{\text{d}}{a}^3=1.256\times {10}^{-11}\left(\frac{\mu {\text{m}}^3}{a^3}\right)\text{\hspace{0.17em}grams}.$$(A.6)

Since m_d > > m_n we obtain $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}\simeq 3623.2{\left(\frac{1}{Z_{\text{d}}}\right)}^2\left(\frac{n_{\text{e}}}{n_{\text{d}}}\right)\left(\frac{m_{\text{n}}}{m_{\text{e}}}\right)a_{\mu \text{m}}^2.$$(A.7)

The average number of dust grains in the disk is $$n_{\text{d}}=\frac{2.2\text{\hspace{0.17em}amu}{n}_{\langle \text{H}\rangle }}{(4/3)\pi {\rho }_{\text{d}}{a}^{3}gd},$$(A.8) $$n_{\text{d}}\approx 2.907\times {10}^{-15}{n}_{\langle \text{H}\rangle }\left(\frac{\mu {\text{m}}^3}{a^3}\right)\left(\frac{100}{gd}\right),$$(A.9)

where gd is the gas-to-dust mass ratio. We have assumed a silicate mass density of 3.0 g cm⁻³. Replacing the dust grain number density by the equation above, we obtain $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}\simeq \frac{{10}^{18}}{Z_{\text{d}}^2}\left(\frac{n_{\text{e}}}{n_{\langle \text{H}\rangle }}\right)\left(\frac{m_{\text{H}}}{m_{\text{e}}}\right)\left(\frac{gd}{100}\right)a_{\mu \text{m}}^5$$(A.10)

For a gas-to-dust mass ratio of 100, $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}\simeq \frac{8.36\times {10}^{20}}{Z_{\text{d}}^2}\chi (\text{e})a_{\mu \text{m}}^5,$$(A.11)

where χ(e) is the relative abundance of free electrons. If we assume that grains have charge 0 > Z_d ≥ − Z_max and using the maximum possible negatively charged grains (Eq. (44)), the ratio becomes $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}\ge 5.2\times {10}^{13}\chi (\text{e})a_{\mu \text{m}}^3.$$(A.12)

Thus, for micron-sized dust grains, the free electron contribution to the ohmic Elsasser number always dominates over the dust contribution down to a free electron abundance of ~1.9 × 10⁻¹⁴. For 0.1 micron grains, the charged dust grain contribution to the ohmic Elsasser number becomes important for free electron abundances below ~1.9 × 10⁻¹¹. If we assumes that the dust charge in the UV-obscured region can be approximated by Eq. (49), the ohmic conductivity ratio is $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}\simeq 1.58\times {10}^{20}\left(\frac{100\text{\hspace{0.17em}K}}{T}\right)a_{\mu \text{m}}^4\chi (\text{e}).$$(A.13)

We note that the ratio is temperature-dependent.

Likewise, we can make a comparison between the contribution of the electrons compared to that of the ions $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{ions},\text{O}}}={\left(\frac{Z_{\text{e}}}{Z_{\text{i}}}\right)}^2\left(\frac{n_{\text{e}}}{n_{\text{i}}}\right)\left(\frac{m_{\text{i}}{\gamma }_{\text{i}}}{m_{\text{e}}{\gamma }_{\text{e}}}\right)$$(A.14)

Assuming singly charged ions and that most of the negative charges are carried by the electrons, $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{ions},\text{O}}}\simeq \frac{m_{\text{i}}{\gamma }_{\text{i}}}{m_{\text{e}}{\gamma }_{\text{e}}}\simeq 0.3{\left(\frac{100}{T}\right)}^{1/2}\frac{m_{\text{n}}}{m_{\text{e}}}\simeq 920{\left(\frac{100}{T}\right)}^{1/2},$$(A.15)

where $$m_{\text{i}}{\gamma }_{\text{i}}\simeq 1.9\times {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 Eq. (36). Equation (A.11) reads $$\frac{{\sigma }_{\text{e},\text{O}}}{{\sigma }_{\text{dust},\text{O}}}\simeq 5.5\times {10}^{13}{a}_{\mu \text{m}}^4\chi (\text{e}).$$(A.17)

Appendix B Effective turbulence coefficient

Figure B.1 shows the effective turbulence coefficient α_eff as a function of the total ionization fraction for a relatively dense gas at density 10¹⁰ cm⁻³ and temperature of 100 K. The four panels correspond to different values of β_mag (10⁴, 10³, 10², 10).

Fig. B.1 Effective turbulence coefficient αeft as a functionof the total ionization fraction for a gas of density n⟨H⟩ = 1010 cm−3 at T = 100 K and four values of the βmag parameters (104, 103, 102, 10). The values for αideal (red short-dashed lines) were derived assuming δ = 0.5. The limits imposed by the ambipolar diffusion resistivity are also shown (dashed lines).

Appendix C Implementing the combined multi-charge grain physics and gas chemistry

We adopteda simplified implementation of the chemical reaction network. In principle, reactions (e.g., charge exchange and recombination) have to be included explicitly for each grain charge (gr^Z_d = gr^Z_min,…, gr³⁻, gr²⁻, gr⁻¹, gr⁰, gr⁺, gr²⁺,…, gr^Z_max). However, with Z_max > 10, the computation cost becomes prohibitive. Instead, the abundance of grains with a charge Z_d is defined by [gr^Z_d] = n_df(Z_d), where Z_d is the charge of the grain (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-charge species alone.

The grain charge distribution f(Z_d) is $${\sum }_{Z_{\text{min}}}^{Z_{\text{max}}}f(Z_{\text{d}})={\sum }_{Z_{\text{min}}}^{-1}f(Z_{\text{d}})+f(0)+{\sum }_1^{Z_{\text{max}}}f(Z_{\text{d}})=1.$$(C.1)

The total abundance is $${\sum }_{Z_{\text{min}}}^{Z_{\text{max}}}[{\text{gr}}^{Z_{\text{d}}}]=n_{\text{d}}.$$(C.2)

We defined the pseudo-charge species (Z⁻, Z⁻, Z) such that $$[Z^{-}]+[Z]+[Z^{+}]=Z_{\text{max}}{n}_{\text{d}},$$(C.3)

where we have split the grain population into three categories: the negatively charged grains $$[Z^{-}]=-n_{\text{d}}{\sum }_{Z_{\text{min}}}^{-1}{Z}_{\text{d}}f(Z_{\text{d}});$$(C.4)

the positively charged grains $$[Z^{+}]=n_{\text{d}}{\sum }_1^{Z_{\text{max}}}{Z}_{\text{d}}f(Z_{\text{d}});$$(C.5)

and the neutral grains, which are modeled by the pseudo-charge species Z of abundance $$[Z]=n_{\text{d}}{Z}_{\text{max}}-[Z^{-}]-[Z^{+}].$$(C.6)