© ESO, 2016

1. Introduction

The discovery of more than 2000 extrasolar planets by the radial velocity and transiting techniques (e.g., Udry & Santos 2007; Borucki et al. 2011; Batalha et al. 2013; Fischer et al. 2014) has signaled the next phase in exoplanet research: the characterization of their atmospheres. Simple molecules such as CO, H₂O and perhaps CO₂ and CH₄ are being detected in a growing number of exoplanet atmospheres (e.g., Seager & Deming 2010; Snellen et al. 2010; Birkby et al. 2013; Fraine et al. 2014; Crossfield 2015; Sing et al. 2016). These atmospheres are thought to be built up largely by the accretion of pebbles and planetesimals in the natal protoplanetary disk (see Johansen et al. 2014 and Benz et al. 2014 for reviews), hence the atmospheres should reflect the chemical composition of the disk. There are two different views on how to treat the chemistry in the midplanes of disks, depending on the scientific focus and heritage.

Planet formation and population synthesis models (e.g., Ida & Lin 2004, 2008; Alibert et al. 2013) consider multiple physical effects taking place in a protoplanetary disk, such as gravitational interactions between bodies, orbital excitation and eccentricity damping, gas drag, accretion of material onto planets, and planet migration in the gaseous disk. Hence, there is a high degree of physical complexity and detail to planet formation processes in these simulations. However, these models do not contain any detailed chemistry. Either they simply use the observed chemical abundances in interstellar ices and assume that these abundances are preserved during disk evolution, or they assume that thermodynamic equilibrium is attained so that chemical abundances are controlled by temperature and pressure only (e.g., Mousis et al. 2010; Johnson et al. 2012; Moses et al. 2013; Marboeuf et al. 2014; Thiabaud et al. 2015). The main observational test is through statistical comparisons with the observed populations of exoplanets and their predicted compositions.

The alternative view starts from detailed physico-chemical models of protoplanetary disks which are closely linked to, and tested by, a wide variety of astronomical observations (see reviews by Williams & Cieza 2011 and Armitage 2011). Starting from an assumed (static) surface density distribution, scale height and disk flaring, such models first determine the temperature structure of dust and gas heated by the central star through calculation of the full radiative transfer of the dust and the thermal balance of the gas (e.g., Dullemond et al. 2007; Nomura & Millar 2005; Woitke et al. 2009; Bruderer 2013). This physical model is then coupled with an extensive gas-grain chemistry network to solve the kinetic chemistry equations at each point in the disk and compute the chemical composition of the gas and ice as a function of time (e.g., Bergin et al. 2007; Furuya & Aikawa 2014; Cleeves et al. 2014b; Reboussin et al. 2015; Walsh et al. 2012, 2015). Since planets are formed in the midplanes of disks, it is particularly important to consider the composition and evolution in the midplanes. To what extent is the initial chemical composition of material that is accreted onto a protoplanetary disk preserved, and what happens to the material after it reaches the midplane of the disk, i.e., to what extent is it reset (Visser et al. 2009; Pontoppidan et al. 2014)? Does planetesimal formation happen so fast that ices are incorporated into large bodies early on in the evolution, preventing further chemical processing (Zhang et al. 2015)?

Additional clues to the chemical evolution in disks come from the observations of comets in our own solar system (Mumma & Charnley 2011). Cometary records suggest that the chemical composition of the pre-solar nebula has been at least partially preserved in the comet-forming zone throughout its lifetime, pointing to little or no chemical processing. However, the original composition of the material that was present in the protoplanetary disk around the Sun when it formed remains unknown, and studies of other disks are needed to provide a framework for our own solar system. Particularly interesting are the recent results from the ESA Rosetta mission finding significant amounts of O₂ in comet 67P/C-G (see Bieler et al. 2015), with similarly high O₂ abundances inferred for comet Halley from a re-analysis of the Giotto data (Rubin et al. 2015). Abundances as high as a few % of solid O₂ with respect to solid H₂O are not yet fully understood. Lastly, the deuteration of water and organics also provides insight into the history of the pristine material from the interstellar medium (ISM; see Ceccarelli et al. 2014).

In their planet population synthesis models, Marboeuf et al. (2014) assumed the initial chemical abundances to be inherited directly from the interstellar ices observed in dense interstellar clouds. A set of eight volatile molecules (H₂O, CO, CO₂, CH₄, H₂S, CH₃OH, N₂, and NH₃, species also considered in this work) were homogeneously distributed in their model disk midplane with relative ratios consistent with interstellar ice observations (Gibb et al. 2004; Öberg et al. 2011a; Boogert et al. 2015). Depending on the physical conditions in different parts of the midplane, as well as the sublimation temperatures of the species, these molecules could then either be assigned to the gas or ice, with the threshold set by the icelines of the species. Icelines (or snowlines) mark the radius in the disk midplane beyond which species exist solely in ice form and are thus depleted from the gas. This occurs at the radius where the accretion rate onto grain surfaces (or freezeout) exceeds the desorption rate from grain surfaces due to the negative temperature gradient in the midplane. The relative rates of these processes are very strong functions of temperature leading to a narrow transition region from gas to ice (moving outwards in radius). The position of the midplane iceline for a particular species will depend on its volatility (i.e., its binding energy). Marboeuf et al. (2014) do not consider any chemical reactions in their models, besides freezeout and desorption.

The positions of the icelines are important because they determine which species are gas and ice at any location in the disk, and thus which material is available to build larger bodies (solids only). If, for example, a giant planet is forming in the disk, the composition of its core will reflect the ice compositions at the different positions in the disk through which the forming planet has moved. The composition of the planet’s atmosphere, on the other hand, will reflect the gas composition at the position where the planet becomes massive enough to accrete an atmosphere onto its surface from the surrounding gas in the disk. Moreover, accretion of icy bodies may still pollute the atmosphere. These pebbles and planetesimals migrate through the disk due to radial drift and may therefore have originated at larger radii. Depending on the pebble and planetesimal sizes, the migration of these objects also affects the location of the icelines (see, e.g., Piso et al. 2015).

Particularly important is the C/O ratio of the solid and gaseous material in the disk (Öberg et al. 2011b). The ratio depends not only on the different volatilities of the chemical species but also on their production or destruction as a consequence of chemical processing. Since H₂O and CO₂ (which are both O rich) freeze out at higher temperatures than species that are more C rich, such as CH₄ and CO, the C/O ratio depends on both the physical structure and chemistry in the disk. Ultimately, the chemical composition of a planet’s core and its atmosphere may thus differ depending on the history of the disk, the formation location of the planet, and any subsequent migration.

To address these questions, we use a physical disk model, in particular its midplane temperature and density, which is the same as that considered in the Marboeuf et al. (2014) population synthesis models. We compute the abundances of chemical species with time using a comprehensive chemical network and different sets of assumptions (see below) to investigate the degree to which chemical evolution/processing affects the resulting abundances of key volatiles in the disk midplane. The sensitivity of our results to the choice of (i) initial chemical abundances (parent cloud inheritance or chemical reset); (ii) the physical conditions (in particular ionisation level); and (iii) the types of chemical reactions included in the model, are also investigated, with details provided in Sect. 2. This generates eight different simulations, the results of which are presented in Sect. 3. Section 4 discusses the validity of the inheritance and reset scenarios, the implications for planet formation, and the extent to which the results hold for other disk models. Section 5 summarises the conclusions from this work.

2. Methods

