Chemical network reduction in protoplanetary disks

TL;DR

This paper demonstrates a species-based network reduction method that simplifies complex chemical networks in protoplanetary disks, enabling efficient modeling of major species' abundances in key disk regions.

Contribution

The study introduces and evaluates a species-based reduction technique that accurately reproduces major chemical abundances with significantly smaller networks in protoplanetary disks.

Findings

Reduced networks with 20-30 species effectively reproduce major species abundances.

The method works well in disk midplane and surface regions.

Most major carriers of elements have similar abundances in reduced and full networks.

Abstract

Protoplanetary disks (PPDs) are characterized by different kinds of gas dynamics and chemistry, which are coupled via ionization, heating and cooling processes, as well as advective and turbulent transport. However, directly coupling gas dynamics with time-dependent chemistry is prohibitively computationally expensive when using comprehensive chemical reaction networks. In this paper, we evaluate the utility of a species-based network reduction method in different disk environments to produce small chemical networks that reproduce the abundances of major species found in large gas-phase chemistry networks. We find that the method works very well in disk midplane and surfaces regions, where approximately 20-30 gas phase species, connected by $\sim$50-60 gas phase reactions, are sufficient to reproduce the targeted ionization fraction and chemical abundances. Most species of the reduced networks, including major carriers of oxygen, carbon and nitrogen, also have similar abundances in the reduced and complete network models. Our results may serve as an initial effort for future hydrodynamic/magnetohydrodynamic simulations of PPDs incorporating time-dependent chemistry in appropriate regions. Accurately modeling the abundances of major species at intermediate disk heights, however, will require much more extended network incorporating gas-grain chemistry and are left for future studies.

Figures

Tables

Table 1. Table 1 “molecular cloud” like initial abundance

H2	$5.00\times {10}^{-1}$	He	$1.40\times {10}^{-1}$
HCN	$2.00\times {10}^{-8}$	CS	$4.00\times {10}^{-9}$
CO	$4.33\times {10}^{-5}$	HCO+	$9.00\times {10}^{-9}$
SO	$5.00\times {10}^{-9}$	C+	$1.00\times {10}^{-9}$
N2	$1.00\times {10}^{-6}$	NH3	$8.00\times {10}^{-8}$
H${}_{}{}^{+}{}_3{}^{}$	$1.00\times {10}^{-8}$	C2H	$8.00\times {10}^{-9}$
H2O(gr)	$2.40\times {10}^{-4}$	CH4(gr)	$5.70\times {10}^{-6}$
CH3OH(gr)	$1.80\times {10}^{-5}$	CO2(gr)	$3.30\times {10}^{-5}$
N	$2.25\times {10}^{-5}$	CN	$6.00\times {10}^{-8}$
Si(gr)	$1.00\times {10}^{-9}$	Mg(gr)	$1.00\times {10}^{-9}$
Fe(gr)	$1.00\times {10}^{-9}$

Table 2. Table 2 Reduced network in disk midplane

Neutrals	H H2 H2O He N2 CO CH4 CO2 O2 NH3
Ions	e- H H${}_{}{}^{+}{}_2{}^{}$ H2O+ He+ N${}_{}{}^{+}{}_2{}^{}$ CO+ CH${}_{}{}^{+}{}_2{}^{}$ CH${}_{}{}^{+}{}_3{}^{}$ CH${}_{}{}^{+}{}_4{}^{}$
	CH${}_{}{}^{+}{}_5{}^{}$ H${}_{}{}^{+}{}_3{}^{}$ H3O+ NH${}_{}{}^{+}{}_4{}^{}$ C+ OH+ HCO+ N2H+

Table 3. Table 3 Reduced network in the disk surface

Neutrals	H C O OH H2O CO Si S H2
Ions	e- H+ H${}_{}{}^{+}{}_2{}^{}$ C+ O+ OH+ CO+ H2O+ Si+ S+
	CH+ H3O+ SiO+ HCO+

Table 4. Table 4 Gas phase reaction list at the disk midplane

Charge exchange reactions	H+ + CH4 $\to$ CH${}_{}{}^{+}{}_4{}^{}$ + H	H+ + H2O $\to$ H2O+ + H	He+ + N2 $\to$ N${}_{}{}^{+}{}_2{}^{}$ + He
	H${}_{}{}^{+}{}_2{}^{}$ + CH4 $\to$ CH${}_{}{}^{+}{}_4{}^{}$ + H2	H${}_{}{}^{+}{}_2{}^{}$ + CO $\to$ CO+ + H2	H2 + He+ $\to$ He + H${}_{}{}^{+}{}_2{}^{}$
	He+ + CH4 $\to$ CH${}_{}{}^{+}{}_4{}^{}$ + He	He+ + H2O $\to$ H2O+ + He
Dissociative recombinations	CH${}_{}{}^{+}{}_5{}^{}$ + e- $\to$ CH4 + H	H${}_{}{}^{+}{}_3{}^{}$ + e- $\to$ H2 + H	H${}_{}{}^{+}{}_3{}^{}$ + e- $\to$ H + H + H
	H3O+ + e- $\to$ H2O + H	HCO+ + e- $\to$ CO + H	N2H+ + e- $\to$ N2 + H
	NH${}_{}{}^{+}{}_4{}^{}$ + e- $\to$ NH3 + H
CRs and X-ray induced reactions	H2 $\to$ H+ + H + e- (0.021)	H2 $\to$ H${}_{}{}^{+}{}_2{}^{}$ + e- (0.883)	H2 $\to$ H + H (0.096)
	He $\to$ He+ + e- (0.478)	H $\to$ H+ + e- (0.440)
