Free Access
Volume 613, May 2018
Article Number A14
Number of page(s) 16
Section Planets and planetary systems
Published online 21 May 2018

© ESO 2018

1 Introduction

Over the past 20 years exoplanetary science has evolved from purely detection (e.g. Udry & Santos 2007; Borucki et al. 2011; Batalha et al. 2013; Fischer et al. 2014) to exoplanet characterisation. One example of characterisation is the deciphering of the molecular content of exoplanet atmospheres which has revealed simple molecules such as CO and H2O, and possibly CO2 and CH4 (e.g. Seager & Deming 2010; Snellen et al. 2010; Birkby et al. 2013; Fraine et al. 2014; Crossfield 2015; Sing et al. 2016).

The C/O ratio of an exoplanet atmosphere (usually derived from relative abundances of the simple molecules listed above) has been proposed to be a possible tool to link an exoplanet to its formation site in the natal protoplanetary disk. This is because disk midplanes have a radial gradient of C/O ratios in gas and ice (see e.g. Öberg et al. 2011). This approach, to couple observed C/O ratios in exoplanets with their possible formation sites, typically assumes that the radial dependence of the midplane C/O ratio is defined solely by the iceline (snowline) positions of the main volatile species in a disk of fixed chemical composition (see e.g. Marboeuf et al. 2014). The position of the iceline for a particular species depends on the temperature structure of the disk midplane, and the binding energy (or desorption energy) of the molecule to ice mantles on dust grains. Icelines thus affect the relative ratios of heavy elements in the gas and ice. For example, beyond the H2O iceline, the gas is depleted in oxygen whilst the ice is enriched in oxygen. Similar changes occur across the CO2, CH4, and CO icelines thereby also altering the local C/O ratio in both the ice and gas. An example of a variable C/O ratio for a disk with a fixed chemical composition (as per the conventional assumption) is shown in Fig. 1 of Öberg et al. (2011). Due to the step-like features of that figure, caused by icelines, the figure will be referred to as the traditional stepfunction.

Protoplanetary disks are not static objects: over the disk lifetime, from formation to dispersal, typically ~0.1–10 Myr, the disk generally cools, spreads, loses material (via accretion onto the star or dissipation) and thus becomes less dense over time (see the review by Williams & Cieza 2011). Hence, the icelines of volatile species also change in time, generally migrating inwards towards the star (Davis 2005; Min et al. 2011; Harsono et al. 2015). To date, most models exploring the location and evolution of icelines with the explicit goal to link to the composition of planet-building material have assumed static midplane conditions (Öberg et al. 2011; Madhusudhan et al. 2014; Thiabaud et al. 2015b). Static models are more straightforward to couple with chemical kinetics. Recently Piso et al. (2015, 2016) explored the evolution of elemental ratios across icelines in both a static and viscously evolving disk model including the radial drift of ice-coated grains of different sizes. The temperature structure is kept fixed in time in these models. A similar approach has been adopted by Ali-Dib et al. (2014) and Ali-Dib (2017). Both of these models include dustand pebble accretion, but do not include kinetic chemistry. Those models which do allow evolving disk midplane conditions also neglect kinetic chemistry for simplicity (see, e.g. Marboeuf et al. 2014; Mordasini et al. 2016).

The impact of kinetic chemistry in planet formation models is yet to be fully explored due to the increased computational demands in coupling a chemical model with, for example, a planet population synthesis model. This is despite chemical kinetics being a necessity in models explaining protoplanetary disk observations for some 20 years since the first detection of molecules, other than CO, in disk atmospheres (see review by Henning & Semenov 2013). This is now changing, as is the traditional picture of a static disk. Recent examples of kinetic gas-grain chemistry in disks can be found in Helling et al. (2014), Walsh et al. (2015) and Cridland et al. (2016), with different levels of treatment of the chemistry and grain growth. Also, there is growing observational evidence for low gaseous CO abundances which can be attributed to either rapid planetesimal growth or chemical evolution (see e.g. Favre et al. 2013; Bruderer et al. 2012; Reboussin et al. 2015; Kama et al. 2016; Miotello et al. 2017).

The disk midplanes suffer from a lack of observational constraints on their chemical compositions. This due to the observational difficulty of probing the midplane through the optically thick higher-lying layers of the disk and to the freeze-out of species in the outer disk midplane. Most molecular emission, whether at millimeter or infrared wavelengths, therefore arises from the warm layer at intermediate disk heights (see Bergin et al. 2007, for a review). Even the N2H+ emission observed with ALMA, which probes the part of the outer disk where CO is frozen out, arises mostly from layers just above the midplane (Qi et al. 2013; van ’t Hoff et al. 2017). Molecules expected to be abundant in midplane ices such as H2O (Hogerheijde et al. 2011; Du et al. 2017), NH3 (Salinas et al. 2016), H2CO (Loomis et al. 2015) and CH3OH (Walsh et al. 2016) have been detected in the gas-phase but the observations either lack spatial resolution or reflect subsequent gas-phase chemistry, making it difficult to infer the outer disk midplane ice abundances. Ultimately, comet abundances may provide the best constraints on this region, at least for our solar nebula disk.

In contrast, the midplane of the warm inner disk where CO is not frozen out has recently been imaged for the first time through CO isotopologue observations of 13C18O by Zhang et al. (2017) for one disk. Such observations are promising, and are perhaps also feasible for other common species, but they require large telescope time investments. Also, the dust continuum emission becomes optically thick in the inner disk, hiding our view of much of the planet-forming midplane. Thus, perhaps ultimately the only information on the midplane chemical composition is to be inferred from the atmospheric compositions of the planets that form out of the midplane, and which remain after the disk is dispersed.

As emphasized in Mordasini et al. (2016), the road from a gaseous disk midplane to an observable exoplanet atmospheric composition is long, requiring many steps, physical (e.g. migration) as well as chemical. A key issue is whether the planets’ heavy element enrichment is controlled by the accreted gas or by the accreted (icy) planetesimals. Mousis et al. (2009), Thiabaud et al. (2015a), and Mordasini et al. (2016) found that the dominant delivery mechanism for heavy elements to a forming gas-giant planet atmosphere is accretion of icy planetesimals rather than direct accretion from the gas, at least for planets with masses less than a few times that of Jupiter. The ultimate C/O ratio in the observable part of the planet’s atmosphere also depends on how well the sublimated planetesimals are mixed in the envelope and on the core-envelope partitioning. In any case, it is clear from models such as those of Mordasini et al. (2016), which explicitly include planetary formation and planetesimal accretion processes, that the compositionof both the gas and the ice of the evolving parent disk is needed to determine the planetary atmosphere composition.

Eistrup et al. (2016, hereafter Paper I) investigated how kinetic chemical evolution affects the composition of volatiles (ice and gas) in a static protoplanetary disk midplane. If interstellar ices are inherited from the collapsing cloud (without alteration) by the protoplanetary disk, then the ice composition is preserved only for the case where there is no source of external ionisation (i.e. no cosmic rays). An important conclusion from Paper I is that if chemistry is efficient on a timescale of ~ 1 Myr, then the C/O ratio of the gas and ice remains less than one always within the CO iceline. Only for the case of low ionisation and full inheritance does the C/O ratio approach 1, and this only occurs for the gas: the ice remains oxygen-rich always.

The effects of chemical kinetics on the composition of gas and ice volatiles in evolving protoplanetary disk midplanes remains to be quantified. We here expand upon the work presented in Paper I to investigate how changing temperature and density in the disk midplane with time modifies the above conclusions. We also considered how disk mass may affect the composition of planet-building material. We furthermore studied the timescale for chemical changes by extending the lifetime of our disk models to ~ 7 Myr at which point we assumed that planet formation has halted and disk dispersal is advanced.

We investigate a wide parameter space, with somewhat simplified physical structures, and two sets of extreme initial abundances. The two scenarios of initial chemical abundances are: cloud inheritance (hereafter “Inheritance” scenario) and chemical reset (hereafter “Reset” scenario). In the inheritance scenario the initial abundances are molecular, and taken to be those found in ices in the parent molecular cloud. In the reset scenario, the initial abundances are atomic, because all molecules are assumed dissociated by high temperatures close to the young star. Both scenarios include H, H2 and He initially, with abundance ratios: N(H)/N(H2) = 10−4, and N(He)/N(H2) = 0.17. Generally, the inner disk midplane should be susceptible to luminosity outbursts that can reset the chemistry, while the outer disk midplane should be expected to inherit the parent cloud abundances, so both scenarios may apply but in different parts of the disk (Pontoppidan et al. 2014). The intention is to investigate chemical trends expected for these two sets of initial conditions and using a more realistic disk, which is becoming cooler and more diffuse over time.

