© ESO, 2015

1. Introduction

Protoplanetary disks provide the ingredients – dust, gas, and ice – for planets and planetesimals such as comets (for a review, see, e.g., Williams & Cieza 2011). In disks around low-mass stars (≲2 M_⊙), planetary systems are thought to form relatively close to the parent star (≲10 AU); hence, the chemical composition of the inner disk region sets the initial conditions and available elemental components for planetary systems. The molecular material within ≈10 AU is generally dense (≳10⁸ cm^-3) and can reach high temperatures (≳100 K) allowing (ro)vibrational excitation of molecules which emit radiation at near- to mid-infrared (IR) wavelengths. The physical conditions within this region drive the chemistry towards the formation of small, simple, stable molecules most of which, fortunately, also have strong (ro)vibrational transitions.

Near- to mid-IR spectroscopy of nearby primordial protoplanetary disks has demonstrated that the inner planet-forming regions are rich in organic volatiles. The Spitzer Space Telescope allowed the first detection of simple organic molecules in protoplanetary disks at IR wavelengths. Lahuis et al. (2006) reported absorption bands of C₂H₂, HCN, and CO₂ in the spectrum of the low-mass young stellar object, IRS 46, attributed to absorption from hot molecules in a disk within a few AU of the embedded star. The following year, Gibb et al. (2007) detected absorption bands from CO, C₂H₂, and HCN in the disk around a member of the binary system, GV Tau, using NIRSPEC on Keck. Ongoing efforts have detected additional molecules in either emission or absorption, including, OH, H₂O, and CH₄, in several nearby disks using both Spitzer and ground-based facilities (see, e.g., Carr & Najita 2008; Salyk et al. 2008; Pascucci et al. 2009; Pontoppidan et al. 2010; Carr & Najita 2011; Fedele et al. 2011; Salyk et al. 2011; Mandell et al. 2012; Bast et al. 2013; Najita et al. 2013; Gibb & Horne 2013; Pascucci et al. 2013). Because of dust opacity, such observations probe the composition of the disk atmosphere only. It remains unclear whether the atmospheric composition is representative of that of the disk midplane within which planetesimals sweep up the bulk of their material.

Several interesting trends have been noticed in the IR observations. Pascucci et al. (2008, 2009) presented results from a low-resolution (R ≈ 64−128) Spitzer/IRS survey of more than 60 sources ranging from brown dwarfs (T_eff ≈ 3000 K) to Sun-like stars (T_eff ≈ 5000 K). The observations demonstrated an underabundance of HCN relative to C₂H₂ in disks around M dwarfs compared with those around T Tauri stars. The authors postulate this could be due to the difference in far-ultraviolet (FUV) luminosity: M dwarf stars have insufficient FUV photons to dissociate N₂, the main nitrogen reservoir, thereby trapping elemental nitrogen which would otherwise be available to form other nitrogen-containing species, e.g., HCN. The authors have since published additional observations at higher spectral resolution (R ≈ 600) and find the same result: the ratio of C₂H₂/HCN line flux and relative column density decreases with increasing spectral type (Pascucci et al. 2013). They also find that the HCN/H₂O line flux ratios are higher for M dwarf stars than for T Tauri stars leading the authors to conclude that the C/O ratio in the inner regions of M dwarf and brown dwarf disks is higher (≈1) than that for disks around T Tauri stars (<1). Within the subset of T Tauri disks, Najita et al. (2013) postulated that a second trend was present. They found a correlation between the HCN/H₂O line flux ratio and disk mass. Their hypothesis is that planetesimal formation is more efficient in higher mass disks and is able to lock up a significant fraction of oxygen (in the form of water ice) thereby increasing the C/O ratio in the inner region of the disk. In summary, the C/O ratio in the inner regions of protoplanetary disks appears to increase with decreasing spectral type, and within a particular sub class of star-disk systems, to increase with increasing disk mass.

Pontoppidan et al. (2010) conducted a similar survey in the 10−36 μm wavelength range, with a source list also covering Herbig Ae/Be stars (T_eff ≳ 10 000 K). They detected strong H₂O mid-IR line emission from 22 T Tauri stars in their sample, with a detection rate on the order of 2/3; however, to their surprise, no disks in their sample of 25 Herbig stars exhibited water (nor OH) line emission. At near-IR wavelengths, Fedele et al. (2011) conducted a high-resolution spectroscopic survey (L-band) of Herbig Ae/Be disks. OH emission was detected in only four objects, mainly flared disks, and, similar to that found by Pontoppidan et al. (2010), H₂O was not detected. Both sets of authors suggested that the stronger FUV flux from the Herbig Ae/Be stars dissociates molecules in the unshielded inner disk region thereby lowering the line flux in relation to the stellar luminosity. Observations from Herschel support this hypothesis: most of the Herbig Ae/Be disks in the GASPS and DIGIT key programmes exhibit strong OI emission 63 μm with a sub sample also showing OH emission (Meeus et al. 2012; Fedele et al. 2013). Many of these sources also have strong OI emission at 6300 Å (see e.g., Acke et al. 2005). In contrast, only HD 163296 (which is a “flat” or settled disk) has a robust (>3σ) water detection between 50 and 220 μm (Fedele et al. 2012, 2013).

One of the first chemical models which concentrated solely on the inner disk (≲ 10 AU) of a T Tauri star is presented in Markwick et al. (2002); however, the authors neglected the influence of the stellar and interstellar radiation fields on the disk physics and chemistry using the reasoning that viscous heating will dominate the inner disk structure. The initial detections of volatiles in the inner regions of protoplanetary disks (Lahuis et al. 2006; Gibb et al. 2007; Carr & Najita 2008) prompted a flurry of astrochemical models mainly focussed on T Tauri disks (e.g, Agúndez et al. 2008; Willacy & Woods 2009; Woods & Willacy 2009; Walsh et al. 2010; Najita et al. 2011). These models differed somewhat in their level of complexity. Agúndez et al. (2008) adapted a model used for photon-dominated regions (PDRs) but neglected heating due to UV excess emission and X-rays and also assumed that the dust and gas temperatures were equal. Willacy & Woods (2009) and Woods & Willacy (2009) calculated the gas temperature by solving the equation of thermal balance; however, the dust temperature and density were assumed and heating by UV excess emission from the star was neglected (D’Alessio et al. 2006). Najita et al. (2011) adopted a similar approach albeit using a reduced chemical network. Walsh et al. (2010) used a self-consistent model for the protoplanetary disk structure including the effects of heating by UV excess emission and X-rays (Nomura & Millar 2005; Nomura et al. 2007) and a chemical network similar in complexity to the work of Willacy & Woods (2009), excluding deuterium. Several works have also concentrated solely on water production in the inner regions of protoplanetary disks (Glassgold et al. 2009; Bethell & Bergin 2009; Meijerink et al. 2009; Woitke et al. 2009; Ádámkovics et al. 2014; Du & Bergin 2014); however, to date, there has been no detailed study of chemistry in the inner regions of M dwarf nor Herbig Ae stars, especially to address the trends seen in the C₂H₂/HCN line ratios.

In this work, we compute the physical and chemical structure of the planet-forming region (≲10 AU) of a protoplanetary disk around stars of different spectral types: (i) an M dwarf star, (ii) a T Tauri star, and (ii) a Herbig Ae star. Our aim is to investigate whether the stellar radiation field plays a role in the observed trends in C₂H₂/HCN and OH/H₂O ratios derived from the mid-IR observations. Given the proposed importance of N₂ photodissociation for the formation of HCN, we investigate the effect of N₂ self shielding using recently computed shielding functions (Li et al. 2013a) on the subsequent nitrogen chemistry. We also use these models to probe the connection between the “observable” gas emission from the disk atmosphere (at near- to mid-IR wavelengths) with the chemical composition of the midplane within which forming planets and planetesimals sweep up the bulk of their elemental building blocks.

The remainder of the paper is structured as described. In Sect. 2, we outline the methods for computing the disk physical and chemical structure, in Sect. 3 we present our results, and in Sects. 4 and 5, we discuss the implications of this work and summarise the main results, respectively.

2. Protoplanetary disk models

2.1. Physical model

The physical structure of each disk model is calculated using the methods outlined in Nomura & Millar (2005) with the addition of X-ray heating as described in Nomura et al. (2007). Because the methodology is covered in detail in a series of previous publications (see, e.g., Nomura & Millar 2005; Nomura et al. 2007; Walsh et al. 2010, 2012, 2014), here, we highlight the important parameters only.

For each disk model, we assume the disk is steady, axisymmetric, and in Keplerian rotation about the central star. We parametrise the kinematic viscosity via the dimensionless α parameter which scales the maximum size of turbulent eddies by the disk scale height, H, and the sound speed of the gas, c_s, i.e., ν ≈ αHc_s. For protoplanetary disks, α ≈ 0.01. We model the structure of a disk surrounding a star of three different spectral types: an M dwarf star, a T Tauri star, and a Herbig Ae star. We list the stellar mass, M_⋆, stellar radius, R_⋆, and effective temperature, T_⋆, of each host star in Table 1. We also list the adopted mass accretion rate of each star-disk system, Ṁ, and the gas mass surface density at 10 AU, Σ_{10 AU}. We assume the T Tauri and Herbig Ae systems have a mass accretion rate typical for pre-main-sequence stars, ~10^-8 M_⊙ yr^-1. Because the accretion signatures from lower-mass stars are not as strong as for the higher-mass systems, we assume a lower accretion rate for the M dwarf system, ~10^-9 M_⊙ yr^-1 (see, e.g., Herczeg & Hillenbrand 2009). We calculate the dust temperature assuming local radiative equilibrium between the absorption and reemission of radiation by dust grains. For the calculation of the FUV extinction by dust grains, we adopt a dust-grain size distribution which replicates the extinction curve observed in dense clouds (Weingartner & Draine 2001). We calculate the gas temperature assuming detailed thermal balance between the heating and cooling of the gas. We include heating via photoelectric emission from dust grains induced by FUV photons and heating due to the X-ray ionisation of H atoms. The cooling mechanisms included are gas-grain collisions and line transitions. The radiation field of each star is simulated as a black body at the stellar effective temperature (as listed in Table 1) with UV excess emission scaled to the relative mass accretion rates for the M dwarf and T Tauri disks. The stellar FUV (912−2100 Å) radiation field at 1 AU is shown in Fig. 1 for all three central stars. The UV excess emission has two components derived from a best-fit model of the observed TW Hya spectrum: a diluted black-body spectrum to simulate bremsstrahlung emission (T_br ~ 25 000 K) and Lyman-α line emission. The Lyman-α line is modelled as a Gaussian with a width ≈2 Å and the peak flux is determined assuming a line/continuum luminosity ratio of 10³ (see Nomura & Millar 2005, and references therein). Because young stars also exhibit strong X-ray emission, we use a TW Hya-like X-ray spectrum (generated by fitting the observed XMM-Newton spectrum) for the M dwarf and T Tauri stars, with a total luminosity, L_x~ 10³⁰ erg s^-1 (Preibisch et al. 2005) and assume L_x ≈ 3 × 10²⁹ erg s^-1 and T_x ≈ 1.0 keV for the X-ray spectrum of the Herbig Ae star (see, e.g., Zinnecker & Preibisch 1994; Hamaguchi et al. 2005).