2.1. Physical disk model

The protoplanetary disk is taken from the models of Alibert et al. (2013), Marboeuf et al. (2014), and Thiabaud et al. (2015), which provide the midplane temperature T(R), pressure p(R), and surface density profiles Σ(R) with radius, R. The disk has a parameterised surface density profile¹, $\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{·}{\left(\frac{\mathit{R}}{\mathrm{5.2}\mathrm{AU}}\right)}^{\mathrm{-}\mathit{\gamma }}\mathrm{·}\exp {\left(\frac{\mathrm{-}\mathit{R}}{{\mathit{a}}_{\mathrm{C}}}\right)}^{\mathrm{2}\mathrm{-}\mathit{\gamma }}\mathit{,}$(1)where, a_C = 20 AU, γ = 0.8 (see the prescription in Alibert et al. 2013, Table 1), and a_C and γ are constrained by observations (Andrews et al. 2010). The surface density is Σ(5.2 AU)= 16 g cm^-2 at R = 5.2 AU and the disk is truncated at R_out = 30 AU. The total mass of the disk is $\begin{array}{ccc}{\mathit{M}}_{\mathrm{disk}}\mathrm{=}{\mathrm{\int }}_{{\mathit{R}}_{\mathrm{0}}}^{{\mathit{R}}_{\max }}\mathrm{2}\mathit{\pi R}\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{R}\mathrm{=}\mathrm{1.3}\mathrm{\times }{\mathrm{10}}^{-3}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\approx }\mathrm{0.13}\mathrm{MMSN}\mathit{,}& & \end{array}$with R₀ = 0.05 AU, R_max defined as the radius at which the cumulative disk mass calculated from inside out reaches the total disk mass within 1%, Σ(R) taken from Eq. (1), and MMSN = 1 × 10^-2M_⊙ from Weidenschilling (1977).

The focus here is on the midplane of the disk, where the gas and dust temperatures are assumed to be coupled. Physical disk models have been developed to explain a wide variety of observations where the emission usually arises from higher up in the disk atmosphere, and where gas/grain decoupling for temperature is significant. However, since the disk vertical structure is not relevant for planet formation in the midplane, a midplane-only model is used here.

The radial grid used here consists of 119 points from R = 0.2 AU to R = 30 AU, with radial step sizes of ΔR = 0.1 AU and 1 AU, inside and outside of 10 AU, respectively.

The disk model includes irradiation from a central star with a spectral type similar to the Sun. The temperature profile from Alibert et al. (2013) is slightly adjusted to remove two features: (i) a physically unrealistic drop at 5.2 AU; and (ii) an imposed lower temperature limit of 20 K in the outer disk. The profile used here follows their profile in the range 0.2 AU ≤ R< 5.3 AU, and uses a power-law function T ∝ R^-0.6 to extrapolate 5.3 AU ≤ R ≤ 30 AU, in agreement with full 2D radiative transfer models of protoplanetary disks (see, e.g., Bruderer et al. 2014). This adjustment is shown in Fig. 1a. In this work, the solid blue temperature profile in Fig. 1a is used throughout the disk, whereas the original Alibert et al. (2013) profile in the outer disk is given by the blue dashed profile. No adjustments are made to the pressure profile from Alibert et al. (2013).

The temperature in the model decreases from T = 434 K at R = 0.2 AU to T = 19.5 K at R = 30 AU in the outer disk. The number density over this radial range spans about 5 orders of magnitude, reaching almost n = 10¹⁵ cm^-3 close to the star, and dropping to about n = 10¹⁰ cm^-3 at 30 AU.

Two different levels for the ionisation rate throughout the disk are considered, a low level and a high level. In particular, recent models have shown that cosmic rays can be excluded from disks by the combined effects of stellar winds and magnetic field structures (Cleeves et al. 2013a, 2014a). The purpose of these two levels of ionisation is to investigate the effect on chemical reactions that are driven by ionisation. Figure 1b shows the ionisation level profiles used in the simulations.

The case with low ionisation level considers ionisation originating from the decay products of short-lived radionuclides (henceforth referred to as SLRs) only. Low ionisation is labeled as “SLR”. The implementation of this ionisation source into the simulations is done using a simplified version of the prescription given in Eq. (30) in Cleeves et al. (2013b), ${\mathit{\zeta }}_{\mathrm{SLR}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{2.5}\mathrm{\times }{\mathrm{10}}^{-19}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{)}\left(\frac{\mathrm{1}}{\mathrm{2}}\right){\left(\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathrm{g}\mathrm{cm}}^{-2}}\right)}^{\mathrm{0.27}}\mathit{,}$(2)where ζ_SLR is the SLR ionisation rate per H₂ molecule in s^-1 and Σ(R) is the surface density of a disk at a given radius R (see Eq. (1)). The simplification ignores the original time dependence of the ionisation rate, as given in Cleeves et al. (2013b); however, including the time dependence will change the ionisation rate by no more than a factor of 2. With this low ionisation level, a higher degree of ionisation is obtained in the inner, denser disk midplane compared with the outer disk (see the red dashed curve in Fig. 1b). This is because the SLR ionisation rate scales with the midplane number density of SLRs, which is assumed to be homogeneous and thus scales with Σ(R), see Eq. (1).