Neutral-neutral reactions	H + H $\to$ H2 (grain surface)
Ion-neutral reactions	C+ + H2O $\to$ HCO+ + H	CH4 + N2H+ $\to$ N2 + CH${}_{}{}^{+}{}_5{}^{}$	CH${}_{}{}^{+}{}_5{}^{}$ + CO $\to$ HCO+ + CH4
	CH${}_{}{}^{+}{}_5{}^{}$ + H2O $\to$ H3O+ + CH4	CO + N2H+ $\to$ HCO+ + N2	H+ + CH4 $\to$ CH${}_{}{}^{+}{}_3{}^{}$ + H2
	H${}_{}{}^{+}{}_2{}^{}$ + H2 $\to$ H${}_{}{}^{+}{}_3{}^{}$ + H	H2 + CH${}_{}{}^{+}{}_2{}^{}$ $\to$ CH${}_{}{}^{+}{}_3{}^{}$ + H	H2 + CH${}_{}{}^{+}{}_4{}^{}$ $\to$ CH${}_{}{}^{+}{}_5{}^{}$ + H
	H2 + CO+ $\to$ HCO+ + H	H2 + H2O+ $\to$ H3O+ + H	H2 + He+ $\to$ He + H+ + H
	H2 + N${}_{}{}^{+}{}_2{}^{}$ $\to$ N2H+ + H	H2 + OH+ $\to$ H2O+ + H	H2O + HCO+ $\to$ CO + H3O+
	H2O + N2H+ $\to$ N2 + H3O+	H${}_{}{}^{+}{}_3{}^{}$ + CH4 $\to$ CH${}_{}{}^{+}{}_5{}^{}$ + H2	H${}_{}{}^{+}{}_3{}^{}$ + CO $\to$ HCO+ + H2
	H${}_{}{}^{+}{}_3{}^{}$ + H2O $\to$ H3O+ + H2	H${}_{}{}^{+}{}_3{}^{}$ + N2 $\to$ N2H+ + H2	H${}_{}{}^{+}{}_3{}^{}$ + NH3 $\to$ NH${}_{}{}^{+}{}_4{}^{}$ + H2
	H + CH${}_{}{}^{+}{}_4{}^{}$ $\to$ CH${}_{}{}^{+}{}_3{}^{}$ + H2	H + CH${}_{}{}^{+}{}_5{}^{}$ $\to$ CH${}_{}{}^{+}{}_4{}^{}$ + H2	He+ + CH4 $\to$ CH${}_{}{}^{+}{}_2{}^{}$ + He + H2
	He+ + CH4 $\to$ CH${}_{}{}^{+}{}_3{}^{}$ + He + H	He+ + CO2 $\to$ O2 + C+ + He	He+ H2O $\to$ OH+ + He + H
	NH3 + H3O+ $\to$ NH${}_{}{}^{+}{}_4{}^{}$ + H2O	C+ + CO2 $\to$ CO+ + CO	NH3 + HCO+ $\to$ CO + NH${}_{}{}^{+}{}_4{}^{}$
Radiative associations	H2 + CH${}_{}{}^{+}{}_3{}^{}$ $\to$ CH${}_{}{}^{+}{}_5{}^{}$	H2 + C+ $\to$ CH${}_{}{}^{+}{}_2{}^{}$

Table 5. Table 5 Gas phase reaction list at the disk surface

Charge exchange reactions	S + C+ $\to$ C + S+	C+ + Si $\to$ Si+ + C	H+ + H2O $\to$ H2O+ + H
	H+ + O $\to$ O+ + H	H+ + OH $\to$ OH+ + H	H+ + S $\to$ S+ + H
	H + CO+ $\to$ CO + H+	H + H${}_{}{}^{+}{}_2{}^{}$ $\to$ H2 + H+	H + O+ $\to$ O + H+
Dissociative recombinations	CH+ + e- $\to$ C + H	CO+ + e- $\to$ O + C	H${}_{}{}^{+}{}_2{}^{}$ + e- $\to$ H + H
	H2O+ + e- $\to$ O + H2	H2O+ + e- $\to$ O + H + H	H2O+ + e- $\to$ OH + H
	H3O+ + e- $\to$ H2O + H	H3O+ + e- $\to$ O + H2 + H	H3O+ + e- $\to$ OH + H2
	H3O+ + e- $\to$ OH + H + H	HCO+ + e- $\to$ CO + H	OH+ + e- + O + H
	SiO+ + e- $\to$ Si + O
CRs and X-ray induced reactions	H2 $\to$ H+ + H + e- (0.021)	H2 $\to$ H${}_{}{}^{+}{}_2{}^{}$ + e- (0.883)	H2 $\to$ H + H (0.096)
	H $\to$ H+ + e- (0.440)
Neutral-neutral reactions	H + H $\to$ H2 (grain surface)	C + OH $\to$ CO + H	H2 + O $\to$ OH + H
	H2 + OH $\to$ H2O + H	H + H2O $\to$ OH + H2	H + OH $\to$ O + H2
Ion-neutral reactions	C+ + H2O + HCO+ + H	C+ + OH $\to$ CO+ + H	H${}_{}{}^{+}{}_2{}^{}$ + O $\to$ OH+ + H
	H2 + C+ $\to$ CH+ + H	H2 + CO+ $\to$ HCO+ + H	H2 + H2O+ $\to$ H3O+ + H
	H2 + O+ $\to$ OH+ + H	H2 + OH+ $\to$ H2O+ + H	H + CH+ $\to$ C+ + H2
	OH + Si+ $\to$ SiO+ + H
Radiative associations	H+ + H $\to$ H${}_{}{}^{+}{}_2{}^{}$	H + C+ $\to$ CH+	H + O $\to$ OH
	H + OH $\to$ H2O
Radiative recombinations	C+ + e- $\to$ C	H+ + e- $\to$ H	S + e- $\to$ S
	Si+ + e- $\to$ Si
Photo reactions	CO $\to$ C + O	H2O $\to$ H2O+ + e-	H2O $\to$ OH + H
	OH $\to$ O + H	OH $\to$ OH+ + e-	S $\to$ S+ + e-
	Si $\to$ Si+ + e-	C $\to$ C+ + e-

Equations

$$R_H=\frac{1}{2}{s}_H{n}_H{v}_H{n}_d{\sigma }_d\eta ,$$

$$F_{\mathrm{pd}}=YG(r,\Sigma )f_s$$

$$\rho (z)={\rho }_0\exp \left(-\frac{z^2}{2H^2}\right)$$