2.2. Chemical model

The network we use to calculate the disk chemical evolution includes gas-phase reactions, gas-grain interactions (freezeout and desorption), and grain-surface chemistry.

2.2.1. Gas-phase network

The basis for the gas-phase chemistry is the complete network from the recent release of the UMIST Database for Astrochemistry (UDfA) termed “Rate12” which is publicly available¹ (McElroy et al. 2013).Rate12 includes gas-phase two-body reactions, photodissociation and photoionisation, direct cosmic-ray ionisation, and cosmic-ray-induced photodissociation and ionisation. In this work, the photodissociation and photoionisation rates are calculated by integrating over the specific reaction cross section for each species and the calculated FUV spectrum at each point in the disk (Walsh et al. 2012) using the cross sections from van Dishoeck et al. (2006). Similar to previous work, we have supplemented this gas-phase network with direct X-ray ionisation reactions and X-ray-induced ionisation and dissociation processes (as described in Walsh et al. 2012, and see Sect. 2.2.5). We have also added a set of three-body reactions compiled for use in combustion chemistry models² (see, e.g., Baulch et al. 2005) which are necessary because three-body processes become increasingly important in the inner disk region where the density and temperature are sufficiently high (≳10¹⁰ cm^-3 and ≳1000 K). We have also included at least one collisional dissociation reaction (AB + M → A + B + M) for all neutral species expected to be abundant in the inner disk. We include a small chemical network involving vibrationally excited, or “hot”, H₂. For each gas-phase neutral-neutral reaction involving H₂ which also possesses an activation barrier, we include a duplicate reaction involving hot H₂ with a barrier reduced by the internal energy of the excited H₂ (≈30 163 K, see, e.g., Bruderer et al. 2012, and references therein).

2.2.2. Gas-grain interactions

We allow the freezeout (adsorption) of molecules on dust grains forming ice mantles and the desorption (sublimation) of ices via thermal desorption and photodesorption (Tielens & Hagen 1982; Hasegawa et al. 1992; Walsh et al. 2010, 2012). We adopt the set of molecular binding energies compiled for use in conjunction with Rate12 (McElroy et al. 2013). We have updated the binding energies in light of recent measurements for HCN (Noble et al. 2012). To simplify the calculation of the gas-grain interaction rates, we assume compact spherical grains with a radius of 0.1 μm and a fixed density of ~10^-12 relative to the gas number density. Each grain thus has ~10⁶ surface binding sites. We include photodesorption by both external photons and photons generated internally via the interaction of cosmic rays with H₂ molecules. We use experimentally determined photodesorption yields where available (see, e.g., Öberg et al. 2009a,b,c). For all other species we use a yield of 10^-3 molecules photon^-1. In the calculation of the freezeout rates, we assume a sticking coefficient, S ~ 1, for all species except H, for which we use a temperature-dependent expression which takes into account both physisorption and chemisorption and describes the decreased sticking probability at higher temperatures (Sha et al. 2005; Cuppen et al. 2010b). We assume the rate of H₂ formation equates to half the rate of arrival of H atoms on dust grain surfaces.

2.2.3. Grain-surface network

For completeness, we also supplement our reaction scheme with grain-surface association reactions extracted from the publicly available Ohio State University (OSU) network³ (Garrod et al. 2008). For those species important in grain-surface chemical reaction schemes, e.g., the CH₃O radical, which are not included in Rate12, we also extract the corresponding gas-phase formation and destruction reactions from the OSU network. The grain-surface network has been further updated to include all studied routes to water formation under interstellar and circumstellar conditions (Cuppen et al. 2010a; Lamberts et al. 2013). The grain-surface reaction rates are calculated assuming the Langmuir-Hinshelwood mechanism only, and using the rate-equation method as described in Hasegawa et al. (1992). We limit the chemically “active” zone to the top two monolayers of the ice mantle. We assume the size of the barrier to surface diffusion is 0.3× the binding energy; in this way, volatile species diffuse at a faster rate than strongly bound species. This value lies at the optimistic end of the range determined by recent off-lattice kinetic Monte Carlo simulations of surface diffusion of CO and CO₂ on crystalline water ice (Karssemeijer & Cuppen 2014). This allows the efficient formation of complex organic molecules via radical-radical association reactions at ≳20 K (see, e.g., Vasyunin & Herbst 2013; Walsh et al. 2014). For the lightest reactants, H and H₂, we use either the classical diffusion rate or the quantum tunnelling rate depending on which is fastest (see, e.g., Tielens & Hagen 1982; Hasegawa et al. 1992). For the latter rates, we follow Garrod & Pauly (2011) and adopt a rectangular barrier of width 1.5 Å. We also include reaction-diffusion competition in which the reaction probability is determined by the relative rates between the barrier-mediated reaction and thermal diffusion (see, e.g., Chang et al. 2007; Garrod & Pauly 2011). Although still relatively simplistic, this takes into account the increased probability of reaction in the limit where the thermal diffusion of the reactants away from a common binding site is slow compared with the barrier-mediated reaction rate.

2.2.4. Photodissociation

One further important process now included is a more robust description of the photodissociation rate of N₂ which includes the effects of self shielding and mutual shielding by H and H₂. N₂ is important because it is thought to be the main nitrogen-bearing molecule in interstellar and circumstellar media. Self (and mutual) shielding occurs predominantly for those species which dissociate via line transitions and occurs when foreground material removes photons necessary for dissociation deeper into the cloud. In this way, the photodissociation rate of molecules which can self shield is reduced relative to the dissociation rates for those species which dissociate via the absorption of continuum photons only. H₂ and CO are famous examples of molecules which can self shield (see, e.g., Federman et al. 1979; Glassgold et al. 1985; van Dishoeck & Black 1988; Lee et al. 1996). H₂ also possesses line transitions which overlap with dissociative states of CO and N₂ leading to shielding of CO and N₂ by foreground H₂, a process termed “mutual” shielding.

