© C. Eistrup and T. Henning 2022

1 Introduction

Exoplanets are thought to form in protoplanetary disks around young stars. Exoplanet atmospheres can be observed, chemically characterised, and molecular constituents of the atmospheres can be identified. More generally, carbon-to-oxygen (C/O) ratios of exoplanet atmospheres can be retrieved from directly imaged exoplanets with existing facilities. These facilities include the gravity experiment at ESO's Very Large Telescope Interferometer (e.g. the exoplanets β Pic b and HR 8799e, see GRAVITY Collaboration 2020; Mollière et al. 2020) and the IGRINS instrument at the Gemini-South Observatory, which was used to retrieve the C/O ratio of the day-side hemisphere of the transiting exoplanet WASP-77Ab (see Line et al. 2021). Such chemical characterisation of exoplanet atmospheres will be accelerated by the James Webb Space Telescope (JWST), both in terms of precision (lower uncertainty on C/O ratios) and in terms of the number of exoplanets for which atmospheres will be chemically characterised (~88)¹.

Insights into atmospheric compositions of exoplanets is a window into understanding what the atmospheres, and the exo-planets as a whole, formed from (see e.g. Mordasini et al. 2016). The idea is that the elemental chemical compositions of exoplanet atmospheres should reflect the combined elemental compositions of the material(s) that went into forming the exoplanet atmospheres. These are the materials that exist in pro-toplanetary disks, and which go into forming planets: gas and solids, where solids have both a refractory (rocky, with high desorption temperatures and molecular binding energies) and a volatile (icy, with low desorption temperatures and molecular binding energies) component. The region of protoplanetary disk midplanes where volatile species start to freeze out as ices usually have temperatures below 150K. As the volatile molecule with the highest molecular binding energy, H₂O freezes out below this temperature (see, e.g. Fraser et al. 2001; Penteado et al. 2017, from the former of which E_bind = 5770K was used for this study).

Planetary atmospheres may be formed by run-away accretion of gas surrounding a solid planetary core (as modelled by Madhusudhan et al. 2014; Notsu et al. 2020, with and without protoplanetary disk chemical evolution, respectively), and accretion of solid bodies moving in disks. There is growing evidence suggesting that micron-sized grains grow during the process of planet formation. When these bodies reach sizes of pebbles, they drift inwards in the disk (see, e.g. Lambrechts & Johansen 2012; Bitsch et al. 2019; Trapman et al. 2019). These pebbles are likely to be accreted by forming planets as they drift. Such accretion could lead to the material contained in the pebbles becoming part of the forming planet, and its atmospheres. In other words, the accretion of pebbles by exoplanets could affect the chemical makeup of the planetary atmospheres.

In order to establish a link between the chemical composition of a planetary atmosphere, and the formation history of the planet and its atmosphere, it is important to constrain the chemical composition of the volatile ices carried on these pebbles, as they drift from the outer, colder regions of the disk, inwards. The volatile ices are of special importance, since these ice components are the main carriers of carbon, oxygen, nitrogen and sulphur, which go into forming the gas-phase phase species constituting exoplanet atmospheres. Pebble drift in protoplanetary disks is an active area of research (e.g. Brauer et al. 2008; Birnstiel et al. 2012, 2015), including the accretion of pebbles onto planets (Mordasini et al. 2016; Bitsch et al. 2015, 2020; Bitsch & Johansen 2016; Bitsch & Battistini 2020; Cridland et al. 2016, 2017), and how this accretion affects the resulting compositions of the exoplanet atmospheres (see, e.g. Madhusudhan et al. 2014; Drummond et al. 2019).

However, most of these efforts have, thus far, assumed a predefined ice abundance for the pebbles, which remains unchanged from the beginning of simulations until the pebbles are accreted onto the forming planet. An exception is the desorption of volatile ices, which is assumed to happen when the disk temperature around the pebbles reaches a certain level, corresponding to an icy species desorbing into the gas-phase. The work by Booth & Ilee (2019), did utilise a chemical kinetics network through the KROME-package (Grassi et al. 2014) to model chemical evolution on drifting solids, although they only considered hydrogenation reactions for ices, and not two-body ice reactions via the Langmuir-Hinchelwood-mechanism (see, e.g. Hasegawa & Herbst 1993; Herbst 2021). Lastly, the work by Soto et al. (2022) coupled disk gas and dust evolution with chemical evolution, including pebble drift. Their results will be discussed in Sect. 4.5.

This paper investigates how the chemical evolution in the volatile ices on growing grains, and the subsequent chemical evolution of pebbles moving through a disk midplane changes the chemical ice composition of the pebbles. The volatile ices will be the focus, as these are the dominant carriers of elemental carbon and oxygen, apart from carbonaceous grains and silicates. The investigation will aid a better understanding of the actual icy compositions of these pebbles as they get accreted onto a forming planet, and will, in turn, provide a better basis of making the connection between characterised exoplanet atmosphere and protoplanetary disk midplane evolution. This will improve the framework that is needed when JWST results of characterised exoplanet atmospheres will be put into the context of how these exoplanets originally formed in their natal protoplanetary disks.

2 Methods