On the other hand, the case of high ionisation level (green solid profile in Fig. 1b) considers contributions from the decay products of SLRs and from cosmic rays (henceforth referred to as CRs), originating externally to the disk (high ionisation is labeled as SLR + CR). Such a high ionisation level has been used in many disk models in the midplane, (see, e.g., Semenov et al. 2004). These CRs are able to penetrate the disk to induce UV photons in the disk midplane via ${\mathrm{H}}_{\mathrm{2}}\stackrel{\mathrm{CR}}{\mathrm{-}0\mathrm{\to }}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\stackrel{{\mathrm{e}}^{\mathrm{\ast }}\mathrm{-}}{\mathrm{-}0\mathrm{\to }}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{\to }{\mathrm{H}}_{\mathrm{2}}\mathit{,}$(3)with the resulting photons generated by radiative decay of ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}$ (see Prasad & Tarafdar 1983). The CR-ionisation rate contribution, ζ_CR, is treated by assuming the following parameterised prescription: ${\mathit{\zeta }}_{\mathrm{CR}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{\approx }{\mathit{\zeta }}_{\mathrm{0}}\mathrm{·}\exp \left(\frac{\mathrm{-}\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathrm{96}{\mathrm{g}\mathrm{cm}}^{-2}}\right)\mathit{,}$(4)where ζ_CR is the CR-ionisation rate per H₂ molecule as function of radius R, ζ₀ = 10^-17 s^-1 is the assumed upper limit to the ionisation rate, and the exponential term represents an attenuation effect for high surface densities (see, e.g., Umebayashi & Nakano 2009). Hence, a higher CR-ionisation rate is reached in the outer disk midplane than in the inner disk (opposite to the case with SLRs only), as represented by the blue dotted profile in Fig. 1b. Mutual neutralisation following collisions of ions with grains can lower the ionisation level, as shown by Willacy et al. (1998). This effect is taken into account in this work. Deuterium chemistry, however, is not considered here, but will be addressed in a separate paper with overlapping authors (Furuya et al. 2016).

2.2. Chemical model

A detailed network is used to compute the chemical evolution of all species, in which many different reactions and pathways are included. The gas-phase chemistry is from the latest release of the UMIST Database for Astrochemistry (see McElroy et al. 2013) termed Rate12. Gas-grain interactions and grain-surface chemistry are included (as described in Walsh et al. 2015, and references therein). The chemistry is solved time-dependently at each radial grid point in the disk. The chemical evolution is assumed to be isolated at these grid points with no exchange of material between the grid points during the evolution. The differences in chemistry at the different points are therefore dictated only by the differences in physical conditions T(R), n(R), and ζ(R). Different chemical species have different volatilities, and thus different temperatures below which the phase change from gas to ice occurs. For each species a binding energy, E_b (K), is adopted. These values are given in Table 1.

Two different versions of the chemical network were utilised to obtain better insight into the importance of different processes: a reduced chemical network, and a full chemical network. The types of chemical reactions considered in the network versions are outlined in Table 2. The reduced chemical network comprises gas-phase chemical reactions, freezeout of gas-phase species onto grain surfaces to form ices, and desorption of ice species off grain surfaces back into the gas phase. The full chemical network also contains grain-surface chemistry in addition to the reactions included in the reduced chemical network. For the reduced chemistry this means that a chemical species becomes non-reactive as soon as it freezes out onto the surface of a grain, and that freeze-out and desorption of a species depend on the accretion and thermal desorption rates only. Photodesorption is excluded from the reduced chemical network to enable direct comparison to model results of Marboeuf et al. (2014), whose only desorption mechanism is thermal desorption. For the full chemistry, on the other hand, photodesorption is included, and chemical processing can continue after freezeout. The motivation behind using a reduced and full chemical network, respectively, is to quantify the effects of gas-phase chemistry, and grain-surface chemistry, respectively.

The simulations are run with two different sets of initial abundances: atomic species or molecular species (see Table 1 for an overview of the initial species in each type of input). All abundances in this paper are with respect to the total number of H nuclei. The molecular abundances in Table 1 are values representative of interstellar ices (see Öberg et al. 2011a; Boogert et al. 2015, and Tables 1 and 2 in Marboeuf et al. 2014). For both sets of initial abundances (atomic and molecular) the elemental ratios are consistent.

The choice of these initial abundances is motivated by the following two scenarios about the history of the midplane material. The first scenario is that the material going into protostellar systems is inherited from the cloud out of which the protoplanetary disk collapsed and formed. This scenario is denoted “inheritance”, and it implies that the material has the same composition as found in dark clouds, especially their ices (see Marboeuf et al. 2014; Mumma & Charnley 2011). The second scenario is the case where the material coming from the dark cloud experiences heating events from the protostar (i.e., accretion bursts or regular stellar irradiation). These heating events are assumed to alter the chemistry in disks significantly (see, e.g., Visser et al. 2015). In the extreme case, the chemistry is reset, meaning that the molecules are assumed to be dissociated into atoms out to R = 30 AU, which can then reform molecules and solids in a condensation sequence, as traditionally assumed for the inner solar nebula (e.g. Grossman 1972). Hence, the scenario considering atomic initial abundances is denoted “reset”. Early chemical models of protoplanetary disks often assumed a set atomic initial abundances (e.g., Willacy et al. 1998; Aikawa et al. 1999; Semenov et al. 2004; Vasyunin et al. 2008; Walsh et al. 2010); however, these early models also did not typically include a comprehensive grain-surface network and focussed solely on the gas-phase chemistry. Some recent physico-chemical models do include surface formation of e.g., H₂O ice, via the Eley-Rideal mechanism (e.g., Kamp et al. 2013; Bruderer et al. 2014; Helling et al. 2014). Early models that did include grain-surface chemistry (e.g. Willacy 2007; Walsh et al. 2010) were limited to simple atom-addition stemming from Tielens & Hagen (1982) and Hasegawa & Herbst (1993). Here, a more comprehensive grain-surface network is used which includes radical-radical recombination, atom addition, and also ice processing (Garrod et al. 2008). This work differs from earlier protoplanetary disk models in that we directly compare and quantify the effects of the inheritance versus reset scenarios for a single disk model using a comprehensive gas-grain chemical network.

Important for investigating chemical evolution is also the size of the grains in the disk midplane, which affects the rates of grain-surface chemical reactions, as well as gas-grain interactions such as freezeout. Here, spherical grains are considered. Fixed sizes of r_grain = 0.1 μm = 10^-5 cm and fixed grain number density 10^-12 with respect to H nuclei are assumed. Dust settling from the upper layers of a disk onto the disk midplane is assumed to happen on a timescale shorter than 1 Myr (Dullemond & Dominik 2004; Aikawa et al. 1999), and this settling will increase the dust density relative to that for the gas. On the other hand, dust coagulation may decrease the dust surface area compared with that of standard ISM dust (Dullemond & Dominik 2005). The consequences of assuming fixed values are briefly discussed in Sect. 4.

All simulations are run for t = 1 Myr. The lifetime of a disk before gas dispersal is found to be in the range <1 to 10 Myr based on observations (see Williams & Cieza 2011; Fedele et al. 2011). The time evolution of the results is briefly discussed in Sects. 3 and 4.

3. Results

This section presents the results from the chemical evolution simulations for key volatiles which contribute significantly to the C/O ratio in the gas and ice: CO, CO₂, H₂O, CH₄, O₂, O, HCN and NO (O₂, O, HCN and NO are presented only where relevant). Chemical evolution results for the key nitrogen bearing species N₂ and NH₃ are also presented. Figure 2 presents the results for the reduced chemistry, and Fig. 3 those for the full chemistry. All figures show midplane abundances of chemical species with respect to H_nuc as a function of radial distance (in AU) from the central star. The arrows to the right of each panel indicate the initial abundance levels assumed for the inheritance scenario, see Table 1. Colour coding of arrows matches that of the plots. CO and CO₂ have the same assumed initial abundances.

In Fig. 2a, the icelines for each key volatile considered here have positions at 0.7, 2.6, 16 and 26 AU, with temperatures of 177, 88, 28 and 21 K for H₂O, CO₂, CH₄, and CO, respectively (see icelines marked as dashed vertical lines in Fig. 2a). The CO iceline is furthest out, because CO is more volatile than the other species, due to its low binding energy, see Table 1. The positions of the icelines of these four volatiles are indicated with arrows in the top parts of each plot presented here.

3.1. Reduced chemical network

Figures 2a and c show the final abundances for the key volatiles when assuming molecular initial abundances, for the low ionisation and high ionisation case, respectively. This is the inheritance scenario. Figure 2a, assuming a low ionisation level, shows negligible changes in the ice abundances of each species (i.e., the initial molecular cloud abundances are preserved), and only minor changes to the gas-phase abundances within each iceline (<10% for CO and <30% for CH₄).

For the higher ionisation rate (see Fig. 2c), the picture looks slightly different. Larger changes in gas-phase abundances are seen here for all the species inside their respective icelines: up to 43% for CO₂, up to 62% for CO, and several orders of magnitude for CH₄. However, the H₂O, CO₂, and CO ice abundances outside their respective icelines are preserved with their initial assumed abundances. This is due to the almost instantaneous freezeout onto grains of these species under the cold and dense physical conditions found in the outer disk, in conjunction with the assumed chemical non-reactivity of these species upon freezeout in the reduced form of the chemical network.

The largest difference in the gas-phase abundances in Fig. 2c when compared with those in Fig. 2a is the destruction of CH₄ gas between 1 <R< 15 AU, as well as the production of CH₄ ice, reaching a peak abundance higher than the initial abundance level between 20 and 25 AU, and also the production of O₂ ice from 20 to 30 AU. The differences in the outer disk can be ascribed to the increasing ionisation level here, as seen in Fig. 1b. The higher ionisation rate creates a destruction pathway for abundant gas-phase species, such as CO, which releases a small proportion of free C and O into the gas-phase for incorporation into other C- and O-bearing molecules, such as CH₄ (via the initiating ionisation-dependent reaction between CO and He⁺, see Aikawa & Herbst 1999) and O₂, via gas-phase reactions (also discussed in Walsh et al. 2015). Beyond their respective icelines, CH₄ and O₂ freeze out onto grain surfaces almost instantaneously whereby they become depleted from the gas and are thus protected from further chemical modification. The gas-phase abundances of both species are not preserved within their respective icelines (16 AU for CH₄ and 21 AU for O₂) because these species are also destroyed in the gas-phase by CR-induced photons (see reaction sequence 3).

Figures 2b and d show the results for the low and high ionisation cases (top and bottom panels, respectively) using atomic initial abundances. In this reset scenario, it is assumed that the gas has undergone an extreme heating event which has fully erased the prestellar composition (see Sect. 2.2). For both ionisation levels, H₂O and CO are the dominant gas-phase species in the very innermost region of the disk midplane, R< 0.3 AU. For R ≥ 0.3 AU, gas-phase CO and O₂ are efficiently produced and reach similar abundance levels (~1.6 × 10^-4 for O₂ and ~1.8 × 10^-4 for CO) out to their respective icelines at 21 AU for O₂ and 26 AU for CO. Gas-phase O₂ is much more abundant (by at least 2 orders of magnitude) in this reset scenario than in the inheritance scenario. In the reset scenario O₂ is formed is the gas-phase via O + OH → O₂ + H (see Walsh et al. 2015), and remains in the gas-phase because it is very volatile (binding energy of E_b = 1000 K), and only freezes out at 24 K, see Table 1. Outside the CO iceline, the CO ice abundance reaches the same level as that for the gas inside the iceline. This CO abundance is about half the H₂O initial abundance in the inheritance scenario. The remainder of the available oxygen (63%) in the outer disk remains in atomic form. These are oxygen atoms that have not had sufficient time to form molecules before freezeout. 30 AU is just around the atom oxygen iceline (which is at 29 AU), so both O gas and ice account for the total O abundance of 3.3 × 10^-4, making atomic oxygen the most abundant O carrier here, almost twice as abundant as CO ice. This result is particular to the reduced chemistry case because all species are rendered chemically inert upon freezeout; in reality, atomic oxygen on and within ice mantles is highly reactive (see, e.g., Linnartz et al. 2015).

Between R = 0.2 and 1.6 AU, gas-phase CO₂ is produced for both ionisation levels. The abundance reaches 9.0 × 10^-5 for the high ionisation case at R = 1 AU (Fig. 2d), which is about 5 times higher than the abundance peak for low ionisation at R = 1 AU (see Fig. 2b). For both low and high ionisation rates, CO₂ ice reaches an abundance of only about 10% of the value assumed in the inheritance scenario outside the CO₂ iceline. In addition, H₂O ice is not efficiently produced in the outer disk in this reset scenario. The reset scenario H₂O abundance levels resemble abundance levels of H₂O naturally produced in the gas phase of 10^-7−10^-6 (see Hollenbach et al. 2009). The abundance of H₂O ice is about an order of magnitude larger in the high ionisation case, showing that ion-molecule reactions in the gas-phase are contributing to the formation of water in the absence of grain-surface chemistry. However, the peak H₂O ice abundance of 2.3 × 10^-5 (reached at 30 AU in Fig. 2d) remains more than an order of magnitude lower than the initial H₂O abundance assumed in the inheritance scenario. The high ionisation level in Fig. 2d also aids the formation of CH₄ ice beyond its iceline, reaching a peak abundances of 3.0 × 10^-6. However, as also found for H₂O and CO₂ ice, the maximum abundance reached for CH₄ is only 17% of the assumed initial abundance for the inheritance scenario.

3.2. Full chemical network

Figure 3a shows the chemical evolution results when assuming the inheritance scenario and a low ionisation rate (SLRs only) for the full chemical network. The abundance behaviour in this figure is practically indistinguishable from that for the case using the reduced chemical network (see Fig. 2a), that is, the initial assumed abundances are preserved in both the gas and the ice. For the higher ionisation rate case (SLRs and CRs, Fig. 3c) the picture is different. Chemical processing by cosmic-ray-induced reactions in both gas and ice occurs, and is most noticeable within the icelines of each species. A large reduction in the gas-phase CO abundance beyond 2 AU is seen, in contrast with the reduced chemistry results. This decrease in CO gas (and ice beyond the iceline) coincides with an overall enhancement in CO₂ ice in the outer disk. CO molecules accreting onto the grain surfaces can react with OH radicals produced in the ice via photodissociation of H₂O ice by CR-induced photons in the full chemistry. This produces CO₂ ice in situ on the grain surfaces via the reaction $\mathrm{C}{\mathrm{O}}_{\mathrm{surf}}\mathrm{+}\mathrm{O}{\mathrm{H}}_{\mathrm{surf}}\mathrm{-}\mathrm{\to }\mathrm{C}{\mathrm{O}}_{\mathrm{2}\mathit{,}\mathrm{surf}}\mathrm{+}{\mathrm{H}}_{\mathrm{surf}}\mathit{,}$(5)which is responsible for the rise in the CO₂ ice abundance between 3 and 15 AU (doubling the initially assumed CO₂ abundance). Within this radial region the dust temperature is between 27 K and 78 K, and because H atoms are very volatile (E_b = 600 K) they can rapidly thermally desorb from the grain surface, impeding further grain-surface reactions involving atomic H. In the outermost disk between 20 and 30 AU, the temperature drops to T ≈ 19.5 K, enabling the more efficient retention of H atoms arriving from the gas or produced in situ within the grain mantles, and thereby increasing the relative rate of H₂O ice production via the reaction, $\mathrm{O}{\mathrm{H}}_{\mathrm{surf}}\mathrm{+}{\mathrm{H}}_{\mathrm{surf}}\underset{{\mathit{\gamma }}_{\mathrm{CR}}}{\rightleftharpoons }{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{surf}}\mathit{.}$(6)This reaction is more efficient than Reaction (5) under these colder conditions. Although balanced by destruction due to photodissociation, the abundance of H₂O ice is increased at the expense of CO₂ ice beyond 20 AU.

CH₄ gas is destroyed between 1 and 16 AU in the full chemistry, similar to the case for the reduced chemistry (see Fig. 2c). The reactions responsible for this destruction of CH₄ at 10 AU, are $\begin{array}{ccc}& & \mathrm{C}{\mathrm{H}}_{\mathrm{4}}\stackrel{{\mathit{\gamma }}_{\mathrm{CR}}}{\mathrm{-}0\mathrm{\to }}\mathrm{C}{\mathrm{H}}_{\mathrm{2}}\stackrel{\mathrm{H}}{\mathrm{-}0\mathrm{\to }}\mathrm{C}{\mathrm{H}}_{\mathrm{3}}\stackrel{\mathrm{O}}{\mathrm{-}0\mathrm{\to }}{\mathrm{H}}_{\mathrm{2}}\mathrm{CO}\mathit{,}\\ & & \mathrm{C}{\mathrm{H}}_{\mathrm{2}}\stackrel{\mathrm{NO}}{\mathrm{-}0\mathrm{\to }}\mathrm{HCNO}\stackrel{\mathrm{O}}{\mathrm{-}0\mathrm{\to }}\mathrm{CO}\stackrel{\mathrm{OH}}{\mathrm{-}0\mathrm{\to }}\mathrm{C}{\mathrm{O}}_{\mathrm{2}}\mathit{,}\end{array}$where γ_CR indicates photodissociation by cosmic-ray induced photons. The carbon is thus converted from CH₄ into CO₂ and H₂CO, but the conversions require CRs (i.e. high ionisation). However, H₂CO does not reach an abundance of more than 10^-12, so CO₂ is the main reservoir for the carbon converted from the CH₄ gas, see Fig. 4. Some carbon from CH₄ is also converted into unsaturated hydrocarbons in the gas phase (as also seen in Aikawa et al. 1999), eventually forming gas-phase C₂H₄. This is the reservoir for the converted carbon in the reduced chemistry, because CO₂ ice requires gran-surface chemistry to form, and hence CO₂ ice does not increase in abundance in the reduced chemistry case. For the full chemistry, however, this C₂H₄ subsequently freezes out onto grain surfaces, undergoes hydrogenation, and forms C₂H₆ which reaches a final abundance of 8.6 × 10^-6, approximately an order of magnitude lower than that of CO₂ (1.2 × 10^-4).

The reaction CO + He⁺ reported by Aikawa et al. (1999), turning carbon from CO into CH₄ and other hydrocarbons, is found here to be a relatively minor channel as compared to the reactions between CH₄ and He⁺, and H₃⁺, respectively, turning CH₄ into CO. The result of this is an efficient conversion of carbon from CH₄ into CO. Outside the CO₂ iceline, this conversion is enhanced. This enhancement is due to the accretion of CO onto grain surfaces, and subsequent rapid reaction of CO with OH to form CO₂, a reaction not included in the models of Aikawa et al. (1999).

The destruction of CO gas due to the conversion into CO₂ has been suggested by e.g. Nomura et al. (2016) to explain the depletion of CO gas inside the CO iceline of the protoplanetary disk TW Hya. Kama et al. (2016a) reported overall carbon depletion in the disk atmospheres of some protoplanetary disks using results from a single-dish survey of CO (J = 6−5) and [CI](³P₁ − ³P₀) line emission with APEX. Schwarz et al. (2016) reported CO depletion in TW Hya from CO isotopologue emission observations with ALMA. Detailed modelling confirmed that carbon is likely depleted in the disk around TW Hya by a factor of ≈100 (Kama et al. 2016b; see also Du et al. 2015). This is in good agreement with the CO gas abundance just inside the CO iceline in Fig. 3c, and the results presented here therefore provide a possible explanation for the presence of this inner-disk CO depletion, found inside of, but somewhat mimicking, the actual CO iceline in TW Hya.

At a temperature of 37.7 K (at R = 10 AU), CO₂ freezes out immediately after production. The abundance increase in CO₂ ice follows the destruction of CH₄ gas, as seen in Fig. 4 which presents the time evolution. It is interesting that this effect is seen under these specific physical conditions, with temperatures ranging from 40 to 150 K. This indicates that the CH₄ destruction happens at radii between the H₂O and CH₄ icelines, where H₂O is not chemically active in gas-phase reactions. More generally, for the high ionisation level, the chemical processing becomes significant after a few times 10⁵ yr. That is shown in Appendix A, where jumps in abundance levels are presented for our four key volatiles between 100 kyr and 500 kyr.

When considering initial atomic abundances (Figs. 3b and 3d) a similar radial behaviour (although not identical) is seen when comparing the results for reduced and full chemistry. Inside the icelines, mainly gas-phase H₂O, O₂ and CO are produced, reaching peak abundances of 3.6 × 10^-4, 1.6 × 10^-4, and 1.8 × 10^-4, respectively. Not as much gas-phase CO₂ is produced within the iceline. Comparing the results for low and high ionisation, the peak abundances of gaseous CO₂ are 2.5 × 10^-5 and 6.6 × 10^-5, respectively, and a negligible amount of gas-phase CH₄ is formed. For the high ionisation rate, Fig. 3d, a larger amount of H₂O ice (an order of magnitude at 2 AU) is produced between 1 and 10 AU, than with a low ionisation rate (Fig. 3b), as is also seen in the model using the reduced chemical network. For H₂O ice, a dip in the abundance is seen around 7 AU in Fig. 3d. This is an effect of the competing productions of H₂O ice and SO₂ ice. At 4 AU, SO₂ ice is produced more slowly than H₂O ice, and reaches a final abundance of 2.7 × 10^-6. At 7 AU, however, the lower temperature favours a more efficient and fast production of SO₂ ice, which is able to lock up atomic O, thereby impeding the simultaneous production of H₂O ice.

Can these models produce abundant O₂ in gas or ice? Significant gaseous O₂ is formed in the inner disk starting from atomic abundances, as discussed above. Using the low ionisation rate (see Fig. 3b), O₂ ice reaches an appreciable abundance of 1–10% that of H₂O ice between 15 and 30 AU, similar to that seen in the results using the reduced chemical network. However, little O₂ ice is found when using the higher ionisation rate, indicating that O₂ ice is susceptible to chemical processing by CR-induced photons. In Fig. 3c using molecular initial abundances, O₂ gas is also produced from 1 AU out to 15 AU but in smaller amounts than starting from atomic initial abundances. It reaches a peak abundance of a few percent of that of H₂O ice at 2 AU.

3.3. Main nitrogen reservoirs

Figures 5 shows the final abundances for N₂, NH₃, HCN and NO (HCN and NO only where relevant) for the full chemistry. Results for both low and high ionisation levels, and both gas and ice species, are plotted in each figure. Figure 5a are the results when assuming molecular initial abundances (i.e. the inheritance scenario), whereas Fig. 5b are those assuming atomic initial abundances (i.e. the reset scenario). For the inheritance-scenario with low ionisation, the initial abundances are largely preserved, as for the C- and O-bearing species, and the icelines are of NH₃ and N₂ are nicely outlined at 2.5 and 30 AU, at temperatures of 90 and 20 K, respectively. These icelines are marked with dashed vertical lines in Figs. 5a. However, for high ionisation there is a destruction of NH₃ and production of N₂ at ≈1.5 AU. The reaction pathways responsible for these features are as follows: $\begin{array}{ccc}& & \\ & & \\ & & \end{array}$NH₃ is converted into N₂ both through ion-molecule reactions (reaction sequence 9), as well as through CR-induced photoreactions (reaction sequences 10 and 11). The timescale of this conversion is a few times 10⁵ yr, as was found for the cases of H₂O, CO, CO₂ and CH₄ and discussed earlier in Sect. 3.2.

When considering atomic initial abundances in Fig. 5b, atomic N is seen to form N₂ gas quickly (as also shown and discussed in Schwarz & Bergin 2014), which is the dominant bearer inside 10 AU. This holds regardless of the assumed ionisation level. In the outer disk, HCN ice is the main reservoir, being 5–10 times more abundant than NH₃ ice outside 10 AU. For low ionisation around 1 AU, NO and HCN are the second and third most abundant N-bearing species, although more than an order of magnitude less abundant than N₂. In the outer disk, NH₃, N₂, and NO all reach ice abundances factors of a few to ten times less than the HCN ice abundance. HCN is produced on very short timescales (<1 yr) in the gas-phase via the following reaction sequence, $\mathrm{C}\stackrel{{\mathrm{H}}_{\mathrm{2}}}{\mathrm{-}0\mathrm{\to }}\mathrm{C}{\mathrm{H}}_{\mathrm{2}}\stackrel{{\mathrm{H}}_{\mathrm{2}}}{\mathrm{-}0\mathrm{\to }}\mathrm{C}{\mathrm{H}}_{\mathrm{3}}\stackrel{\mathrm{O}}{\mathrm{-}0\mathrm{\to }}\mathrm{HCO}\stackrel{\mathrm{N}}{\mathrm{-}0\mathrm{\to }}\mathrm{HCN}\mathit{,}$(12)with subsequent freeze-out of HCN onto grains. NO is produced mainly via the well-known gas-phase reaction between N and OH, whereas NH₃ ice is formed in situ on the grain surfaces via hydrogenation of atomic N (see also Walsh et al. 2015). Even with grain-surface chemistry included, when beginning with atomic initial abundances, NH₃ ice formation is less efficient than that for those species more reliant on gas-phase formation followed by freezeout.

4. Discussion

4.1. Compositional diversity at different radii: inheritance vs. reset

The third column of Table 3 shows the final abundances for the full chemical model with low ionisation level and molecular initial abundances at three different radii throughout the disk midplane (see Fig. 3a). This is the model for which the results show best preservation of the initial abundances (listed in Table 1 for the inheritance case). The remaining columns show the fractional abundances of each species at each radial distance for all other models using full chemistry relative to this reference model. The fractional abundances have been calculated using the formula given in the footnote of Table 3. Fractional abundance values smaller than 1 means abundances lowered relative to the reference model whereas values larger than 1 means increased abundances. Values of 0 indicate abundances more than two orders of magnitude lower than for the reference model, so as to avoid large numbers in the table. A value of 1 means no change from the reference model. Three different radii are considered which span the radius of the disk and probe three distinct regions: (i) H₂O-ice rich only (1 AU); (ii) H₂O-ice rich and CO₂-ice rich (10 AU); and (iii) volatile-gas poor (30 AU, i.e., gas fully depleted of volatiles via freezeout).

At R = 1 AU, CO gas is enhanced in all three models by factors of 1.5 to 2.5, with a greater enhancement seen in the atomic case. This is generally at the expense of all other considered species: H₂O ice and gas, CH₄ gas, and CO₂ gas. Both CH₄ gas and H₂O ice show extreme depletion for the atomic cases, for the reasons discussed in Sect. 3.

Moving to 10 AU, i.e., beyond the CO₂ iceline, a different behaviour is seen. CO₂ ice is enhanced by factors of 1.5 to 1.8 in the other three models, this time generally at the expense of CO and CH₄ gas and H₂O ice. The extreme depletion of water ice in the atomic cases also extends into this region. However, an enhancement of a factor 1.3 is seen in CO gas for the atomic case and high ionisation, showing that the higher ionisation rate can also facilitate gas-phase formation as well increased ice processing.

Moving outwards to 30 AU, where most volatiles have accreted onto grain surfaces, different behaviour is seen yet again. This time CH₄ ice is enhanced by a factor of 1.7 to 2.4 at the expense of CO gas and ice (with the latter depleted by almost an order of magnitude) for the high ionisation cases. H₂O and CO₂ ice are enhanced and depleted by up to 22% and 17% respectively, with high ionisation only. Little change is seen for low ionisation, despite beginning the calculation with atomic abundances. For this scenario, it appears that in the outer disk (30 AU), the resulting abundances of CO₂ and H₂O ice reproduce the inherited values to within a few percent.

Finally, it is important to recognize that at all radii and in all scenarios, including the full reset case, the chemistry is not in thermodynamic equilibrium, i.e., the abundances are not simply set by the overall elemental abundances and pressure (see also the discussion in Henning & Semenov 2013).

4.2. C/O ratio

The results in Table 3 show that for different sets of assumptions or models setups, a large diversity is seen in the resulting computed abundances of dominant C- and O-bearing volatiles. The deviations are also radially dependent, and align with the positions of icelines. The fractional deviations are sufficiently large to affect the C/O ratio in the gas and ice which will go into forming the building blocks of planets.

Figure 6 shows the C/O ratios for gas (dashed) and ice (solid) species in the four different model setups using the full chemistry. In the calculation of the ratio, all dominant C- and O-bearing species were taken into account. In addition to the main considered volatiles (CO, H₂O, CO₂ and CH₄), these also include CH₃OH, O₂, HCN, NO, and atomic O where appropriate. The horizontal lines indicate the canonical C/O ratio (at 0.43) and C/O = 1.

For the inheritance scenario and low ionisation (Fig. 6a), the C/O ratio for the gas resembles a step function. Moving outwards in radius, the steps coincide with the icelines of H₂O, CO₂, and CH₄, respectively. This profile is very similar to that in Fig. 4 in Öberg et al. (2011b) in which the C/O ratios in the gas and ice were assumed to be dictated solely by the positions of icelines. Only for the case of inheritance and low ionisation level is there a region in the disk, between 3 and 16 AU, in which the gas is carbon rich, i.e., C/O > 1. The ice ratio for the same model shows that overall the ice remains carbon poor (or oxygen rich), but does become relatively more carbon rich moving outwards in the disk as the icelines for CO₂ and CH₄ are surpassed.

Considering the inheritance scenario with the high ionisation level (Fig. 6c), the chemical processing induced by CRs has a noticeable affect on the C/O ratio in the gas and ice. The C/O gas-phase ratio remains less than 1 everywhere due to the destruction of CH₄ gas and the production of O₂ gas. The gas-phase C/O ratio appears to increase significantly beyond 26 AU; however, this is beyond the CO iceline where most molecules (except H₂) are depleted from the gas. The C/O ice ratio looks similar to that for the case of low ionisation. The ratio is a bit higher for the low ionisation case within 3 AU (73% at 1.6 AU), and slightly lower beyond 3 AU (15% at 10 AU) reflecting the repartitioning of atomic carbon and oxygen from H₂O ice into CO₂, CO, and O₂ (see Fig. 3). Note that Helling et al. (2014) conducted a similar study and found that chemistry induced by cosmic rays is important for the C/O ratio in a disk midplane in both gas and ice; however, because they only included limited grain-surface chemistry, they observed stronger effects in the gas phase C/O ratio (C/O → 1 between the H₂O and CO₂ icelines) than in this work (see Figs. 6b and d).

Turning attention to the reset scenarios in Figs. 6b and 6d, the picture changes significantly. Between 3 and 16 AU the C/O ratio in the ice is higher than that in the gas (although both ice and gas remain carbon poor). This is opposite to both inheritance cases; the dominant carriers of gas-phase C and O in these models in this region are CO and O₂, and the dominant ice component is CO₂ ice as opposed to H₂O ice. Thus the C/O ice ratios tends towards ≈0.5 whereas in the gas it tends towards ≈0.3.

The peak at ≈2.5 AU for the profile in Fig. 6d is due to HCN ice being produced and freezing out at a slightly higher temperature than CO₂. This causes an increase in the C/O ice ratio in this local region (see Fig. 5b). It is only seen in the low ionisation case because N-bearing species formed via gas-phase chemistry in the inner disk (<3 AU) are able to survive to 10⁶ yr. The presence or otherwise of this large peak is dependent upon the relative binding energies of HCN and CO₂ assumed in the model; recent measurements of thermal desorption of pure HCN ice do derive a higher binding energy than both CO₂ and NH₃ (e.g., Noble et al. 2013), but in a mixed ice they may be more similar.

The C/O ratio plots in Fig. 6 highlight that the different model setups considered here (especially the inheritance versus reset scenarios) result in very different C/O ratio profiles for the material in the planet-forming regions of the disk midplane.

4.3. Implications for planet formation and comets

4.3.1. Giant planet atmospheres

Giant planets can accrete their atmospheres either directly from the surrounding gas or through accretion of icy planetesimals, or both. In addition, radial migration of a forming planet can influence the makeup of the resulting atmosphere as the planet moves through a gradient in gas and icy planetesimal composition. Figure 6a suggests that if only accretion from the gas is considered, a gas-giant planet forming between 3 and 16 AU may be able to form a carbon-rich atmosphere in the case of inheritance with low ionisation rates. For high ionisation rates, the atmosphere would still be relatively carbon rich with respect to the canonical ratio (horizontal dashed lines at C/O ratio 0.34 in Figs. 6), yet with a ratio which remains <1. In both cases, the atmosphere can be polluted by accreting icy planetesimals which are oxygen rich, lowering C/O in the planet’s atmosphere.

For the reset scenario, a gas-giant planet forming between 3 and 16 AU would accrete C-poor gas; however, if the volatile component accreted by the planet were dominated by icy planetesimals, the resulting atmosphere may become carbon rich relative to the canonical value.

The low abundance of ice (relative to the gas) inside the CO₂ iceline for the reset scenario is also interesting from the perspective of the overall core-envelope partitioning. If the core of the planet was to form from the solid material available at 1 AU in Fig. 3d, then the bulk of the planet would not be built up of volatile ices, but rather of more refractory components (rocks). The composition of the forming atmosphere will then be set solely by the composition of gas accreted onto the forming planet.

The CO₂-rich ice mantles formed beyond 16 AU in the reset scenario also have an interesting implication for the first stage of planet formation, the formation of pebble-sized objects. The sticking efficiency of 100 μm-sized CO₂ ice particles was recently determined to be an order of magnitude lower than that for similarly sized H₂O ice particles (Musiolik et al. 2016). Hence, under these particular conditions, the first steps of planet formation may be impeded.

4.3.2. Cometary composition

With molecular initial abundances (i.e., inheritance), the composition is generally preserved (with a few already mentioned exceptions, see Sect. 3). Hence, planetesimals forming under these specific conditions will be composed of “inherited” material. The overlap in the composition of comets, considered pristine remnants of the solar nebula, and interstellar ices certainly supports the hypothesis that a significant fraction of disk midplane material may be inherited from the molecular cloud (see, e.g., Mumma & Charnley 2011; Öberg et al. 2011a). For the reset case, cometary and interstellar abundances would be significantly different, at least for comets formed inside 30 AU.

It is interesting to consider whether a protoplanetary disk midplane, which inherits only water ice, is able to synthesis O₂ at a level similar to that seen in comets because this could be an observational constraint on the inheritance scenario. Our models can produce significant amounts of O₂ gas in the inner disk, as was also seen and discussed in full 2D protoplanetary disk models of Walsh et al. (2015). Bieler et al. (2015) find an O₂-to-H₂O ice ratio of ≈4% in comet 67P/C-G, with O₂ closely associated with the H₂O ice. They concluded O₂ was likely primordial in origin, i.e., originating from the parent molecular cloud and/or protostellar envelope. This result is surprising due to the high reactivity of O₂ ice. The confirmation of O₂ in comet 1P/Halley at a level similar to that in 67P/C-G suggests that O₂ is a dominant cometary component (Rubin et al. 2015). In Figs. 3a and 3c, molecular oxygen is produced only for the high ionisation case and reaches a peak abundance of <0.1% that of water ice. This reflects the efficient conversion of O₂ into water ice, once formed. The reset case with low ionisation produces more O₂ gas and ice with respect to H₂O ice, but it remains difficult to envisage a scenario in which H₂O and O₂ ice would be well mixed. Nevertheless, these results are sufficiently interesting to warrant a detailed study exploring a broader range of parameter space, which will be conducted in future work.

HCN has been detected in numerous comets with an abundance ratio of the order of a few × 0.1% relative to water (see, e.g., Mumma & Charnley 2011; Wirström et al. 2016). The results presented here show that HCN ice is unlikely to be produced in situ via the processing of ices containing only N₂ and NH₃, which is opposite to the case for O₂.

4.4. Caveats of model assumptions

Several assumptions have been made in our models, in addition to the set of initial abundances, the ionisation rate, and the chemical network used. We have considered only a single temperature and density profile, as well as single dust-grain size and density. Although the disk surface density used is typical of those derived for nearby protoplanetary disks (see, e.g., Williams & Cieza 2011), in reality, disk midplane temperature profiles depend on a multitude of parameters including, disk mass, stellar spectral type, dust distribution and opacity, and disk flaring index (which influences the amount of UV absorbed and/or scattered by the disk surface layer). The iceline locations determined here are consistent with those determined in protoplanetary disk models which are considered to be representative of the pre-solar nebula (see, e.g., Öberg et al. 2011b, and references therein). For disks around warmer stars, the icelines would move outwards in radius, with the opposite result for disks around colder stars. All else being equal, this will only affect the radial extent over which we see each chemical effect. A significantly warmer inner disk will also allow gas-phase reactions with reactions to occur more readily which may perturb the chemistry differently to that seen here. Also, our adopted disk structure is static in time. This was chosen to keep focus on chemical evolution for fixed physical conditions. Future work will utilise an physical disk model that is evolving in time.

There is observational evidence for both dust-grain settling and growth in protoplanetary disks (e.g., Dullemond & Dominik 2004, 2005; Williams & Cieza 2011). Settling towards the midplane will increase the dust-to-gas mass ratio, whereas dust growth to ≈mm sizes, will act in the opposite sense to decrease the total available surface area for freezeout and grain-surface reactions. The combination of both effects can either decrease or increase the total surface area of dust grains per unit volume. A smaller surface area (i.e., generally larger grains) will decrease the rate of freezeout (although freezeout is found to be almost instantaneous at the high densities of disk midplanes), and decrease the rate of grain-surface processing of ice mantles. Hence, this will lessen the effects of grain-surface chemistry, particularly for the case of high ionisation. Larger dust grains are also able to drift inwards as they become decoupled from the gas and affect the locations of icelines, by as much as 60% (as argued by Piso et al. 2015). The effects of different grain sizes will be investigated and addressed in future work.

Only two ionisation levels have been assumed here, to probe what are considered to be two extreme cases; however, there is also speculation in the literature that the pre-solar nebula may have been exposed to an increased cosmic-ray rate due to nearby supernovae early in the Sun’s lifetime (see, e.g., the review by Adams 2010). We have demonstrated here that a relatively high level of ionisation causes more chemical processing than a lower level. We would expect this trend to continue for ionisation levels ≫10^-17 s^-1. A lower level of ionisation would likely not have any effect, as the low level considered here is already seen not to significantly alter the chemistry. An increase or decrease in total disk mass (or surface density) will have the same effect as decreasing and increasing the ionisation rate due to CRs (especially in the outer disk), and increasing and decreasing the ionisation rate due to SLRs (especially in the inner disk).

The choice of initial abundances in the inheritance-scenario assumes an instant transition from molecular cloud phase to the disk midplane for the gas, with no chemical processing taking place “on the way”. Likewise, the reset-scenario assumes all gas content to be dissociated into atoms at the time it reaches the midplane. These two sets of abundances are extremes. Recent models indicate that chemical processing of the material does take place on its way to the midplane (Drozdovskaya et al. 2016), so a realistic initial composition of the material is more likely to be a mix of atoms and molecules. Considering the extreme cases, however, provides insight into the degree of chemical processing taking place.

Finally, we extracted our chemical abundances at a single time step only, at 10⁶ yr, which is considered to be representative of the lifetime of the pre-solar nebula and nearby young protoplanetary disks (e.g., Williams & Cieza 2011; Henning & Semenov 2013). A shorter time would lead to less extreme effects due to ion-molecule reactions and CR-induced photoreactions, because these reactions typically have long timescales (≳10⁵ yr; see Fig. 4); however, this strongly depends on the assumed cosmic-ray ionisation rate. A longer time would lead to the opposite case. The ALMA detection of gaps and rings in the disk around the young (<10⁶ yr) protostar, HL Tau, suggests that grain growth and planet(esimal) formation in protoplanetary disks may occur much earlier than heretofore considered (ALMA Partnership et al. 2015); hence, shorter chemical timescales are worthy of further investigation in future work. Specifically, in our models, assuming the inheritance-scenario and high level of ionisation, significant changes to the abundances of key volatiles take place after a few times 10⁵ yr (for more details, see Appendix A).

5. Conclusions

The models presented in this work have examined the importance of kinetic chemistry on the molecular composition (gas and ice) in protoplanetary disk midplanes. The main conclusions are listed below.

• The disk midplane composition reflects that of interstellar ices only for the case of low ionisation (SLRs only) in the inheritance scenario. The partitioning between gas and ice is determined solely by iceline positions, as is assumed in the planet population synthesis models. The inclusion of grain-surface chemistry has a negligible effect.

• Assuming a higher rate of ionisation (SLRs plus CRs) and inheritance leads to an increase in the abundance of CH₄ ice beyond its iceline, and a significant depletion of gas-phase CH₄ in the critical region between 1 and 15 AU. Cosmic-ray-induced chemistry enables the release of free carbon from CO in the outer disk (>15 AU) which is incorporated into gas-phase methane which freezes out. This naturally leads to low gas-phase CO abundances as is observed in some disks. On the other hand, cosmic-ray-induced chemistry efficiently destroys methane gas in the CH₄-poor region. The conclusion holds whether grain-surface chemistry is included or not.

• When grain-surface chemistry is considered, the CO₂ ice to H₂O ice ratio is increased, with CO and CH₄ gas destroyed at the expense of an increase in CO₂ ice. The critical reactions are the photodissociation of H₂O ice to form OH radicals within the ice mantle, which subsequently react with CO to form CO₂ ice. This reaction is able to proceed faster than H + OH recombination within ≈20 AU because the very volatile H atoms are quickly lost to the gas phase for grains warmer than 20 K. Beyond this radius, the reformation of water ice wins. The partitioning between N₂ and NH₃ is similarly affected.

• For the extreme reset scenario in which all elements are initially in atomic form, the picture changes significantly. Without grain-surface chemistry, gas-phase CO, O₂, and atomic oxygen are the main carbon and oxygen-bearing species beyond 1 AU. The chemistry does not have sufficient time to incorporate all available initial elemental oxygen into molecules by 10⁶ yr. Gas-phase water and CO₂ do form and subsequently freeze out, albeit achieving much lower abundances than in the inheritance scenario. A higher ionisation level helps to increase the production of H₂O and CH₄.

• With grain-surface chemistry, the abundances of H₂O and CO₂ ice increase significantly, demonstrating the absolute necessity of grain-surface chemistry for the synthesis of these two dominant ice components. The final abundance ratios reached in the very outer disk (30 AU) for H₂O and CO₂ are similar to those for the inheritance scenario, regardless of the ionisation level; however, the higher ionisation level does impede the abundance of O₂ ice at the levels seen in comets and enables a conversion from CO to CH₄.

• In the reset scenario, species other than the main considered volatiles are also produced in non-negligible quantities: O₂, HCN, and NO. The higher ionisation level generally helps the production of HCN and NO (at the expense of N₂ and NH₃) and impedes the survival of O₂ ice.

• The inclusion of chemistry has a significant impact on the C/O ratio of both gas and ice in the planet-forming region which is expected to influence the resulting composition of forming planet(esimal)s. The ices remain, on the whole, dominated by oxygen (i.e., C/O < 1). For both inheritance cases, the gas is carbon rich relative to the canonical value; however only for the low ionisation case is there a reservoir of gas-phase material with C/O > 1. For both reset scenarios, the ice becomes more carbon rich than the gas, which is opposite to the inheritance case.

The results presented here show that under certain conditions, highlighted above, chemistry can have a profound effect on the composition of the planet-forming material in disk midplanes. Chemistry influences the partitioning of elemental carbon, oxygen, and nitrogen, into molecules of differing volatilities, such that the positions of ice lines alone, are not necessarily adequate for determining the ratio of C/O in neither the gas, nor the ice. Only under the extreme case of full inheritance and low ionisation, are the elemental ratios determined solely by the positions of icelines. This conclusion is also similar to that for the assumption that ices are already locked up in larger bodies by ≈10⁵ yr.

The work presented here follows the time evolution of chemistry in a static protoplanetary disk which is the simplest physical case. In reality, disk conditions evolve with time, at the same time as planetesimals are forming and migrating within the midplane. Future plans include determining the influence of chemistry in an evolving protoplanetary disk (where the density, temperature, and ionisation rate also vary with time), and to couple the outputs of these models with planet formation tracks to determine, in a quantitative manner, the influence on the resulting composition of gas-giant planetary atmospheres.

---

Acknowledgments

The authors thank Amaury Thiabaud and Ulysse Marboeuf for making their disk density and temperature profile available and for many inspiring discussions. The authors also thank Karin Öberg for a fruitful discussion about CO gas depletion inside the CO iceline. Astrochemistry in Leiden is supported by the European Union A-ERC grant 291141 CHEMPLAN, by the Netherlands Research School for Astronomy (NOVA), and by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize. CW also acknowledges the Netherlands Organisation for Scientific Research (NWO, grant 639.041.335).

References

Appendix A: Timescales

Figure A.1 shows the abundances of the four key volatile, taken at four different evolutionary times. Note the different y-axes ranges. It is seen that significant changes occur between 100 and 500 kyr, that is after a few times 10⁵ yr. This goes for all the four species, gas as well as ice, out to about 25 AU.

It is interesting to see that H₂O and CO₂ are respectively produced and destroyed in the gas phase in the inner disk, and vice versa in the outer disk. For CO and CH₄ there are turn-over radii in the inner disk, inside of which they are both being produced, and outside of which they are both being destroyed. These turn-over radii correspond to the iceline positions of H₂O and CO₂ for the case of CH₄ and CO respectively. This nicely illustrates the dependance on gas phase H₂O and CO₂ for the production/destruction of CO and CH₄ in the gas phase.

All Tables

All Figures