The inclusion of the self (and mutual) shielding of N₂ is now possible due to the work by Li et al. (2013a) in which they present parametrised shielding functions calculated using a high-resolution model spectrum of N₂. These shielding functions are publicly available for download⁴ to use in astrochemical models. To use the computed shielding functions which are parametrised in temperature and H, H₂, and N₂ column density, one needs to a priori calculate the foreground column densities of each species. In protoplanetary disk models this is not trivial because the dissociating radiation can have multiple sources (e.g., stellar photons and interstellar photons). In addition, the chemistry needs to be calculated in series (from the inside outwards and from the surface downwards) rather than in parallel which can significantly increase computation time. For these reasons, a more pragmatic approach is adopted here, similar to that used in Visser et al. (2011): an “effective” shielding column is calculated based on the calculated FUV integrated flux relative to that assuming no intervening material. The FUV extinction, τ_UV, is given by ${\mathit{\tau }}_{\mathrm{UV}}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{=}\mathrm{-}\ln \left[\frac{{\mathit{G}}_{\mathrm{UV}}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}}{{\mathit{G}}_{\mathit{\star }}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{+}{\mathit{G}}_{\mathrm{ext}}}\right]\mathit{,}$(1)where G_UV(R,Z) and G_⋆(R,Z) are the calculated and the geometrically diluted unattenuated stellar FUV integrated fluxes at a grid point (R,Z), and G_ext is the external unattenuated FUV flux. The effective visual extinction is related to the UV extinction via the empirical relation, ${\mathit{A}}_{\mathrm{v}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}$ ≈ τ_UV(R,Z) / 3.02 mag, and the “effective” H₂ column density is calculated using, ${\mathit{N}}_{{\mathit{H}}_{\mathrm{2}}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}$ ≈ 0.5 × (1.59 × 10²¹) ×${\mathit{A}}_{\mathrm{v}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}$ cm^-2 (Bohlin et al. 1978). For simplicity, the effective shielding column densities for N₂ at each point in the disk are estimated by assuming N₂ has a fixed (rather conservative) fractional abundance, ~10^-5, with respect to H₂. Finally, the photodissociation rate for N₂ is given by, $\begin{array}{ccc}{\mathit{k}}_{\mathrm{ph}}^{{\mathit{N}}_{\mathrm{2}}}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{0}}^{{\mathit{N}}_{\mathrm{2}}}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{\times }{\mathit{\theta }}_{{\mathit{N}}_{\mathrm{2}}}\mathrm{[}{\mathit{N}}_{{\mathit{H}}_{\mathrm{2}}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathit{,}{\mathit{N}}_{{\mathit{N}}_{\mathrm{2}}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathit{,T}\mathrm{(}\mathit{R,Z}{\mathrm{)}}^{\mathrm{]}} {\mathrm{s}}^{-1}\mathit{,}& & \end{array}$(2)where ${\mathit{k}}_{\mathrm{0}}^{{\mathit{N}}_{\mathrm{2}}}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathit{\lambda }}{\mathit{G}}_{\mathrm{UV}}\mathrm{\left(}\mathit{R,Z,\lambda }\mathrm{\right)}\hspace{0.17em}{\mathit{\sigma }}_{{\mathit{N}}_{\mathrm{2}}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\lambda }$ is the unshielded photodissociation rate for N₂ and θ_N2 is the shielding function which is a function of H₂ and N₂ column density and temperature, T.

We adopt a similar method for the computation of the H₂ and CO photodissociation rates using shielding functions calculated by Lee et al. (1996) and Visser et al. (2009), respectively.

2.2.5. X-ray-induced reactions

We include a set of X-ray-induced reactions which we duplicate from the existing set of cosmic-ray-induced reactions contained in Rate12 (McElroy et al. 2013). The reaction rates are estimated by scaling the cosmic-ray-induced reaction rate by the ratio of the local X-ray and cosmic-ray ionisation rates, i.e., k_XR ≈ k_CR × (ζ_XR/ζ_CR). This is a common assumption in chemical models of X-ray irradiated environments (see, e.g., Maloney et al. 1996; Stäuber et al. 2005). The “secondary” X-ray ionisation rate, ζ_XR, is calculated at each grid point in the disk by taking into account the local X-ray spectrum and the explicit elemental composition of the gas (Glassgold et al. 1997). We also include a set of “primary” X-ray ionisation reactions (for further details see Walsh et al. 2012).

2.2.6. Initial Abundances

To generate a set of initial abundances for input into the disk model, we run a dark cloud model (T_gas = T_dust = 10 K, n = 10⁴ cm^-3, and A_V = 10 mag). We use the low-metal elemental abundances from Graedel et al. (1982) supplemented with updated values for O, C, and N based on diffuse cloud observations: 3.2 × 10^-4, 1.4 × 10^-4, and 7.5 × 10^-5 relative to total hydrogen nuclei density, respectively (Cardelli et al. 1991, 1996; Meyer et al. 1998). In this way, we begin the disk calculations with an appreciable ice reservoir on the grain mantle built up over the lifetime of the pre-stellar core prior to disk formation.

In Table 2 we list the abundances (with respect to total H nuclei density) of abundant C-, N-, and O-bearing species at times of 1.0, 3.2, and 10.0 × 10⁵ years. We limit the listed species to those which have an abundance ≳1% that of water ice, the dominant O-bearing species at late times (>10⁵ years). Over the relatively large time steps listed, the trend from atomic to molecular gas can be seen, as can the freezeout of volatile species formed in the gas phase, such as CO and N₂. Ice species formed in situ via hydrogenation of atoms on and within the grain mantle (CH₄, NH₃, and H₂O) show a general trend of increasing abundance towards late times (~10⁶ years) as does the abundance of CO₂ ice which is formed primarily via the reaction between CO and OH. The behaviour of H₂CO and CH₃OH is more complex: although both are formed via the hydrogenation of CO ice, over time, processing of the ice by the cosmic-ray-induced UV field causes a depletion in CH₃OH at late times in favour of H₂CO.

The median relative ice abundances in dense, quiescent clouds measured in absorption against background stars is 100:31:38:4 for H₂O:CO:CO₂:CH₃OH (Öberg et al. 2011). These values are also in line with those measured in low-mass protostellar envelopes which also include measurements for NH₃ and CH₄ (H₂O:NH₃:CH₄ = 100:5:5, Öberg et al. 2011). We opt to use initial abundances at a time of 3.2 × 10⁵ years which corresponds to an ice ratio of 100:15:6:1 for H₂O:CO:CO₂:CH₃OH and 100:7:22 for H₂O:NH₃:CH₄. This is a compromise between the set of abundances at early and late times: at the former, the CO abundance is low compared with observations, whereas at the latter, the CH₃OH abundance is low compared with observations. However, for the physical conditions in the inner disk, we expect the initial abundances to be important only in the midplane of the disk where the chemical timescales can be long compared with the disk lifetime.

3. Results

3.1. Disk physical structure

In Fig. 2 we display the physical structure of each disk as a function of radius and height (scaled by the radius) for the M dwarf disk (left-hand column), T Tauri disk (middle column), and Herbig Ae disk (right-hand column). The lower density in the atmosphere of the Herbig Ae disk is because the scale height of the Herbig disk is smaller than that for the disks around the lower-mass stars ($\mathit{H}\mathrm{=}{\mathit{c}}_{\mathrm{s}}\mathit{/}\mathrm{\Omega }\mathrm{\propto }{\mathit{M}}_{\mathit{\star }}^{-0.5}$, where c_s and Ω are the sound speed and Keplerian angular velocity, respectively). The surface density of the M dwarf disk is around an order of magnitude lower than the other two objects; hence, the lower number density of gas throughout. In all three disks, the gas and dust temperatures decouple in the disk atmosphere such that the gas is significantly hotter than the dust. There are several general trends with increasing spectral type: (i) the gas and dust temperatures increase; (ii) the strength of the FUV flux in the disk surface increases; (iii) the strength of the X-ray flux in the disk surface decreases. For the M dwarf and T Tauri disks, the X-rays penetrate deeper into the disk atmosphere than the FUV photons. For all three disks, the midplane is effectively shielded from all sources of external radiation, including the central star and the interstellar radiation field. The increasing importance of viscous heating in the midplane is indicated by a temperature inversion below which the temperature begins to increase with depth (see the second and third rows of Fig. 2).

3.2. C₂H₂ and HCN

3.2.1. Chemical structure

In Fig. 3 we display the fractional abundance with respect to gas number density of C₂H₂ and HCN as a function of disk radius, R, and disk height divided by radius, Z/R. The dotted and dot-dashed lines represent the τ = 1 surface at 3 and 14 μm, respectively. These are determined using the dust opacity table adopted in the computation of the disk physical structure. Throughout the remainder of the paper, we adopt the ad hoc definition of the disk atmosphere as the material above the τ(14 μm) = 1 surface.

As the stellar effective temperature increases, several trends are evident: (i) the disk molecular layer is pushed deeper into the disk atmosphere; (ii) the fractional abundances of C₂H₂ and HCN decrease in the atmosphere; and (iii) the extent over which both species reach a significant abundance increases in the disk midplane. C₂H₂ and HCN are relatively abundant in the molecular layer of the M dwarf disk reaching maximum fractional abundances of ≈5 × 10^-6 and ≈5 × 10^-5, respectively. The peak C₂H₂ fractional abundance in the molecular layer of the T Tauri disk is ≈1 × 10^-7, whereas that for the Herbig Ae disk is negligible (<10^-11). The peak HCN fractional abundance in the molecular layer for the T Tauri and Herbig Ae disks is around two (~10^-7) and four (~10^-9) orders of magnitude lower than that for the M dwarf disk. Thus, the model results suggest that the relative molecular complexity in the disk atmosphere decreases with increasing stellar effective temperature in line with increased photodestruction.

C₂H₂ reaches a relatively large fractional abundance (≳10^-6) only in specific regions in the midplane of the T Tauri and Herbig Ae disks whereas in the M dwarf midplane it is much lower (≲10^-9). In contrast, HCN reaches a relatively high fractional abundance (≳10^-6) over a greater spatial extent when compared with C₂H₂. An apparent HCN “snow line” moves outwards as the disk midplane temperature increases (see Fig. 3). However, emission from midplane HCN at near- to mid-IR wavelengths is likely obscured by dust in the disk atmosphere as shown by the locations of the τ = 1 surface at 3 and 14 μm in Fig. 3.

3.2.2. Chemistry of C₂H₂ and HCN

A snapshot of the dominant chemical reactions controlling the abundance of C₂H₂ and HCN (and related species) in the disk atmosphere is given in Figs. 4 and 5, respectively. The reactions shown are those which contribute ≳ 10% to the formation and destruction rates at the position of peak fractional abundance in the atmosphere at a radius of 1 AU. The networks are similar to those presented in Agúndez et al. (2008) and Bast et al. (2013) except that we also include the dominant destruction mechanisms.

Free carbon and nitrogen (necessary for incorporation into molecules and radicals such as CH, NH, CN, and C₂) are released from the main gas-phase reservoirs (CO and N₂) via photodissociation for the Herbig Ae disk. For the two cooler disks, both X-ray-induced dissociation and reactions with He⁺ (which itself is produced by direct X-ray ionisation) are dominant. Many reactions involving He⁺ have measured rate coefficients (Adams & Smith 1979; Anicich et al. 1977).

Whether or not a molecule survives in the disk atmosphere with an appreciable abundance requires a delicate balance between formation (via neutral-neutral chemistry or ion-molecule chemistry) and destruction (via photodissociation or X-ray-induced dissociation). The abundance and distribution of HCN in the atmosphere of all three disks is primarily controlled by formation via H₂ + CN (Baulch et al. 1994), and destruction via photodissociation (van Dishoeck et al. 2006). In the Herbig Ae disk, destruction via reaction with atomic hydrogen also plays a minor role. This reaction has a large reaction barrier (12 500 K, Tsang & Herron 1991) and is only significant in very hot gas (≳1000 K). CN has numerous formation routes via neutral neutral reactions: N + C₂ (Smith et al. 2004), NH + C (Brownsword et al. 1996), and NO + C (Chastaing et al. 2000) where NO is formed via the reaction, N + OH (Wakelam et al. 2012). Only the formation of NH (via N + H₂) possesses a substantial reaction barrier (18 095 K, McElroy et al. 2013).

For the M dwarf and T Tauri disks, photodissociation at Lyman-α is more significant than that by the FUV continuum background because around 70−80% of the FUV flux is contained in the Lyman-α line (see also Fogel et al. 2011). This percentage is in line with that determined towards observations of M dwarf and classical T Tauri stars (see, e.g., France et al. 2013, 2014). HCN is treated as though it photodissociates via line transitions (see, e.g., Lee 1994) and it also has a non-negligible photodissociation cross section at 1216 Å (σ = 3× 10^-17 cm², van Dishoeck et al. 2006).

HCN is relatively more abundant in the M Dwarf disk due to the weaker FUV flux leading to decreased destruction via photodissociation. Ion-molecule reactions also contribute to HCN formation in the M dwarf disk (via HCNH⁺ + e⁻) and destruction (via HCN + X⁺, where X⁺ is H⁺, C⁺, H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, H₃O⁺, and He⁺) as shown in Fig. 5. These reactions all have measured rate coefficients (Huntress 1977; Clary et al. 1985; Anicich et al. 1993; Semaniak et al. 2001). X-ray-induced photodissociation also contributes to the destruction of HCN at the level of 15−20%. The ion-molecule formation route is triggered by the formation of CN⁺ and HCN⁺ via the reactions, N + CH⁺, NH + C⁺, and N + CH$\begin{array}{c}\mathrm{+}\\ \mathrm{2}\end{array}$ (Viggiano et al. 1980; Prasad & Huntress 1980). HCN⁺ and HCNH⁺ are then formed via CN⁺ + H₂ and HCN⁺ + H₂ (Raksit et al. 1984; Huntress 1977).

In all cases, midplane gas-phase HCN is synthesised via neutral-neutral chemistry: it is not related to the desorption of HCN ice as is the case for traditional snow lines. The abundance is mediated by formation via H₂ and CN and destruction via collisional dissociation (Baulch et al. 1994; Tsang & Herron 1991). This formation route has a reaction barrier of 820 K and thus requires warm temperatures for activation (>200 K). These temperatures are surpassed in the midplane of each disk due to the inclusion of heating via viscous dissipation (for details see Nomura & Millar 2005, and references therein). Viscous heating dominates over stellar heating in the midplane within radii of ≈0.40, 2.5, and 1.1 AU for the M dwarf, T Tauri, and Herbig Ae disks, respectively. The Herbig Ae transition radius is less than that for the T Tauri disk because the former has significantly stronger stellar heating.

In Fig. 4, a snapshot of the chemistry of C₂H₂ is presented. In all three disks, the formation and destruction of C₂H₂ in the atmosphere is dominated by the neutral-neutral reaction H₂ + C₂H (Laufer & Fahr 2004) and photodissociation (van Dishoeck et al. 2006), respectively. C₂H₂ is also preferentially photodissociated at Lyman-α wavelengths in both cases (σ≥ 4 × 10^-17 cm², van Dishoeck et al. 2006), for similar reasons as discussed above for HCN. In the M dwarf disk, C₂H₂ is also destroyed via reactions with C, H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, and He⁺ (Kim & Huntress 1975; Chastaing et al. 1999; Laufer & Fahr 2004), as well as via X-ray-induced photodissociation (at the level of ≈15−20%). There are also formation routes via ion-molecule chemistry (C₂H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ + e⁻) which contribute at the level of a few percent and are barrierless.

C₂H₂ is significantly less abundant in the atmosphere of the Herbig Ae disk than in the other two disks. In the absence of efficient ion-molecule pathways (triggered by the formation of C⁺ via CO + He⁺), an important first step in the formation of carbon-chain molecules is the formation of CH via C + H₂ which has a large reaction barrier (11 700 K). CH can then react barrierlessly with atomic C to give C₂. In the atmospheres of protoplanetary disks, the gas temperature is controlled by the strength of the UV field. Carbon-chain growth is impeded in the atmosphere of the Herbig Ae disk because the increased photodissociation counteracts the temperature-activated gas-phase chemistry. In Fig. 4, we split the carbon chemistry into that dominated by UV radiation (and a higher gas temperature) highlighted in red on the left-hand side, and that dominated by X-ray radiation (and a lower gas temperature) highlighted in blue on the right-hand side. The left-hand side represents the chemistry more dominant in Herbig Ae disks and the right-hand side represents the chemistry more important in M dwarf disks. The chemistry tends towards the middle of this reaction scheme in environments where the UV radiation is too strong for the survival of species other than CO, C, and C⁺.

In summary, the gas-phase chemistry depends, not only on the strength of the FUV radiation (which controls the gas temperature), but also on the adopted ionisation sources and spectra and corresponding rates propagated throughout the disk. This is especially true for the M dwarf disk in which ion-molecule chemistry also plays a role in the formation and destruction of HCN and C₂H₂ as shown in Figs. 4 and 5. The liberation of free carbon and nitrogen is the principle underlying reason for the importance of X-ray chemistry.

3.2.3. Column densities

A higher C₂H₂ and HCN fractional abundance is seen in the M dwarf disk atmosphere compared with the others; however, this does not necessarily translate into an observable column density, especially given the more tenuous nature of the M dwarf disk and the fact that the molecules peak in fractional abundance higher in the disk atmosphere (where the density is also lower). Figure 6 displays the vertically integrated column densities of both species as a function of disk radius over the entire vertical extent (left-hand column) and down to the τ = 1 surface at 14 μm (right-hand column). There is a general trend that the column density peaks at inner radii, then falls off sharply at a particular radius which moves outwards with increasing stellar spectral type. This fall off occurs at a smaller radius for C₂H₂ than for HCN. This behaviour is also seen in the T Tauri model presented in Agúndez et al. (2008).

For C₂H₂, the T Tauri and Herbig Ae disks achieve a similar peak column density ~10²⁰ cm^-2 at radii of ≈0.2 and 0.3 AU respectively. The column density then falls with increasing radius to ≲10¹⁷ cm^-2 beyond 0.5 AU for the T Tauri star and beyond 1 AU for the Herbig Ae star. The total column density for the M dwarf disk lies orders of magnitude lower, reaching a peak value of ~10¹⁷ cm^-2 at 0.1 AU and falling to ≲10¹⁵ cm^-2 beyond 1 AU. Comparing with the column densities calculated down to the τ(14 μm) = 1 surface, for the two warmer disks, the values are lower by between three and four orders of magnitude because the total column density is dominated by midplane C₂H₂. In contrast, the values for the M dwarf disk remain comparable beyond ≈1 AU because in this case, the total column is dominated by atmospheric C₂H₂. The M dwarf disk has a larger column density in the atmosphere than the T Tauri disk beyond 0.2 AU. Although C₂H₂ does not appear to be abundant in the atmosphere of the Herbig Ae disk (according to Fig. 3), there is still a significant column density: this is because there is a thin layer (only one to two grid cells wide) of relatively abundant C₂H₂ which overlaps with the τ = 1 surfaces. This is likely a feature of the grid resolution of our model and thus is a numerical artefact.

For HCN, the T Tauri and Herbig Ae disks reach a peak column density of ~10²¹ cm^-2 at similar radii to those for C₂H₂. The radial behaviour of the column density then follows the spatial extent of HCN in the disk midplane (see Fig. 3). Similar to C₂H₂, the peak value for the M dwarf disk is lower (~10²⁰ cm^-2 at 0.1 AU), and the column density remains constant beyond ≈0.8 AU at a few ×10¹⁵ cm^-2. For the corresponding values down to the τ(14 μm) = 1 surface, the M dwarf disk has a higher column density than the T Tauri and Herbig Ae disks beyond 0.6 AU and 4 AU, respectively. Although the fractional abundance of HCN is lower in the Herbig Ae disk, the molecular layer is located deeper in the disk atmosphere where the density is higher leading to the apparently large calculated column density of HCN (for R< 4 AU).

3.2.4. Comparison with observed trends

The preferred method for comparing model results with observations is the simulation of the line emission; however, this is a non-trivial matter involving careful consideration of collisional and radiative excitation in the IR (see, e.g., Pontoppidan et al. 2009; Meijerink et al. 2009; Thi et al. 2013; Bruderer et al. 2015). Moreover, the dust properties and size distribution (including, e.g., grain growth) become crucial for the calculation of the emitted spectrum (see, e.g., Meijerink et al. 2009). This is beyond the scope of the work presented here which is focussed on the chemistry; however, this is planned future work. Here, we compare the calculated column densities (and ratios) with those derived from observations to investigate if the chemical models are at least able to reproduce the observed abundances and related trends.

The C₂H₂/HCN column density ratio derived for disks around cool stars ranges from 0.87 to 4.3 (see Table 6 in Pascucci et al. 2013). In contrast, the C₂H₂/HCN column density ratios derived for T Tauri disks lies between 0.006 and 0.43 (see Table 4 in Carr & Najita 2011). Salyk et al. (2011) derive a range from 0.13 to 20 for their sample of T Tauri stars. In Fig. 7, the model ratios in the disk atmosphere are plotted as a function of radius for the T Tauri disk (red dashed lines) and M dwarf disk (gold dotted lines), overlaid with the observed range. The fine lines are the equivalent ratios for the column densities integrated down to the τ = 1 surface at 3 μm. The ratio is flat for the M dwarf disk (~0.1) and lies roughly one order of magnitude lower than the observed range (gold hatched zone). On the other hand, the ratio for the T Tauri disk increases with radius reaching a peak of ≈4. Also, the T Tauri ratios beyond 0.7 AU lie well within the wide range of observed values (red shaded zone).

For HCN, Carr & Najita (2011) derive best-fit column densities for their sample of T Tauri disks between 1.8 and 6.5 × 10¹⁶ cm^-2; however, assuming extremes in the optical depth of the HCN emission expands this to between 0.2 and 31 × 10¹⁶ cm^-2. For C₂H₂, the derived range is 0.02 to 1.6 × 10¹⁶ cm^-2. Salyk et al. (2011) derive column densities between 0.05 and 0.63 × 10¹⁶ cm^-2 and between 0.05 and 1.0 × 10¹⁶ cm^-2 for HCN and C₂H₂ respectively. The observed column density ranges are indicated by the solid colour regions in Fig. 6. The T Tauri model results give good agreement with both the absolute column densities and the ratio of C₂H₂/HCN in the vicinity of the radius within which the observed emission originates.

In the inner region of the disk (≲2 AU), the models qualitatively reproduce the trend that the C₂H₂/HCN ratio is higher in M dwarf disks than in T Tauri disks; however, the ratio predicted in the model lies lower than that observed. Hence, the M dwarf model is either overpredicting the HCN column density or underpredicting the C₂H₂ column density. Pascucci et al. (2013) derive column densities ranging from 4.5 to 50.1 × 10¹⁶ cm^-2 for C₂H₂ and from 1.1 to 18.1 × 10¹⁶ cm^-2 for HCN. The range of observed column densities are indicted by the gold hatched regions in Fig. 6. The model values for C₂H₂ are more than one order of magnitude lower whereas the values for HCN lie within a factor of a few of the observed column densities. However, the absolute values will depend to a degree on the adopted model parameters, such as disk surface density. It is also possible that the excitation mechanisms for C₂H₂ and HCN differ in cool stars relative to T Tauri stars which means that the relative line emission no longer traces the relative abundances nor column densities. The excitation of HCN in protoplanetary disks was investigated in detail by Bruderer et al. (2015) who conclude that HCN abundances derived assuming LTE should differ by no more than a factor of three from those derived assuming non-LTE. They also conclude that C₂H₂ will behave similarly to HCN because of the presence of similar infrared bands through which the excitation can be pumped. We plan to test this via simulations of the molecular emission in future work.

Pascucci et al. (2013) speculate that the self-shielding of N₂ may play a role in determining the C₂H₂/HCN ratio: this can potentially lock up a greater fraction of atomic nitrogen thereby impeding the production of other nitrogen-bearing species, such as HCN. Here, we have also shown that X-ray-induced chemistry may play an important role in both the production and destruction of molecules in the atmosphere of disks around cool stars. We further investigate the role of N₂ self-shielding and X-ray-induced chemistry in Sects. 4.1 and 4.2, respectively.

3.3. OH and H₂O

3.3.1. Chemical structure

Another interesting trend seen in the IR data is the lack of water detections in disks around hotter stars. In Fig. 8 we display the fractional abundance of OH (top row) and H₂O (bottom row) as a function of radius (R) and height divided by the radius (Z/R). OH is more extended and resides in a layer slightly higher in the disk than HCN, C₂H₂, and H₂O, in line with the hypothesis that the OH chemistry is driven by photodissociation. The peak fractional abundance of OH is 8, 5, and 4 × 10^-5 for the M dwarf, T Tauri, and Herbig Ae disks, respectively. The corresponding values for H₂O are 5, 5, and 4 × 10^-4, again showing a general (albeit very shallow) decline in molecular complexity with increasing spectral type. The fractional abundance of OH is negligible (<10^-11) in the midplane of all three disks. The extent over which H₂O is abundant in the disk midplane is related to the thermal desorption of H₂O ice (we assume a binding energy of 5570 K, Fraser et al. 2001). The H₂O “snow line” shifts outwards in radius with increasing stellar spectral type with positions at 0.35, 1.5, and 6.1 AU for the M dwarf, T Tauri, and Herbig Ae disk, respectively. Beyond these radii, H₂O is frozen out on dust grains in the midplane.

3.3.2. Chemistry of OH and H₂O

Figure 9 shows a snapshot of the dominant reactions contributing to the formation and destruction of OH and H₂O (and related species) at the position of peak abundance in the disk atmosphere at a radius of 1 AU. Gas-phase H₂O is predominantly produced in the disk atmosphere via the neutral-neutral reaction, H₂ + OH (using the rate coefficient from Oldenborg et al. 1992). In the cooler M dwarf disk, thermal desorption of water ice dominates at small radii (≈1 AU). Destruction is typically via photodissociation: for the M dwarf and T Tauri disks, this is again dominated by Lyman-α photons (van Dishoeck et al. 2006). In the M dwarf disk, there are additional destruction routes via ion-molecule reactions with H⁺ and H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ (Kim et al. 1974; Smith et al. 1992) and X-ray-induced photodissociation contributes at the level of 12−15% (similar to that found for C₂H₂ and HCN). Similarly, the gas-phase abundance of OH is mediated by production via the photodissociation of H₂O and the reaction between H₂ and O (Baulch et al. 1992), with destruction via photodissociation (van Dishoeck et al. 2006) and reactions with H₂ and H (Tsang et al. 1986; Oldenborg et al. 1992). As in the case for H₂O, dissociation by Lyman-α dominates over that by the background FUV continuum. This can be summarised in the following, rather succinct, reaction scheme, $\mathit{OH,}\hspace{0.17em}\mathit{h\nu }\mathrm{\left]}{\mathit{H}}_{\mathrm{2}}\mathit{OH}\mathrm{\right[}\mathit{H,}\hspace{0.17em}\mathit{h\nu }\mathrm{\left]}{\mathit{H}}_{\mathrm{2}}{\mathit{H}}_{\mathrm{2}}\mathit{O}\mathrm{\right[}\mathrm{desorption}\mathrm{]}\mathrm{freezeout}{\mathit{H}}_{\mathrm{2}}{\mathit{O}}_{\mathrm{ice}}$(3)(see also Bethell & Bergin 2009; Woitke et al. 2009, and Fig. 9). Thus, the water chemistry in the disk atmosphere is rather simple with the relative abundances of O, OH, and H₂O controlled primarily by the relative abundances of H and H₂, the gas temperature (which needs to be sufficiently high to activate the neutral-neutral chemistry) and the strength of the FUV radiation field (necessary for photodissociation). We find that reactions with vibrationally excited or “hot” H₂ (see Sect. 2.2.1) increase in importance towards the upper atmosphere. However, these reactions do not contribute significantly at the position of peak abundance and deeper and hence, do not influence the total column density of neither OH nor H₂O. For a detailed review of water chemistry see van Dishoeck et al. (2013).