This paper utilises the Walsh chemical kinetics-code (used in Walsh et al. 2015; Eistrup et al. 2016, 2018), which is based on the UMIST RATE 12 network (McElroy et al. 2013), and includes gas-phase reactions, gas-grain interactions, grain-surface reactions (including the Langmuir-Hinchelwood mechanism), and accounts for grain growth (following Eistrup et al. 2022). The focus here is the midplane of a protoplanetary disk. This region is considered well-shielded from radiation from the central star. However, cosmic rays (CRs) are assumed to be able to penetrate through the upper layers of the disk, and, with surface density-dependent attenuation length of ${\lambda }_{\text{CR}}\text{=96}\frac{\text{g}}{{\text{cm}}^2}$ (see Umebayashi & Nakano 1981, 2009), these CRs are able to ionise H₂ and He, which subsequently produces a secondary UV-field.

These UV photons can then photo-dissociate species in the gas phase and in the ice, with reaction rate coefficients k_pd in both phases calculated using:

$$k_{\text{pd}}=\alpha {\left(\frac{T}{300\text{ K}}\right)}^{\beta }\frac{\gamma }{1-\omega }{s}^{-1},$$(1)

where α is the CR ionisation rate (surface-density-attenuated) seen in the bottom plot of Fig. 1, β accounts for the fact that some species are photo-dissociated by discrete line absorption (see Gredel et al. 1987; Woodall et al. 2007; Chaparro Molano & Kamp 2012), γ is the efficiency of the cosmic-ray ionisation events as defined in Eq. (8) of Gredel et al. (1989), and ω is the albedo of dust grains in the far-UV. For two-body, grain-surface reactions, a ratio of diffusion-to-binding-energies $\frac{E_{diff}}{E_{bind}}=0.5$ is assumed. Two monolayers of molecules are assumed to be chemically active in the grain surface ices. Furthermore, for grain-surface hydrogenation reaction, the barrier width for quantum tunneling between surface sites is set to q_qt = 1Å. Both these assumption are motivated by Cuppen et al. (2017).

The Walsh code is used to track the chemical evolution of pebbles in a protoplanetary disk midplane. The chemical code considers two different physical effects, as described below.

Grain growth is assumed to take place locally in the disk midplane, at two specific radii (128 and 200 AU), with constant local temperatures of 25 and 19.5K. The two radii are selected to represent different starting locations of pebbles migration, with associated differences in starting temperatures (at 25 K CO is in the gas phase, whereas at 19.5 K, it is frozen out as ice). The grain growth follows the models by Krijt et al. (2020), which tracks the total available grain surface area of the grain population, along with the number of grains of different sizes. For each time step in the evolution of this grain population, the chemical code is run to track the chemical evolution, as it depends on temperature, density, ionisation level and available grain surface area.

The chemical evolution is run alongside grain growth for 0.1 Myr, locally at both radii. This represents the time it takes to grow grains from micron-sizes, where the grains are wellcoupled with the local gas motion, to mm-sizes, at which point the grains are assumed to dynamically decouple from the gas, lose orbital angular momentum to the gas, and hence start drifting inwards in the disk (see, e.g. Krijt et al. 2020).

After 0.1 Myr of grain growth and chemical evolution at each radius, all solid material is assigned the ice composition resulting from the local chemical evolution (see Table 1), and all solids are assumed to be compact, spherical pebbles of radius R = 0.6 mm, with uniform chemical compositions. These pebbles are subsequently assumed to drift inwards through the midplane at constant speed. In reality, it is predominantly pebbles that drift (solid bodies with Stokes numbers ~ 1), and neither smaller grains (that are dynamically coupled to the gas), nor larger solids (which are dynamically independent of the gas) are expected to drift as efficiently as pebbles. However, in this work, it is assumed that upon reaching 0.1 Myr of local chemical evolution and grain growth at a given radius, all grain mass in the grain population is carried in 0.6 mm sized spherical pebbles.

Following Eq. (2) from Birnstiel et al. (2012), a silicate volume density for the pebbles (3 g cm⁻³), and a MMSN surface density structure² from Eq. (2.6) in Hayashi (1981), a pebble size of 0.6 mm corresponds to Stokes numbers St = 0.24 (at disk mid-plane radius R = 128 AU) and St = 0.47 (at R = 200 AU). Using these Stokes numbers, the radial drift speed v_r of the pebble drifting from 128 AU or 200 AU can be estimated using Eqs. (17) and (18) from Brauer et al. (2008):

$$v_{\left(\text{r}\right)}=-\frac{2v_{\text{n}}}{\text{St}+\frac{1}{\text{St}}}=-\frac{2\frac{c_{\text{s}}^2}{2V_{\text{K}}}\left(\delta +7/4\right)}{\text{St}+\frac{1}{\text{St}}},$$(2)

where $c_{\text{s}}=\frac{k{T}_{\text{gas}}}{\mu }$ is the isothermal sound speed for a given midplane gas temperature, T_gas, and mean molecular weight, μ (taken here to be μ = 2.3m_proton for a gas mainly composed of H₂ and He), and $V_K={\Omega }_{\text{K}}r=\sqrt{\frac{G{M}_{★}}{r^3}}$ is the Keplarian velocity, with the gravitational constant G, the central stellar mass M_★ = M_⊙, and r being the disk midplane radius. δ in Eq. (2) is the power-law index in the assumed disk surface density profile (Hayashi 1981, from Eq. (2.6)), taken to be 1.5 in this work.