2 Methods

The evolving physical disk model from Alibert et al. (2013) was used. This disk model was adopted in planet population synthesis studies. It features time evolution of disk midplane temperature and pressure as a function of radius (0.2–30 AU) over a timescale of 7 Myr. The approach for calculating the ionisation rate and chemical evolution follows that presented in Paper I: we summarise the main aspects here.

2.1 Evolving disk model

An evolving disk here refers to a disk that is cooling in time whilst also losing mass. In this work power laws were fitted to the temperature structure evolving in time as originally computed by Alibert et al. (2013).

The disk structure adopted from Alibert et al. (2013) and used in Paper I will be referred to as the 0.1 MMSN disk, having an initial mass of 1.3 × 10−3 M or 0.13 times the minimum mass solar nebula, or MMSN (Weidenschilling 1977). The reason for fitting power laws was that the initial time-dependent temperature structure from Alibert et al. (2013) featured a step, that is, an unphysical radial drop, at around 5 AU, as well as an artificial lower temperature cut-off at 20 K. Both of these features were smoothed out with a power law profile in Paper I. The same unphysical drop in temperature was observed for later time steps. In order to smooth out all timesteps, whilst ensuring a temperature profile that is monotonically decreasing in time for all radii, power laws were fitted to all temperature profiles from Alibert et al. (2013). For timesteps up to 400 kyr of evolution, the power laws were fitted using model data from 1 to 2.1 AU, and for all later timesteps from 0.2 to 0.6 AU. The resulting power law fits were then extrapolated to the entire radial range of the disk. This procedure gave the temperature profiles that can be seen for selected timesteps (as labelled) in Fig. 1a.

Using these temperature profiles, the midplane density profiles for all timesteps were computed using the midplane pressure profiles from Alibert et al. (2013), assuming H2 to be the main gas constituent. The resulting midplane number density radial profiles at selected timesteps are shown in Fig. 1b.

The surface density profiles at the first timestep were computed using the prescription and values for Table 1, Disk “8” in Alibert et al. (2013). Taking into account disk evolution, the surface density at later timesteps at a radial distance R was computed by scaling the initial surface density Σt=0(R) to the evolution in midplane number density n(R) as seen in Fig. 1b. Hence, the surface density at point R at time t is given by (1)

To drive chemical evolution, ionisation of the disk midplane is needed, a conclusion reached by our and other models (see e.g. Helling et al. 2014; Aikawa et al. 1996; Willacy et al. 1998; Walsh et al. 2012; Cleeves et al. 2013b). Paper I found that chemical evolution was observed to be significant within 1 Myr evolution time only when short-lived radionuclides (SLRs) and cosmic rays (CRs) were included. In this paper, the contribution to ionization by SLRs and CRs were treated in a similar way to that in Paper I: the ionisation rate from decay of SLRs is computed using the prescription from Eq. (30) in Cleeves et al. (2014b), including the time-dependent term, which was omitted in Paper I. (2)

The contribution from CRs were calculated assuming attenuation depending on how much material the CRs have to pass through on their way to the midplane. The attenuation factor was calculated using the (now) evolving surface density of the disk in the vertical direction giving the following ionisation rate: (3)

where Σ(R, t) is the evolving surface density profile, taken from Eq. (1), and ζ0 = 10−17 s−1 (as in Paper I).

The SLR ionisation rate profiles (“low” ionisation) and the SLR + CR ionisation contribution profiles (“high” ionisation) are plotted at selected timesteps in Fig. 1c and d. We note that the SLR ionisation level decreases rapidly with time and becomes very low, of order 10−22 s−1 by ~ 7 Myr. In contrast, the CR ionisation level increases with time, especially in the inner disk, due to the decreasing surface density.

The timestep grid of Alibert et al. (2013) was adopted, with physical conditions changing for every timestep of their models. However, timesteps were added at early evolutionary stages (1 104 yr) in order to account for short-term chemical effects. The first timestep in the time grid in this work is 1 yr, and the physical conditions at 1 yr fromAlibert et al. (2013) were adopted for this timestep. Four logarithmic decades were applied ranging from 1 to 104 yr, with physical conditions during this time fixed at the conditions at 1 yr. After 104 yr, the time grid and the associated, continuously changing physical conditions of Alibert et al. (2013) were adopted. The time goes up to 7 Myr, in order to also investigate long-term chemical evolution. Until 7 Myr, this disk model goes through 401 timesteps.

2.2 Higher disk mass: 0.55 MMSN disk

For investigating the effect of a higher disk mass, a more massive evolving disk physical structure from Alibert et al. (2013) was adopted (Table 1, “Disk 12” of Alibert et al. 2013). The initial mass of this disk is 0.55 MMSN = 0.0055 M, determined in the same manner as was the case for the less massive disk in Paper I. Since this disk has the same steep temperaturedrop as was the case with the 0.1 MMSN disk, and also features the 20 K lower temperature limit, power laws were also fitted to the temperature profiles at each timestep. In order to have the temperature profiles decreasing in time monotonically at each radial point, for this disk, the power law fitting was done for the first timestep (1 yr) and the last timestep (7 Myr) for the evolvingdisk temperature profiles between 1.1 and 1.5 AU. For all intermediate timesteps, power law indices intermediate to the two fitted indices were adopted, so that the power law index is continuously and linearly decreasing in time, with the step-wise decrease in index value being constant for each timestep. The original disk temperature structure from Alibert et al. (2013) was then adopted inside 1.1 AU, whereas the fitted (for first and last timesteps) and imposed power law structures were applied outside 1.1 AU. For the outer parts of the disk, at late timesteps, the temperaturegets very low (below 8 K), which is unphysical. A lower limit of 8 K to the temperature was therefore imposed.

The midplane number density and ionization level structures for this 0.55 MMSN disk were computed in the same way as described for the 0.1 MMSN disk in Sect. 2.1. The resulting midplane temperature, density and ionisation profiles for the midplane of the 0.55 MMSN disk are presented in Fig. A.1. The evolving 0.55 MMSN disk goes through 438 time steps during the 7 Myr modelled evolution. The first timestep is at 1 yr, and from 1 yr through 104 yr the time resolution is the same as for the same early timesteps for the 0.1 MMSN disk model, as described in Sect. 2.1. After 104 yr, the timesteps from Alibert et al. (2013) were adopted.

The disk structures differ for the two different disk masses. Since the ionisation level is dominated by the contribution from CRs at all times throughout the disk (see Figs. 1 and A.1, panels c and d), the decreasing surface density means the degree of ionisation in the inner disk is increasing in time. Because of the higher densities, the degree of ionisation in the inner disk is increasing faster in time for the 0.55 MMSN disk than for the 0.1 MMSN disk.

The temperature structure for the 0.55 MMSN disk is hotter in the inner disk, and has a steeper radial gradient, than that for the 0.1 MMSN disk. The low temperature in the outer disk for the 0.55 MMSN disk is due to lesser heating from irradiation from the central star. This is because more material shields the radiation for the 0.55 MMSN disk than for the 0.1 MMSN disk.

Figure 2 shows the evolution of the disk masses for the two disks, normalised to their initial masses. Both disks experience significant mass loss during the evolution, with the 0.1 MMSN disk losing 95% of its initial mass, and the 0.55 MMSN disk losing 98% of its initial mass, after 7 Myr.

thumbnail Fig. 1

Midplane temperature and density (top) and ionisation rates (bottom) as functions of disk radius for different evolutionary times for the 0.1 MMSN evolving disk model. For this model, the starting conditions at 1 yr are given by the orange profiles shown here, which feature slightly different physical conditions than that for the static disk structure adopted by Paper I (see text for details).

2.3 Chemical model

The chemical model included gas-phase chemistry, gas-grain interactions and grain-surface chemistry. The gas-phase chemistry was from the latest release of the UMIST Database for Astrochemistry (McElroy et al. 2013) termed RATE12. The rates for gas-grain interactions and grain-surface chemistry were calculated as described in Walsh et al. (2015, and references therein). A gas-to-dust mass ratio of 100 was adopted. The chemical model used here was identical to the “Full chemistry” setup in Paper I. Photodissociation and photodesorption due to cosmic-ray induced UV photons were included.

