Open Access
Issue
A&A
Volume 632, December 2019
Article Number A44
Number of page(s) 28
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201732187
Published online 29 November 2019

© W. F. Thi et al. 2019

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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−2 M or even higher, and disappear after a few million years (Alexander et al. 2014). The disk material can evaporate, form giant planets, or fall onto the star. Gas can accrete from the inner rim of disks to the star at a rate of 10−9–10−7 M yr−1 (Calvet et al. 2000; Muzerolle et al. 2004; Johns-Krull et al. 2000; Hartmann et al. 1998), although episodic accretion rates of up to 10−5 M yr−1 can occur for young disks (Audard et al. 2014). Accretion shocks onto the stellar surfaces result in ultraviolet excess emission. Gas accretion is also traced by optical emission lines (Mendigutía et al. 2012, 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 ν = αcsh, where cs 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)cs in the TW Hya disk, while Flaherty et al. (2018) constrained α to lie between 0.04 cs and 0.13 cs. Flaherty et al. (2015) provided a low upper limit of vturb < 0.03cs for the turbulence in the HD 163296 disk. Earlier studies with the IRAM Plateau de Bure interferometer provided limits of vturb ≤ (0.3−0.5)cs (Dartois et al. 2003; Piétu et al. 2007). Hughes et al. (2011) constrained the value vturb ≤ 0.1cs for the TW Hya disk and vturb ≤ 0.4cs 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 H2, 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 Trp 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 = TgasTdust (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, (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−4. Collisional ionization of metals can also occur (2)

where M is a neutral metal (e.g., Fe, Mg, S, Na, Si) and N is H, H2, 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 fPAH = 10−1–10−3 compared to their interstellar abundance of 3 × 10−7 (Tielens 2008). We chose to use the peri-condensed circumcoronene (C54 H18) 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+, PAH2+, PAH3+; see Table 1). The effective radius of a PAH is computed as (Weingartner & Draine 2001b) (3)

where NC is the number of carbon atoms in the PAH. The radius for the circumcoronene is aPAH (C54H18) = 4.686 × 10−8 cm. The PAH ionization potential can either be taken from the literature when they are measured or estimated using (Weingartner & Draine 2001b) (4)

where W0 is the work function assumed to be 4.4 eV (7.05 × 10−12 erg) and ZPAH 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).

Table 1

Circumcoronene electron affinity and ionization potential.

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 = πa2. The charge of the grain is Zd with a discrete charge distribution function f(Zd) with a minimum charge Zmin and maximum charge Zmax, such that .

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 W0 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 Zd. The effective work function Weff, equivalent to the ionization potential for an atom or molecule, becomes (Weingartner & Draine 2001b) (5)

where Wc is the extra work function defined by (6)

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

where W0 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 (8)

where x0 = 2.5 and λ0 = 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 Zd reads (9)

When a grain is negatively charged the photodetachment threshold energy is (Weingartner & Draine 2001b) (10)

where the electron affinity is (11)

The band gap Ebg = 0 for metals and semimetals, while it can be several eVs for other materials. Weingartner & Draine (2001b) assumed (W0Ebg) = 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 Emin follows the definition of Weingartner & Draine (2001b) (12)

For a singly negatively charged grain, assuming Zd = −1 and a = 1 μm, one gets pd(Zd, a) ≃ 2.9 eV or λpd ≃ 0.4 μm for a pure silicate grains; and pd(Zd, 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 (13)

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

thumbnail 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.

Open with DEXTER

2.4.3 Electron attachment

The electron attachment rate coefficient onto neutral grains is (14)

where Te is the electron temperature assumed equal to the gas kinetic temperature, Se 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 me 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 (Zd > 0) the electron recombination rate is enhanced by Coulomb attraction (15)

where Wc is the extra work function (see Eq. (6)). On the contrary, the recombination rate is lower if the grain is negatively charged (Zd < 0): (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) (17)

where Td 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 Wtherm is equal to Weff for positively charged grains and to EA + Emin 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. (18)

where the work Wcd is equal to Weff for positively charged grains and to EA(Zd + 1) + Emin 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): (19) (20)

where and are gas-phase and surface ions, respectively, or (21)