Based on these assumption, T_gas = 25 K at r = 128 AU and T_gas = 20 K at r = 200 AU, the drift speeds for the 0.6 mm-sized pebble is 25 m s⁻¹ at 128 AU and 43 m s⁻¹ at 200 AU. Taking these drift speeds to be constant throughout drift, the drift time scales τ_d to get from 128 AU to 1 AU is τ_d (128 AU) = 2.4 × 10⁴ yr, and τ_d(200AU) = 2.2 × 10⁴ yr. If the size of the pebble corresponded to at Stokes number St= 1 at r = 200AU, then the drift speed would be v_r = 55 ms⁻¹ and the drift timescale would be τ_d = 1.7 × 10⁴ yr. This is therefore taken to be the minimum drift timescale. Several factors, including pressure bumps, may cause the actual timescale to be longer than this minimum. As such, the short drift timescale calculated above is derived assuming a smooth disk. This is further discussed in Sect. 4.6.

From these drift speeds and drift timescales, three different bracket choices are made for the drift speeds of the pebbles. These choices capture the minimum drift timescale, as well as two longer drift timescales, and are motivated by Brauer et al. (2008); Birnstiel et al. (2012, 2015); Birnstiel & Andrews (2014):

1. Slow, where the drift time from origin (128 AU or 200 U) to 1 AU is 1 Myr.

2. Medium, where the drift time from origin to 1 AU is 100 kyr.

3. Fast, where the drift time from origin to 1 AU is 10 kyr.

Pebbles drifting from both 128 AU and 200 AU were assumed to be of size R = 0.6 mm (and each starting location features a unique ice composition for the pebbles, see Table 1). This was the grain size assigned to all grains when running the chemical code during radial drift. In order to conserve the total grain mass (which was assumed constant), the grain number density with respect to H nuclei, N_grain, is scaled down in accordance with how much mass is carried in a 0.1 μπι (spherical) grain versus how much is carried in a R_adopt = 0.6 mm pebbles, giving us an adopted grain number density per H nuclei:

$$N_{\text{grain,adopt}}\text{=}{\left(\frac{0\text{.1μm}}{0\text{.6mm}}\right)}^3\text{×1}{\text{.3×10}}^{\text{-12}}{\text{H}}^{\text{-1}}{\text{=6×10}}^{\text{-24}}{\text{H}}^{\text{-1}},$$(3)

where 1.3 × 10⁻¹² is the grain number density with respect to H nuclei, assuming that all grains are made of silicates, are 0.1 μm in radius, and a gas-to-dust ratio of 100.

To ease the implementation of the grain drift into the chemical code, the treatment of the chemical evolution of the drifting solids is simplified. While in reality drifting solid bodies move through different physical environments (with increasing temperature and density, which is included here), the drifting solids likely also move through gas and icy solids in the disk midplane that have experienced different types or degrees of chemical evolution than have the solids drifting from further out in the disk. For this reason, it is reasonable to assume that an ice-covered drifting solid body will experience, not only increasing temperatures and densities, but also different chemical compositions of the surrounding gas than the drifting body experienced further out. Furthermore, when crossing icelines of different molecules, it is likely that the gas density inside the iceline will be higher and have a higher abundance of the molecule associated with the iceline, stemming from the desorption of that given molecule off the grain surfaces of solids, which have drifted across the iceline and deposited this molecule as gas, prior to the drifting of the solid bodies that are tracked here.

Here, this situation is simplified by assuming the following: the material that the models start out with stays the same throughout the chemical evolution. That is to say, the material is considered a closed parcel, where no material is added or lost. Within this parcel, solids and gas coexist, chemical reactions can take place between volatile chemical species, and molecules can freeze out onto solid surfaces, or desorb off the surfaces. However, the sum total of elemental C, O, N and S, and the total mass of solids stay constant (see Table 2). In order to mimic the effects of drift on the chemical composition of the solids in the parcel, the parcel of gas and solids is exposed to changes in physical environment (by changing temperature, particle density and ionisation level) according to the different physical conditions found further in in the disk midplane, see Fig. 1. The temperature and density structures in the disk midplane are calculated following the prescriptions by Umebayashi & Nakano (2009). The ionisation structure is calculated assuming a cosmic ray ionisation rate of ξ = 10⁻¹⁷ s⁻¹, and then accounting for surface density-dependent attenuation of this ionisation rate at each radius in the disk midplane, following the prescription from Umebayashi & Nakano (2009).

This simplified treatment of pebbles drifting in a disk mid-plane is justified by the fact that increasing temperatures likely lead to less ice chemistry, and accretion of gas molecules on the solid body surfaces, than would be the case assuming constant or decreasing temperatures. Furthermore, some of the timescales for drift considered here are short (10kyr and 100 kyr) in comparison with timescales for ice chemistry on grains (see, e.g. Eistrup et al. 2018; Eistrup & Walsh 2019), so the effect of the solid bodies drifting through gas of varying compositions in reality may be negligible for the composition of the ices of these solids. However, modelling assuming ice-covered pebbles to be moving through varying compositions of gas corresponding to different regions of the disk midplane, should be pursued in future work. This may be especially interesting if assuming that ice species on pebbles that drift early desorb into the gas phase in the inner disk, thereby moving material from ice in the outer disk to the gas in the inner disk (see, Soto et al. 2022, and discussion hereof in Sect. 4.5). Pebbles that drift later may, therefore, in reality, drift through regions with higher gas densities than early-drifting pebbles, and this may increase the efficiency of gas-grain interactions for gas-phase molecules that collide with, and stick to, the surface of a grain.

Fig. 1 Temperature T (top), number density ρ (middle), and ionisation level ξ (bottom) structures in the disk midplane, as functions of radius. Temperature and density structures are calculated from the prescription for the Minimum Mass Solar Nebula structure in Hayashi (1981); Aikawa & Herbst (1999). The ionisation structure is calculated assuming a cosmic ray ionisation rate of ξ = 10−17 s−1, and accounting for surface density-dependent attenuation of this rate at each radius in the disk midplane, following the prescription from Umebayashi & Nakano (2009).