The adopted grain size is 0.1 μm. As discussed in Paper I, this assumption is certainly not correct, but the effects of grain growth to larger sizes and decreasing midplane gas-to-dust ratios due to settling tend to compensate each other. The important dust parameter for the chemistry is the amount of surface area available for freezeout and subsequent grain-surface chemistry. This is often parameterised as a combinationof grain size and number density (set by the assumed gas-to-dust mass ratio). The assumption on the constant grain size will be further discussed in Sect. 3.6.3. More generally, all trends found in this paper are robust but the exact timescale for chemical changes and the magnitude of the chemical changes will depend on the assumed grain size and the gas-to-dust mass ratio.

For details about the two sets of initial abundances (the inheritance and reset scenarios), see Table 1 in Paper I. The overall C/O ratio utilised is 0.34. In Paper I it was found that if the disk was exposed to ionising cosmic rays, efficient production of CO2 ice and O2 gas and ice takes place at intermediatedisk radii, whereas CH4 gas was converted to CO and CO2. Also, if full chemical reset into atoms occurs en route into the disk midplane (e.g. via a shock or a luminosity outburst), then the abundances of dominant volatiles can differ from interstellar values by more than an order of magnitude (see also Drozdovskaya et al. 2016). Most striking was the depletion of H2O ice advancing enhanced gas-phase CO and O2. This difference was larger for the case where cosmic rays were excluded from the reset scenario. Helling et al. (2014) also found cosmic rays to be of importance to chemistry that alters the C/O ratio in the midplane. Their results, however, were different from those in Paper I, due to differences in the grain-surface chemistry utilised.

Throughout this work, whenever distinguishing between gas-phase and icy species is necessary, an “i” in front of a molecule involved in a reaction will denote that the molecule is icy, and a “g” in front will denote a gas-phase molecule. For reference, an iceline is the midplane radius inside of which a given chemical species is in the gas phase, and outside of which the species is in the ice.

thumbnail Fig. 2

Evolving disk masses for the two disks, normalised to the initial disk masses of 0.1 and 0.55 MMSN. For both disks, more than half of the initial mass is lost before 0.5 Myr, and by 1–1.5 Myr evolution, only 20% of the initial masses remain.

3 Results

Figure 3 shows abundances as a function of radius for selected volatile species in gas and ice after 7 Myr evolution, for four different model setups. As indicated in the boxes in each panel, the model setups explore a range of static and evolving disks, for the inheritance versus reset scenarios, and at both low and high ionisation levels.

3.1 Static versus evolving disk: shifting icelines

In Paper I it was found that given low ionisation and inherited initial chemical abundances, chemical evolution was insignificant on the timescale of 1 Myr. Figure 3a shows the same case, but now for a physically evolving disk structure (the 0.1 MMSN disk), as described in Sect. 2, and for up to 7 Myr. It is evident from Fig. 3a that the chemistry inthe midplane remains defined by the iceline positions, and that chemical processing has caused no significant chemical changes. This conclusion holds for the 0.55 MMSN evolving disk structure as well. Hence, the conclusion from Paper I, that in the case of inherited abundances and low ionisation level, chemical evolution is insignificant, can be extended to a timescale of 7 Myr for these cases.

The main effect of utilising an evolving disk structure instead of a static one is that the icelines shift inwards with time. This in evident when comparing Figs. 3c and d, which both feature the inheritance scenario at high ionisation, but for a static and an evolving disk structure, respectively. For Fig. 3d the icelines are shifted inwards in time, whereas for 3c the icelines are static corresponding to the static temperature structure. The additional effects of adopting an evolving disk structure instead of a static structure are the retention of H2O ice around 0.7 AU in panel d and one to two orders of magnitude lower CH4 gas abundance in the inner disk. The iceline positions for the disks are determined from Fig. 3 as the maximum distance R from the star at which a volatile is more abundant in gas than in ice. Figure 4 features the locations of the thermal icelines of selected volatile species as functions of time. The solid profiles indicate the movement of the icelines over time. We note that the axes scales in Fig. 4 are linear.

For H2O, NH3, and CO2 in the 0.1 MMSN disk, the icelines shift by 30–35% between 0.5 and 7 Myr of evolution (e.g. from 0.7 to 0.5 AU for H2O). For CH4, CO, and N2 the icelines are shifted inwards by ~60%. The larger shifts for the latter, more volatile species are due to the larger decreases in temperature in the outer disk, due to the decreasing slopes of the temperature profiles with time (see Fig. 1a). For the 0.55 MMSN disk, the icelines shift farther than for the 0.1 MMSN disk. This is due to the steeper temperature profile for the 0.55 MMSN disk, as can be seen when comparing top left panels in Figs. 1 and A.1. The temperature drops by ~50% from 0 to 2 Myr at 10 AU in the 0.55 MMSN disk, versus a temperature drop by ~20% in the 0.1 MMSN disk over the same timeframe and radius.

3.2 Timescales of chemical changes

Figures 5 and 6 show abundances for five key volatiles as function of midplane radius at six timesteps: 0.1, 0.5, 1, 2, 5 and 7 Myr, for the inheritance scenario model at high ionisation. The left panels are for the 0.1 MMSN disk, and the right panels are for the 0.55 MMSN disk. A time of 0.1 Myr was chosen as the first timestep because it was found in Paper I thatfor the inheritance scenario no chemical changes had taken place yet at this evolution time. An age of 0.5 Myr is representative of the timescale corresponding to the end of the embedded phase of star formation, when the protoplanetary disk is revealed. This was also found in Paper I to be the chemical evolution timescale for the high ionisation case. In that work, 1 Myr was considered the full evolution time: here, we consider additional and later timesteps at 2, 5, and 7 Myr. The first is chosen as it is the median age of T Tauri stars and represents the lifetime of the gas-rich protoplanetary disk. A time of 5 Myr represents an aged disk like TW Hya. By 7 Myr the gaseous disk is thought to have fully dissipated resulting in a gas-poor debris disk (see e.g. Fedele et al. 2010, and Fig. 2).