$$\begin{array}{rl}{\zeta }_X& =\frac{L_X}{10^{29}erg/s}{r}_{\mathrm{au}}^{-2.2}[{\zeta }_1\left(e^{-({\mathrm{N}}_{H1}/{\mathrm{N}}_1{)}^{\alpha }}+e^{-({\mathrm{N}}_{H2}/{\mathrm{N}}_1{)}^{\alpha }}\right)\\ & +{\zeta }_2\left(e^{-({\mathrm{N}}_{H1}/{\mathrm{N}}_2{)}^{\beta }}+e^{-({\mathrm{N}}_{H2}/{\mathrm{N}}_2{)}^{\beta }}\right)]\end{array}$$

$$T_g(r,\Sigma )=T_{\mathrm{mid}}(r)\log \left[10+{\left(\frac{0.1\mathrm{g}\cdot \mathrm{c}{\mathrm{m}}^{-2}}{\Sigma }\right)}^3\right].$$

$$G(r,\Sigma )=\frac{10^7}{r_{\mathrm{au}}^2}\exp \left[-{\left(\frac{\Sigma }{{\Sigma }_{\mathrm{FUV}}}\right)}^2\right],$$

$$B[i]=\sum _{j=1}^{N^{\prime }}{\left(\frac{n_i}{g_j}\frac{\partial f_j}{\partial n_i}\right)}^2$$

Full text

1. Introduction

The formation of planets, and perhaps the existence of extraterrestrial life, is intrinsically linked to the chemistry of planet formation. Chemical processes in protoplanetary disks (PPDs) offer rich observational diagnostics that help constrain disk properties (e.g. Williams & Cieza, 2011; Henning & Semenov, 2013). The volatile content in PPDs is directly relevant to the abundance and composition of solids (Öberg et al., 2011), which affects grain surface properties and dust coagulation (Wada et al., 2009; Okuzumi et al., 2016), and will eventually be reflected in the composition of planets (Madhusudhan et al., 2011). Formation of complex organic molecules in disks may serve as initial seeds of extraterrestrial life (Öberg et al., 2015).

Chemical processes in PPDs are coupled with many aspects of the PPD gas dynamics. The level of disk ionization (including charge distribution), resulting from ionization-recombination chemistry, largely determines the coupling between gas and magnetic field (e.g. Umebayashi & Nakano, 1988; Sano et al., 2000; Fromang et al., 2002; Ilgner & Nelson, 2006; Okuzumi, 2009; Aikawa et al., 2015). Such coupling further governs the overall gas dynamics and angular momentum transport, and eventually determines the global disk structure and evolution (e.g., Wardle, 2007; Turner et al., 2014; Xu & Bai, 2016; Bai, 2017). In the disk atmosphere, heating and cooling processes are also directly coupled with chemistry (Glassgold et al., 2004; Nomura & Millar, 2005; Kamp et al., 2010), which determines the thermal structure in the disk surface layer. The heating and cooling processes also strongly affect disk mass loss through (magneto-)photoevaporation (Gorti et al., 2009; Owen et al., 2010; Bai et al., 2016; Bai, 2016), and contribute to disk dispersal. Moreover, detailed knowledge of thermodynamics and chemistry is essential in interpreting observations of emission lines from the warm disk surface layer (Carr et al., 2004; Simon et al., 2016). More recently, we showed in Xu et al. (2017) the impact that gas dynamics can have on CO depletion in the upper disk atmospheres.

Despite the evidence above that chemistry and gas dynamics in PPDs affect each other, they are generally modeled independently. Most studies of disk gas dynamics assume a steady-state chemistry. For instance, tabulated magnetic diffusivities based on steady-state ionization chemistry are commonly used in non-ideal magnetohydrodynamic (MHD) simulations of PPDs (e.g., Bai & Stone, 2013; Gressel et al., 2015; Bai, 2017). Disk chemistry studies have revealed that in most disk regions steady-state is never reached, however, but rather the chemistry continues to evolve time dependently throughout the disk life time.

It is highly desirable to couple realistic gas dynamics and chemistry in time-dependent simulations of PPDs. In recent years, the growing computational power has encouraged more efforts towards this direction. Examples include star-formation in molecular clouds (Glover et al., 2010), and on the simulation of disk ionization (e.g. Nomura et al., 2009; Semenov & Wiebe, 2011; Heinzeller et al., 2011; Furuya et al., 2013; Furuya & Aikawa, 2014)111There are other works focusing on the chemistry from core collapse to early stages of disks that follow the chemistry along either analytical trajectories or post-processing of simulation trajectories without feedback (Visser et al., 2009, 2011; Yoneda et al., 2016). and photoevaporation (e.g. Wang & Goodman, 2017b; Nakatani et al., 2017). However, the stiffness of chemical reaction networks requires implicit solvers that constantly invert matrices, where the cost scales as the number of species cubed. Therefore, evolving (magneto-)hydrodynamics with time-dependent chemistry is prohibitively expensive if a large chemical network is adopted. Existing works of this kind all rely on substantially reduced chemical network to minimize computational cost, and it is the goal of this paper to test the validity of reducing network size in the context of PPDs.

Small chemical networks can be constructed by either selecting important species and reactions from large, comprehensive chemical networks (Ruffle et al., 2002; Wiebe et al., 2003), or by building it up from scratch using chemical intuition (Keto et al., 2014). In the case of PPDs, we favor the first reductive approach since the chemistry is not always intuitive and the latter approach therefore risk omitting key reactions and reactants. A species-based approach was originally described by Tomlin et al. (1992) and Pilling (1997) in the context of combustion chemistry, and was later introduced to astrochemistry by Ruffle et al. (2002) and Rae et al. (2002). An alternative, reaction-based reduction method was presented by Wiebe et al. (2003) based on selecting important reactions rather than species. Using both methods, Semenov et al. (2004) (hereafter SWH04) studied chemical network reduction for PPDs, aiming at accurately reproducing the ionization fractions. They found good agreement between ionization fractions in the reduced and full networks, but had to develop a new reduced network for each disk region (defined by a disk radius and height), which somewhat limits the utility of the presented networks.