3 Results

Six models were computed, for each of the six permutations of the model setup:

• –

  Drifting starting at 128 AU or 200 AU.

• –

  Drifting to 1 AU taking 10 kyr, 100 kyr or 1 Myr.

The model results are show in Fig. 2. The left-hand panel of the figure shows drifting from 128 AU, whereas the right-hand panel is drifting from 200 AU. Ten chemical species are tracked from 128 AU (H₂O, CO₂ CH₃OH, NH₃, HCN, CH₄, SO, SO₂, H₂S and S₂), whereas 11 are tracked from 200 AU (CO additionally). The reason for excluding CO ice from the tracking starting at 128 AU is that the temperature there increases to 25 K, which is above the iceline temperature of CO, so throughout the tracking of solids drifting from 128 AU, CO is mostly in the gas phase, and not ice on the drifting solids.

The solid profiles in both panels are for slow drift (1 Myr), the dashed profiles are for medium drift speed (100 kyr), and the dotted profiles are for fast drift (10 kyr). Following the profile for an ice species from right-to-left in each panel corresponds to following the drifting pebbles as they move inwards at constant radial speed through the disk midplane. It is seen in both top panels that the 4 most abundant ice species (H₂O, CO₂ CH₃OH, NH₃, and CO in the right-hand panel) do not change abundances much when comparing the slow drift to the fast drift. For ice species less abundant than 10⁻⁵ (SO, SO₂, H₂S, S₂ and HCN), the changes in abundances across drifting speeds are larger, which reflects the fact that these species are inherently less abundant than the other four species, and are therefore more affected by small changes to the abundances of the other four species.

For all species, it is seen that encountering their icelines, they do not transfer from gas to ice instantly, but instead a gradual decrease in ice abundance is seen going inwards. This is a consequence of icelines being treated here, not as a fixed transition radius (with temperature), but as a balance between the effects of freeze-out (for gas molecules sticking onto grain surfaces as ice) and desorption (for ice molecule that desorb into the gas phase). These two effects are at play for all physical conditions, but under colder conditions the freeze-out term in the chemical kinetics code will dominate and most molecules of a given species will exist as ice on grain surfaces, and under warmer condition the desorption term will dominate, and most of the molecules will exist in the gas-phase. Around the iceline of a species, the terms will be equal, and the molecules of a given species will co-exist, half in the gas, half in the ice. See, e.g. Eistrup et al. (2018) for the behaviour of the gas abundances in relation to the gradient of the ice abundances as shown here.

In order to gain a better impression of the difference that it makes to drift from 128 AU versus 200 AU, and the difference between slow and fast drift, comparisons between the initial abundances of the solids, and the abundances of the solids just before the drift across the iceline of different ice species are shown in Fig. 3. The y-axes for the top panels in this figure is linear, whereas it is logarithmic for the bottom panels. The values of the y-axes are the ratios of final abundances of each molecule (taken before the solid has crossed the iceline of given ice species, see fourth column in Table 1) to the initial abundances of each molecule (which are the initial abundance values at the beginning of the drifting). The initial abundances of the ice species are different, depending on the initial drift location, see Table 1. For each of the eight ice species featured in Fig. 3, the fourth column in Table 1 indicates which radii (same for both drifting speeds) the final abundances of these species are being adopted from.

The x-axis categories are the same for all nine species. The first three bars for each species are for fast drift (10 kyr), medium drift (100 kyr) and slow drift (1 Myr). The first and third, and fifth bars are for drifting from 128 AU, and the second, fourth and sixth bars are for drifting from 200 AU. For CH₃OH ice in the top row, and H₂S in the bottom row, a horizontal line indicates a y-axis-value of 1, corresponding to a net zero change in abundance of the species during the drifting.

For H₂O ice it is seen that drifting, no matter the speed or origin, causes a decrease in the abundance. The decrease is smaller (<0.4%) for faster drift than for medium drift (<2%) and slow drift (~ 10–20%). For CO₂, the opposite is seen: a small increase (<1%) for fast drift, a larger increase (<10%) for medium drift and a more significant increase for low drift of 50% when drifting from 200 AU, and 65% when drifting from 128 AU. This trend, either larger increase/decrease in abundances with slowing drift speed is visible for all species except CH₃OH and SO₂ ice, suggesting that the faster the drift, the less changes to the initial abundances. This is logical, since faster drift allows less time for chemical processing in the ices, than does slower drift. For the species in the top row, the changes in ice abundances for fast drift is <5%. For medium drift (100kyr), the changes in abundances range from <2% for H₂O ice and <5% for NH₃ ice, to <8% for CO₂ ice and less than 15% for HCN ice and CH₃OH ice. H₂O ice and NH₃ ice both decrease in abundance for all drifting scenarios, as compared with the initial abundances, whereas HCN ice and CO₂ ice increase in abundance for all scenarios. CH₃OH ice decreases in abundance for all scenarios, except for slow drift originating at 200 AU, in which case there is a ~5% increase in abundance.

For sulphur-bearing species in the bottom row, it is generally seen that SO ice and SO₂ ice increase in abundances for all drift speeds (and S₂ ice for fast and medium drifts), but most for slow drift (albeit from low initial abundances). The elemental oxygen for this increase is sourced from H₂O, which decreases in abundance for all drifts, except fast drift from 200 AU (in which case it increases in abundance by ~10%).