As seen in Figs. 5 and 6, chemistry is continuously ongoing throughout the evolutionary time of 7 Myr. The chemical abundance changes over 7 Myr vary from factors of a few (e.g. a 3 times decrease in H2O ice abundance at 10 AU in Fig. 5a to orders of magnitude (e.g. a CH4 gas depletion at 5 AU in Fig. 6a, resulting in a significantly different chemical “picture” from the beginning to the end of the evolution. For all species considered in these plots, the shifting icelines can be seen as an inward shift with time of the radius where the gas and ice profiles cross each other.

Overall, with the exceptions of H2O, CO and CO2 inside 0.4 AU to be discussed in Sect. 3.3, the trends are similar for the 0.1 MMSN and 0.55 MMSN disks. Thus, in the following discussion only the chemical evolution of the 0.1 MMSN disk case is considered.

To visualise the main chemical changes more clearly, Fig. 7 highlights the evolution of key species with time at each of four chosen radial positions: 1, 5, 10 and 30 AU. At 1 AU large changes happen between 1 and 7 Myr. The abundances of CO and O2 gas seem to converge at the end of the evolution, indicating that a steady state situation may be reached where these two species are the main carbon and oxygen reservoirs. H2O ice is still dropping in abundance, thus, if evolution were to continue for a longer time, a further decrease in H2O ice abundance would likely be seen.

At 5 AU, no steady state seems to have been reached after 7 Myr. CO and O2 gas and H2O and CO2 ice have near equal abundances to within a factor of two, and a steady increase in C2 H6 ice abundance is seen throughout the evolution reaching a final abundance of ~ 10−5, making it the third most abundant carrier of elemental carbon, accounting for ~ 15% of the carbon. CH4 gas is depleted due to reactions with He+, which is a consequence of the high ionisation level, as also mentioned in Paper I. A “reverse” evolution of CO gas and CO2 ice is also evident, taking place at about the same time (between ~1.6–5 Myr), indicating interdependency between production and destruction of CO and CO2. The decreasing abundance of CO gas and increasing abundance of CO2 ice from ~0.1 Myr to ~2 Myr can be explained by CO gas colliding with the grains, and subsequently reacting fast with OH (faster than CO can desorp) to produce CO2 on grain surfaces, as also described in Paper I. The reverse evolution for CO and CO2, after 2 to 3 Myr, is due to the increasing midplane ionisation level over time (see Fig. 1d). For both a static and an evolving disk, the continuous CR ionisation of the midplane causes an increasing CR-induced photon fluence there (we note that the fluence increases faster for the evolving disk due to the evolving ionisation level). This increasing fluence amplifies the photodissociation of CO2 ice into CO and O on the grains, which starts to dominate over the CO2 production reaction iCO + iOH → iCO2 + iH. The result is a decreasein CO2 ice and an increase in CO gas after 2–3 Myr. Thus, for an evolving disk, and the physical conditions considered here, the CO2 ice abundance reaches a maximum. The more diffuse the midplane becomes over time, the higher the ionisation contribution from CRs stays accordingly, and the shorter time it takes for theCO and CO2 evolutions to invert, in what will now be referred to as an “abundance inversion”.

Going out to 10 AU, the evolution is very similar to that at 5 AU. The main difference is the rapid decrease in O2 and CO gasabundances at the end of the evolution, which is due to the shifting of the O2 iceline from outside to inside of 10 AU (see Fig. 4), and the destruction of CO gas. The latter is caused by the rapid conversion of CO into CO2 ice mentioned above. These changes to O2 and CO gas at 10 AU remove them from the gas phase. An effect of this freeze-out is also the extra availability of O2 and CO on the grain surfaces which both react to form H2O and CO2 ices, which are seen to increase in abundances simultaneously with the freeze-out of CO and O2.

Out at 30 AU it is cold enough for all volatiles, except CO, to be frozen out initially. As seen in Fig. 7, within the abundance limits, CO gas and ice are co-existing because R = 30 AU is just inside the initial thermal CO iceline, with the ice abundance accounting for a few percent of the total CO abundance. The CO gas abundance decreases after ~0.4 Myr, at which point the midplane has cooled enough for the CO iceline to move inwards of 30 AU, after which the CO ice becomes more abundant than the gas. However, the CO ice abundance after about 1 Myr is only a few percent that of the initial CO gas abundance, indicating that CO is processed on grain surfaces upon freeze-out. CO2 ice is also destroyed after 1 Myr of evolution, but CH4 ice is steadily produced to become the dominant carbon-bearing species after 7 Myr.

This production of CH4 ice at 30 AU is linked to the destruction of CH3OH ice through photodissociation by CR-induced photons which takes place simultaneously after 1 Myr. CH3OH ice is destroyed on the grains inside and outside of the location of CO iceline. It is not reproduced from hydrogenation of CO because CO is being converted into CO2. Outside the CO iceline CO and CH3OH are both being photodissociated on the grains caused by CR-induced photons, through the reactions,

This means low abundances of oxygen-containing complex organic molecules. In the cold outer midplane with temperatures down to around 15 K, the mobility of molecules on the grains is low, so hydrogenation reactions become dominant. iCH3 + iH → iCH4 is causing the increase in CH4 ice abundance, with CH4 being the dominant carbon carrier after 7 Myr of evolution at 30 AU. Likewise, the hydrogenation of C2 H4 ice and C2 H5 ice is responsible for the production of C2 H6 ice. H2O ice is formed from hydrogenation of OH leftover from the photodissociation of CH3OH, including the reaction:

The abundance of C2 H6 ice increases with increasing radial distance. In addition to C2 H6 ice, HCN ice is also abundantly produced at 30 AU, and these two species jointly account for 40% of the elemental carbon, with CH4 accounting for an additional 50%, and the remaining 10% is in other molecules, such as CH3OH (further discussed in Sect. 4.2) by 7 Myr. H2O ice accounts for more than 99% of the elemental oxygen.

Overall, it is evident from these figures and discussion that chemical evolution continues up to (at least) 7 Myr. That was also a conclusion of Kamp et al. (2013). Under the conditions and assumptions presented here, the majority of changes to the chemical abundances take place between 0.5 and 7 Myr, and as in evident from Figs. 5 and 6, the changes are significant (from factors of a few to orders of magnitude changes for all species).

thumbnail Fig. 3

Abundances after 7 Myr evolution for four different model setups. Panel a: inheritance scenario with low ionisation and evolving disk structure. Panel b: reset scenario, high ionisation and evolving disk structure. Panel c: inheritance scenario with high ionisation and static (unchanging) disk structure. Panel d: inheritance scenario with high ionisation and evolving disk structure. The blue, red-and-black and green arrows to the right of each panel indicate the initial abundances assumed for H2O, CO/CO2 and CH4, respectively, in the inheritance scenario.

thumbnail Fig. 4

Evolving iceline positions for main volatiles, for 0.1 MMSN (left) and 0.55 MMSN (right) evolving disks.

thumbnail Fig. 5

Abundances of H2O, CO and CO2 gas (dashed) and ice (solid) as a function of midplane radius at five timesteps: 0.5, 1, 2, 5, and 7 Myr. The left-hand and right-hand panels show the 0.1 MMSN and 0.55 MMSN disk results, respectively.

thumbnail Fig. 6

Same as Fig 5, but for CH4 and O2.

thumbnail Fig. 7

Time evolution of selected volatiles, in gas as ice at 1, 5, 10 and 30 AU, for the 0.1 MMSN disk with high ionisation and inheritance.

3.3 Varying disk masses

Figures 5 and 6 present radial abundance distributions at different evolutionary timesteps for CO, CO2, CH4, H2O and O2. Left panels are for the low mass evolving disk structure. Right panels are for the high mass evolving disk structure.

Overall, when comparing the two disk masses for all five species it is seen that there are no significant differences in the abundance evolutions. However, the temperature structures are different, thus the icelines are located and shifted by different radial distances (see Sect. 3.1). Only in the case for CO2 gas, is the evolution in the innermost disk (inside 0.4 AU) very different (by up to 3 orders of magnitude). This is due to the high temperatures and densities here, causing the abundance to be dictated by chemical equilibrium. This is more pronounced for the 0.55 MMSN disk, which is warmer in the inner tenth of an AU, than the 0.1 MMSN disk.

3.4 Importance of initial abundances: inheritance vs reset

The choice of initial abundances was shown in Paper I to have a large effect on the final abundances after 1 Myr of evolution. In Fig. 3, panels b and d show abundances after 7 Myr evolution for an evolving disk with a high level of ionisation, with panel b assuming a reset scenario for the initial chemical abundance, and panel d assuming the initial abundances to be inherited. Since both cases have assumed the same evolving physical structure, the icelines are identically located in both panels.

Inside the H2O iceline at ~0.5 AU, where all species are in the gas phase, similar abundance levels are seen for H2O, CO and CO2, in both scenarios. CH4 gas, on the contrary, is absent from the inner disk in the reset scenario. In relation to this absence, a destruction of CH4 gas is evident in the inheritance scenario as seen at 1 AU in Fig. 7a, which indicates that CH4 is being continuously destroyed in the gas phase, a destruction happening more slowly inside the H2O iceline in Fig. 3, than it does outside that iceline.

Redirecting focus from the inner to the outer disk, the abundances of the ices outside the O2 iceline at around 10 AU are very similar (at 30 AU CO2 is a factor of two different, CH4 is ~10% different, and H2O is only 1% different) when comparing the reset and the inheritance scenarios in Figs. 3b and d for chemical evolution by 7 Myr. Indeed, outside 10 AU the choice of initial abundances seems irrelevant for the resulting abundances, given a long enough evolution time. In other words: in this region, which is very relevant for comet formation, chemical steady state is close, and the fingerprint of the interstellar chemical abundances is erased if the gas disk lasts for 7 Myr.

The inner (mostly) gaseous midplane, and the outer (mostly) icy midplane produce similar abundance ratios independent of the choice of initial abundances. The intermediate region in Figs. 3b and d, roughly between the H2O and O2 icelines at 0.5 AU and 10 AU, is quite different. The dominant gas phase species are CO and O2, both at approximately the same abundance level within a factor of 2, in both the reset and the inheritance scenarios. H2O and CO2 ices, however, are one-to-two orders of magnitude lower for the reset scenario than for the inheritance scenario. Looking into the time evolution of CO, O2, H2O and CO2 in the inheritance scenario with high ionisation at 5 AU shown in Fig. 7b, it is evident that none of these species have yet reached steady state (although it is close) by 7 Myr. Additional evolution time would therefore likely bring the abundances of H2O and CO2 ices in the inheritance scenario further down, and CO and O2 gases further up, similar to the plot for the reset scenario in Fig. 3b.

3.5 O2 as a significant disk midplane molecule

Molecular oxygen is often not considered to be a dominant molecule in planet formation models. The detection of a high abundance of O2 ice to H2O ice of 3.8 ± 0.85% on Comet 67P reported by Bieler et al. (2015) has triggered models to explain the origin of this molecular oxygen (see Mousis et al. 2016; Taquet et al. 2016).

In disk chemistry models by Walsh et al. (2015) O2 is abundantly produced in the intermediate layers of a disk. Here, the O2 evolution in the midplane is traced up to 7 Myr. Figure 6 shows abundance profiles for O2 at different timesteps for the inheritance scenario. It is evident that between the H2O and O2 icelines, the abundance of O2 gas continuously increases with time. O2 is produced on the grain surfaces through the following reaction pathway:

O2 ice subsequently desorbs (if inside the O2 iceline), and steadily builds up a high abundance of ~10−4 wrt H2O over 7 Myr evolution. The initial abundance of H2S is 6 × 10−6 (20 times lower than the O2 gas abundance by 7 Myr), and the HS molecule left over from the last step in the reaction chain above is recycled for conversion of more oxygen atoms into O2 via the second step above. Henceforth, the initial abundance of sulphur is seemingly irrelevant for the final abundance of O2, since the sulphur merely acts as a catalyst for the conversion of oxygen atoms into molecular oxygen.

3.6 Caveats of model assumptions

3.6.1 Choice of grain-surface reaction parameters

An important part of the chemical evolution presented here is caused by grain-surface chemistry. The importance of grain-surface chemistry was also highlighted in Paper I, as well as in Walsh et al. (2014). However, the grain-surface reactions are sensitive to the choice of grain surface parameters. In this work, the barrier height for quantum tunnelling bqt between grain surface sites has been set to bqt = 1 Å, and the ratio of diffusion energy Ediff between grain surface sites and molecular binding energy Edes has been set to EdiffEdes = 0.5. This is a conservative take where most grain-surface chemistry happens at moderate temperatures (≳ 25 K). If bqt was to be increased, the hydrogenation reactions would be slowed down, and vice versa for decreasing bqt. If EdiffEdes was to be decreased, then grain-surface radical-radical recombination could happen at lower temperatures. Such changes would affect chemical evolution, but an exploration of this parameter space is left for future work.

3.6.2 Dependence on ionisation

As highlighted in several previous studies, including Paper I, Aikawa et al. (1996), Helling et al. (2014), Cleeves et al. (2014a) and Walsh et al. (2015), the chemical evolution taking place in protoplanetary disk midplanes is highly dependent on the degree of ionisation. The timescale of chemical evolution has been found in this work to scale inversely with the ionisation level. It is found that an ionisation level 10 times higher than that presented here will shorten the chemical evolution time scale 10 times, and vice versa.

In the situation where the midplane is significantly shielded from CR ionisation, for example by a high column density attenuating the CR particles, or by magnetic field effects (see Cleeves et al. 2013a), then no significant chemical evolution takes place over the disk lifetime of 7 Myr. In the case of complete exclusion of CRs, this is seen in Fig. 3a, where the ionisation level is solely based on SLR decay products in the midplane. This ionisation level starts out between 10−19 and 10−18 s−1, and drops orders of magnitude over the evolution time, as seen in Fig. 1c. Thus, for midplane ionisation levels of a few times 10−19 s−1 or lower, 7 Myr of evolution is insufficient time for the ionisation to cause significant chemical evolution. Along the same lines, if chemical steady state is to be achieved within 1 Myr of evolution, in the outer disk, then the ionisation level has be ~ 10−16 s−1 (about ten times higher than assumed for high ionisation in the outer disk here) or higher. The effects of X-ray ionisation on disk chemistry in planet-forming midplanes will be explored in future work.

3.6.3 Choice of grain sizes and distribution

Like in Paper I, the physical grain parameters have been kept constant throughout the evolution time. However, in real disks several processes take place for grains during the lifetime of a disk, all of which could change the grain parameters adopted here. These processes include grain settling to the midplane, agglomeration and fragmentation, as well as grain radial drift. The physics of grain evolution has been studied extensively (see e.g. Birnstiel et al. 2016; Kataoka et al. 2014; Okuzumi et al. 2012, and references therein).

One way of theoretically investigating the effects of changing dust parameters, is to simply adopt a different dust size as input for the models. Increasing the dust size will reduce the total surface area-to-mass ratio, prolonging the freeze-out and desorption timescales, and vice versa for decreasing dust size. However, the grain-surface chemical evolution of the ices will likely simply be slowed down for a reduced surface area-to-mass ratio, and sped up for an increased surface area-to-mass ratio. Thus, the same chemical effects might be seen, just on different timescales (see Sect. 2.3). Recent work by Facchini et al. (2017) found that the dust growth timescale in the disk midplane is significantly longer (3 orders of magnitude) than the timescale of freeze out and desorption, for example for CO.

The dust growth timescale is, however, still shorter than the timescale for chemical evolution. Although this means that the grain size utilised for grain-surface reactions here in fact is not static (as assumed) over the chemical evolution time, the grain growth might not affect the chemistry that differently from what has been presented in this work: grain growth likely does not accumulate small compact spherical grains into larger and compact grains with a predictable change in grain surface areas. Rather the growth may lead to fluffy aggregates or grains with uneven surfaces (Blum & Wurm 2008). Hence, the change to the grain surface area-to-mass ratio may not be significant. If grains are subsequently compactified resulting in a reduced grain surface area-to-mass ratio, it may simply lead to the aforementioned longer timescale for grain-surface reactions.

4 Discussion

4.1 Iceline shifts: comparison to other studies

The positions and movements of volatile icelines have been subject to much investigation recently, see Piso et al. (2015), Ali-Dib (2017), Öberg & Bergin (2016), and references therein. Piso et al. (2015) modelled shifting iceline positions over 3 Myr assuming both radial drift of particles, and a viscous disk with gas accretion. In these models, they found shifted icelines of volatiles compared to a static disk situation. These shifts, however, are not caused by evolving temperature structures as is the case in the present work, but by radial drift of different sized particles. Thus, radial drift can also shift the icelines independently by setting different radii for desorption of icy grain mantles depending on the sizes of the grains. Despite these differences, both evolving temperatures and radial drift of particles act to shift the icelines inwards. Piso et al. (2015) found that the H2O iceline was shifted inwards by up to 40%, whereas the icelines of CO2 and CO were shifted inwards by up to 60% and 50%, respectively. These shifts are comparable in size to those found in this work, but not caused by the same mechanism. Thus, the combination of different physical disk evolution effects could possibly result in even larger iceline shifts.

4.2 What happens to gaseous CO?

The dominant ices in the outer disk midplane at the end of evolution are all stable molecules, for both choices of initial abundances. For the reset scenario this results from hydrogenation of carbon and oxygen atoms on the grain surfaces, whereas for the inheritance scenario, the oxygen and carbon atoms are mainly supplied from the CR-induced photodissociation of CO and CO2.

CO gas has been studied extensively through observations of protoplanetary disks. As noted in the introduction several observational studies have found surprisingly weak CO emission from disks, more than can be accounted for by simple freeze-out and photodissociation. In particular for TW Hya, Schwarz et al. (2016) reported depletion by two orders of magnitude through CO isotopologue observations, and Zhang et al. (2017) recently determined the CO midplane iceline with an associated freeze-out temperature of K, higher than expected from the binding energy of pure CO ice (~ 21 K). Zhang et al. (2017) attributed this higher freeze-out temperature to CO being frozen out onto a mix of H2O and CO2 ices. In relation to these reported depletions, the findings for CO in this work are interesting. The outer icy disk here is not the only region lacking CO. A decrease in CO gas abundance is also seen to take place for the inheritance scenario for the region between the CO2 and the CO icelines in Fig. 5c around ~3 AU, until 2 Myr of evolution. The effect of CO freezing out and subsequently reacting into CO2 and chemically converting into other species such as HCN and C2 H6 means that CO gas is effectively reduced in abundance by up to 2 orders of magnitude, within the CO iceline. This reduction in CO gas, however, is only happening so long as the CR-induced photon fluence is low enough to not dissociate CO2 into CO on the grains, as discussed in Sect. 3.2. The models show that chemical CO gas depletion in a disk midplane could therefore be a temporary effect. As discussed earlier, the timescale of CO depletion, which lasts until the abundance inversion takes place, will depend on the ionisation conditions, and thus the density, in the midplane, with a lower ionisation level prolonging the CO depletion timescale.

An additional large drop in CO gas abundance is seen starting by 5 Myr at ~10 AU for the 0.1 MMSN disk, see Fig. 5c (the inheritance scenario). Thus, if one was to determine the CO iceline from depletion of gas-phase CO isotopologues in an old disk like TW Hya, the desorption temperature for CO in this model could appear to be at ~29 K (similar to the freeze-out temperature for CO found by Zhang et al. 2017) at ~10 AU, rather than at ~21 K at 16 AU, for the 0.1 MMSN disk. This appearance of the CO iceline to be at a higher temperature and smaller radius than prescribed by the CO binding energy is due to the chemical depletion of CO gas, happening outside the O2 iceline, resulting in a “fake” CO iceline. A fake CO iceline was also seen in Paper I, although that was different from the one here: in Paper I the fake iceline was seen outside the CO2 iceline by 1 Myr, thus at a smaller radius and at an earlier evolutionary stage. The fake iceline seen here coincides with the thermal iceline of O2, and the reduction in CO gas abundance outside 10 AU is due to the freeze-out of O2. O2 ice is hydrogenated into HO2 ice on the grains, and subsequently dissociated into OH ice. This OH, in turn, reacts with colliding CO gas to produce CO2 ice, thereby reducing the CO gas abundance. This effect will be further investigated in future work.

In Fig. 8 the abundance evolution of the main carbon-bearing volatiles, CO, CO2, CH4, HCN, C2 H6, CH3OH are plotted at 10, 20 and 30 AU, alongside “Other”, which accounts for the carbon that is not contained in the plotted species. For all three radii, the main contributors to Other are HCO and HNC, as well as complex molecules and hydrocarbons such as C3 H4, CH2CO, H2 CCO, CH3OCH3 and CH3COCH3, which all show abundances >10−7 with respect to H nuclei. It is evident at 10 AU in panel a that CO2 locks up the majority of the carbon, through the grain-surface reaction with CO and OH, and CO gas is seen to decrease in abundance. Here, C2 H6 ice is accounting for about 10% of the elemental carbon, with less than 10% of the carbon locked up in Other.

At 20 AU in panel b the picture looks different after 7 Myr. The other species take up about 30% of the elemental carbon, which is more than any single species alone. At 20 AU, the temperature decreases from about 30 to 18 K during the evolution time. That means that at the end of evolution, all molecules are frozen out onto the grains, and due to the temperature regime, most speciesare mobile on the grain surfaces, and able to react to form larger complex molecules, like those included in the Other category.

In panel c at 30 AU, the temperature goes down to 15 K after 7 Myr, meaning all species, including H atoms, can stick to the grains. H atoms, along with lowered mobility of molecules on the grain surfaces, means that hydrogenation reactions become dominant. Thus carbon is locked up in CH4, HCN and C2 H6 and CH3OH, with only around 10% elemental carbon in other species.

Recent work by Yu et al. (2016, 2017) also find this chemical conversion of CO into more complex ices using an evolving disk model in the midplane. Their disks have temperatures of ~35–40 K at 20–30 AU. They also found CO destruction in the disk midplane, by modelling chemical evolution using the older RATE06 network. However, their models produce the C2 H5 radical in the ice, which in this work is hydrogenated into C2 H6 ice. Also, they find only about 5% increase in their CO2 ice abundance, indicating that the CO gas destruction is not feeding the production of CO2 ice on the grains. Reboussin et al. (2015) did see conversion of CO gas to amongst other CO2 and C2 H6 ices on the grains using the NAUTILUS model, in agreement with this work. Besides these studies, models focussing on analysis of chemical evolution in disk midplanes are limited. Cridland et al. (2017) used a chemical network to evolve chemistry in the midplane before forming planets out of the material, but did not focus on the chemical response to the (varying) physical conditions and reaction types dominating in the midplane. Helling et al. (2014) found a different effect of chemistry on the C/O ratios that what is seen in this work, because their treatment of grain-surface chemistry was simpler.

4.3 Evolving elemental ratios

An important aspect of chemical evolution in a disk midplane is the relation to the composition of exoplanets, particularly their atmospheres. In Paper I it was shown that unless the inheritance scenario with low ionisation was considered, the C/O ratio of the chosen model setup changed over the course of 1 Myr of evolution. It was demonstrated that the C/O ratio radial stepfunction in Fig. 1 of Öberg et al. (2011) was modified. In this work, the evolution of the C/O ratios are traced up to 7 Myr instead of 1 Myr, and the effects of both chemical evolution and the evolving temperature structure are shown.

The longer timescale causes larger changes to the dominant volatile abundances and significant changes to the traditional C/O ratio stepfunction, even more than what was seen Paper I. The larger abundance changes are apparent in Fig. 7, where changes of up to two orders of magnitude are seen from 1–7 Myr, thus significantly larger than the changes happening until 1 Myr (the case from Fig. 4 in Paper I). Considering the more realistic scenario of an evolving, cooling disk (see Fig. 3 panels b and d) rather than a static one (Fig. 3 panel c), similar trends in chemical abundance changes to the static case are seen. The main effect is that volatile icelines are shifted inwards with time, because the disks cool. The physical conditions at each timestep of the simulations determine the chemical composition at that timestep. This holds for all scenarios and tells us that the chemical timescale is faster than the physical timescale in these models, andthat a static disk model may be sufficient for modelling chemistry. However, a planet forming in the crossing-zone of an iceline may experience a drastic change to the composition of accreting material over time.

In Fig. 9a–d, the evolving C/O, C/H, and O/H and N/H ratios, respectively, in gas and ice are presented for the evolving0.1 MMSN disk. The yellow profiles are at 0.1 Myr only, where as yet no significant chemical evolution has taken place, and the ratios are therefore the canonical values corresponding to the input abundances. Thus, for C/O for example, this yellow profile resembles the stepfunction from Öberg et al. (2011).

For the C/O ratio evolution in panel a, it is clear at later times that the evolving chemistry has caused the C/O ratios in gas and ice to change. At 5 AU, the C/O ratio in the gas drops from 1.2 to 0.3 over the course of 7 Myr. For the ice at the same distance, the C/O ratio increases from 0.2 to almost 0.5, thus making the ice between 3 and 10 AU more carbon rich than the gas. By 7 Myr the C/O ratios in both gas and ice are below 0.5, and thereby the C/O gas and ice ratio profiles are significantly different than the Öberg et al. (2011) stepfunction, and the profiles correlate with the C/O ratio profiles for the reset scenario from Paper I. That indicates that steady state is close after 7 Myr of evolution with inherited abundances. It is expected that so long as the initial C/O < 1, the same trends would be seen in C/O ratio evolutions, if a different set of elemental ratios were chosen.

The C/H evolution in panel b shows a significant turnover at the CO2 iceline at 2 AU, with more carbon being contained in gas than in ice inside the iceline, and vice versa outside of it. This iceline is thus critical for the metallicity of the material, and has important implications for the carbon content of the forming planets. Besides this distinct iceline, other features of the C/H ratios are the decreasing level in ice between the CH3OH (at 1.3 AU by 7 Myr) and the CO2 icelines. This is due to the destruction of CH3OH on the grain surfaces, as well as the CO and CO2 dependent interconversion of C/H ratio in gas and ice between the CO2 and the O2 icelines (see Sect. 3.2).

Panel c features the O/H ratio evolution. Here the H2O, CO2 and O2 icelines at ~0.8, 2 and 10 AU, outline the radii of main changes. Between the H2O and CO2 icelines, the ice starts as two times more oxygen rich than the gas, but after 7 Myr, the gas can become 3 times more oxygen rich than the ice. Between the CO2 and O2 icelines, the ice starts an order of magnitude more oxygen rich than the gas, but after 7 Myr, the gas and ice are equally oxygen rich. This is due to the increasing abundance of gaseous O2.

Although this work so far has focussed on carbon and oxygen-bearing species, for completeness the evolving nitrogen content in the material is shown in panel d. The decreasing nitrogen ice content between 2 and 8 AU is due to the destruction of NH3 ice (iceline at ~2.5 AU) in this region, which causes an increased nitrogen level in the gas, with N2 (iceline at ~20 AU) becoming the main nitrogen gas carrier (see Paper I for further discussion). Schwarz & Bergin (2014) also modelled the evolution of nitrogen-bearing species in a disk, and found different partitionings of the nitrogen into HCN, NH3 and N2, depending on their model assumptions.

thumbnail Fig. 8

Time evolution of the main volatile species at 10, 20 and 30 AU, plotted on linear scales to ease comparison to Fig. 7 of Yu et al. (2016) and Fig. 2 of Yu et al. (2017). The H2O ice abundance in these plots is at least 80. A rapid decrease of gaseous CO is seen within 1 Myr. The CO2 ice is plotted at half its abundance to emphasize the steep drop for CO gas before 1 Myr at all three radii.

thumbnail Fig. 9

Elemental ratios of carbon, oxygen and nitrogen in gas and ice as a function of radius, at six timesteps: 0.1, 0.5, 1, 2, 5 and 7 Myr. Panel a: C/O ratios. Panel b: C/H ratios. Panel c: O/H ratios. Panel d: N/H ratios. Gas profiles are dashed, ice profiles are solid. The y-axes have different scales, in particular in panel a the y-axis is linear. The grey boxes indicate the radial regions that the iceline of a volatile move over during the evolution.

4.4 Impact for planet formation and planet atmospheric composition

Recent years has seen enormous effort put into the field of planet formation and planet population synthesis models. Physical effects like grain growth, migration and radial drift, and a planet’s acquisition of an atmosphere from a disk have been included in increasingly sophisticated models (see e.g. Lambrechts & Johansen 2012; Alibert et al. 2013; Ali-Dib et al. 2014; Madhusudhan et al. 2014; Benz et al. 2014; Bitsch & Johansen 2016, and references therein). Planet population synthesis models however have neglected chemical evolution, although Cridland et al. (2016, 2017) did consider kinetic chemistry to model the material that go into forming planetary atmospheres.

Mordasini et al. (2016) conducted planet formation synthesis models introducing a “chain” of models to self-consistently predict the atmospheric transmission spectra of a hot Jupiter exoplanet from modelling the formation of the planet as it forms in the disk. One prediction from their work was that for hot Jupiters less massive than a few MJup, heavy element enrichment in the envelope from impacting planetesimal polluters dominate the final atmospheric composition over the enrichment obtained from gas accretion. This is interesting in the context of the findings in this work, that the C/O, C/H, O/H and N/H ratios in both gas and ice evolve over time. For hot Jupiters, it is thus not necessarily the heavy element content of the gas in the disk that will determine the planet’s resulting atmospheric elemental ratios.

The findings of this work, that chemistry matters in disk midplanes under certain assumptions, may complicate the challenge of relating observed exoplanet atmospheric elemental ratios to the formation history of the planet. Not only is it shown in this paper that chemical evolution can change the radial elemental ratios in gas and ice over time, but as Mordasini et al. (2016) stated, the disk material – gas or ice – that the atmosphere of the planet is formed from, may have a mixed origin.

The result of these two effects can, however, be explored through modelling. By having a planet atmosphere form from different relative contributions of gaseous and icy material in the disk, and by assuming (as a start) one of the time evolved elemental ratio abundances presented here, the final elemental ratio in the atmosphere, and the resulting transmission spectrum, can be generated. It can then be assessed whether or not the effects are observable, with current or future observational facilities. This would also be a first step to better implement disk chemistry into the “chain” of models, as requested by Mordasini et al. (2016). First attempts at this were indeed carried out by Cridland et al. (2016, 2017). However, they only considered planets formed inside 4 AU. An expansion of their work would thus be welcomed.

Lastly, the outer disk volatile ice abundances presented here obviously inspire a comparison to observed abundances in comets in the Solar System, in order to constrain the models with Solar System conditions. Such a comparison is currently in progress and will be published in an upcoming paper.

5 Summary

In this paper, the implications of chemistry in evolving disk midplanes on the composition of planet-forming material have been investigated.

Chemical evolution is an ionisation-driven long-term effect, which for the disk structures considered here becomes significant after a few times 105 yr of evolution, and goes on until at least 7 Myr, in the region between the H2O and O2 icelines (~0.5–10 AU, for the 0.1 MMSN disk). Considering the more realistic scenario of an evolving, cooling disk rather than a static one, similar trends in chemical abundance changes to the static case are seen. The main effect is that volatile icelines are shifted inwards with time, because the disks cool. Hence, using a static disk structure representative of the conditions in the disk during the accretion epoch of a forming planet may be sufficient for modelling chemistry. Changes of similar magnitudes are seen for all scenarios for the 0.1 MMSN disk and the 0.55 MMSN disk. That means that assuming different disk masses does not significantly affect the chemical evolution. However, higher mass disks than those studied here should be explored.

After 7 Myr of evolution, inside the H2O iceline and outside the O2 iceline, the fingerprint of interstellar abundances has been erased, for the high ionisation level, and chemical steady state has been achieved. That is to say, for ice in the outer disk, it becomes irrelevant whether the inheritance or the reset scenario is assumed. In the inner tenths of an AU, chemical equilibrium conditions are likely achieved. An abundance inversion is seen for CO and CO2 between the CO2 and CO icelines, where the abundance evolutions invert at a time which depends on the ionisation level. For the high ionisation level, CO gas reaches its minimum abundance by ~1.6 Myr, and CO2 reaches its maximum abundance at ~2.5 Myr.

By 7 Myr CO gas is seen to deplete inside the CO thermal iceline, from outside the O2 iceline, which could mimicthe iceline being further in and at warmer temperatures than the CO binding energy prescribes. Ice species such as HCN and C2 H6 become important carbon-bearing species in the outer disk by a few Myr. Between the H2O and O2 icelines, O2 gas is the dominant oxygen bearing species. Just outside the O2 iceline, it is warm enough for molecules to move on the grain surfaces, but not cold enough for H atoms stick on the grains. This results in the production of larger, more complex molecules on the grain surfaces.

Chemical evolution changes the elemental ratios in the gas and ice significantly over time. This results in a region where the ice has a higher C/O ratio than the gas. In order to relate observed atmospheric C/O ratios in exoplanets to formation locations and conditions in the disk midplane, chemical evolution is important to consider, independenton whether an atmosphere is built primarily from gas accretion or from planetesimal impacts. Generally, chemical evolution acts to significantly lower the C/O ratio of gas, such that neither the ratio in gas or in ice is at or above unity, but rather the C/O ratios in both gas and ice evolve to be below 0.5 (given an overall C/O ratio of 0.34).

If planet formation lasts for longer than a few times 105 yr in the disk midplane and cosmic ray ionisation is present, then the chemical abundances and elemental ratios in gas and ice will change alongside with planet formation. Henceforth, chemical evolution needs to be addressed when predicting the chemical makeup of planets and their atmospheres, as well as when inferring formation locations of exoplanets based on their observed atmospheric elemental ratios.


The authors thank Amaury Thiabaud and Ulysse Marboeuf for making their disk density and temperature profile available and for many useful discussions. The authors also thank Nikku Madhusudhan for a helpful discussion regarding the C/O ratio evolution and its effect on the metallicity, and the referee for useful comments that improved the presentation of the paper. 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 prise. CW also acknowledges the Netherlands Organisation for Scientific Research (NWO, grant 639.041.335) and the University of Leeds for financial support.

Appendix A Physical structure for 0.55 MMSN disk

Figure A.1 show the evolving physical structure for the 0.55 MMSN disk, featuring temperature, density and ionisation level evolution.

thumbnail Fig. A.1

Same as for Fig. 1, but for the 0.55 MMSN disk.


  1. Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ali-Dib, M. 2017, MNRAS, 467, 2845 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ali-Dib, M., Mousis, O., Petit, J.-M., & Lunine, J. I. 2014, ApJ, 785, 125 [NASA ADS] [CrossRef] [Google Scholar]
  4. Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24 [NASA ADS] [CrossRef] [Google Scholar]
  6. Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C. 2014, in Protostars and Planets VI, eds. H. Beuther, R. F. Klessen, C. P. Dullemond, & T. Henning (The University of Arizona Press), 691 [Google Scholar]
  7. Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil (The University of Arizona Press), 751 [Google Scholar]
  8. Bieler, A., Altwegg, K., Balsiger, H., et al. 2015, Nature, 526, 678 [Google Scholar]
  9. Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35 [NASA ADS] [CrossRef] [Google Scholar]
  10. Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bitsch, B., & Johansen, A. 2016, A&A, 590, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013a, ApJ, 772, 5 [NASA ADS] [CrossRef] [Google Scholar]
  16. Cleeves, L. I., Adams, F. C., Bergin, E. A., & Visser, R. 2013b, ApJ, 777, 28 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014a, ApJ, 794, 123 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cleeves, L. I., Bergin, E. A., Alexander, C. M. O., et al. 2014b, Science, 345, 1590 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cridland, A. J., Pudritz, R. E., & Alessi, M. 2016, MNRAS, 461, 3274 [Google Scholar]
  20. Cridland, A. J., Pudritz, R. E., Birnstiel, T., Cleeves, L. I., & Bergin, E. A. 2017, MNRAS, submitted, [arXiv:1705.02381] [Google Scholar]
  21. Crossfield, I. J. M. 2015, PASP, 127, 941 [NASA ADS] [CrossRef] [Google Scholar]
  22. Davis, S. S. 2005, ApJ, 620, 994 [Google Scholar]
  23. Drozdovskaya, M. N., Walsh, C., van Dishoeck, E. F., et al. 2016, MNRAS, 462, 977 [NASA ADS] [CrossRef] [Google Scholar]
  24. Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98 [NASA ADS] [CrossRef] [Google Scholar]
  25. Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38 [NASA ADS] [CrossRef] [Google Scholar]
  28. Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Fischer, D. A., Howard, A. W., Laughlin, G. P., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. F. Klessen, C. P. Dullemond, & T. Henning (The University of Arizona Press), 715 [Google Scholar]
  30. Fraine, J., Deming, D., Benneke, B., et al. 2014, Nature, 513, 526 [NASA ADS] [CrossRef] [Google Scholar]
  31. Harsono, D., Bruderer, S., & van Dishoeck, E. F. 2015, A&A, 582, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Helling, C., Woitke, P., Rimmer, P. B., et al. 2014, Life, 4, 142 [NASA ADS] [CrossRef] [Google Scholar]
  33. Henning, T., & Semenov, D. 2013, Chemical Reviews, 113, 9016 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  34. Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  35. Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Kamp, I., Thi, W.-F., Meeus, G., et al. 2013, A&A, 559, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
  39. Loomis, R. A., Cleeves, L. I., Öberg, K. I., Guzman, V. V., & Andrews, S. M. 2015, ApJ, 809, L25 [NASA ADS] [CrossRef] [Google Scholar]
  40. Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJ, 794, L12 [Google Scholar]
  41. Marboeuf, U., Thiabaud, A., Alibert, Y., Cabral, N., & Benz, W. 2014, A&A, 570, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 [NASA ADS] [CrossRef] [Google Scholar]
  44. Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41 [NASA ADS] [CrossRef] [Google Scholar]
  46. Mousis, O., Marboeuf, U., Lunine, J. I., et al. 2009, ApJ, 696, 1348 [NASA ADS] [CrossRef] [Google Scholar]
  47. Mousis, O., Ronnet, T., Brugger, B., et al. 2016, ApJ, 823, L41 [NASA ADS] [CrossRef] [Google Scholar]
  48. Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19 [NASA ADS] [CrossRef] [Google Scholar]
  49. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [NASA ADS] [CrossRef] [Google Scholar]
  50. Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
  51. Piso, A.-M. A., Öberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Piso, A.-M. A., Pegues, J., & Öberg, K. I. 2016, ApJ, 833, 203 [NASA ADS] [CrossRef] [Google Scholar]
  53. Pontoppidan, K. M., Salyk, C., Bergin, E. A., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. F. Klessen, C. P. Dullemond, & T. Henning (The University of Arizona Press), 363 [Google Scholar]
  54. Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Science, 341, 630 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  55. Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Salinas, V. N., Hogerheijde, M. R., Bergin, E. A., et al. 2016, A&A, 591, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Schwarz, K. R. & Bergin, E. A. 2014, ApJ, 797, 113 [NASA ADS] [CrossRef] [Google Scholar]
  58. Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91 [NASA ADS] [CrossRef] [Google Scholar]
  59. Seager, S. & Deming, D. 2010, ARA&A, 48, 631 [NASA ADS] [CrossRef] [Google Scholar]
  60. Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59 [NASA ADS] [CrossRef] [Google Scholar]
  61. Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  62. Taquet, V., Furuya, K., Walsh, C., & van Dishoeck, E. F. 2016, MNRAS, 462, S99 [NASA ADS] [CrossRef] [Google Scholar]
  63. Thiabaud, A., Marboeuf, U., Alibert, Y., Leya, I., & Mezger, K. 2015a, A&A, 580, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Thiabaud, A., Marboeuf, U., Alibert, Y., Leya, I., & Mezger, K. 2015b, A&A, 574, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Udry, S., & Santos, N. C. 2007, ARA&A, 45, 397 [NASA ADS] [CrossRef] [Google Scholar]
  66. van ’t Hoff, M. L. R., Walsh, C., Kama, M., Facchini, S., & van Dishoeck, E. F. 2017, A&A, 599, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114 [NASA ADS] [CrossRef] [Google Scholar]
  68. Walsh, C., Millar, T. J., Nomura, H., et al. 2014, A&A, 563, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Walsh, C., Loomis, R. A., Öberg, K. I., et al. 2016, ApJ, 823, L10 [NASA ADS] [CrossRef] [Google Scholar]
  71. Weidenschilling, S. J. 1977, Ap&SS, 51, 153 [Google Scholar]
  72. Willacy, K., Klahr, H. H., Millar, T. J., & Henning, T. 1998, A&A, 338, 995 [NASA ADS] [Google Scholar]
  73. Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
  74. Yu, M., Willacy, K., Dodson-Robinson, S. E., Turner, N. J., & Evans, II, N. J. 2016, ApJ, 822, 53 [NASA ADS] [CrossRef] [Google Scholar]
  75. Yu, M., Evans, II, N. J., Dodson-Robinson, S. E., Willacy, K., & Turner, N. J. 2017, ApJ, 841, 39 [NASA ADS] [CrossRef] [Google Scholar]
  76. Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nat. Astron., 1, 0130 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

Midplane temperature and density (top) and ionisation rates (bottom) as functions of disk radius for different evolutionary times for the 0.1 MMSN evolving disk model. For this model, the starting conditions at 1 yr are given by the orange profiles shown here, which feature slightly different physical conditions than that for the static disk structure adopted by Paper I (see text for details).

In the text
thumbnail Fig. 2

Evolving disk masses for the two disks, normalised to the initial disk masses of 0.1 and 0.55 MMSN. For both disks, more than half of the initial mass is lost before 0.5 Myr, and by 1–1.5 Myr evolution, only 20% of the initial masses remain.

In the text
thumbnail Fig. 3

Abundances after 7 Myr evolution for four different model setups. Panel a: inheritance scenario with low ionisation and evolving disk structure. Panel b: reset scenario, high ionisation and evolving disk structure. Panel c: inheritance scenario with high ionisation and static (unchanging) disk structure. Panel d: inheritance scenario with high ionisation and evolving disk structure. The blue, red-and-black and green arrows to the right of each panel indicate the initial abundances assumed for H2O, CO/CO2 and CH4, respectively, in the inheritance scenario.

In the text
thumbnail Fig. 4

Evolving iceline positions for main volatiles, for 0.1 MMSN (left) and 0.55 MMSN (right) evolving disks.

In the text
thumbnail Fig. 5

Abundances of H2O, CO and CO2 gas (dashed) and ice (solid) as a function of midplane radius at five timesteps: 0.5, 1, 2, 5, and 7 Myr. The left-hand and right-hand panels show the 0.1 MMSN and 0.55 MMSN disk results, respectively.

In the text
thumbnail Fig. 6

Same as Fig 5, but for CH4 and O2.

In the text
thumbnail Fig. 7

Time evolution of selected volatiles, in gas as ice at 1, 5, 10 and 30 AU, for the 0.1 MMSN disk with high ionisation and inheritance.

In the text
thumbnail Fig. 8

Time evolution of the main volatile species at 10, 20 and 30 AU, plotted on linear scales to ease comparison to Fig. 7 of Yu et al. (2016) and Fig. 2 of Yu et al. (2017). The H2O ice abundance in these plots is at least 80. A rapid decrease of gaseous CO is seen within 1 Myr. The CO2 ice is plotted at half its abundance to emphasize the steep drop for CO gas before 1 Myr at all three radii.

In the text
thumbnail Fig. 9

Elemental ratios of carbon, oxygen and nitrogen in gas and ice as a function of radius, at six timesteps: 0.1, 0.5, 1, 2, 5 and 7 Myr. Panel a: C/O ratios. Panel b: C/H ratios. Panel c: O/H ratios. Panel d: N/H ratios. Gas profiles are dashed, ice profiles are solid. The y-axes have different scales, in particular in panel a the y-axis is linear. The grey boxes indicate the radial regions that the iceline of a volatile move over during the evolution.

In the text
thumbnail Fig. A.1

Same as for Fig. 1, but for the 0.55 MMSN disk.

In the text

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

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

Initial download of the metrics may take a while.