3.3.3. Column densities

Figure 10 shows the vertically integrated column densities of OH (top row) and H₂O (bottom row) as a function of radius over the entire vertical extent of the disk (left-hand columns) and down to the τ = 1 surface at 14 μm (right-hand column). The column densities show less structure compared with those for C₂H₂ and HCN. Because the majority of OH is in the disk atmosphere, the column densities for both cases are similar and range between a few times 10¹⁵ to 10¹⁷ cm^-2. The OH column density also does not vary greatly with radius. The Herbig Ae disk generally has the largest column density and the M dwarf has the lowest, although the values for the T Tauri and M dwarf disks are similar beyond ≈1.6 AU. All three disks show very high total column densities of H₂O (≳10²¹ cm^-2) within each respective snow line. The H₂O column densities also show a similar trend to that for C₂H₂ and HCN. The column densities start high in the inner disk then fall off sharply at a distinct radius for each disk (which moves outwards with increasing spectral type). This radius generally lies beyond that for C₂H₂ and HCN. Similar to the case for HCN, gas-phase H₂O in the disk midplane is likely obscured by dust at near- to mid-IR wavelengths. The bottom right-hand panel shows that the column density of “visible” H₂O is on the order of a few times 10¹⁹ cm^-2 within the snow line for the T Tauri and Herbig Ae disks. The column density drops to ~10¹⁶ cm^-2 beyond 1.5 AU for the T Tauri disk. The value for the M dwarf model remains constant over most of the radial extent of the disk at ~10¹⁷ cm^-2.