Fig. 2 Abundances of volatile ices on a drifting pebble as a function of radius (time), where the pebbles moves inwards with time. Left is drift from 128 AU, right is drift from 200 AU. Top panels are for species with abundances larger than 10−5 wrt H (H2O, CO2 CH3OH, NH3 and CH4). Bottom panels are for species with abundances ranging from 10−8–10−5 (HCN, SO, SO2, H2S and S2). Solid profiles are for slow drift (1 Myr), dashed profiles are for medium drift (100 kyr) and dotted profiles are for fast drift (10 kyr). Colours correspond to the same species in both panels. We note that the top and bottom panels have different dynamical y-axis ranges.

Fig. 3 Bar plots corresponding to ratios of final abundances of different species relative to initial abundance at the onset of drifting. Initial abundances are seen in columns two and three in Table 1, and the radii at which final abundances are read out of the simulation are seen in column four in Table 1. The y-axis indicates the ratios of these values, for the six different permutations of modelling setup (on the x-axis). Drift time increases (drift becomes slower) going from left to right in each panel. We note that species which encounter their icelines further out (see Fig. 2), will experience less drifting time to change their initial abundances than will species with icelines further in. Top and middle rows have linear j/-axes, bottom row has logarithmic y-axis. We note that all y-axes feature different dynamical scales.

4 Discussion

This section puts the results from the previous section into the context of the astrophysical process of planet formation. This section will also discuss the shortcomings of the method used, and propose an improved framework to be developed in the future.

4.1 The matter of drift speed

This paper has utilised assumptions for drift speeds that vary by factors of up to 100. While realistic drift speeds have not been derived analytically or numerically from first principles, the speeds used are motivated by a range of drift speeds found in the literature (Birnstiel et al. 2012, 2015). The aim was thus not to provide a precise prediction to a specific pebble drift scenario, but rather to explore and compare the chemical changes that happens when assuming three different representative drift speeds.

It is clear from Fig. 2 that the changes to all abundances as a function of radius (time) is monotonous with changing drift speeds. In addition to this, it is also clear from the panels in Fig. 3 that the longer (drift) time is assumed, the larger changes to volatile abundances can be caused by chemical evolution. That means that assumption of either slower or faster drifts than assumed here can be, at least qualitatively extrapolated from the results presented here. However, it is noted that for the four most abundant species (H₂O, CO₂, CH₃OH and NH₃), the abundance changes amount to <70% of their initial abundances overall, and for fast and medium drift speeds, these changes are all <20%. These changes are thus smaller than the relative ratios of the different species individually. In other words, and as is also seen in Fig. 2, the above mentioned four species which are top four in the initial abundance hierarchy, do not switch places individually in this hierarchy.

Lastly, it is noted that the difference in abundances across drift scenarios for the four most abundant volatile ices is significantly larger for the slow drift scenarios, relative to the medium and fast drift scenarios. This is explained by the timescale of chemical evolution in ice in disk midplanes, which is usually >100 kyr. Thereby, it is only in the case of slow drift that significant abundance changes are seen, whereas the differences between fast and medium drifts are comparatively smaller.

4.2 Changes for main volatiles

As seen in the top left panel in Fig. 3, H₂O changes only by less than 20% for all drifting scenarios. It is, however, seen that all drifting scenarios cause a decrease in the abundance of H₂O ice, which is caused by the photodissociation of H₂O ice by CR-induced photons via the reaction

$${\text{iH}}_2\text{O }\underleftarrow{\text{γCR}}\text{ iOH + iH.}$$(4)

Here, "i" denotes an ice species. A generally decreasing trend (except for medium drift speeds originating at 200 AU) is seen for CH₃OH ice, which is also caused by CR-induced photodissociation:

$${\text{iCH}}_3\text{OH }\underleftarrow{\text{γCR}}{\text{ iOH + iCH}}_3\text{,}$$(5)

and the elemental oxygen and carbon lost from these two species are mainly processed into CO₂ ice, which sees an increase in abundance for all drifting scenarios, via the two-body grain-surface reaction:

$${\text{iOH}}_3\text{+ICO}\to {\text{iCO}}_2\text{ + iH.}$$(6)

Ice chemistry during drift therefore processes elemental oxygen carried in OH from H₂O ice and CH₃OH ice into CO₂ ice, with an interaction with the sticking CO from the surrounding gas being the source of elemental carbon in the form of CO that goes into CO₂ ice. The iCH₄ produced via reaction (5) further undergoes hydrogenation to form CH₄ on the grain surface, which, due to its low binding energy of E_bind = 1090K (see Table 1), desorbs into the gas, when the midplane gas temperature is T_gas > 25 K. These effects have been previously modelled and described in, for example, Eistrup et al. (2016, 2018); Bosman et al. (2018). An implication of this is that the drifting of any one pebble may lead to decreasing CO gas abundances in the midplane, not just at one radius, but along the entire radial drift trajectory that the pebble follows, as long as it remains outside the icelines of H₂O, CH₃OH and CO₂. Furthermore, the dissociation of CH₃OH ice leading to the production of CH₄ gas may enhance the CH₄ gas abundance along the pebble's drift trajectory, although the amount of CH₄ gas produced depends of the initial assumption for the abundance of CH₃OH ice.