In this paper, we adopt the species-based method, which is simpler and gives more definitive criterion in selecting species and reactions. Compared to SWH04, we expand the scope of the inquiry to also include major coolants and some abundant molecules of interest to astrochemists and planet formation theorists. We also investigate the cost of adding key photoprocesses to the disk atmosphere reduced network (e.g. self-shielding) and the cost (in terms of network size) of combining reduced networks at different disk radii to enable broader applications. Our network reduction focuses on the the disk midplane (within $\sim 2$ scale heights) and disk surface layers (where FUV can penetrate). The chemistry in these two regions are relatively simple, yet has wide applications in the disk dynamics community.

In the midplane region, the ionization-recombination chemistry can strongly affect the gas dynamics via non-ideal MHD effects (e.g., Turner et al., 2014). While recombination timescales in disks are typically shorter than dynamical timescales, justifying steady-state approximations (Bai, 2011), they become comparable in the absence of small grains. For this latter case, Turner et al. (2007) showed that incorporating time-dependent chemistry in a grain-free network can help revive the Ohmic dead zone in the inner PPDs in the conventional picture of layered accretion, and the free-electron abundance depends on the chemistry. Coupling gas dynamics with chemistry in this region is also important when studying advective and turbulent transport of chemical species, especially across snowlines, where dynamical timescales are often comparable to chemical timescales in the disk midplane (Stevenson & Lunine, 1988; Kretke & Lin, 2007; Krijt & Ciesla, 2016).

In the disk surface, the mass loss process via (magneto-)photoevaporative wind is sensitive to the heating/cooling processes, which are intrinsically coupled to chemistry. In this region, chemical species can experience rapid varying environments in terms of density, temperature and radiation field due to advection in the outflow, invalidating the assumption of chemical equilibrium (Wang & Goodman, 2017b). Coupling chemistry with dynamics in the disk surface layer is also crucial to provide reliable observational diagnostics to constrain disk and wind properties.

Chemical networks containing a few tens (rather than hundreds) of species are typically manageable in simulations where the computational cost spent on chemistry calculations is comparable or smaller than that spent on hydrodynamics/MHD (e.g., Wang & Goodman (2017a)), and this is the goal of our network reduction in the disk midplane and surface layers for the above contexts. On the other hand, chemistry at the intermediate layer is known to be much more complex (SWH04), and may not be amendable to substantial network reduction.

The paper is organized as follows. In §2, we describe our original chemical network, disk model and the network reduction method. In §3, we present the results of our reduction in disk midplane and disk surface, and compare it with previous works. We summarize and discuss the applications in §4.

2. The Methodology of Network Reduction

In this section, we first describe the chemistry model and disk model that are used for our calculations in §2.1 and §2.2. Because the physical conditions in disks vary substantially in the vertical direction, we treat the models separately at the disk midplane and surface layers (neglecting the intermediate layer). In §2.3, we describe the reduction method. For each layer, we conduct network reduction at different radii, and then combine the results.

2.1. Chemistry Model

The full chemical reaction network is a standard astrochemical nework from UMIST database which includes 9 elements (H, He, C, N, O, S, Mg, Si, Fe) and 450 gas-phase species (no isotopologues and species with more than ten carbon atoms) (McElroy et al., 2013). Based on these species, 3940 gas-phase reactions are extracted from the latest version of the UMIST database. Photo-reactions only included in the disk surface models since the disk midplane region is almost completely shielded from UV photons (see last paragraph of this subsection for how photo-reactions are treated in the disk surface models).

A single population of dust grains is included with maximum grain charge of $\pm$ 4. Grains are assumed to be spherical with fixed size a=0.1$\mu$m and density $\rho$=3g cm*-3*. Dust to gas mass ratio is fixed to $f_{g}=10^{-4}$ to obtain the total surface area found in more realistic grain coagulation/fragmentation simulations in PPDs (Birnstiel et al., 2011)222Note that from the point of view of the chemistry, what matters the most is the surface area per volume element.. The grain-related reactions included in this otherwise gas-phase chemistry model are adsorption and thermal desorption, recombination with electrons and ions, as well as neutralization reactions among charged grain species. They are given in Table 3 and 4 of Ilgner & Nelson (2006), with binding energies taken from the UMIST database. In this paper, we refer to “mantle species” as counterparts of gas-phase species adsorbed onto grain surfaces. Formation of H2 on dust grains is included, where we assume a constant probability $\eta$ for a pair of adsorbed hydrogen atoms to form H2 molecules on grain surfaces. The effective reaction rate $R_{H}$ for H + H $\to$ H2 is given by (Cazaux & Tielens, 2002, 2004)

[TABLE]

where $s_{H}$ is the sticking coefficient of atomic H on grains (taken to be 1), $v_{H}$ is the thermal velocity of atomic H, $\sigma_{d}$ is the collisional cross section between H atoms and dust grains, which is approximately geometric, $n_{H}$ and $n_{d}$ are the number densities of atomic hydrogen and dust grains. We adopt $\eta=0.2$ as standard value (Cazaux & Tielens, 2004).

In the disk surface models, we further include far-UV (FUV) radiation, which drives a series of photo-ionization and photo-dissociateion reactions. Accordingly, we extract additional 136 photo-reactions from the UMIST database for disk surface chemistry. We also include photodesorption of CO and H2O. Using Eq. (6) of Hollenbach et al. (2009), the flux $F_{pd}$ of adsorbed particles leaving the surface due to photodesorption is given by

[TABLE]

where $Y=0.001$ for H2O and 0.01 for CO is the yield in molecules per photon. $G(r,\Sigma)$, which is defined in Section 2.2, is the photon flux incident on the unit surface area of grain surface and is parameterized in units of the Habing flux $G_{0}$ (Habing, 1968). $f_{s}$ is the surface fraction covered by the molecule. It is however fixed to unity, for simplicity.

We also implement a simplified model to test self-shielding effect for H2 where shielding factor is only a function of surface number density and estimated as $f=\rm min[1,(n/{10^{14}cm^{-2}})^{-3/4}]$ (Wolcott-Green et al., 2011). Shielding for CO and CI cannot be approximated so easily and is not included in this proof-of-concept study. Photoreactions rates for CO and CI are directly adopted from UMIST database. More discussions on the effect of shielding can be found in section 3. Apart from the grain process listed above, no grain surface chemistry is considered in our models.