3.3.4. Comparison with observed trends

Here, we compare the column densities calculated by the models with those derived from the observations to test whether the chemical calculations are able to reproduce the observed abundances and trends. For their sample of classical T Tauri stars in which both OH and H₂O were detected, Salyk et al. (2011) derive OH column densities from 0.04 to 6 × 10¹⁶ cm^-2. The H₂O column densities range from 4 × 10¹⁷ to 1.8 × 10¹⁸ cm^-2. Carr & Najita (2011) derive a wider range of values (4 × 10¹⁷ to 7.9 × 10²⁰ cm^-2). The observed ranges are indicated by the solid colour shaded regions in Figure 10. Despite the model T Tauri disk not being representative of any particular source, there is significant overlap between the calculated column densities in the disk atmosphere and observed values for both species.

For the Herbig Ae disks, only OH has been robustly detected at IR wavelengths (Fedele et al. 2012, with one exception discussed below). The range of observed column densities (derived from VLT/CRIRES data) lie between 1.3 and 20 × 10¹⁶ cm^-2 and the model values lie within a factor of a few of this range (indicated by the striped region in Fig. 10, Fedele et al. 2011). The maximum model H₂O column density in the disk atmosphere is more than an order of magnitude larger than the maximum upper limits derived by both Salyk et al. (2011) and Fedele et al. (2011). OH and H₂O have been detected in the disk around HD 163296 in Herschel/PACS data at far-IR wavelengths with a relative column density of OH/H₂O ≈ 1 and an excitation temperature ranging from 200−500 K (Fedele et al. 2012). However, this emission originates further out in the disk atmosphere (15−20 AU) than that expected at shorter wavelengths (see also the analysis in Fedele et al. 2013).

The Herbig Ae model appears to overproduce gas-phase water in the inner disk atmosphere relative to observations; however, the derived upper limits are dependent on the chosen emitting radius and gas temperature. Figure 11 shows the ratio of OH to H₂O column density as a function of disk radius for the T Tauri model (blue dashed lines) and the Herbig Ae model (green dot-dashed lines). The fine lines represent the equivalent ratios for the columns integrated down to the τ(3 μm) = 1 surface. The model T Tauri ratios lie within the observed range for radii less than ≈1 AU. The ratio for the Herbig Ae model has a value ~10^-3−10^-2 within ≈ 8 AU, beyond which it increases and quickly tends to ≈1 at larger radii. The analysis of the non-detections of OH and H₂O in Salyk et al. (2011) suggest a lower limit to the OH/H₂O ratio of a few × 10^-3 assuming an emitting radius which ranges from a few AU to several 10’s of AU. Fedele et al. (2011) suggest a ratio of > 1−25 within an emitting radius as far out as 30 AU based on near-IR data. The model values are on the cusp of the Salyk et al. (2011) lower limit in the very inner region and tend towards the Fedele et al. (2011) lower limit at 10 AU.

The lack of hot water emission from the innermost regions of disks around Herbig Ae stars remains a puzzle, especially given the hypothesis that water can shield itself from photodissociation by the stellar radiation field at column densities, N(H₂O) ≳ 2 × 10¹⁷ cm^-2 (assuming a photodissociation cross section, σ_H2O ≈ 5 × 10^-18 cm²; Bethell & Bergin 2009; Du & Bergin 2014; Ádámkovics et al. 2014). The column densities predicted here for the Herbig Ae disk are in line with those by other work (see, e.g., Woitke et al. 2009). One possible explanation is that turbulent mixing within the planet-forming region can help sequester water in the midplane where, if the temperature is sufficiently low, it can become trapped as ice on dust grains. This would require that as the dust grain grow, either via ice mantle growth or via coagulation, they become decoupled from the gas and remain in the shielded (and cold) midplane (see, e.g., Stevenson & Lunine 1988). There are two caveats to this theory: the first is that the temperature of the midplane would need to be below ≈150 K within ≈10 AU. The midplane temperature in our Herbig Ae model is generally too high for water to reside as ice on dust grains except beyond ≈6 AU (see Fig. 2); however, the exact temperature profile of the disk is sensitive to numerous factors including, for example, the adopted disk surface density, dust-grain size distribution, and degree of flaring (see later). It also remains to be demonstrated whether large-scale radial mixing between the warm inner midplane and cold outer midplane is a viable mechanism in disks around hotter stars. The second caveat is that one would expect the abundance of OH to also be affected by vertical mixing in the atmosphere since OH and H₂O are chemically coupled. For a given radiation field, the rate coefficients for the primary routes to formation and destruction are known; hence, there may be additional destruction routes for gas-phase water, not yet included in the networks. One potential explanation which remains to be explored in protoplanetary disks is the effect of rotationally excited OH which can be produced via photodissociation of H₂O at Lyman-α wavelengths (see, e.g., Fillion et al. 2001). The rate coefficient for the OH(v,j) + H reaction can be enhanced by several orders of magnitude relative to that for ground state OH (Li et al. 2013b). These reactions may be important for shifting the ratio of O/OH/H₂O in the disk atmosphere (see Eq. (3)).

Pontoppidan et al. (2010) also extensively discuss several hypotheses for the lack of hot H₂O in Herbig Ae disks; including (i) an intrinsically lower abundance by an, as yet, unknown physical or chemical mechanism; (ii) veiling of the molecular features by the strong mid-IR background (if the water line luminosity is a weaker function of stellar spectral type than the continuum); and (iii) a well-mixed disk atmosphere with the canonical gas-to-dust mass ratio of 100. The results presented here show that for a well-mixed atmosphere, the observable column density of water vapour remains high. However, a second outcome of grain growth in protoplanetary disks (in addition to that discussed above) is the greater penetration of FUV radiation which can push the molecular layer deeper into the disk atmosphere (see, e.g., Aikawa & Nomura 2006); hence, a higher fraction of the gas-phase water may be “hidden” from view. This remains to be confirmed specifically for the inner regions of disks. Another factor to consider, and mentioned previously, is the disk gas and dust structure. Disks around Herbig Ae/Be stars have been classified into groups based on the shapes of their SED at mid-IR wavelengths (Meeus et al. 2001). Group I disks are postulated to have a flared structure which allows the disk to capture more FUV photons which increases the gas and dust temperature throughout the disk, whereas group II disks are “flatter”, capture less FUV and are thus much colder. Dust grain growth and settling has been postulated as the reason behind the apparent dichotomy of Herbig Ae/Be disks (see, e.g., Dullemond & Dominik 2004) and it has also been suggested that group I disks may be transitional in nature, i.e., they have evidence (usually confirmed by spatially resolved imaging) of a significant gap in the inner disk (see, e.g., Grady et al. 2015). The disk model we have used assumes a flared disk in hydrostatic equilibrium without a gap: whether significant grain growth and the presence of an inner gap affects the abundance and distribution of gas-phase water remains to be explored.