For nitrogen, the abundance of NH₃ ice decreases for all drifting scenarios, due to CR-induced photodissociation producing NH₂ and NH. The elemental nitrogen in these species is then processed into N₂ which desorbs into the gas-phase for _Tgas > 19 K, and HCN ice, which sees an increase in abundance for all drift scenarios. For sulphur, H₂S decreases for all drift scenarios, and SO, SO₂ and S₂ all increase. Interestingly, SO ice is the dominant carrier of sulphur for almost its entire drift time towards its iceline (with the exception of the outermost part of the disk in the beginning of the drift, where H₂S ice dominates). Inside the SO iceline (R ~ 15 AU), SO₂ ice dominates, and between 20 and 40 AU S₂ ice peaks (for medium and fast drifts, see Fig. 2) with an abundances of ~10⁻⁶ wrt H, and is the second-most abundant sulphur-carrying species. It is noted that observational evidence by, e.g. Kama et al. (2019); Le Gal et al. (2021) suggests that most sulphur in disks may be carried in refractory material. In that context, it is possible that the volatile sulphur species modelled here only account for minor amounts of elemental S.

Figure 4 shows the ice compositions of the pebbles for the modelling scenarios with drift starting at 200AU. Each row represents a chemical element, and each pie chart shows how much of the total amount of an element is carried in different ice species. The first column shows the compositions at the beginning of drift. The second, third and fourth rows represents the final ice compositions for fast (10 kyr), medium (100 kyr) and slow drift (1Myr), respectively. These final compositions are taken at a midplane radius of 18.9 AU, which is the innermost radius at which none of the considered ices have desorbed yet (except for CO ice, which is also featured in the top left chart).

Reading the pie charts for each element, from left to right, corresponds to considering longer and longer evolution times for the chemistry. Since the pebbles in all scenarios drift at constant radial speeds, the "Fast" column features the pebble ice composition after ~9kyr (R = 18.9 AU is ~90% of the way from 200 AU to 1 AU on a 10 kyr timescale). Similarly, the "medium" and "Slow" columns features pebbles compositions after ~90 kyr and ~0.9 Myr of chemical evolution, respectively.

These pie charts show which ice species carry more of one element than other species, and how this changes with drift speed (evolution time). For carbon -and oxygen-carrying species in the top two rows, for example, there is a clear trend from left-to-right that CO₂ ice is increasing in abundance. In the top left chart, CO₂ ice accounts for slightly less than half of total elemental carbon at the beginning of pebble drift. The fraction of elemental carbon in CO₂ then grows, and for slow drift, CO₂ ice accounts for 70% of all carbon (primarily at the expense of destruction of CH₃OH ice). For elemental oxygen in the second row, the picture is similar for the CO₂ ice fraction of elemental oxygen (the abundance of CO₂ ice with respect to H₂ was counted double in these oxygen pie charts, as CO₂ carries two oxygen atoms per molecule), which grows with drift time. As such, the fraction H₂O:CO₂~5:1 at the start of drift (first column), but changes to H₂O:CO₂~3:2 for slow drift in the last column.

For nitrogen, NH₃ ice remains dominant for different drift timescales, although HCN ice is steadily increasing in abundance with evolution time. H₂S ice is dominant as a sulphur carrier at the start of drift, and for fast and medium drift. However, for slow drift, SO ice takes over as dominant sulphur carrier, with a tiny fraction of sulphur carried in SO₂ ice. Less abundant volatile species, such as C₂H₆, C₂H₄, C₂H₂ and HNC ices are not shown in this work. Out of these, C₂H₆ ice reaches an abundance of 4 × 10⁻⁶ wrt H at ~ 20AU for slow drift, but for all other drift timescales, these species remain < 10⁻⁶ wrt H (HNC ice is consistently less abundant than HCN ice, for all drift timescales).

Fig. 4 Pie charts showing the percentage of each chemical element (from top to bottom: carbon, oxygen, nitrogen and sulphur) that is carried in different volatile ices (the abundances of CO2 and SO2 were doubled when accounting for the oxygen budget in the second row, and the abundance of S2 was likewise doubled when accounting for the sulphur budget in the fourth row). Left column features the starting abundances for pebble drift starting at 200 AU. The second, third and fourth columns represent fast (10 kyr), medium (100 kyr) and slow drift (1 Myr), respectively. Going from left to right therefore corresponds to increasing the chemical timescales. The pie charts in these columns represent pebble ice compositions at 18.9 AU, which is the innermost midplane radius at which none of the shown volatile molecules have desorbed yet (see Table 1), except for CO ice, which is only present in the ices at the beginning of drifting (see "Carbon" and "Oxygen" pie charts in the first column).

4.3 Use for modelling planet formation via pebble accretion

One key goal for the study was to improve the understanding of chemistry in drifting solids, and to provide the planet formation modelling community with insights into these effects. Several efforts over the years have assumed that the chemical abundances assumed initially in models of planet formation would be conserved, even when solids in protoplanetary disks evolve, grow and move. The exception to this is the effects of freeze-out and desorption of molecules, which accounts for the transfers of material between the gas and ice phases.

An important insight from this study is that chemical evolution during drift does not significantly alter the volatile abundances of the dominant carriers of carbon and oxygen. This is the case, especially for slow and medium drift timescales. Quantitatively, the changes to the dominant volatile species, for slow and medium drift is <20%. Following this, when modelling the volatile composition of pebbles that drift from the icy, colder disk and inwards (to be accreted by a forming planet) it is reasonable to assume that the individual abundances of the ices on the pebbles remain largely unchanged during drift, with the exception of desorption of ices into the gas-phase, if icelines are crossed.

4.4 The impact of drift origin