Some caution should be exercised for ignoring grain-surface chemistry. In the disk surface layers, strong radiation fields are expected to maintain all molecules in the gas-phase and grain surface chemistry is thought to only play a minor role, though explicit tests are needed. In the midplane region, most volatiles except molecular hydrogen exist in the form of icy grain mantles, and hence grain surface chemistry could be important. In particular, grain-surface reactions offer important pathways for forming complex organic molecules (e.g., Walsh et al. 2014). The general outcome depends on a combination of initial chemical conditions, as well as ionization and radiation characteristics. On the other hand, the disk midplane models we consider start with all O, C and N already in molecular form, making some species less subjective to grain surface chemistry, including water, which mainly form in low temperature (i.e. molecular clouds; Cleeves et al., 2014).

2.1.1 Initial Abundances

Two types of initial abundances are considered in this paper. The first is “molecular cloud” like initial abundances, with abundances shown in Table 1. Compared to Table 2 of Cleeves et al. (2016), the only change we made is that we assume metal species are initially adsorbed on grain surface. The second type is the atomic initial abundances, where we convert all species in the “molecular cloud” initial abundances into their atomic phase.

For chemistry in the disk midplane, we use the “molecular cloud” initial abundances since comet observations indicate a large degree of inheritance in disk midplanes (Visser et al., 2011; Taquet et al., 2016; Alexander et al., 2018). At the disk surface, we find that both initial abundances give essentially the same results and we use the atomic initial abundances for simplicity.

2.2. Disk model

We adopt a thin disk model with surface density profile $\Sigma(r)=500\times r_{\rm{au}}^{-1}\rm{g\cdot cm^{-2}}$ and midplane radial temperature profile $T_{\mathrm{mid}}(r)=200\times r_{\rm{au}}^{-0.6}\rm{K}$ (Andrews & Williams, 2007; Andrews et al., 2009), where $r_{\rm{au}}$ is disk radius measured in AU. Disk is assumed to be vertically isothermal and the vertical density profile is simply determined from hydrostatic equilibrium

[TABLE]

where $\rho$ is gas density, $H$ is scale height and defined as $H=c_{s}/\Omega_{K}$, where $\Omega_{K}$ is the midplane Keplerian frequency, $c_{s}=\sqrt{k_{B}T_{\rm mid}(r)/\mu m_{p}}$ is the midplane sound speed, with mean molecular weight taken to be $\mu=2.34$.

Multiple ionization sources are considered in disk midplane. Cosmic-rays give an effective ionization rate of $\zeta_{\rm CR}=10^{-17}$ s*-1* exponentially attenuated with column density to the disk surface normalized by 96 g cm*-2* (Umebayashi & Nakano, 1981). A spatially uniform ionization rate by radioactive decay is set to $7.0\times 10^{-19}$ s*-1*. For X-ray ionization $\zeta_{X}$, a fitting formula from Bai & Goodman (2009) to Igea & Glassgold (1999) is used by assuming X-ray temperature $T_{X}$ = 3 KeV and X-ray luminosity $L_{X}$ = 1030 erg s*-1*.

[TABLE]

where NH1,2 is the column density of hydrogen nucleus vertically above and below the point of interest, $\zeta_{1}=6.0\times 10^{-12}s^{-1}$, N1 = $1.5\times 10^{21}{\rm cm^{-2}}$, $\alpha$ = 0.4, $\zeta_{2}=1.0\times 10^{-15}s^{-1}$, N2 = $7.0\times 10^{23}{\rm cm^{-2}}$, and $\beta=0.65$.

2.2.1 Disk model in the surface layer

In this paper, the surface layer is the defined as the region where FUV penetration is high (attenuation is less than $e^{-1}$). This layer is characterized by a photon-dominated chemistry and strong gas-heating. We prescribe a simple gas temperature profile that increases logarithmically with decreasing surface density $\Sigma$ (i.e., column density vertically above the point of interest) as

[TABLE]

Density at disk surface always deviates from Gaussian due to the rapidly varying vertical temperature profile, magnetic pressure support and disk wind. For our purposes, it suffices to simply approximate the gas density by $\rho=\Sigma/H$333At the surface, mean molecular weight varies as more molecules are destroyed, but for our purpose it suffices to approximate it as a constant.. With photochemical reactions being the key ingredients in disk surface chemistry, these approximations should suffice for the purpose of network reduction.

FUV radiation is characterized by a small penetration depth of $\sim 0.01-0.1$ g cm*-2* followed by rapid attenuation (Perez-Becker & Chiang, 2011), and hence is neglected in disk midplane network. Attenuation of FUV results mostly from very small dust grains (e.g. PAHs), whose properties, especially abundances, are highly uncertain (Geers et al., 2006, 2009; Keller et al., 2008). We parameterize the FUV fluxes at the disk surface as a function of radius and column density to disk surface $\Sigma$ to mimic the presence of tiny dust grains

[TABLE]

where the value $10^{7}$ at 1 AU is considered typical for T Tauri stars Bergin et al. (2004), and we adopt $\Sigma_{\rm FUV}=0.01\mathrm{g\cdot cm^{2}}$ as the characteristic FUV penetration depth (assumed to be a constant). While not rigorous, this prescription simply aims to mimic the attenuation of FUV radiation over height from more realistic calculations (e.g., Nomura & Millar, 2005; Walsh et al., 2010), and $\Sigma_{\rm FUV}$ could be larger for lower tiny grain/PAH abundances. The general properties at a representative radius of $R=10$ AU in the surface region are summarized in Figure 1.

2.3. Reduction method

We evolve the full network for a certain time (i.e., 1 Myr), and the reduction is performed based on the species abundances and reaction rates at this time. While in principle the reduction should be performed at different times to capture different stages of chemical evolution, we show in Appendix A that the reduced network constructed at the last time step agrees well with full network throughout almost the entire evolutionary process. This time dependence can be, to some extent (though not self-consistent), considered as to mimic the effect of advection when coupled with dynamics (Semenov et al., 2004), where fluid advection keeps mixing new species into the system to bring it out of equilibrium.