4. Discussion

4.1. On the importance of N₂ self shielding

To quantify the importance of N₂ shielding in the disk atmosphere, in Fig. 12, we show the fractional abundance of N, N₂, and HCN as a function of Z/R at R = 1 AU. An equivalent figure for R = 10 AU is shown in Fig. A.1. Results from three different models are presented: (i) the fiducial model (black dashed lines, with N₂ shielding and X-rays included); (ii) a model without N₂ shielding (purple lines); and (iii) a model without X-ray-induced chemistry (green lines). The inclusion of shielding has an affect on the relative abundances of N and N₂ in a narrow region of the disk only at both 1 and 10 AU for all disks. The ratio of N/N₂ is generally more affected as the central stellar effective temperature increases. Although the fractional abundance of atomic nitrogen can vary by more than one order of magnitude, this translates to a difference on the order of a factor of a few only for the HCN abundance. The results demonstrate that N₂ shielding alone is not able to account for the change in C₂H₂/HCN column densities and line flux ratios seen from M dwarf to T Tauri stars.

The above conclusion holds for models in which the dust and gas are assumed to be well mixed. If a significant fraction of the dust has grown and settled to the midplane (see, e.g., Dominik et al. 2005), this can lead to a relatively dust-poor disk atmosphere and allow greater penetration of FUV radiation. In that case, the importance of molecular (or self) shielding increases relative to dust shielding (Visser et al. 2009; Li et al. 2013a). Whether the same conclusion holds for M dwarf disks with advanced grain growth and settling remains to be confirmed.

4.2. On the importance of X-ray-induced chemistry

In Fig. 12 (and Fig. A.1) models with and without X-ray-induced chemistry are also shown (black dashed lines versus green solid lines). X-ray-induced chemistry is significantly more important for the M dwarf and T Tauri disks than for the Herbig Ae disk. In both cases and at both radii, the abundance of N₂ is significantly increased in the disk atmosphere when X-ray chemistry is neglected which leads to a decrease in the ratio of N/N₂. The abundance of HCN is also significantly perturbed by the exclusion of X-ray chemistry. The abundance of HCN is increased in the disk atmosphere despite the reduction in available free nitrogen indicating that X-rays are important for HCN destruction even at heights where the FUV field is strong. The higher penetration depth of X-rays versus FUV photons allows an increase in N/N₂ deeper into each disk and hence leads to an increase in HCN in the molecular layer (see also Aikawa & Herbst 1999). For the T Tauri disk at 1 AU, X-rays also help to destroy HCN deeper down towards the disk midplane.

The results show that the inclusion or exclusion of X-ray-induced chemistry can have a profound affect on the position and the value of peak fractional abundance of molecules such as HCN. To investigate whether carbon-bearing species are similarly affected, in Fig. 13, the fractional abundance of C₂H₂ is shown as a function of Z/R at a radius of 1 AU for both the M dwarf and T Tauri disks with and without X-ray-induced chemistry (black dashed lines and orange solid lines, respectively). The results show that X-rays are important in the disk molecular layer for releasing free carbon into the gas phase for incorporation into species such as C₂H₂. The inclusion of X-rays in the M dwarf disk increases the peak abundance of C₂H₂ by more than two orders of magnitude and increases the extent over which C₂H₂ is relatively abundant ≳10^-9 with respect to gas number density. The results for the T Tauri disk are even more extreme with a three orders of magnitude increase in the peak fractional abundance from ≈ 10^-11 to ≳10^-8 when X-ray-induced chemistry is included.

As discussed in Sect. 3.2, the efficacy of the X-ray-induced chemistry is because of the generation of He⁺ which in turn reacts with those molecules robust to photodissociation, CO and N₂. This creates a steady supply of free and reactive atomic atoms and ions for incorporation into other molecules via more traditional ion-molecule chemistry (see, e.g., Herbst 1995). C₂H₂ is more reliant on efficient ion-molecule chemistry for its formation than HCN because neutral-neutral pathways to carbon-chain growth have significant activation barriers (see Fig. 4). The barriers en route to HCN, in comparison, are lower (as discussed in Sect. 3.2). Thus, switching off X-ray chemistry has a larger effect on the magnitude of the peak abundance reached by C₂H₂ in the atmosphere as both ion-molecule and neutral-neutral pathways are inhibited in the cooler disks.

4.3. On the influence of initial nitrogen reservoirs

The results have so far suggested that X-ray-induced chemistry is crucial for efficient molecular synthesis in the disk atmosphere in the planet-forming regions of protoplanetary disks around cool stars. An additional scenario to consider is the effect of the assumed initial abundances at the beginning of the calculation. In the results presented thus far, a set of initial abundances from the output of a dark cloud model were used (see Table 2). The calculations begin with a ratio of N:N₂:N_2ice:NH_3ice equal to 1.0:0.21:0.39:0.26. Taking inspiration from the recent work by Schwarz & Bergin (2014), we run an additional set of calculations in which we assume all species are in atomic form and (i) nitrogen is also in atomic form; (ii) nitrogen begins as N₂ gas; (iii) nitrogen begins as N₂ ice; and (iv) nitrogen begins as NH₃ ice. In this way, we investigate the degree of chemical processing in the disk atmosphere for each initial nitrogen reservoir.

In Fig. 14, the fractional abundances of gas-phase N₂ (top row), NH₃ (middle row), and HCN (bottom row) are shown as a function of Z/R at 1 AU. An equivalent plot for R = 10 AU is presented in Fig. A.2. The results show that the calculated abundances in the disk atmosphere (Z/R≳ 0.1) are independent of the form of the initial nitrogen reservoir at a radius of both 1 and 10 AU for all models. The chemistry in the disk atmosphere has achieved steady state by 10⁶ years and has “forgotten” its origins. On the other hand, the abundances in the disk midplane are very sensitive to the initial nitrogen reservoir. The results for cases (ii) and (iii) show that similar abundances are achieved regardless of whether N₂ begins in gas or ice form. At 1 AU the dominant factor is whether nitrogen begins as NH₃ ice (case (iv)). Here, the abundance of N₂ in the midplane is lower relative to the fiducial model with the difference increasing with increasing spectral type. Correspondingly, the abundance of gas-phase NH₃ is significantly higher increasing by around one order of magnitude compared with the fiducial model. At the other extreme, beginning with all species in atomic form generates the lowest abundance of NH₃ in the disk midplane, especially in the case of the M dwarf and T Tauri disks with differences between 1 and 6 orders of magnitude when compared with the model in which nitrogen begins as ammonia ice. NH₃ ice is thought to be produced in situ on the surfaces of grains or on or within the ice mantle; hence, the efficiency of the conversion from atomic N to NH₃ is very sensitive to temperature and a sufficient flux of both atomic N and H must reside on the surface for the reaction to proceed.

The initial nitrogen reservoir also has an effect on the midplane HCN abundance. For the T Tauri and Herbig Ae disks, at 1 AU, HCN reaches the highest abundance for the atomic model, followed by the model in which N begins as ammonia ice. Both models which have nitrogen initially in the form of N₂ produce the lowest abundance of HCN in the disk midplane. The results at 10 AU show less spread except in the case of the Herbig Ae disk. NH₃ and HCN are both significantly enhanced in the midplane for the model in which N begins as NH₃ ice.

At elevated temperatures, HCN is formed via the reaction between H₂ and CN (Baulch et al. 1994), the latter of which has various routes to formation via atomic nitrogen and nitrogen hydrides, e.g., N + CH or C + NH. Both examples are barrierless reactions (Brownsword et al. 1996; Smith et al. 2004; Daranlot et al. 2013). Conversely, reactions which directly produce CN or HCN from N₂ (e.g., C + N₂ or CH + N₂) have large reaction barriers (≫10 000 K, Baulch et al. 1994; Rodgers & Smith 1996). In the dense, warm midplane, NH₃ is more easily broken apart by cosmic-ray-induced photodissociation than N₂ thus releasing nitrogen hydrides into the gas phase (Gredel et al. 1989; Heays et al. 2014). There is also a direct (and barrierless) route to HCN from NH₃, $\mathit{N}{\mathit{H}}_{\mathrm{3}}\mathrm{+}\mathit{CN}\mathrm{-}\mathrm{\to }\mathit{HCN}\mathrm{+}\mathit{NH}2$(4)(Sims et al. 1994). Hence, the nitrogen chemistry in the disk midplane (and resultant abundance and distribution of N-bearing molecules) is sensitive to whether nitrogen begins in the form of N₂ or NH₃ with the latter resulting in an increase in N-bearing species in the midplane.

4.4. Comparison with previous models

Astrochemical models of the inner regions of protoplanetary disks have often been neglected in favour of the outer disk (≫10 AU), motivated by the larger molecular inventory observed via emission at (sub)mm wavelengths (see, e.g., the recent review by Dutrey et al. 2014). However, focus will return to the planet-forming regions of protoplanetary disks driven by spatially-resolved observations of molecules at (sub)mm wavelengths with the ALMA Large Millimeter/Submillimeter Array (ALMA, see, e.g., Qi et al. 2013) and the launch of the James Webb Space Telescope (JWST) in 2018 (see, e.g., Gardner et al. 2006).

Early models of the planet-forming region spanned a wide range of complexity in both physics and chemistry and have primarily focussed on disks around T Tauri stars. Willacy et al. (1998) used a one-dimensional dynamical model to determine the chemical composition of the midplane from 0.1 to 100 AU. They used a chemical network which included gas-phase chemistry and gas-grain interactions (freezeout and desorption). They conclude that neutral-neutral chemistry is more important than ion-molecule chemistry for controlling the abundances in the midplane, similar to that found here. They also concluded that the chemical composition of the inner midplane was not dependent on the initial molecular abundances. This is a different conclusion to that presented here where we find that the composition of the initial ice reservoir plays a crucial role in the subsequent gas-phase chemistry; however, the model presented here uses a chemical network which is significantly more expansive than that adopted in Willacy et al. (1998) and includes grain-surface chemistry and ice mantle processing.

Markwick et al. (2002) calculated the two-dimensional chemical composition of the inner region assuming a disk heated internally by viscous heating only, and using a similar chemical network to that from Willacy et al. (1998). They also included a simple prescription for the X-ray ionisation rate throughout the vertical extent. In general, Markwick et al. (2002) compute significantly higher column densities for C₂H₂ than found in this work and they also do not produce OH in the disk surface layer. This is because they neglected photodissociation and photoionisation by photons originating from the central star and the external interstellar radiation field, processes that we find are important for governing the chemistry in the “observable” molecular layer in the inner regions.