This study has considered drift from two distinct midplane radii. Pebbles originating at each of the two started out with different ice compositions (see Table 1), corresponding to the abundance as its origin radius, as it has been evolved chemically for 100 kyr before drifting started. As seen from Fig. 3, the differences in chemical evolution between the two radii is generally between 30% and a factor of 2 (with the exception of SO and H2S). Because the actual drift times for drifting scenarios originating at both radii are the same (the fast drift speed from 128 AU is slower than the fast drift speed from 200 AU, because the latter needs to drifts almost twice as far in the same time), the changes between the chemical evolution happening from the two radii of origin is likely to be caused by the differences in initial ice abundances for the two cases. These initial differences are, however, not large (see Table 1).

4.5 Comparison to other work

Recently, Soto et al. (2022) modelled chemical evolution in an evolving disk with drifting pebbles, using the same chemical kinetics code as is used here (Walsh et al. 2015). They accounted for pebbles drifting through new material (gas and dust), as they moved through the disk, which is a more realistic approach that the closed parcel-approach that was used here. Their drifting time scales are on the shorter side of what is used here: for their models with stellar accretion rates of m = 10^-8 M_⊙ yr⁻¹, they find that only pebbles drifting from radii <10 AU make it to their inner disk region within 100 kyr. A radius of 10 AU in their Fig. 1 corresponds to a midplane temperature of 60–70 K, which is below their CO₂ midplane freeze-out temperature of ~80 K. For lower stellar accretion rates of m = 10⁻⁹ M₀ yr⁻¹, their pebble drift speeds are faster (drift times shorter), and so pebbles from the outer disk (> 10 AU) can influence their inner disk on at 100 kyr timescale.

Soto et al. (2022) predicts that pebble drift may transport a significant amount of volatiles from the outer disk to the inner disk, with the results that the outer disk loses volatile ices, and the gas abundances around volatile icelines in the inner disk increase. Figure 2, panels d–f (including pebble drift) in that paper also shows that the abundances of H₂O ice and CO₂ ice increase with decreasing miplane radius, which is taken to result from chemical evolution in ices during pebble drift. However, a direct comparison between their results for pebble ice compositions during drift, and the results presented in this paper, is not straightforward. This is because their work tracked the abundance of a given ice species at a given midplane radius at a fixed time, where the abundance of an ice species at that radius was made up of the contributions from different pebbles that had drifted from different radii (and had different sizes). In this work, in contrast, the ice composition at a given radius is made up of one (average) pebble of one size drifting from one of two fixed radii, which was tracked over time, but assuming different drift time scales.

4.6 Caveats of model assumptions

The modelling setup presented in this paper is a simplification of reality. Most importantly, the pebble drift chemistry modelled here assumes a close parcel of solids and gas, which is being heated up. This means that the drifting pebble is not drifting through new material (material residing further in, in the disk midplane), so there is neither new solids material and new ice compositions being encountered, nor is there new surrounding gas, and possibly higher gas densities around icelines (when a given species desorbs), and increasing gas densities generally moving inwards, as the ices from more and more drifting pebbles desorb and are added to the gas.

It is not straightforward to predict how the effects of a more complete framework including these effects would differ from the effects seen in this paper. However, given that JWST is soon to commence operations, it is timely and prudent for the planet formation community, and the astrochemistry community, to come together and develop appropriate modelling frameworks. These need to account for ice-covered pebbles drifting inwards through a disk, where it will encounter gas and solid material that has experienced a different chemical evolution from itself. The frameworks also need to track the effects that pebble drift chemistry has on the ice composition of the pebbles, at it is accreted onto a forming planet.

The opacity in the disk midplane is not evolving along with the dust growth or pebble drift in this work, and the ice mantles on the pebbles are not assumed so thick that parts of them are shielded from UV. Both these aspects should be addressed in future work, in order to make the treatment of the physical and chemical effects more realistic.

The modelling setup presented here only considered one set of initial abundances (per starting radius) for the drifting pebbles. These sets of initial abundances, in turn, were based on global elemental ratios of O:C:N:S of 520:180:60:6.2 (as used in Eistrup et al. 2016, 2018). Different choices of global elemental ratios may lead ices of different volatile to evolve differently than seen here. However, following results from Eistrup et al. (2018), one could predict that the same trends for increase and decrease of different volatile would be seen, as the ones seen in this paper, even if starting abundance levels and magnitudes of changes might be different. Future work will explore the impact of different chemical starting conditions.

The drift time scales assumed here do capture the minimum drift timescale as calculated in Sect. 2, but are generally longer than this. The 100kyr and 1 Myr timescales were chosen to explore chemical evolution in pebbles, assuming that pebbles do not drift inwards at a constant, maximum speed, but also accounting for such effects as dust traps (see, e.g. van der Marel et al. 2013; Baruteau et al. 2019; Pinilla et al. 2020, and reference therein). Such traps are generally assumed to either slow down, or halt the drift of solids in a disk. If assuming that a realistic disk does feature more radial structure in the gas and refractory component than what is given by the smooth MMSN power-law surface density profile used in this work, then the drift timescale becomes longer than the minimum timescale calculated earlier. In the context of such longer timescales, this work offers some predictions about the chemical evolution in pebbles, and how this evolution may be different for such a case, as compared to the case for fast drift.