Our network reduction method is iterative and is described as follows. First, we specify a number of species whose abundances we are most interested in. They serve as initial species of the reduced network and new species will be added iteratively following our procedures. The criteria for adding new species is controlled by the sensitivity parameter B[$i$], which measures the importance of the $i$th species to species already selected in the reduced network. The definition of $B[i]$ reads (Tomlin et al., 1992; Ruffle et al., 2002)

[TABLE]

where the summation over $j$ is restricted to all species selected in the current iteration of the reduced network (containing $N$’ species), $n_{i}$ is the number density of the $i$th species, $f_{j}$ is the net formation or destruction rate of the $j$th species, and $g_{j}$ is the larger one between total formation and total destruction rates of the $j$th species (without adsorption and desorption). The term inside the bracket represents the sensitivity of the net production/destruction rate of the $j$th species to the abundance of $i$th species. After the summation, $B[i]$ indicates the sensitivity of the reduced network to the abundance of $i$th species. Larger $B[i]$ means that the $i$th species is more important for determining the abundances of the species in the reduced network. Sensitivity $B[i]$ is calculated only for gas phase species and dust grains but not for mantle species, because mantle species desorption rate can be very temperature sensitive and may overwhelm the $B[i]$ parameter and adversely affects the reduction process. On the other hand, mantle species will automatically be included if their gas-phase counterparts are selected.

After computing $B[i]$ for each species, we sort all species according to their $B[i]$. If the ratio between the smallest value of $B[i]$ within the reduced network and the largest value of $B[i]$ outside the reduced network is higher than some threshold (see below), we consider that the currently selected set of species is sufficiently self-contained, and the identification of necessary species is complete. Otherwise, we incorporate one new species with the largest $B[i]$ into the network (i.e., the most necessary new species), and proceed to the next iteration.

To illustrate this procedure, we show in Figure 2 the sorted distribution of sensitivity B[$i$] in disk midplane at 1AU after 5 iterations (top panel) and after 22 steps when the iteration is completed (bottom panel, see §3.1 for more detail). As more species are included, the difference in $B[i]$ between the selected species and the remaining species becomes more and more prominent, and eventually, a clear cutoff is developed. Quantitatively, we find that 2 to 3 orders of magnitude drop in sensitivity between the least necessary species within the reduced network and the most necessary species outside the reduced network is a good criterion for terminating the iteration.

2.3.1 Application to Our Disk Model

We run network reduction at 1, 10 and 100 AU at both the midplane, and the disk surface with $\Sigma=\Sigma_{\rm FUV}=0.01$ g cm*-2*. We expect these locations to cover broad disk environments in the midplane region (within $\sim 2H$) as well as the surface layer ($\Sigma\lesssim\Sigma_{\rm FUV}$). Note that the intermediate layer in between is neglected in this work. The resulting six reduced networks will then be combined into two. One for the midplane and one for the surface.

The initial reduced network we choose in both disk midplane and disk surface contain electrons, because ionization fraction is an important factor for gas dynamics. In the midplane region, we also incorporate CO in the initial reduced network. Besides being an important gas tracer, CO is also directly related to HCO+, which is a dominant charge carrier and is crucial for disk ionization balance. In the disk surface, we also incorporate major coolants (i.e. CO, H2O, Si, S), whose abundances determines the heating/cooling rates that crucially affect disk mass loss and observational diagnostics. The chemistry of Si and S in the disk surface is dominated by photoreactions and is mostly isolated from other species. Removing them from initial reduced network would in general only affect species directly associated with them, e.g. their ionized counterparts.

Once the species are selected at different radii and combined together, a reduced network can be extracted from full network by keeping all reactions that only involve these species. The next step is a sensitivity analysis based on reaction rates to eliminate redundant reactions (see Tomlin et al. (1992) for more detail). Here we only keep important formation and destruction pathways based on their reaction rates. For each species in the combined reduced network, we only choose the dominant formation and destruction reaction channels with their total reaction rates exceeding 99% percentage of total formation/destruction rates. Dominant reactions are selected at different radii and then combined together for disk midplane and surface regions, respectively.

3. Results

We apply the technique described in previous section and construct reduced chemical networks. In §3.1, we discuss network reduction in the disk midplane regions, and in §3.2, we apply the same technique at disk surface dominated by photochemical reactions. In §3.4, we compare the results with previous works.

3.1. Network reduction in disk midplane

Electron and CO are used as our initial reduced network. In disk midplane, 22, 23, and 10 gas-phase species are selected at 1AU, 10AU, 100AU separately. Adding these necessary species together and keeping only dominant reactions involved by these species, the final reduced network has 28 species and 53 gas-phase reactions (see Table 2 for the species list and Table 4 in Appendix B for the reaction list).

Primary carriers for charges, hydrogen, oxygen, carbon and nitrogen are all selected in the reduced network. In Figure 3, we show the disk ionization fraction from the full network (top panel) and log(ne,reduced/ne,full) (bottom panel) as a function of spatial locations after 1 Myr’ evolution, where $n_{\rm e,reduced}$ are $n_{\rm e,full}$ are the electron number densities from the reduced and full network, respectively. The ionization fraction from our chemical network is consistent with previous works by Ilgner & Nelson (2006); Bai & Goodman (2009) (e.g. See Figure 4 in Bai & Goodman, 2009). The ionization fraction in the reduced network agrees with full network very well within 2 scale heights above/below the midplane. Differences are well within a factor of 2, and are within $\sim 10\%$ in most regions. The largest deviation occurs at higher altitude and within a radius of a few and tens of au. This is not surprising given that original network reduction was conducted at the disk midplane. While the deviations can be reduced by incorporating more species, we also note that the chemistry in these upper layers (i.e., the intermediate layer) is intrinsically more complex (SWH04), which we leave for future work for more detailed investigation that include grain surface chemistry in the original network, and potentially a different reduction scheme.