More recently, Agúndez et al. (2008) explored the chemistry in the inner regions of a disk around a T Tauri star using a model similar to that adopted for photon-dominated regions (PDRs). They find similar conclusions to here: temperature-activated neutral-neutral chemistry helps to build chemical complexity in the disk atmosphere. However, we also find that X-ray driven chemistry is potentially very important for building additional complexity, by releasing free atomic (or ionic) carbon and nitrogen into the gas deeper into the disk atmosphere. Qualitatively, we see the same behaviour in the column densities in the disk atmosphere with radius: the column density starts higher then decreases at a radius specific to each molecule. C₂H₂ decreases first followed by HCN, then H₂O. OH, on the other hand remains flat. The peak column densities calculated for C₂H₂ and HCN by Agúndez et al. (2008) (~10¹⁶ cm^-2) are somewhat lower than those computed here for the T Tauri disk. We also find much larger column densities for both OH and H₂O in the atmosphere. The reason for the particularly low OH column density is unclear but Agúndez et al. (2008) do neglect heating via UV excess emission from the star and also assume that the gas and dust temperatures are equal; whereas we find that the gas is significantly hotter than the dust in the region of the inner disk atmosphere where the molecules reside. The column densities that we calculate for OH and H₂O for the T Tauri disk are also in line with those determined in the work by Bethell & Bergin (2009) and Glassgold et al. (2009); however, unlike Bethell & Bergin (2009), we do not consider dust grain settling and so do not need to invoke H₂O self shielding as a mechanism for explaining the survival of gas-phase water in the disk atmosphere. It is interesting that Bethell & Bergin (2009) (and follow up work by Du & Bergin 2014) also predict an increase in column density of H₂O in the disk atmosphere with increasing stellar FUV luminosity (and corresponding hardening of the radiation field).

4.5. A connection between the disk atmosphere and the planet-forming midplane?

Currently, the only means to probe the inner planet-forming regions of protoplanetary disks is via near- to mid-infrared observations. Future high-angular-resolution observations at (sub)mm wavelengths with ALMA may also elucidate the molecular composition of the inner region; however, due to the high column densities in the inner disk, dust opacity begins to affect even (sub)mm line emission. The higher sensitivity and spectral resolution of JWST will also allow the measurement of absorption features by gases and ices other than water in nearby edge-on protoplanetary disk systems.

Given that forming planets sweep up material primarily from the disk midplane, it is worth exploring to what degree the composition of the disk atmosphere reflects that of the midplane. During the planet-formation process, the molecules accreted (whether gas or ice) are ultimately reprocessed in the planet atmosphere. Recent population synthesis models suggest that the main contribution to the heavy element content in the atmospheres of forming planets are ices accreted during the formation of the planetary embryo and icy planetesimals which are captured during the gas accretion stage and vaporised in the atmosphere (see, e.g., Thiabaud et al. 2014).

In Fig. 15, the percentage contribution of each of the dominant molecular carriers of oxygen (top row), carbon (middle row), and nitrogen (bottom row) are plotted as a function of radius. In Fig. 16, we present the equivalent data for the disk atmosphere only. The dominant carriers are identified as those species which contribute most to the total column density of each element (for the data plotted in Fig. 15) and to the column density in the disk atmosphere, down to the τ(14 μm) = 1 surface (for the case of the data plotted in Fig. 16). The data identified as “Other” refers to the summation over all other species which contribute to the elemental abundance, but which individually do not contribute significantly. The identity of “Other” depends on the location in the disk and the particular disk model but are typically complex organic ices (e.g., H₂CO, CH₃OH, etc.) in the outer regions of the two cooler disks, and gas-phase hydrocarbons (e.g., C_nH_m) in the inner regions of the two warmer disks. For the most part, the plotted species contribute ≈90−100% to the total elemental abundance.

The dominant carriers throughout the vertical extent is set by the locations of snow lines: gas-phase H₂O in the inner regions is superseded by H₂O ice once the snow line is surpassed. An interesting result for the M dwarf disk is that CO₂ ice carries most of the oxygen beyond ≈3 AU indicating an efficient conversion from H₂O ice to CO₂ ice on and within ice mantles. This is not seen in the two warmer disks because the higher temperatures make it difficult for CO to reside on the grains sufficiently long for reaction with OH. Beyond the H₂O snowline, the dominant gas-phase carrier of elemental oxygen is CO: thus the gas-phase C/O ratio seemingly tends towards 1, whereas the ice mantle remains oxygen-rich. However, in each disk, there are regions where molecular oxygen has a non-negligible contribution to the total oxygen column, on the order of 10%. The two-dimensional fractional abundances of important gas-phase oxygen-bearing molecules not discussed in detail in the text, such as CO, O₂, and CO₂, are presented in Fig. A.3.

The dominant carbon carrier is CO; in the outer disk of the cooler M dwarf, CO₂ ice and other complex organic ices such as CH₃OH take over. In the inner region of the T Tauri and Herbig Ae disks, gas-phase CH₄ and other hydrocarbons begin to contribute at the level of 30−40%: the physical conditions in the inner regions of the hotter disks are such that the chemistry tends towards thermochemical equilibrium in which hydrocarbons dominate over CO. In the outer regions of the two hotter disks, again either CO₂ ice or gas also contributes at the level of 30%. For the disk around TW Hya, Bergin et al. (2014) suggest that gas-phase CO is an order of magnitude lower in abundance than that expected if 100% of the freely available (i.e., not contained in refractory material) elemental carbon were locked up in gas-phase CO. Our results for the M dwarf disk show a depletion in CO in the outermost regions consistent with this level of depletion; however, we find that significantly more carbon is contained in CO gas in the warmer disks (≈40−50%, at least within a radius of 10 AU).

The picture for nitrogen is more simple: in all cases, gas-phase N₂ is the primary carrier. In the M dwarf disk, NH₃ contributes at the level of 10% whereas in the two hotter disks, there is a region where the contribution from gas-phase NH₃ approaches 50% which is temperature dependent. This is again because the conditions in the inner disk midplane approach thermochemical equilibrium. The low abundance of NH₃ ice in this region suggests that planetary atmospheres which gain the bulk of their heavy elements from planetesimal accretion will be depleted in nitrogen relative to carbon and oxygen. The two-dimensional fractional abundances of important gas-phase nitrogen-bearing molecules not discussed in detail in the text, such as N₂ and NH₃, are presented in Fig. A.6.

In Fig. 16, the equivalent values for the disk atmosphere are plotted. For the oxygen carriers, the main difference is the increased contribution of gas-phase O₂ to the oxygen budget in the atmosphere (up to ≈50%). Where O₂ dominates the oxygen budget depends on the disk model, moving from the innermost regions of the M dwarf disk to the outermost region of the Herbig Ae disk. O₂ is formed in the atmosphere via reaction between atomic oxygen and OH and destroyed via reaction with atomic hydrogen and carbon, with photodissociation increasing in importance as the spectral type of the star increases. Figure 9 shows the dominant formation and destruction mechanisms for O₂ in the disk atmosphere, and the two-dimensional abundance distributions and column densities of atomic and molecular oxygen are given in the Appendix. Where O₂ dominates, CO₂ also makes a non-negligible contribution in the atmosphere (up to ≈30%).

The story for carbon is more simple: CO dominates in the disk atmosphere over much of the radial extent of all disks, with some contribution from CO₂, as found for the full column values. In the two warmer disks, the gas-phase carbon carriers in the disk atmosphere switch to “Other”, in this case, various gas-phase hydrocarbons, C_nH_m, which is not representative of the total column. Gas-phase HCN also has a non-negligible contribution to both the carbon and nitrogen budget in the atmosphere (up to 30% and 60% respectively).

Similarly, gas-phase N₂ dominates the nitrogen budget throughout most of the disk atmosphere; however, in the innermost regions of the two warmer disks, gas-phase NH₃ and HCN take over. Again, this is not representative of the nitrogen budget throughout the disk vertical extent which is dominated by N₂ (≈90%). In Fig. 5 we also show those reactions responsible for the formation and destruction of N₂ and NH₃ in the disk atmosphere.

In Fig. 17, the C/O ratio is plotted as a function of radius for three cases. Also plotted is the assumed underlying elemental ratio (C/O = 0.44, grey dashed lines). The left-hand panel shows the ratio when summed over the dominant ice reservoirs (H₂O and CO₂ ice). The ice is more oxygen rich than the underlying elemental ratio; however, for the two cooler disks, as CO₂ begins to freezeout and/or form, the C/O ratio of the ice tends towards the underlying ratio. On the other hand, the gas is either representative of the elemental ratio or more carbon rich (as shown in the middle panel). The behaviour of the C/O ratio depends on the locations of the H₂O and CO₂ snow lines which move outwards with increasing stellar spectral type. When H₂O ice freezes out and/or forms, the C/O ratio in the gas increases as oxygen is removed. Once CO₂ ice begins to freeze out and/or form, the gas C/O ratio begins to decrease again and the ice correspondingly becomes relatively more carbon rich. This is most clearly seen for the M dwarf disk. In the final panel, the C/O ratio calculated using “observable” tracers only is shown (CO, CO₂, H₂O, C₂H₂, CH₄, and HCN gas). Comparing the final panel with the middle panel, for the two cooler disks, the gas in the disk atmosphere appears significantly more carbon rich when calculated using just the listed tracers. Notably, the C/O ratio for the M dwarf disk appears ≳ 1 throughout most of the disk atmosphere (≳1 AU). This is primarily due to the presence of gas-phase O₂ in the atmosphere which is another “hidden” reservoir of atomic oxygen (in addition to H₂O ice and CO₂ ice).

The chemical model results suggest that disks around cooler stars might appear more carbon rich without the need for additional sinks (neither chemical nor physical) to account for the depletion of oxygen. However, we stress that detailed radiative transfer calculations are necessary to confirm definitively whether the chemical models replicate the trends seen in the observations. Because this is a non-trivial calculation for a two-dimensional disk structure requiring careful consideration of the dust structure and size distribution, this is beyond the scope of the work presented here but is planned future work. Assuming that the trends in the chemical models do translate into trends in the simulated emission, the derivation of the C/O ratio using CO, CO₂, H₂O, CH₄, C₂H₂, and HCN alone, may overestimate the underlying C/O ratio by up to a factor of two. Similarly, the C/N ratio may be overestimated by between one and several orders of magnitude if observations of HCN and NH₃ in the disk atmosphere alone are used to determine the underlying C/N ratio. In the case that the dominant source of heavy elements in a planetary atmosphere is icy planetesimals rather than the gas, then this overestimation increases to a factor of 10. Although less extreme than the estimation of the C/N ratio, a factor of two is sufficient to incorrectly assume the formation of carbon-rich versus oxygen-rich planetary atmospheres, and the dominant carbon and oxygen carriers subsequently found therein (see, e.g., Madhusudhan et al. 2013).