The gas-phase species acquires the approach energy, which is the sum of the adsorption energy (ΔEads(Xg)) and the Coulomb energy ΔECoulomb. The parameter ΔErec 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 (ΔErelax < 0). It can also be used to break the weak bond (ΔEads(Xs) ~ 0.1–0.2 eV) between the species Xs 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: (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, (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 (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) (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 (26)

Both recombination and charge exchange reactions proceed with the rate (27)

where Etherm is an energy equal to the endothermicity of the reaction. For an exothermic charge exchange, the energy is null (Etherm = 0). The ion temperature Tion is equal to the gas thermal temperature Tgas and we assumed Sion = 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 (28)

for neutral and positively charged grains. We assume that the adsorption and Coulomb energy are negligible. For negatively charged grains (29)

where IPneu is the ionization energy of the neutral species. Again for Zd = −1 and a = 1 μm, EA (Zd + 1, a) ≃ 2.99 eV and Emin(Zd < 0, a) ≃ 2.9 eV, thus (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, (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 (32)

where the reduced mass is basically that of the gas-phase species μ =mn and the polarizability is αpol = 10−24, where the barrier term becomes (33)

in practically all cases, the geometrical rate dominates, (34)

Anions can also react with positively charged grains (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 H3 O+ and NH 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, W0 = 8 eV, a maximum photon energy of max = 13.6 eV, an elementary charge e in cgs of 4.803206815 × 10−10, and the conversion 1 eV = 1.60217657 × 10−12 erg, the maximum charge (assuming that photoejection is the dominant electron ejection process) reads (36)

For the minimum grain charge, we adopted again the formalism of Weingartner & Draine (2001b) (37)

where (38)

For a silicategrain 1 micron in radius, Uait ≃ 702 V and Zmin ≃−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 (kRD(Zd < 1) = ke,-), we can derive an upper limit to the amount of negative charges on a grain from Eqs. (14) and (17), (39)

where (40)

Assuming that Se = 1, r = 0.5, and Td = Tg = T, (41)

and (42)

For dense regions with ne = 1012 cm−3 and T = 100 K, the term B is greater than unity. The term kTlnB ≃ 2.37 × 10−5Tln(5.27 × 10−8T−3∕2ne) eV is always negligible compared to W0 (=8 eV), thus a strong lower limit to the grain charge can be derived: (43)

Assuming elementary charge in cgs of e = 4.80321 × 10−10 statC and W0 = 8 eV, Ebg = 5 eV, and f ~1, we obtain (44)

In this study, we did not include photoejection due to X-ray photons, and we adopted in our numerical models the limits (45)

Table 2

Examples of work function and ionization potentials (≤13.6 eV).

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 Se = Si (Evans 1994) (46)

where we can set x = Zde2akT. Using the medium global neutrality nelec = nion + Zdnd, the equation becomes (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 nion ~ nelec an approximate solution to the equation is (48)

where nnucl is the average number of nucleons of the cations. If the main ion is a molecular ion with 19 nucleons (HCO+), the approximation becomes (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 (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 (51)

This last equation together with Fig. 2 shows that, unless the ionization fraction is below 10−13, 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 Zd when photons are absent has a standard deviation of δZd = 0.5|Zd|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.

thumbnail 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.

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

Open with DEXTER

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 Ptherm(r, z) to magnetic pressures Pmag(r, z), by the function (52)

where δ is a parameter between 0.5 and 1. The βmag (beta magnetic) parameter is related to other gas parameters (53)

where Eth and Emag are the thermal and magnetic energy respectively, ρ is the gas mass density, is the sound speed with mean molecular mass μn = 2.2 amu (atomic mass unit), Bz(r) is the vertical component of the magnetic field, and vA is the Alfvèn speed in the disk vertical direction defined by (54)

The values for αideal are bound between 2 × 10−3 and 0.1 for βmag between 1 and 106. 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), (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) (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−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−4 when all the carbon is ionized in the upper atmospheres (weakly charged plasma). Hydrogen is predominantly molecular. The ionization fraction at disk surfaces can be much 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 mj, charge Zj, number density nj, and drift velocity relative to neutral vj. The current density J including all the charge-bearing species j (electrons, ions, PAHs, dust grains) with charge Zj is (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 (58)

where (59)

is the drag coefficient (cm3 g−1 s−1) with 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 ∑jnjZj = 0 we obtain (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) (61)

where and 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−1. The knowledge of the conductivities (or equivalently the resistivities) is central to the efficiency of MRI-driven turbulence in disks. The conductivities depend on the charge carriers and hence on the disk chemistry. Conversely, the decay of 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−3) at gas temperature Tgas (K) with a mean mass mn (in grams) and mass density ρn (in grams cm−3). We define for each charged species j (j = e, atomic and molecular ions, PAH, gr, where gr stands for dust grains) the unit-free plasma βj parameter (which should not be confused with βmag) as (62)

where Zj is the charge of species j of mass mj (in grams); B is the magnetic field strength, for which we only consider the vertical component Bz (in units of gauss, Gs, with 1 Gs = 1 cm−1∕2 g1∕2 s−1); and e is the elementary charge (in units of statC). Finally, is the drag frequencies of the positively or negatively charged species j of number density nj 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 (63)

The collisionrates , 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 (64)

The atomic and molecular ions have average collision rates with neutral species following the Langevin rate (65)

Specifically for the HCO+-H2 system, Flower (2000) provided the equation (66)

which shows a weak temperature dependency. At 100 K, the rate becomes 2.5 × 10−9. For PAH ions, we use the largest value between the temperature-independent Langevin and the geometrical rate (67)

where aPAH is the radius of the PAH. The average collision rate between the negatively charged grains with the neutral gas species is approximated by (68)

We assume that the grains are basically immobile with respect to the gas. For a fully molecular gas mn ≈ 2.2 amu = 3.65317 × 10−24 grams, (69)

The ohmic diffusivity is characterized by the dimensionless Elsasser number, which is defined as the ratio of the Lorentz to the Coriolis force, (70)

where (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 × 1013 cm. The ohmic resistivity is ηO = (c2∕4π)∕σO. The ohmic Elsasser number can be rewritten as the sum of the contribution of each charged species to the ohmic conductivity: (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: (73)

This criterion translates to (74)

In the approximation of vertical isothermal disk (h = cs∕Ω) and in the condition of vA = cs (βmag = 2, αideal = 1), the criterionbecomes (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 (~B2) compared to the energy dissipated by ohmic resistivity.

thumbnail Fig. 4

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

Open with DEXTER

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 (76)

where the ambipolar resistivity ηA is (Wardle & Ng 1999) (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 (78)

where the Pedersen conductivity is (79)

The Hall conductivity σH is defined by (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 (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: (82)

The maximum possible value for α is shown in Fig. 4. For Am = 1, max(α) ≃ 10−2. 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: (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 (84)

when , and αeff ≃ 0 otherwise. It has been shown that, when , 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 (∝ krecnionne). 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⟩ = 1010 cm−3 and 1012 cm−3). 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−3 down to a total ionization fraction of ~10−8–10−10 depending on the gas density. At disk surfaces, where the ionization fraction is high and the value of βmag low, the efficiency of MRI turbulence is probably limited by the ambipolar diffusion resistivity to a maximum value of ~0.25.

thumbnail 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.

Open with DEXTER

5 Disk model

We chose to model a typical protoplanetary disk with a total gas+solid mass of 10−2 M (Woitke et al. 2016). The disk extends from rin = 1 au to rout = 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 amin = 0.1 μm to amin = 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 fPAH is set to 0.01 (fPAH = 1 corresponds to an abundance of 3 × 10−7). The disk flares with an index of 1.15. The gas scale-height is 1 au at 10 au (10%). The cosmic ray flux is assumed to be constant throughout the disk at 1.7×10−17 s−1. 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 102 to 106. For comparison, we also ran a model with no (MRI) turbulence, i.e., the line broadening is purely thermal (vturb = 0 km s−1) and other models with a constant turbulence width vturb = 0.05, 0.1, 0.15, and 0.2 km s−1 and no turbulence heating throughout the disk (passive disk models). The central object is either a T Tauri star (M* = 0.7 M, Teff = 4000 K, L* = 1.0 L) or a Herbig Ae star (M* = 2.3 M, Teff = 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) (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 (86)

Adopting δ = 0.5, (87)

Assuming that the pressure term can be described by a barotropic law (88)

the gas heating rate is (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 . 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 cs by vturb = cs (Hughes et al. 2011). Another relationship such as vturb = αeff cs has been proposed (Simon et al. 2013). The local line width becomes , where vth 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 = 104. The figures show the different disk structures with radius r in au in the horizontal axis and zr 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−4) from atomic carbon ionization stick to the grains, rendering them negatively charged. The total charge reaches the value of the elemental abundance of carbon (1 to 3 × 10−4). As carbon becomes neutral at lower 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 Tdust well below the sublimation temperatures (Tsubl. ~ 1500 K) even though the gas temperature can reach up to ~10 000 K (Fig. F.1). We have not considered the destruction of dust grains in warm gas due to chemisputtering.

In the midplane inner disk region (r < 1 au with total charge fraction <10−12), the radiation field is too weak to eject the 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 fPAH 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 = 104 (left panels) and βmid = 106 (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 = 104 and βmid = 106. In unobscured regions, the electron fractional abundance χ(e) is relatively high in the range 10−8–10−4. Due to the limit on negative charges on grains (Eq. (50)), in a gas with total ionization fraction greater than a few 10−12, the free electrons will be the main contributor to the 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−13: (90)

The approximated ohmic Elsasser distribution for a disk model with βmid = 104 is shown in Fig. G.1.

The disk distribution for the effective turbulent parameter αeff is shown in Fig. 9 for βmid of 104 (left panels) and 106 (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 = 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 = 104, CO will be emitted significantly in the turbulent warm molecular disk layers with vturbcs between 0.25 and 0.4. For the βmid = 106 model, vturbcs is between 0.1 and 0.3. For completeness, the sound speed cs for the two models are plotted in Fig. E.1. The CO emission comes mostly from gas with sound speed ~0.3–0.5 km s−1. Figure 10 shows the main heating and cooling processes for the disk models with βmid = 104 (left panels) and βmid = 106 (right panels). In the dead zone where , 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 (zr ~ 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 H2 formation on dust grains and by the chemical reactions. The disk atmosphere region where PAHs are not frozen onto dust grains in the βmid = 104 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 > 106). 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 H2 O lines can trace the mechanically heated gas.

The actual CO line profiles computed in non-LTE for CO, 13CO, and C18O 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 = 106 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 vturb, we can constrain vturb down to 0.05 km s−1 (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 13CO transitions (J = 3–2 and J = 2–1) show the strongest sensitivity to βmid or vturb. A low peak-to-trough ratio indicates a high level of turbulence. The 12CO lines are emitted in the region where vturb/cs reaches its maximum. The optically thinner 13CO lines are emitted in the layers closer to the midplane where vturb/cs is lower. The peak flux is mostly emitted in the outer disk, where the C18O lines can be optically thin. Therefore, the peak fluxes for the C18O lines are low and the peak-to-trough ratios for the C18O lines are less sensitive to the turbulence broadening than the optically thick 12CO and 13CO lines.

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

thumbnail 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.

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

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

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

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

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

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

Open with DEXTER

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 cs for a βmid = 104 in the midplane. The computed turbulent line widths assuming βmid = 104 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 (vturb < 0.03 cs). 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 vturb = 0.2 cs. 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 12CO J = 3–2 line. In our MRI and fixed vturb 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 12CO J = 3–2 and J = 2–1 profiles show higher peak-to-trough ratios than the 13CO profiles (see right panels of Fig. 13) in contrast to our model results. A higher 12CO/13CO 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 N2 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−19 s−1. Their dead zone is defined by the criterion REM ≤ 3300, where REM is the magnetic Reynolds number (Balbus & Terquem 2001). The molecular ions HCO+ and N2 H+ have theoretical abundances that depend sensitively on the treatment of ionization and recombination in the gas (Rab et al. 2017) and on the N2 self-shielding for N2 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 N2 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, N2 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 cs and 0.5 cs. ALMA observations of CS in the DM Tau disk reveals a non-thermal line width between 0.05 and 0.15 km s−1 at 300 au (Semenov et al. 2018). Sulfur chemistry is complex and CS lines are weak, even in massive disks Dutrey et al. (2011).

The magnetic Reynolds number and the ohmic Elsasser number are then related by (91)

assuming cs = hΩ. Flaherty et al. (2017) used DCO+ lines in addition to CO and C18O 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 < 105, i.e., that CO line profiles can be used to constrain α down to ~5 × 10−3. The CO abundances can also vary over the lifetime of the disks (Yu et al. 2017), which may affect the CO line profile.

The CO rotational lines that are seen with large beams in the submillimeter domain probe a large disk volume, but may not be sensitive to the location in the protoplanetary disk s with the highest vturbcs 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 vturbcs 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 Tdust 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 H2 into H by an energetic ionization event (X-ray, cosmic rays, radioactive decay), and subsequently H. 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−4 for αsettle (Pinte et al. 2016). This corresponds to a disk model with βmid ~ 106.

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.

thumbnail 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).

Open with DEXTER

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 H2 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 vturb), 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 (A.1)

where χ(e) is the free electron fractional abundance. The equation can be rewritten as (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: (A.3)

For the electrons, (A.4)

(memn) and for the dust grains (A.5)

where the average mass of a grain is (A.6)

Since md > > mn we obtain (A.7)

The average number of dust grains in the disk is (A.8) (A.9)

where gd is the gas-to-dust mass ratio. We have assumed a silicate mass density of 3.0 g cm−3. Replacing the dust grain number density by the equation above, we obtain (A.10)

For a gas-to-dust mass ratio of 100, (A.11)

where χ(e) is the relative abundance of free electrons. If we assume that grains have charge 0 > Zd ≥ − Zmax and using the maximum possible negatively charged grains (Eq. (44)), the ratio becomes (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−14. 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−11. If we assumes that the dust charge in the UV-obscured region can be approximated by Eq. (49), the ohmic conductivity ratio is (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 (A.14)

Assuming singly charged ions and that most of the negative charges are carried by the electrons, (A.15)

where (A.16)

assuming that mi > > mn. 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 (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 1010 cm−3 and temperature of 100 K. The four panels correspond to different values of βmag (104, 103, 102, 10).

thumbnail 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).

Open with DEXTER

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 (grZd = grZmin,…, gr3−, gr2−, gr−1, gr0, gr+, gr2+,…, grZmax). However, with Zmax > 10, the computation cost becomes prohibitive. Instead, the abundance of grains with a charge Zd is defined by [grZd] = ndf(Zd), where Zd is the charge of the grain (integer) and f(Zd) is the normalized discrete distribution (fractional abundance) of grains with charge Zd. The values for Zmin and Zmax are chosen such that f(Zmin) = 0 and f(Zmax) = 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(Zd) is (C.1)

The total abundance is (C.2)

We defined the pseudo-charge species (Z, Z, Z) such that (C.3)

where we have split the grain population into three categories: the negatively charged grains (C.4)

the positively charged grains (C.5)

and the neutral grains, which are modeled by the pseudo-charge species Z of abundance (C.6)

The averagegrain charge is (C.7)

The pseudo-neutral grain species Z includes the actual neutral grains and the positively charged grains exchanging charges with cations of index j to become more positively charged and absorbing a photon and ejecting an electron: (C.8)

Considering only this reaction, the chemical reaction differential equation is (C.9)

Another way to write the reaction differential equation is (C.10)

where the positive grain-charge-averaged charge-exchange rate is (C.11)

or (C.12)

The absorption of a neutral grain or of a positively charged grain leads to the ejection an electron: (C.13)

The rate for the pseudo-charge species Z is , whose derivation follows (C.14)

We obtain (C.15) (C.16)

Likewise, the pseudo-neutral grains model the thermal emission of electrons actual neutral and positively charged grains: (C.17)

The rate of thermionic emission is (C.18)

The pseudo-neutral species Z is also used to model electron attachment reactions (C.19)

with the rate (C.20)

where (C.21)

The electron recombinations of positively charged grains are modeled by the pseudo-reaction (C.22)

In the case of the electron recombination with a grain with n positive charges and rate kre(n), the sum of all the recombination chemical rate reactions reads (C.23)

We can derive the rate : (C.24)

The pseudo-negative charge species Z models the reactions leading to neutralize the ions A of index j: (C.25)

The generic reaction is written as (C.26)

For example, the pseudo-reaction (C.27)

encompasses the dissociation recombination reactions of HCO+ with all negatively charged grains. (C.28)

where (C.29)

The pseudo-charge species Z also participate to the photodetachment reaction (C.30)

with the rate (C.31)

Finally Z can emit electrons due to the thermionic effect (C.32) (C.33)

To illustratethat our concept is valid, we assume first that the grains can have at most one positive or negative charge, in other words, Zmin = −1 and Zmax = 1. Then the reactions in Table C.1 correspond to the actual reactions presented in Table C.2.

Table C.1

Pseudo-reactions involving the pseudo-grain charge species Z, Z, and Z+.

Table C.2

Reactions involving the grain charge species gr, gr, and gr+.

The definition of the rates above depends on the charge distribution f(Zd), which has to be known or solved simultaneously with the chemical reactions. In a steady-state treatment, the chemistry and the grain chargedistribution are solved alternatively.

Using an estimate of the electron and ion abundances, the grain charge distribution f(Zd) is calculated by iteration. The average grain charge ⟨Zd⟩ is determined by balancing the photoemission (photoejection and photodetachment) rates and the thermionic rate with the recombination and charge exchange rates. Then all the chemical rates between a gas-phase species and the pseudo-charge species Z, Z, and Z+ can be derived. The electron abundance is computed from the global gas neutrality considering the charges on gas-phase species, PAHs, and dust grains. The grain charge distribution f(Zd) is subsequently re-evaluated knowing the abundances of the cations and electrons. This iterative process ensures thatthe grain charges and the gas-charge chemistry are computed self-consistently. In particular, the global neutrality of the medium is enforced. We assume that a dust grain can have charges from Zmin to Zmax, which are numerical free parameters. The results do not depend on the exact values for Zmin or Zmax provided that ZminZd − 10 and ZmaxZd + 10. We assumed in this work for the numerical reason that Zmin = −Zmax. If the chemistry is solved time-dependently, the grain charge fractional distribution has to be solved together with the chemical reaction differential equations. Solving a time-dependent grain charge distribution entails many challenges whose resolution goes beyond the scope of this paper.

Appendix D Disk model parameters

The protoplanetary disk parameters are summarized in Table D.1.

Table D.1

Model parameters and values for the reference model.

Appendix E Disk vertical component of the magnetic field

The vertical component of magnetic field, the distribution of βmag, and the sound speed distribution are shown for the disk model with βmid = 104 (left panels) and βmid = 106 (right panels) in Fig. E.1.

thumbnail Fig. E.1

Resulting vertical component of the magnetic field (upper panels) and of the distribution of βmag (middle panels) for disk model with βmid = 104 (left panel) and βmid = 106 (right panel). We note that βmag = βmid at all radii r at z = 0. The contour with βmag = 1 defines the disk surface. The sound speed structures are shown in the lower panels. The red contours in the lower panels show vturb∕cs.

Open with DEXTER

Appendix F Disk model temperature structure

The disk dust and gas thermal structure are shown in Fig. F.1.

thumbnail Fig. F.1

Temperature structure of disk dust (left panel) and gas (right panel). The top panels correspond to the βmid = 104 model and the lower panels to the βmid = 106 model.

Open with DEXTER

Appendix G Analytical ohmic Elsasser number

Figure G.1 shows the distribution of ΛOhm using an analytical approximation for βmid = 104.

thumbnail Fig. G.1

Analytical approximation to the ohmic Elsasser number for a model with βmid = 104.

Open with DEXTER

References

All Tables

Table 1

Circumcoronene electron affinity and ionization potential.

Table 2

Examples of work function and ionization potentials (≤13.6 eV).

Table C.1

Pseudo-reactions involving the pseudo-grain charge species Z, Z, and Z+.

Table C.2

Reactions involving the grain charge species gr, gr, and gr+.

Table D.1

Model parameters and values for the reference model.

All Figures

thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 4

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

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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 , 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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).

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

Resulting vertical component of the magnetic field (upper panels) and of the distribution of βmag (middle panels) for disk model with βmid = 104 (left panel) and βmid = 106 (right panel). We note that βmag = βmid at all radii r at z = 0. The contour with βmag = 1 defines the disk surface. The sound speed structures are shown in the lower panels. The red contours in the lower panels show vturb∕cs.

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

Temperature structure of disk dust (left panel) and gas (right panel). The top panels correspond to the βmid = 104 model and the lower panels to the βmid = 106 model.

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

Analytical approximation to the ohmic Elsasser number for a model with βmid = 104.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.