Besides electrons, the final reduced network also yields relatively accurate abundances for most other species in the reduced network. Figure 4 shows the comparison between reduced network (thin lines) and full network (thick lines) for several ions (top panel) and neutral species (bottom panel) as a function of radius. In the top panel, calculation from the full network shows that dust grains are dominant charge carriers in inner disk ($\sim$ 1AU). This is replaced by H3O+ at $\sim$2AU. In between a few AU and $\sim 40$ AU (where the CO snow line is located), HCO+ is the dominant ion. Beyond the CO snowline, N2H+ abundance increases rapidly and takes over to become the dominant ion, see also (Qi et al., 2013). Near 100AU, H${}_{3}^{+}$ replaces N2H+ and becomes the dominant charged particle. Metal elements Fe and Mg are not present in the reduced network because they almost entirely freeze out onto the grains in the temperature range explored. In a grain-free network where metals are dominant charge carriers, we tested that they will be selected in the reduced network, as expected.

Accurately capturing the ice lines of volatile species is important in many aspects of planet formation (e.g., Öberg et al., 2011; Madhusudhan et al., 2011), whose exact locations are varing throughout PPDs lifetime (e.g. Piso et al., 2015). In the bottom panel of Figure 4, we see a clear transition from gas phase water to water ice at $\sim$2-3AU, gas phase CO to CO ice at $\sim$40AU both in full network and reduced network. This is not surprising because the transition is primarily determined by the balance between desorption and sublimation.

We also compare the vertical profiles of the dominant neutral and ion species in the reduced network (black lines) with those in the full network (red lines), as shown in Figure 5. In the top three panels, we see that dominant neutrals show very good agreement with full network at 1 au and 100 au. At 10 au, and above one scale height at 1 au we do see differences on the order of a few for several species, and at 10 au the discrepancy is above an order of magnitude for CO at the most elevated disk heights considered. The dominant ions, as shown in bottom three panels, are accurately reproduced by the reduced network at all radii and vertical altitudes.

3.2. Network reduction at disk surface

In the disk surface, with prescribed vertical temperature profile, FUV radiation field and additional photoreactions, we perform network reduction at the FUV front where surface mass density is given by $\Sigma=\Sigma_{\rm FUV}=0.01{\rm g\cdot cm^{-2}}$. Important coolants i.e. CO, Si, S, and H2O are adopted together with electron as initial species in the reduced network (e.g. Wang & Goodman (2017a)). By performing the reduction, 13, 16, 16 gas-phase species are selected at 1 AU, 10 AU and 100 AU separately. Merging the results, the final reduced network contains 23 gas phase species and 58 gas phase reactions including photoreactions (see Table 3 for the species list).

Top three panels in Figure 6 show fractional ionization and dominant ions as a function of surface mass density $\Sigma$ at different radii. We see that the ionization fraction first increases dramatically with decreasing $\Sigma$ and approaches a constant value $\sim 10^{-4}$ when $\Sigma$ is around a few times $10^{-2}\rm{g\cdot cm^{2}}$, as carbon becomes fully ionized. For even lower $\Sigma$ and higher temperature, ionization fraction tends to increase further as H+ replaces C+ to become the dominant ions. At all radii, the reduced network well reproduces the ionization fraction and species abundances when $\Sigma\lesssim\Sigma_{\rm FUV}$, i.e. abundances are within a factor of a few with a handful of exceptions.

At higher densities where UV radiation is shielded, the reduced network does not accurately reproduce the ionization fraction because more reactions involving polyatomic species become progressively more important. The chemistry in this intermediate layer is more complex and a separate network reduction is needed. We also conducted a separate network reduction at $\Sigma=0.1$g$\cdot$cm2 where the FUV field is largely shielded. We find that $\sim$ 80 species are required to accurately capture the bulk chemistry in this region, even if only gas-phase chemistry is considered. Similar conclusion was also reached by SWH04. Besides, gas-grain chemistry and non-thermal desorption beyond photodesorption, which we have ignored, are expected to play an increasingly important role in this layer.

Self and mutual-shielding play an important role in the disk surface. We tested the reduction from both within and without H2 self-shielding implemented and obtain the same network. While CO and CI self-shielding are not implemented, we tested reduction with CO photodissociation and CI photoionization turned off (the limiting case assume very effective self-shielding) and again obtained roughly the same network, with only a few different species444More rigorously, which is not proceeded here, one could run network reduction with simplified prescriptions of self and mutual-shielding (e.g., using varied $A_{v}$ values), and then compare the results between full and reduced networks by implementing more complex shielding functions of, e.g., Visser et al. (2009).. These results suggest that the reduced network is reasonably robust and could also be used in models with more accurate self-shielding prescriptions.

Our reduced network well reproduce the abundances of major coolants such as O, OH, CO, Si, and S in the (UV) photon-dominated surface layer, as shown in the middle and bottom panels of Figure 6. Among these major species, only H2O abundance shows relatively large deviation from the full network results by one order of magnitude at $\Sigma\leq 10^{-2}\rm g\cdot cm^{-2}$. This occurs when H2O is no longer the dominant species in these regions, which then carries less weight based on the reduction method. This layer largely covers the bulk and base of the (magneto-)photoevaporative flow (Perez-Becker & Chiang, 2011; Bai et al., 2016; Wang & Goodman, 2017b; Nakatani et al., 2017). Therefore, we expect our reduced network to provide accurate results for simulating such outflows when coupled with dynamics, and it will also help provide important information for observational diagnostics.

3.3. Necessity of Time-dependent Chemistry

The importance of treating the chemistry time-dependently rather than assuming steady state can be assessed by comparing chemical and dynamical timescales: time-dependent chemistry is necessary in dynamical simulations when relevant chemistry timescales exceed the dynamical timescale.

The chemical timescale can vary dramatically from species to species and from reactions to reactions. Here, for our purposes, we give a simple definition that may be considered as a rough estimate for the bulk chemical evolution as follows. We first evolve the full chemical network for 10 million years. We then define the timescale as the time when half of the species’ abundances reaches within 10% compared to their final abundances. Also, we can separate out the ionization timescale, which can be similarly defined as the time the ionization fraction reaches to within 10% of the final ionization fraction.