4.6. On the abundance of O₂ in protoplanetary disks

The models predict that gas-phase O₂ is a significant oxygen reservoir in the inner regions of protoplanetary disks. O₂ may contain up to 10% of the total available oxygen over the full vertical extent of the disk, and this increases to up to 50% when the disk atmosphere only is considered (see Figs. 15 and 16).

O₂ has proven to be a somewhat elusive molecule in interstellar and circumstellar environments. Upper limits determined by Submillimeter Wave Astronomy Satellite (SWAS) and Odin towards nearby molecular clouds constrain O₂/H₂ ≲ 10^-7 (Goldsmith et al. 2000; Pagani et al. 2003). On the other hand, ISO upper limits for the abundance of O₂ ice in dark clouds are much more conservative (O₂/H₂O < 0.6, Vandenbussche et al. 1999). However, gas-phase O₂ has been successfully detected towards two warmer sources, the dense core, ρ Oph A (Larsson et al. 2007; Liseau et al. 2012), and Orion (Goldsmith et al. 2011). Furthermore, Yildiz et al. (2013) report a deep Herschel-determined upper limit towards the low-mass Class 0 protostar, NGC 1333-IRAS 4A (O₂/H₂ < 6 × 10^-9). The authors find that O₂ is absent in both the outer cold envelope and inner hot core and conclude that the material entering protoplanetary disks is likely poor in molecular oxygen (gas and ice).

Early gas-phase only chemical models routinely overpredicted the abundance of O₂ in dark clouds, which was postulated to form primarily via ion-molecule chemistry (see, e.g., Bergin et al. 2000). More modern and sophisticated gas-grain models are able to reproduce the low gas-phase abundance of O₂ provided the conversion of O₂ ice into H₂O ice is included, and the chemistry is allowed to evolve for sufficiently long time scales (see, e.g., Bergin et al. 2000; Roberts & Herbst 2002; Yildiz et al. 2013). Recent laboratory experiments have shown that the hydrogenation of O₂ ice is rapid at low temperatures (Ioppolo et al. 2008; Miyauchi et al. 2008).

The origin of gas-phase O₂ in the inner regions of protoplanetary disks is different to that expected in dark clouds. In both cases, the vital reaction is formation via O + OH. This reaction has been well studied⁵ across the temperature range of interest for circumstellar environments. However, the origin of OH in dark clouds is via the dissociative recombination of H₃O⁺ which is generated via successive proton-donation reactions originating from O + H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$. In the inner regions of disks, OH is generated via the reaction between O + H₂, and atomic oxygen, in turn is released from CO via photodissociation, X-ray-induced dissociation, or via reaction with He⁺. O₂ is destroyed via photodissociation and reactions with C and H to yield CO and OH, respectively (see Fig. 9). These latter two reactions have well-constrained rate coefficients (Geppert et al. 2000; Miller et al. 2005). Gas-phase O₂ is able to persist in the disk atmosphere for the same reasons as gas-phase H₂O: the gas temperature is sufficiently high to activate the required neutral-neutral chemical reactions. The model T Tauri disk generates column densities of OH and H₂O in good agreement with those observed towards T Tauri stars (see Fig. 10); hence, given that the rate coefficients for the important reactions are well studied, the prediction that gas-phase O₂ may also be relatively abundant in T Tauri disks is substantive.

5. Summary

In this work, the chemistry and resulting molecular composition of the planet-forming regions (<10 AU) of protoplanetary disks has been explored, with the aim to investigate potential reasons for the trends seen in near- to mid-infrared observations. The results demonstrate that, as the effective temperature of the central star increases, the molecular complexity of the disk atmosphere decreases, showing that the FUV luminosity of the host star plays an important role in determining the molecular composition of the disk atmosphere. The weaker FUV flux impinging upon disks hosted by M dwarf stars allows molecules to reach relatively high abundances in the atmosphere: X-ray-induced chemistry can further increase molecular complexity by driving a rich ion-molecule chemistry, helping to qualitatively explain the large column densities of small organic molecules, such as, C₂H₂ and HCN seen in M dwarf disk atmospheres. The key process is the liberation of free carbon and nitrogen from their main molecular reservoirs (CO and N₂, respectively) via X-rays. The results shown here, in conjunction with the mid-IR observations, suggest that these objects are good targets for ALMA which can further help elucidate chemistry in disks towards the low-mass low-luminosity regime.

The chemical models suggest that the gas in the inner regions of M dwarf disks is generally more carbon rich than that in disks around T Tauri stars. When the C/O gas-phase ratio is calculated using only observable tracers in the disk atmosphere, then the ratio appears larger than it actually is (and C/O → 1). This is because the models predict that gas-phase O₂ is a significant reservoir of oxygen in the disk atmosphere beyond the water snowline.

The results also demonstrate a degree of chemical decoupling between the disk atmosphere and the midplane. The gas is generally more carbon-rich than the midplane ices (see Fig. 17). This is further corroborated by our studies on the importance of the initial nitrogen reservoir (whether atomic nitrogen, molecular nitrogen, or ammonia). We find that this does not play a role in determining the resulting composition of the observable molecular layer: the chemistry is at steady state. However, the initial reservoir is important for determining the composition in the disk midplane where, generally, the chemical timescales are longer. For example, icy planetesimals forming in disks where much of the initial nitrogen is locked up in NH₃ ice may be more nitrogen rich than those for which nitrogen was contained primarily in the more volatile N₂ (see also Schwarz & Bergin 2014).

Whether the trends seen in the chemical models can reproduce those derived from observations remains to be confirmed via calculations of the molecular emission and will be conducted in future work. Assuming that the chemical trends do translate into observable trends, near- to mid-IR observations of the dominant tracers in the atmosphere may overestimate the underlying C/O and C/N ratios of the gas and ice in the region in which forming planets sweep up most of their material by up to a factor of 10 and more than an order of magnitude, respectively.

Despite the results qualitatively demonstrating some of the observed trends, several issues remain. The large column densities of water vapour and correspondingly strong water emission lines at near- to mid-IR wavelengths predicted by models of Herbig Ae disks have not been confirmed by observations, indicating that there may be a heretofore unconsidered destruction mechanism for gas-phase water at high temperatures. For example, reactions of rovibrationally excited OH with atomic hydrogen may shift the ratio of O/OH/H₂O in the atmosphere. Alternative hypotheses are that the molecular line emission is veiled by the strong stellar continuum emission (as discussed in Pontoppidan et al. 2010) or that the disk dust structure plays an important role.

A further issue is that the M dwarf disk model predicts a C₂H₂/HCN ratio which is an order of magnitude lower relative to the observations. An increase in the underlying C/O or C/N elemental ratio in the disk atmosphere may help explain the high C₂H₂/HCN ratio in M dwarf disks with the enrichment in carbon relative to oxygen and nitrogen caused by vertical or radial mixing. In this scenario, less volatile species (e.g., H₂O and NH₃) are transported to regions where it is sufficiently cold for freezeout onto dust grains (see, e.g., Stevenson & Lunine 1988; Meijerink et al. 2009). If the ice-coated dust grains are sufficiently large, they become decoupled from the gas and settle to the disk midplane. In this way, the gas can become enriched in more volatile species, e.g., CO, which alters the underlying elemental balance of the atmosphere. However, as shown in the results here, molecular nitrogen (which is volatile) is significantly more abundant in the disk atmosphere than NH₃. Hence, this mechanism is unlikely to enrich the disk atmosphere in carbon relative to nitrogen but perhaps the chemistry of nitrogen-bearing species is perturbed by the depletion of oxygen via this mechanism.

The gas-phase chemical network used for HCN and C₂H₂ is (to our knowledge) relatively complete, with experimentally measured and/or critically reviewed reaction rate coefficients used where available. We have shown that self-shielding of N₂ alone is not sufficient to explain C₂H₂/HCN in the inner regions and has a more significant effect in disks around stars with a higher FUV luminosity. It remains to be confirmed whether this also holds for disks in which grain growth and settling have generated a relatively dust-poor atmosphere for which molecular shielding dominates over dust shielding. Instead, X-ray-induced chemistry is more important for releasing atomic nitrogen from N₂ in the M dwarf disk. Furthermore, we find that the exclusion of X-ray-induced chemistry only increases the discrepancy with observation, because the formation of C₂H₂ is primarily driven by ion-molecule chemistry in the M dwarf disk atmosphere, as shown in Figs. 4 and 13. However, all chemical networks suffer from a degree of uncertainty; hence, a systematic sensitivity study of the relative abundance of these two species over the parameter space of physical conditions in protoplanetary disks around cool stars, is worthy of further exploration. This would also confirm the hypothesis presented here, that is, that X-ray chemistry is responsible for the higher C₂H₂/HCN ratio observed in cool stars.

We finish by stating that the future is bright for near- to mid-IR astronomy with Mid-InfraRed Instrument (MIRI, Wright et al. 2004) on JWST and METIS (Brandl et al. 2014), currently being developed for installation on the European Extremely Large Telescope (E-ELT). Both facilities will have significantly higher spectral resolution (~3000 and 100 000 respectively) across the wavelength range of interest for probing the molecular composition of the planet-forming regions of nearby protoplanetary disks.

Online material

Appendix A: Supplementary figures

In Figs. A.1 and A.2, we show results at a radius of 10 AU from the same models described in Figs. 12 and 14.

In Figs. A.3 to A.6 we show the fractional abundances

as a function of radius, R, and height scaled by the radius, Z/R, for those species not discussed in detail in the text. These include important oxygen- (e.g., O₂, CO, and CO₂), carbon- (e.g., C_nH_m) and nitrogen-bearing species (e.g., N₂ and NH₃).

---

Acknowledgments

The authors thank Drs. Alan Heays, Kenji Furuya, and Mihkel Kama for useful discussions and an anonymous referee for their insightful comments. C.W. acknowledges support from the European Union A-ERC grant 291141 CHEMPLAN and from the Netherlands Organisation for Scientific Research (NWO, program number 639.041.335). H.N. acknowledges the Grant-in-Aid for Scientific Research 23103005 and 25400229. She also acknowledges support from the Astrobiology Project of the CNSI, NINS (grant numbers AB261004, AB261008). A portion of the numerical calculations were carried out on SR16000 at YITP in Kyoto University.

References

All Tables

All Figures