Lastly, this paper has used a single pebble size of R = 0.6 mm for all chemical modelling. This size was motivated by grain growth models by Krijt et al. (2016), which also informed the accompanying paper Eistrup et al. (2022). From pebble drift dynamics there is a general trend that the larger the pebble, the faster it will drift. Adopting one pebble size for three different drift speeds is therefore not necessarily realistic. However, as was shown in Eistrup et al. (2022), the larger the pebble (or solid body, generally), the slower the efficiency of ice chemistry on the surfaces of the solid bodies. Therefore, by extension, it can be derived that if a larger pebble size than 0.6 mm had been used for modelling the chemical evolution in the ices over the drift timescales considered here, then this evolution would had been slower than what has been shown in this paper. The results shown here can therefore be considered as extremes for the degree to which chemical evolution can change the icy composition of pebbles. If larger pebbles had been considered for the modelling here, then smaller abundance changes would likely had been observed, but the general trends in the chemical processing (increasing or decreasing abundance for a given volatile ice as a function of pebble drift) would had remained.

5 Conclusion

This paper has presented a novel framework for modelling the chemical evolution that ice-covered pebbles with a size of R = 0.6 mm will experience as they drift inwards from two different radii of origins in a protoplanetary disk, on three different drift timescales. Below are the main conclusions summarized:

• –

  Drift speed matters. The slower a pebble drifts, and the longer time it therefore takes it to move from the outer disk to 1 AU, the larger the changes to the abundances of the volatile ice species on the pebble.

• –

  Generally, pebble drift chemistry causes CO₂ ice, HCN ice and SO ice to increase in abundance. Elemental oxygen is processed into CO₂ ice (and, to a minor extent, into SO ice), elemental carbon is processed into CO₂ (and, to a minor extent, HCN ice), elemental nitrogen is processed into HCN, and elemental sulphur is processed into SO. For all other volatile ices considered, pebble drift causes a decrease in their abundances.

• –

  For medium and fast drift (where the pebble is drifting for 10kyr or 100kyr), the four dominant volatile ices (H₂O, CO₂, CH₃OH and NH₃) undergo abundance changes of <20%. This means that, assuming such drift time scales, it is reasonable to assume that no chemical evolution happens in the ices, when pebble drift is modelled in relation to planet formation. However, for slower drifts, including possible slowing effects such as pressure bumps, it is possible that chemical evolution during the drift of the pebbles needs to be modelled explicitly.

In the context of understanding the chemical evolution that takes place in ices on drifting pebbles, it may be a reasonable approximation to assume that the chemical compositions in the ices remain unchanged during drift, if the drift takes place on a short time scale (tens of thousands of years). We note that a choice of chemical abundances in the ices on the pebbles still should reflect the chemical evolution that the material may have undergone prior to the drift starting.

Acknowledgements

The authors thank Bertram Bitsch and Dmitry Semenov for interesting discussions that improved the methodology and the quality of the manuscript. Also thank you to the anonymous referee. Research for this paper was supported by the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28.

References

All Tables

All Figures

Fig. 1 Temperature T (top), number density ρ (middle), and ionisation level ξ (bottom) structures in the disk midplane, as functions of radius. Temperature and density structures are calculated from the prescription for the Minimum Mass Solar Nebula structure in Hayashi (1981); Aikawa & Herbst (1999). The ionisation structure is calculated assuming a cosmic ray ionisation rate of ξ = 10−17 s−1, and accounting for surface density-dependent attenuation of this rate at each radius in the disk midplane, following the prescription from Umebayashi & Nakano (2009).

In the text

Fig. 2 Abundances of volatile ices on a drifting pebble as a function of radius (time), where the pebbles moves inwards with time. Left is drift from 128 AU, right is drift from 200 AU. Top panels are for species with abundances larger than 10−5 wrt H (H2O, CO2 CH3OH, NH3 and CH4). Bottom panels are for species with abundances ranging from 10−8–10−5 (HCN, SO, SO2, H2S and S2). Solid profiles are for slow drift (1 Myr), dashed profiles are for medium drift (100 kyr) and dotted profiles are for fast drift (10 kyr). Colours correspond to the same species in both panels. We note that the top and bottom panels have different dynamical y-axis ranges.

In the text

Fig. 3 Bar plots corresponding to ratios of final abundances of different species relative to initial abundance at the onset of drifting. Initial abundances are seen in columns two and three in Table 1, and the radii at which final abundances are read out of the simulation are seen in column four in Table 1. The y-axis indicates the ratios of these values, for the six different permutations of modelling setup (on the x-axis). Drift time increases (drift becomes slower) going from left to right in each panel. We note that species which encounter their icelines further out (see Fig. 2), will experience less drifting time to change their initial abundances than will species with icelines further in. Top and middle rows have linear j/-axes, bottom row has logarithmic y-axis. We note that all y-axes feature different dynamical scales.

In the text

Fig. 4 Pie charts showing the percentage of each chemical element (from top to bottom: carbon, oxygen, nitrogen and sulphur) that is carried in different volatile ices (the abundances of CO2 and SO2 were doubled when accounting for the oxygen budget in the second row, and the abundance of S2 was likewise doubled when accounting for the sulphur budget in the fourth row). Left column features the starting abundances for pebble drift starting at 200 AU. The second, third and fourth columns represent fast (10 kyr), medium (100 kyr) and slow drift (1 Myr), respectively. Going from left to right therefore corresponds to increasing the chemical timescales. The pie charts in these columns represent pebble ice compositions at 18.9 AU, which is the innermost midplane radius at which none of the shown volatile molecules have desorbed yet (see Table 1), except for CO ice, which is only present in the ices at the beginning of drifting (see "Carbon" and "Oxygen" pie charts in the first column).

In the text