In Figure 7, we show the bulk chemistry timescale at different radii at disk midplane (black line) and disk surface (blue line), together with the ionization timescale in the disk midplane (red line). Orbital (dynamical) timescale (cyan line) is also plotted. Chemistry timescale in the disk midplane is dominated by neutral gas reactions, while disk surface chemistry timescale is dominated by photochemistry. Although photochemistry is a rather rapid process, volatile species could be gradually converted to less volatile species and frozen onto grains, which results in longer chemical timescales, especially in the outer disk region555The same calculation for the grain-free case shows that the inner disk chemistry timescale is mostly unaffected while the chemistry timescale at 100AU is reduced by  2 order of magnitude, which is still larger than dynamical timescales.. We see that for both regions, the timescale is either comparable or larger compared to orbital timescales, and hence justifying the importance of coupling gas dynamics with chemistry (the only exception being the the surface at 1 au). Compared to neutral species, ionization fraction reaches steady state in relatively smaller timescales, especially in the inner disk. However, in the outer disk (e.g. 100AU), ionization timescale is also comparable with local orbital timescale.

3.4. Comparison with previous work

Our results are generally consistent with the work done previously by SWH04, who started from UMIST95 database, when comparing our results to theirs at the specific radii and disk heights they modeled. The size of our reduced network is similar with SWH04 if we combine their reduced networks at different radii together. They focused on reproducing the chemistry at individual disk locations through location-specific reductions. However, it is not clear how well such a network would function at the regions in between where the individual reduced networks were derived. Our reduced network was developed to work for either the entire midplane or surface, and we have verified that ours agree with the full network for the entire midplane or surface.

In the disk surface, our reduced network is larger than SWH04 for their X-ray and UV dominated chemistry (see Table 7 and 8 in SWH04). This is primarily because besides ionization fraction, we also want to keep dominant coolants in the reduced network and accurately capture their abundances. In the intermediate region where FUV radiation is shielded, SWH04 also claimed that the chemistry is the most complex and in some cases more than 100 species are needed to accurately reproduce the ionization fraction. Our preliminary studies find similar results and we leave for future work for more comprehensive studies of this region.

4. Summary and Discussion

In this paper, we aim to construct reduced chemical reaction networks for PPDs that both reproduce the abundances of key chemical species of interest, and is simple enough to use time-dependently in hydrodynamic/magnetohydrodynamic simulations of PPD gas dynamics. In this context, the chemical species of prime interest are those related to disk ionization (electrons and major ions), thermodynamics (important coolants) and abundant species that affect planet formation (e.g., CO, H2O). The species based method can also be applied to other species of interest, and can help identify important chemical processes associated with those species. Starting from a complex chemical reaction network with 450 gas-phase species and roughly 4000 gas-phase reactions, we follow species-based network reduction method and conduct network reduction at both the midplane and disk surface regions (with photoreactions) over a wide range of radii.

In the disk midplane region, we find that around 28 gas phase species ($\sim$95% reduction) are adequate to reproduce the ionization fraction and the abundances of major species such as CO and H2O over a wide range of radii. Relative error in ionization fraction is well within a factor of two up to two scale heights above or below disk midplane. Dominant ions and neutral species are also well captured by the reduced network in most disk regions. In the disk surface above the FUV front, roughly 20 gas phase species are found to be sufficient to reproduce the ionization fraction as well as the abundances of major coolants.

However, right below the FUV front (i.e., intermediate layer in between the midplane and surface), the chemistry is much more complex, and we find $\sim$ 80 species are needed to adequately reproduce the bulk chemistry in this region. Furthermore this underestimates the complexity since grain surface chemistry is not considered in this study, and they certainly play an important role at intermediate disk heights. Because of this missing link, our reduced network is not yet appropriate in dealing with all problems regarding coupling chemistry with dynamics (e.g., turbulent mixing of radicals from surface towards the midplane (Semenov & Wiebe, 2011)).

We have ignored grain surface reactions other than H2 formation in this work. Non-thermal desorption and self shielding are considered but are treated in a simplified fashion. The natural next step is to study the application of this method to grain-surface chemistry. In the mean time, grain-surface reactions also important in the intermediate layer that we have not covered in this work. This is again left for future studies which will be essential for applications requiring chemistry information across the entire disk vertical extent.

Overall, our results have demonstrated that small scale chemical networks can reasonably reproduce important chemistry in both the midplane and surface layers of PPDs. They can potentially be implemented to be coupled with hydrodynamics/magnetohydrodynamics at affordable computational cost, allowing for more realistic and comprehensive studies of PPDs that couple chemistry with gas dynamics.

We are indebted to an anonymous referee whose carefulness and rigor leads to significant improvements of our manuscript. We also thank D. Semenov, L.I. Cleeves and Fujun Du for useful discussions. XNB acknowledges support from Tsinghua University.

Appendix A A. Time dependent errors caused by reduced network

In this Appendix, we discuss the errors in the reduced network by comparing the reduced and full networks throughout the duration of chemical evolution. This may be considered as a first step to assess the robustness of the reduced network when coupling to dynamics, as fluid advection constantly stirs chemistry out of equilibrium. While we are not doing network reduction at a series of times (to account for different chemical evolutionary stages) and combining the results together, we find that fortunately, the reduced network constructed at the final time ($t=10^{6}$ yr) works equally well at earlier times (see below). This result has already being discussed in Wiebe et al. (2003).

Figure 8 shows the time evolution of the relative errors $\log\left(|n_{\rm reduced}-n_{\rm full}|/n_{\rm full}\right)$ for a number of major species between our reduced network and the full network. We see that for both the disk midplane (top panels) and disk surface (bottom panels) regions, the relative error is smaller than $\sim$ 10% for majority of species (except CO in the disk midplane at 10AU) during all evolution times. While this is encouraging, we comment that further tests (e.g., in 1D) are needed to fully assess the robustness of the reduced network in a dynamical model.

Appendix B B. Reaction lists in disk midplane and in the disk surface

In the disk midplane, the combined reduced network for 1AU, 10AU and 100AU contains 28 gas phase species and 53 gas phase reactions. In the disk surface, reduced network contains 23 gas phase species and 58 gas phase reactions.