© The Authors 2023

1 Introduction

It is commonly assumed that Solar System bodies have bulk compositions that are representative of the material present in the protosolar nebula (PSN) from which they formed. If the PSN was homogeneous in composition, gas and ice giants would be expected to reflect this homogeneity. However, observations show that giant planets present a range of supersolar metallicities. In Jupiter’s atmosphere, the abundances of volatile elements were found to be ~ 1.5–6.1 times higher than their protosolar values (Atreya et al. 2003; Mousis et al. 2018; Li et al. 2020), with a few exceptions attributed to interior processes. In Uranus and Neptune, volatile abundances can reach up to about 100 times their protosolar values (Lindal et al. 1987, 1990; Baines et al. 1995; Karkoschka & Tomasko 2009; Sromovsky et al. 2014).

Numerical models show that the dynamics of icy pebbles and their vapors around icelines is an efficient mechanism to produce local changes in the composition of the PSN (Booth et al. 2017; Desch et al. 2017). This process efficiently concentrates species around their respective icelines, creating the compositional gradient that may be responsible for the volatile enrichment in giant planets of our Solar System. The radial transport of trace species and the effect of icelines have been investigated using accretion disk models to assess the composition of the PSN. The resulting compositional profiles would then be used to constrain the formation conditions of gas giants (Mousis et al. 2019; Schneider & Bitsch 2021; Aguichine et al. 2022) or ice giants (Owen et al. 1999; Monga & Desch 2015; Mousis et al. 2020). This approach has also been used to explain the diversity among the cometary compositions and to determine their source or the general origins of their building blocks (Mandt et al. 2020; Mousis et al. 2021a).

The form taken by volatiles that have fallen from the interstellar medium (ISM) onto the PSN is still an open question. Icy pebbles that are present at the earliest stages of the PSN may have formed in the very cold environment of the ISM, at temperatures of 10 K or below (Gibb et al. 2004). At such low temperatures, H₂O condenses in an amorphous structure that can efficiently trap other volatile species (Mayer & Pletzer 1986; Jenniskens et al. 1995). At ~ 135 K, the amorphous ice transitions to crystalline ice and releases trapped volatiles (Bar-Nun et al. 2007). This temperature is much higher than the usual sublimation temperature of pure condensates. In the PSN, the heliocentric distance at which this release occurs is called the Amorphous to Crystalline Transition Zone (ACTZ), and is located at approximately 5 au (Mousis et al. 2019). However, it is not clear whether this amorphous ice survives the fall onto the PSN. For instance, it has been proposed that amorphous dust was heated up to crystallization when entering the PSN (Visser et al. 2013), as a consequence of the presence of shockwaves in the accreting PSN (Miura et al. 2017; Burkhardt et al. 2019). As a result, many circumstellar disk models treat volatile species as pure condensates, with sublimation temperatures computed from thermodynamic tables (Dodson-Robinson et al. 2009; Ciesla et al. 2015; Öberg & Wordsworth 2019). Alternatively, volatile species may also be trapped in clathrate form in the PSN, when sufficient amounts of H₂O are available (Lunine & Stevenson 1985; Gautier & Hersant 2005; Mousis et al. 2021b). A species, i, is usually trapped in clathrate at a higher temperature than the one needed to sublimate its pure condensate form, except for the cases of CO₂ and NH₃. Because CO₂ clathrate and NH₃ monohydrate form at lower temperatures than their respective condensates at nebular pressures, they are not considered in our model. In the following, we refer to the heliocentric distance at which volatiles are entrapped in or released from clathrates as the clathration line. For this reason, volatile species can remain adsorbed on amorphous ice or trapped in clathrates, then they are released several au closer to the Sun - than they otherwise would be if they were in the form of pure condensates.

In this study, we aim to quantify the influence of the presence of various icelines in the PSN, including the ACTZ and multiple clathration lines, on the nature of the main volatile reservoirs that were at play during the formation of the first icy grains in the disk. To do so, we use an existing PSN model that already describes the condensation and sublimation of pure ices, as well as the transport of species in solid and vapor forms (Aguichine et al. 2020, 2022), along with prescriptions of amorphous ice destabilization, as well as clathrate crystallization and dissociation added. Eleven key species are considered in our approach, namely: H₂O, CO, CO₂, CH₄, H₂S, N₂, NH₃, Ar, Kr, Xe, and PH₃. Each of these species can exist in the forms of vapor, crystalline ice, clathrate (monohydrate in the case of NH₃), or amorphous ice in the PSN. Two scenarios are investigated, each of them corresponding to a different initial state of the system. In scenario 1: the PSN is filled with volatile species in pure condensates or vapor form, depending on their location in the PSN. In scenario 2: the PSN is filled with volatile species adsorbed on amorphous ice beyond the ACTZ, and in vapor form in regions closer to the Sun. Because the ACTZ is located beyond the snowline, H₂O is found in crystalline form between the snowline and the ACTZ. Figure 1 represents the different forms of volatiles in the PSN in these two cases. In scenario 1, volatile species are in the vapor phase between the snowline and their respective icelines. If enough H₂O is available, these vapors can form clathrates. In scenario 2, icy grains of amorphous ice that drifted inward of the ACTZ release the adsorbed volatile species as vapors. The outward diffusion of these vapors can lead to the formation of clathrates or pure condensates (or both).

Section 2 describes the PSN model, as well as the transport modules used in our calculations. It also depicts the different formalisms that have been added to the disk model to mimic the formation and destabilization of clathrates, as well as the desorption of volatiles from the amorphous ice particles crossing the ACTZ. In Sect. 3, the radial profiles of the abundances of the different species are represented in various forms as a function of time in the PSN and in the individual cases described by the two scenarios. Section 4 is devoted to a discussion of the sensitivity of our results in light of the variation of the disk parameters. Our model is then used to compare the compositions of its resulting volatile reservoirs with those of the H₂O-poor comet C/2016 R2 PanSTARRS (R2), Jupiter, Uranus, and Neptune. This allows us to discuss the formation conditions of those bodies in the context of the two scenarios. Section 5 presents our conclusions.

2 Volatile transport and evolution model

In this section, we describe the protosolar nebula model employed in our simulations, along with the modules calculating the transport of dust particles and vapors within the disk. Source and sink terms related to sublimation and condensation of pure ices as well as to clathrate destabilization and formation within the disk are also depicted.

2.1 Protoplanetary disk model

The disk model used here is the one described in Aguichine et al. (2020) and Mousis et al. (2020). The evolution of the PSN is governed by the following differential equation (Lynden-Bell & Pringle 1974): (1)

which describes the time evolution of a viscous accretion disk of surface density, Σ_g, and viscosity, v, assuming invariance in the orbital direction and hydrostatic equilibrium in the azimuthal direction. This equation can be rewritten as a set of two first-order differential equations coupling the gas surface density, Σ_g, field and mass accretion rate, Ṁ: (2)

The first equation is a mass conservation law and the second one is a diffusion equation. The mass accretion rate is expressed as a function of the gas velocity field, v_g, and the radius, r, as Ṁ = −2πΣ_gv_gr.

The dynamical viscosity is calculated using the prescription of Shakura & Sunyaev (1973): (3)

where α is the viscosity coefficient, c_s is the sound speed in the PSN, and Ω_K is the Keplerian frequency; also, α is estimated to be in the 10⁻⁴−10⁻² range, based on models calibrated on disk observations (Hartmann et al. 1998; Hersant et al. 2004; Gautier & Hersant 2005; Birnstiel et al. 2012; DraZkowska & Alibert 2017; Armitage 2019). The sound speed, c_s, is expressed as follows: (4)

where μ_g is the mean molecular mass of the gas in the PSN, assumed here to be equal to 2.31 g mol⁻¹, T is the midplane temperature, and R is the ideal gas constant.

Two energy sources are considered in our model, namely viscous heating, and the constant irradiation by the local environment of ambient temperature, T_amb = 10 K. Irradiation from the young Sun is neglected because the presence of shadowing is assumed in the outer part of the disk (Ohno & Ueda 2021). This allows the disk temperature to decrease to the condensation temperature of Ar (~20 K), allowing this species to be trapped in Jupiter’s building blocks in a way consistent with the supersolar abundance observed in its envelope (Mousis et al. 2009b, 2012, Mousis et al. 2021b). The temperature profile is computed by summing the energy production rates of both energy sources (Hueso & Guillot 2005): (5)

where σ_sb is the Stefan-Boltzmann constant, while τ_R and τ_P are the Rosseland and Planck optical depth, respectively. Here, we assume τ_P = 2.4τ_R, a case corresponding to the opacity generated by dust grains smaller than 10 μm (Nakamoto & Nakagawa 1994); τ_R is derived from the Rosseland mean opacity, K_R, via the following expression (Hueso & Guillot 2005): (6)

Here, K_R is computed as a sequence of power laws of the form K_R = K₀ρ^aT^b, where ρ denotes the gas density at the mid-plane, and K₀, a, and b are constants that are obtained by fits on observational data in different opacity regimes (Bell & Lin 1994).

The initial state of the model is computed from the self-similar solution derived by Lynden-Bell & Pringle (1974): (7)

By combining Eqs. (7) and (2), and assuming , which corresponds to the case for an early disk (Lynden-Bell & Pringle 1974), the initial profiles of the dust surface density and mass accretion rate are given by: (8)

where r_c is the centrifugal radius and Ṁ_acc,0 is the initial mass accretion rate onto the central star, set to 10^−7.6 M_⊙ yr⁻¹ (Hartmann et al. 1998). The disk mass is related to the surface density profile via the following expression: (9)

where R_min and R_max are the inner and outer bounds of our model. The total disk mass is set to 0.1 M_⊙ and most of it (99%) is encapsulated within ~200 au. The centrifugal radius, r_c, is determined by solving Eqs. (8) with (9) for the chosen values of mass accretion rate and disk mass. Figure 2 represents the thermody-namic profiles of our PSN model assuming α = 10⁻³, and at t =10⁴, 10⁵, and 10⁶ yr of the disk evolution.

Fig. 1 Two outcome scenarios for volatile species explored in this paper. Top panel: the case where volatiles are initially delivered in pure condensate form to the PSN (scenario 1). Bottom panel: the case where volatiles are released in vapor form in the PSN when amorphous grains cross the ACTZ region (scenario 2). Pure condensates, clathrate, and amorphous ice pebbles are represented as blue, brown and red circles, respectively. Vapor is represented as purple dots. The iceline, clathration line, and ACTZ are represented as blue, brown, and red solid lines, respectively. Once delivered to the disk, the phase (solid or gaseous) of each species is determined by the positions of the corresponding condensation, hydration, or clathration lines. Except for the case of CO2, which vaporises at a higher temperature than its clathrate form, hydration, or clathration lines of the volatiles considered are closer to the Sun than their respective icelines. Gaseous volatiles condense or become entrapped (depending on the availability of water ice) when diffusing outward of the locations of their condensation, hydration, or clathration lines. Conversely, volatiles condensed or entrapped in grains or pebbles are released in vapor form when drifting inward of their lines. Peaks of abundances form close to each phase-transition line (see text). Those enrichments are represented by higher solid and vapor concentrations in the panels.

Fig. 2 Profiles of the disk midplane temperature, pressure and surface density calculated at t = 104, 105, and 106 yr as a function of heliocentric distance, assuming α = 10−3, shown from top to bottom.

2.2 Dust dynamics

To determine the size of the dust pebbles, we rely heavily on the two-population algorithm developed by Birnstiel et al. (2012). This algorithm relies on the key idea that the dynamics of dust pebbles of many different sizes can be well approximated by the dynamics of only two populations of particles, in which all particles have the same representative sizes. The first group corresponds to the small population, where grains are of constant size: a₀ = 0.1 μm. The second group represents a large population, where pebbles have a representative size a₁, which depends on the characteristics of the flow.

In the disk, pebbles grow by sticking collisions via the following law: (10)

where τ_growth is the growth timescale, (11)

where Σ_b is the total surface density of solids, and ϵ_g is the dust growth efficiency through mutual sticking set to 0.5 (Lambrechts & Johansen 2014). Then, we compute the Stokes number of pebbles as a function of their sizes (Johansen et al. 2014): (12)

The top and bottom lines of Eq. (12) correspond to the Epstein and Stokes regimes, respectively. The limit between both regimes is fixed by the gas mean free path , computed by equating the two terms in Eq. (12). The term ρ_b is the pebbles mean bulk density: (13)

computed as the average of each species’ bulk density, ρ_b,i, weighted by their solid surface density, Σ_b,i.

Observations indicate that disks are rich in small dust and suggest that fragmentation is a dominant process (Williams & Cieza 2011). Based on this observation, our approach considers fragmentation and radial drift as the growth-limiting mechanisms. These mechanisms set an upper limit on the highest Stokes number that particles can achieve. The first limitation results from the fragmentation occurring when the relative speed between two pebbles due to turbulent motion exceeds the fragmentation velocity, u_f. This upper limit is given by (Birnstiel et al. 2012): (14)

where u_f is set to 10 ms⁻¹ and the factor f_f = 0.37 accounts for the fact that the representative size of the large population is smaller than the biggest size particles can achieve before they fragment.

A second limitation for dust growth is determined by the drift velocities of the different pebbles. When pebbles drift faster than they grow, this sets another upper limit for the Stokes number (Birnstiel et al. 2012): (15)

where P is the disk midplane pressure, v_K the keplerian velocity, and f_d = 0.55 has the same origin as f_f.

When dust grains drift at a high velocity and collide with other particles on their path, they can fragment. This induces a third upper limit for the Stokes number (Birnstiel et al. 2012), given by: (16)

where the factor N = 0.5 accounts for the fact that only larger grains fragment when colliding.

In the algorithm, all limiting Stokes numbers are computed and compared with the Stokes number derived from Eq. (12). At each time step, the smallest Stokes number found in this comparison becomes the reference Stokes number which, in turn, sets the value for the representative size, a₁, of the large population. The representative size of the small population is always a₀, and their Stokes number is always computed in the Epstein regime (top line of Eq. (12)).

Finally, the two-population algorithm of Birnstiel et al. (2012) introduces f_m the fraction of the mass contained in the large population. Among the three size-limiting mechanisms, if particle drift is the most limiting one (St_drift = min (St_frag, St_drift, St_df)), then the fraction of the mass contained in the large population is f_m = 0.97. Otherwise, f_m is set to 0.75 (Birnstiel et al. 2012). The mean grain size ā is then given by: (17)

2.3 Trace species evolution model

Trace species are considered in four distinct forms: vapors, pure condensates, entrapped in clathrates, or forming a monohydrate (case for NH₃ only), and adsorbed in amorphous ice. In our model, a distinct surface density is attributed to each of these forms, with Σ_v,i, Σ_p,i, Σ_c,i, and Σ_a,i corresponding to species i in vapor, pure condensate, clathrate or hydrate, or amorphous ice phases, respectively. Their time and radial evolution is governed by the advection-diffusion equation (Birnstiel et al. 2012; Desch et al. 2017): (18)

where D_i is the diffusion coefficient and is v_i is the radial speed. is a source or sink term that accounts for phase changes, counted positive/negative when some matter is created or lost.

For surface densities of vapors, we assume D_i = D_g and v_i = v_g because vapors are well coupled to the PSN gas and evolve similarly. The gas diffusivity, D_g, is assumed to be equal to the viscosity, v, and the gas velocity is (Shakura & Sunyaev 1973): (19)

At each time and location, we assume that dust particles are formed from a mixture of all available solids. As a consequence, surface densities of solid phases, namely clathrates (and NH₃ monohydrate), amorphous ices, and pure condensates, are evolved with the same diffusion coefficient, D_s, and radial velocity, v_s. For particles of a given size, a, and Stokes number, St, the diffusion coefficient is given by (Birnstiel et al. 2012): (20)

The dust radial velocity is expressed as the sum of gas drag and drift velocities (Birnstiel et al. 2012): (21)

where the drift velocity is (Weidenschilling 1997): (22)

The diffusion coefficient and radial velocity of solids are computed for the small and the large populations, that is, for particles of sizes a₀ and a₁. The diffusion coefficient, D_s, and radial velocity v_s used to evolve surface densities of solids are then given by mass-averaged diffusivities and velocities of the small and large population (Birnstiel et al. 2012): (23)

2.4 Sources and sinks of trace species

We follow the approach of Aguichine et al. (2020) to depict the sources and sinks for both the solid and vapor phases of the different species. A pure condensate of species i undergoes sublimation if its partial pressure is lower than the corresponding equilibrium pressure. Sublimation results in a sink term forpure condensates during the time step, Δt (Drazkowska & Alibert 2017): (24)

where μ_i is the molar mass of species i, is the mean bulk density of grains, P_eq,i is the equilibrium pressure, and ā is the mean size of grains. The second part of the minimum function ensures that no more than the available quantity of pure condensate sublimates. Equilibrium curves of pure condensates are given in Appendix A.

Conversely, a gas of species i forms a pure condensate if its partial pressure is larger than the corresponding equilibrium pressure. The condensation rates results in a source term for pure condensates (Drążkowska & Alibert 2017): (25)

We also added the possibility of NH₃ monohydrate and clathrate crystallization in the PSN. Assuming that enough crystalline water is available, these solids form first during the cooling of the disk because their crystallization temperatures are higher than those of the corresponding pure condensates (see Fig. 3). The only exceptions to that rule are CO₂ and NH₃, namely the only species that condense at a higher temperature than their hydrates at nebular conditions (see Fig. 3). To compute the source and sink terms of trace species in clathrates, we used the same prescription as that used for pure condensates (Eqs. (24) and (25)). The equilibrium pressures of pure condensates are replaced by those of clathrates and NH₃ monohydrate (see Appendix B). The formation of clathrates and NH₃ monohydrate also requires the presence of a minimum amount of crystalline water, resulting in a limit for their source term : (26)

This expression takes the amount of available crystalline water and subtracts the amount of water that is already used to trap currently existing clathrates. In this expression, S_k is the stoichiometric ratio between the species k and water. This ratio is set to 5.75, 5.66, and 1 in the cases of type I clathrate, type II clathrate, and NH₃ monohydrate, respectively. This condition sets an upper limit on the clathrate source term, that can be equal to 0 if all the crystalline water is already holding clathrates. If all conditions are met for a trace species to be able to form pure condensates, clathrates, and monohydrates, we prioritize the solid phase that has the greatest P_i – P_eq,i value. Prioritizing the largest Pi – P_eq,i value is equivalent to prioritizing the highest solidification rate. Such considerations are not taken into account when performing a prioritization test among clathrate formation and condensation of pure condensates, since clathrates only form when pure condensates are not stable.

When grains containing amorphous ice cross the ACTZ, water crystallizes and releases all the adsorbed volatiles irreversibly. The corresponding sink term is then: (27)

with the only condition that must be satisfied is: T > T_ACTZ.

From the various rates of condensation, crystallization and sublimation, we can finally derive the overall sink and source term for the vapor: (28)

Fig. 3 Equilibrium curves of pure condensates (solid lines), clathrates, and NH3 monohydrate (dashed lines) in a pressure-temperature domain relevant to PSN conditions (see the appendix for the relevant data). Two clathrate structures are considered in our model, namely, type I and type II with stoichiometric factors of 5.75 and 5.66, respectively. Partial pressures are calculated by considering the species abundances given in Table 1.

3 Results

Simulations have been performed in the case of two distinct scenarios. In scenario 1, particles initially released in the PSN are only made from pure condensates. During their inward drift, they can sublimate, condense again and/or form various clathrates and a NH₃ monohydrate, depending on the local temperature and pressure conditions of the disk. In scenario 2, particles initially released in the PSN are only made from amorphous ice. During their inward drift, water contained in these particles transitions to a crystalline structure when crossing the ACTZ and leads to the release of adsorbed volatiles in the gas phase. The released vapors can in turn condense into pure ices and/or form various clathrates and NH₃ monohydrate, following the prescription described in the previous section. Both scenarios are explored with the assumption of 0.1 M_⊙ for the mass of the PSN and a viscosity parameter of α = 10⁻³.

The initial surface density of a species, i, is: (29)

where x_0,i is the initial mass fraction of the species, i. At t = 0, the partial pressures of the different species are computed at each point of the grid. In scenario 1, if the partial pressure of a given species is below its equilibrium pressure, then the corresponding location is filled with vapor. Otherwise, this location filled with pure ice. In scenario 2, species are in the vapor phase where T ≥ T_ACTZ and in an amorphous ice phrase otherwise. In both scenarios, thermodynamic equilibrium is established after a few time steps (~1 yr).

The initial PSN composition is derived from the protosolar elemental abundances tabulated by Lodders et al. (2009). We assume that all C is distributed between CO, CO₂, or CH₄, with the remaining O forming H₂O. We have set CO:CO₂:CH₄ = 10:4:1 in the PSN gas phase. The CO:CO₂ ratio is derived from ROSINA measurements of comet 67P/C-G between 2014 August and 2016 September (Mousis et al. 2014). The CO:CH₄ ratio is consistent with the production rates measured in the southern hemisphere of the 67P/C-G in October 2014 by the ROSINA instrument (Le Roy et al. 2015). Sulfur is assumed to be half in H₂S form and half in refractory sulfide components (Pasek et al. 2005). We also assumed N₂:NH₃ = 1:1, a value predicted by thermochemical models that take into account catalytic effects of Fe grains on the kinetics of N₂ to NH₃ conversion of the PSN (Fegley 2000; Mousis et al. 2009a). The molar abundances of the different species are derived from the gas phase abundances given in Table 1.

Figure 4 represents the radial evolution of the water mass abundance in the PSN at different epochs of its evolution in scenario 1 and scenario 2 (top row and bottom row, respectively). Both cases produce very similar water abundances and show a water enrichment peak that is about 10 times its initial abundance at the location of the snowline (~2.8 au). However, the availability of crystalline ice is much more limited in scenario 2 due to H₂O being mostly in an amorphous state. Indeed, in this case, crystalline ice is present only in the ~2.8−4 au region.

3.1 Scenario 1: initial delivery of pure condensates

Figure 5 represents the time and radial evolution of the abundance ratios (with respect to initial abundances) of the different species existing in various phases in the case of scenario 1, respectively. After 10 kyr of PSN evolution, the species considered are present under all their possible forms - except for CO₂ and NH₃, which are never enclathrated. By successive order with progressing heliocentric distance (and decreasing temperature and pressure conditions), we first find the different vapors, then narrow regions corresponding to the presence of clathrates, and, finally, an outer region that is only populated with pure condensates. Regions where clathrates exist expand from 7 to 12 au, depending on the species considered. At this early stage of the PSN evolution, there is no significant enrichment that can be observed for any species.

After 0.1 Myr of disk evolution, Fig. 5 also shows that beyond 6 au, clathrates coexist with pure condensates and their abundances decrease steeply - except for CO₂ and NH₃, which that exist only as pure ices. A few au beyond that, pure condensates become the only solid structures existing in the outer PSN. Depending on the PSN ther modynamic conditions and the availability of crystalline water, when the different species become fully enclathrated, their abundances form unique enrichment peaks located at their clathration lines, reaching up to ~15 times their initial values. On the other hand, if the budget of crystalline water is not high enough, then the species are only partly enclathrated. They then form two narrow enrichment peaks at their condensation and clathration lines, reaching ~5 and ~15 times their initial abundances, respectively. Table 2 displays the heliocentric distance and the value of the enrichment peak (relative to the initial abundance) for each species under consideration at t = 0.1 Myr of the PSN evolution. Depending on the species considered, the enrichment peaks range between 2 and 18 times the protosolar values and are located in the 2–11 au region. The closest and furthest peaks from the Sun are those of the water snow line at 2.8 au and N₂ iceline at 10.8 au, respectively. Two enrichment peaks, located at the clathration and icelines, are found in the cases of CH₄ and Kr.

In Fig. 5, all species, except CO, N₂, and Ar, exhibit a dip in the surface densities of the pebbles of pure condensates around 10.5 au. In this region, the total surface density of icy pebbles is increased due to the combined actions of CO, N₂, and Ar icelines located at 10.4, 10.9, and 10.7 au, respectively. The condensations rates of CO, N₂, and Ar locally increase the surface density of solids, which are in excess compared with their loss rates via inward drift. On the other hand, because the icelines of the other species are located closer to the Sun, their surface densities progressively decrease in this region, as a result of the inward drift of particles.

After 1 Myr of PSN evolution, all peaks are smoothed out in the gas phase. Because of the inward drift of pebbles, most of the species are in vapor form when the solids start to deplete. The total mass of solids has decreased by a factor of four compared to the beginning of the simulation. This effect is more significant in the case of clathrates forming at very low temperatures (less than 50 K), since the fraction of crystalline water used to form higher temperature clathrates increases with time.

Fig. 4 Time evolution of the mass abundance of water, defined as the radial profile of , for both scenarios at t = 104, 105, and 106 yr. Solid lines represent H2O in the gaseous phase (orange line), crystalline phase (blue line), and in amorphous phase (red line).

Fig. 5 Time and radial evolution of species’ mass abundances normalized to their initial values in gaseous phase (orange line), pure condensate form (blue line), and clathrate (green line), at t = 104, 105, and 106 yr in scenario 1.

3.2 Scenario 2: initial delivery of amorphous ices

Figure 6 represents the time and radial evolutions of the abundance ratios (with respect to the initial abundances) of the different species existing in various phases in the case of scenario 2. After 10 kyr of PSN evolution, with progressing heliocentric distance, first the different vapors are found, then an outer region that is only populated with amorphous ice. The budget of solid phase is dominated by species trapped in amorphous ice. Pure condensates only form when the vapors desorb from amorphous ice at the ACTZ (shown at ~5 au), diffuse outward and cross an iceline. Such a process is efficient only if the condensation line is close to the ACTZ, leading to very narrow regions with the presence of pure condensates only in the case of NH₃.

After 0.1 Myr of disk evolution, it is notable that only pure condensates of NH₃, H₂S, CO₂, and PH₃ form. Table 3 displays the heliocentric distance and the value of the enrichment peak (relative to the initial abundance) for each species under consideration at this epoch of the PSN evolution. All peaks are located in a much narrower region, centered at the location of the ACTZ, compared with the scenario 1, with values ranging from 7 to 10 times the initial abundances. The peak locations are influenced by the presence of narrow (less than 1 au) regions filled with pure condensates. In those regions, corresponding to the icelines, the abundance of pure condensates exceed that of amorphous ice, thus influencing the location of the enrichment peaks.

After 1 Myr of PSN evolution, as a result of pebble drift, all peaks have been smoothed and the volatiles trapped in amorphous water ice phase are strongly depleted by factors reaching more than 100, compared with their gaseous abundances in the inner disk. The amount of pure condensates exceeds that of amorphous ice in some narrow regions, except in the cases of H₂S and PH₃. For reasons identical to those invoked at the same epoch of PSN evolution in scenario 1, the surface density of solids is strongly decreased in the region centered at ~7 au.

4 Discussion

In this section, we discuss the implications of our model for various bodies of the solar system. We first investigate the sensitivity of our model to the variation of its input parameters and then provide fits of the volatile composition of comet R2 and those of Uranus and Neptune.

4.1 Sensitivity to parameters

The stability of our results has been tested against the variations of the pebble density, disk's mass, and the viscosity parameter in the 0.1–1 g cm⁻², 10⁻²–10⁻² M_⊙, and 10⁻⁴–10⁻² ranges, respectively. Variation in pebble density leads to results similar to those presented in Sects. 3.1 and 3.2. Lower density pebbles drift over shorter timescales at given size, but are also smaller because of a lower fragmentation limit (see Sect. 2.2). On the other hand, smaller pebbles drift over longer timescales, implying that both effects are (roughly) mutually counterbalanced and produce only minor variations in the abundance profiles. Both the mass and a viscosity parameter of the disk affect its viscous evolution and, thus, the locations of the various icelines. Less massive disks are cooler because of their reduced viscous dissipation (see Eq. (5)) and display their icelines closer to the host star. For example, we find that the condensation and clathration lines of the species investigated in our study are ~2 au closer to the Sun in a 10⁻² M_⊙ disk, compared with a 10⁻¹ M_⊙ disk. As another example, the enrichment peak associated with the CO iceline ranges between 7 and 11 au from the Sun when the α-value is varied between 10⁻⁴ and 10⁻². The magnitude of the abundance peaks is also affected by the variation of a. The abundance of the CO peak ranges between about 10 and 5 times, respectively, its initial PSN abundance when the a value is varied between 10⁻⁴ and 10⁻². In our disk model, the contribution from the irradiation by the Sun to the disk's midplane temperature is not considered because it is assumed that the outer PSN is shadowed by its inner region (Ohno & Ueda 2021). This allows the disk to reach temperatures low enough to enable the condensation of ultravolatiles (CO, N₂, Ar, etc.) at the current locations of the giant planets, assuming the absence of migration during formation. When this contribution is included by considering the formalism depicted in Adams et al. (1988) and Ruden & Pollack (1991), the temperature profile becomes slightly warmer, implying that the icelines and clathra-tion lines are moved outward by 1–2 au. Despite these changes, the general trend of our results is not impacted.

Fig. 6 Time and radial evolution of species' mass abundances normalized to their initial values in gaseous phase (orange line), pure condensate form (blue line), and amorphous form (red line), compared with their initial mass abundances, at t = 104, 105, and 106 yr in scenario 2.

Fig. 7 N2/CO (top panel) and CO/H2O (bottom panel) abundance ratios in pebbles represented as a function of heliocentric distance in the case of scenario 1, and compared with those measured in R2’s coma (blue bar) at 0.1, 0.5, and 1 Myr of the PSN evolution. Both ratios measured in R2 are simultaneously matched by our model at a distance of 7 au and 1 Myr.

4.2 Implications for comet R2

R2 is a long-period comet displaying an unusually high N₂/CO ratio of 0.006−0.008 (Biver et al. 2018; Opitom et al. 2019). Another peculiar characteristic of this comet is its heavy depletion in terms of H₂O, with a CO/H₂O ratio of about 312 (McKay et al. 2019).

Figures 7 and 8 represent the radial profiles of the CO/H₂O and N₂/CO ratios calculated in the PSN pebbles with our nominal model as a function of time in the cases of scenario 1 and scenario 2, respectively. Both ratios are compared with those measured in R2’s coma (blue horizontal bar). Assuming that R2 formed from a unique set of building blocks, we require both ratios to be reproduced by our model at the same heliocentric distance and epoch. Simulations performed in the case of scenario 1 reproduce both ratios at a heliocentric distance of 7 au, after 1 Myr of PSN evolution. Our result is then consistent with those derived from the study of Mousis et al. (2021b), based on a simple approach that does not consider the interplay between the clathrate and water ice reservoirs.

Simulations performed in the case of scenario 2 reproduce the N₂/CO ratio at 8 au after 1 Myr of PSN evolution. However, the CO/H₂O ratio is not matched by our model, even if a peak is seen at 8 au and 1 Myr. We should note that the positions and magnitudes of those peaks can change when the α-parameter and mass of our disk model are varied. This implies that, even if scenario 1 provides a better match of R2’s composition with our nominal model, scenario 2 cannot be excluded. Interestingly, our calculated peaks are within a zone of dynamic instability in the early Solar System, which is more likely to result in ejection of planetesimals than capture by the Oort Cloud. This is a possible explanation for the lack of R2-like comets observed today (Anderson et al. 2022).

Fig. 8 N2/CO (top panel) and CO/H2O (bottom panel) abundance ratios in pebbles represented as a function of heliocentric distance in the case of scenario 2, and compared with those measured in R2’s coma (blue bar) at 0.1, 0.5 and 1 Myr of the PSN evolution. Only the N2/CO ratio is reproduced in R2 at 8 au and 1 Myr. The CO/H2O ratio is however approached by our model at the same location and epoch of the PSN evolution.

4.3 Implications for the composition of Jupiter

One-σ error bar measurements made at Jupiter by the Galileo probe and the Juno spacecraft indicate C, N, O, S, P, Ar, Kr, and Xe abundances that are ~1.5 to 6 times higher than the protosolar values (Atreya et al. 2003; Wong et al. 2004; Mousis et al. 2018; Li et al. 2020). To explain those features, it has been proposed that Jupiter’s atmosphere could reflect the composition of icy planetesimals either made of amorphous ice (Owen et al. 1999) or from pure condensates and/or clathrates (Gautier et al. 2001; Gautier & Hersant 2005; Mousis et al. 2018, 2021b). Alternatively, it has been proposed that this supersolar metallicity could result from the accretion of already pre-enriched PSN gas (Mousis et al. 2019; Aguichine et al. 2022).

Figure 9 represents the time evolution of the sum of the elemental enrichments calculated in vapor and solid phases at the heliocetric distance of 4 au, compared with their protosolar values, and in the cases of our two scenarios. Following the approach of Aguichine et al. (2022), we focused on the composition of the PSN at 4 au, chosen as the location of Jupiter’s formation. This distance is, in the model, beyond the water ice-line but inward of the icelines of all other trace species. The dust-to-gas ratio can easily become greater than 2 to 3 times the protosolar composition in this region and could ease the formation of a proto-Jupiter core via the streaming instability (Yang et al. 2017).

Figure 9 shows that the measured elemental enrichments are all matched by our model after 0.8–1 Myr and 50–100 kyr of the PSN evolution in scenario 1 and scenario 2, respectively. Our models suggest that the consideration of clathrate formation in addition to the crystallization of pure condensates in the PSN (scenario 1) still allow for the formation of a Jupiter-like planet from supersolar gases originating from the disk, compared with models considering the crystallization of pure condensates only, such as that developed by Aguichine et al. (2022). Our model also suggests that the presence of multiple condensation and clathration lines does not alter the formation of supersolar vapors subsequent to their release from amorphous ice (scenario 2). Previous works exploring this possibility did not consider the formation of icelines in their models (Monga & Desch 2015; Mousis et al. 2019). In particular, Saturn’s tropospheric abundances of C, N, S, and P have been measured to be between 3 and 13 times their protosolar abundances (Atreya et al. 2018), suggesting that its metallicity is higher than that of Jupiter. Assuming that Saturn formed in the vicinity of the CO₂ iceline to account for its high carbon enrichment, which corresponds to a heliocentric distance of ~6 au at early epochs in our PSN model, our calculations show that this range of enrichments is reproduced within 0.2–0.3 Myr in scenario 1 and 0.2 Myr in scenario 2. This implies that Saturn could have formed earlier than Jupiter in scenario 1, whereas in scenario 2, it could have formed later. An earlier formation of Saturn could have reduced the inward flux of pebble and vapors, lowering the value of elemental enrichment peaks at Jupiter’s location. However, the height of enrichment peaks is highly sensitive to the disk parameters (e.g., the α-value; see Aguichine et al. 2022). Therefore, if the accretion of pebbles by Saturn reduces the height of the enrichment peaks at Jupiter's location, it is still possible to achieve a metallicity similar to what is measured in its atmosphere.

Fig. 9 Time evolution of the elemental abundances (relatives to the pro-tosolar values) at a heliocentric distance of 4 au in the cases of scenario 1 (top panel) and scenario 2 (bottom panel). The blue area corresponds to the range covered by the elemental abundances derived from spacecraft measurements (see text).

4.4 Implications for the composition of Ice Giants

The composition of the deep atmospheres of Uranus and Neptune is shrouded in mystery since most of the heavy constituents condense at pressures deeper than may readily be probed remotely (Mousis et al. 2020). The only determinations that can be used so far in our model are the C/N and C/S ratios, which have been found equal to or higher than ~175 and ~35, respectively (Asplund et al. 2009; Karkoschka & Tomasko 2009, 2011; Irwin et al. 2018, 2019a,b). Figures 10 and 11 represent the radial profiles of the C/N and C/S elemental ratios in pebbles as a function of time in the PSN, compared with the minimum values measured in the tropospheres of the two ice giants. In both cases, these ratios can be reproduced in the ~7–8 au region after 1 Myr of PSN evolution. The formation distance of solids in the PSN is model-dependent, but our simulations suggest that the giant planets did form in a more compact configuration than the current one - which is also in agreement with several dynamical models (Tsiganis et al. 2005; Lykawka et al. 2010; Guilera et al. 2011). Given their high metallicities, the formation of Uranus and Neptune requires a higher surface density of solids than the values derived here for the outer PSN. This assumption is at odds with the fact that the planetesimal density and collision probability are both low in the outer disk. One way to overcome this difficulty would be to assume the formation of the two giants via the accretion of pebbles directly onto the planetary embryo, which can still work efficiently far from the host star (Helled et al. 2014; Bitsch et al. 2015, 2018; Armitage 2019).

Fig. 10 Radial profiles of the C/N (top panel) and C/S (bottom panel) ratios calculated in pebbles at different epochs of the PSN evolution in the case of scenario 1. The horizontal dashed line represents the minimum ratio measured in the tropospheres of Uranus and Neptune. The orange area encompasses the current locations of the ice giants in the solar system.

4.5 Limitations of the model

Over the recent years, substructures have been largely observed in protoplanetary disks. The ALMA/DSHARP survey showed that protoplanetary disks are not smooth and that substructures are ubiquitous (Andrews et al. 2020; Jennings et al. 2022). Substructures can be produces by planet-disk interactions, in the form of spiral density waves (Muto et al. 2012; Zhang et al. 2021) and gaps (Bae et al. 2017). Radiative hydrodynamic models show that these substructures can be also formed by instabilities in the disk (Lovelace et al. 1999; Lovelace & Romanova 2014; Blanco et al. 2021). Such features translate into local density and pressure variations that act as dust traps which locally enhance the dust surface density. Although our model does not take into account such disk substructures, we show that the icelines, clathration lines, and the ACTZ correspond to vapor and pebble enrichment peaks, which can lead to instabilities.

In this work, giant planets are assumed to be formed in situ. In scenario 1, Jupiter is formed within 0.8–1 Myr. This timescale is compatible with a formation via pebble accretion (Bitsch et al. 2015, 2018; Alibert et al. 2018; Armitage 2019; Venturini & Helled 2020) and gravitational instability (Boss 1997; Zhu et al. 2012; Kratter & Lodato 2016). On the other hand, in scenario 2, Jupiter forms in less than 0.1 Myr. This timescale is consistent with a gravitational instability which can be triggered as early as 0.1 Myr (Zhu et al. 2012). An early formation of Jupiter allows us to account for the observed carbonaceous chondrites dichotomy (Kleine et al. 2020). Although, the disk instability model is consistent with formation timescales found in both scenario 1 and scenario 2, it is important to note that observations and models suggest that, at minimum, amorphous ice should be heated up to its crystallization temperature when falling from the presolar cloud onto the PSN (Visser et al. 2013). Those findings suggest that scenario 1 is the most likely scenario, implying that our results remain consistent with the core accretion model.

One limitation of our model is the fact that planet migration is not considered. During or after formation, planets migrate inward or outward (Masset & Papaloizou 2003; Bitsch et al. 2015; Schneider & Bitsch 2021). During migration, planets accrete material from different parts of the disk with various compositions. Although planetary migration contradicts an in situ formation hypothesis, 3D hydrodynamical simulations indicate that the rates of type I and type II migrations could be several times slower than the prescription usually employed in the literature of planet formation (Chrenko & Nesvorný 2020; Lega et al. 2021; Chametla & Chrenko 2022). Assuming an in situ formation of Jupiter in our model to reproduce its observed metallicity is expected to remain valid in light of these revised migration rates because our derived timescales are still quite short, compared to the PSN evolution.

Fig. 11 Radial profiles of the C/N (top panel) and C/S (bottom panel) ratios calculated in pebbles at different epochs of the PSN evolution in the case of scenario 2. The horizontal dashed line represents the minimum ratio measured in the tropospheres of Uranus and Neptune. The orange area encompasses the current locations of the ice giants in the Solar System.

5 Conclusion

In this work, we investigate the impact of clathrate formation on the radial distribution of volatiles in the PSN, considering two scenarios, each of them corresponding to a distinct initial condition. To do so, we used a 1D protoplanetary disk model coupled with modules describing the evolution of trace species in the vapor phase, as well as the dynamics of dust and pebbles. This model also considers the different sources and sinks for the volatile phases considered (vapors, pure condensates, or clathrates).

In scenario 1, we assume that the volatiles were delivered to the PSN in the form of pure condensate grains. In this case, we show that clathrates can crystallize and form enrichment peaks at about 10 times the value of the initial abundances at their clathration lines, which are closer to the Sun than their corresponding icelines. The amount of clathrates formed in the PSN depends on the local abundance of crystalline water, which in many cases, acts as a limiting factor in our model. In scenario 2, we assumed that the volatiles were delivered to the PSN in the form of amorphous grains. Under those conditions, volatiles are only released from amorphous ice when the icy grains are heated up to ~135 K, namely, at the ACTZ location. An enrichment peak up to about seven times the initial abundances then forms at the ACTZ location. In this case, clathrate formation is not possible because there is no crystalline water ice available beyond the ACTZ in the PSN. All the enrichment peaks of pure condensates are also located close to the ACTZ. Our investigation shows that both scenario 1 and scenario 2 can reproduce the known compositions of comet R2, Jupiter, Uranus, and Neptune. This implies that our model does not allow us to formally rule out the presence of amorphous ice during the early phases of the PSN, assuming that those planetary bodies accreted from pebbles or pebble-made planetesimals. More planetary composition data, such as the in situ measurements of Saturn’s atmosphere, are needed to understand the initial formation conditions of the PSN.

Acknowledgements

The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille Université - A*Midex, a French "Investissements d’Avenir programme" AMX-21-IET-018. This research holds as part of the project FACOM (ANR-22-CE49-0005-01_ACT) and has benefited from a funding provided by l’Agence Nationale de la Recherche (ANR) under the Generic Call for Proposals 2022. O.M. acknowledges support from CNES. J.I.L. was supported by the Juno project. We acknowledge Sarah E. Anderson for helpful discussions about her dynamical simulations regarding the origin of comet R2. We thank the anonymous referee for their useful comments that improved the quality of this paper.

Appendix A Vapor pressures of pure condensates

The vapor pressures of the pure condensates considered in this work, except H₂O and NH₃, follow the law: (A.1)

where the a_k factors have been determined experimentally (Fray & Schmitt 2009) and are summarized in Table A.1.

The vapor pressures of H₂O is given by (Wagner et al. 2011): (A.2)

where P_tp and T_tp are the pressure and the temperature of the triple point of water, respectively.

The phosphine equilibrium pressure is given by Stull (1947): (A.3)

Appendix B Dissociation pressures of NH₃ monohydrate and clathrates

The dissociation pressures of NH₃ monohydrate and clathrates follow the Antoine law: (B.1)

with A and B parameters determined from experiments and summarized in Table B.1.

References

All Tables

All Figures

Fig. 1 Two outcome scenarios for volatile species explored in this paper. Top panel: the case where volatiles are initially delivered in pure condensate form to the PSN (scenario 1). Bottom panel: the case where volatiles are released in vapor form in the PSN when amorphous grains cross the ACTZ region (scenario 2). Pure condensates, clathrate, and amorphous ice pebbles are represented as blue, brown and red circles, respectively. Vapor is represented as purple dots. The iceline, clathration line, and ACTZ are represented as blue, brown, and red solid lines, respectively. Once delivered to the disk, the phase (solid or gaseous) of each species is determined by the positions of the corresponding condensation, hydration, or clathration lines. Except for the case of CO2, which vaporises at a higher temperature than its clathrate form, hydration, or clathration lines of the volatiles considered are closer to the Sun than their respective icelines. Gaseous volatiles condense or become entrapped (depending on the availability of water ice) when diffusing outward of the locations of their condensation, hydration, or clathration lines. Conversely, volatiles condensed or entrapped in grains or pebbles are released in vapor form when drifting inward of their lines. Peaks of abundances form close to each phase-transition line (see text). Those enrichments are represented by higher solid and vapor concentrations in the panels.

In the text

	Fig. 2 Profiles of the disk midplane temperature, pressure and surface density calculated at t = 104, 105, and 106 yr as a function of heliocentric distance, assuming α = 10−3, shown from top to bottom.
In the text

Fig. 3 Equilibrium curves of pure condensates (solid lines), clathrates, and NH3 monohydrate (dashed lines) in a pressure-temperature domain relevant to PSN conditions (see the appendix for the relevant data). Two clathrate structures are considered in our model, namely, type I and type II with stoichiometric factors of 5.75 and 5.66, respectively. Partial pressures are calculated by considering the species abundances given in Table 1.

In the text

	Fig. 4 Time evolution of the mass abundance of water, defined as the radial profile of , for both scenarios at t = 104, 105, and 106 yr. Solid lines represent H2O in the gaseous phase (orange line), crystalline phase (blue line), and in amorphous phase (red line).
In the text

	Fig. 5 Time and radial evolution of species’ mass abundances normalized to their initial values in gaseous phase (orange line), pure condensate form (blue line), and clathrate (green line), at t = 104, 105, and 106 yr in scenario 1.
In the text

	Fig. 6 Time and radial evolution of species' mass abundances normalized to their initial values in gaseous phase (orange line), pure condensate form (blue line), and amorphous form (red line), compared with their initial mass abundances, at t = 104, 105, and 106 yr in scenario 2.
In the text

	Fig. 7 N2/CO (top panel) and CO/H2O (bottom panel) abundance ratios in pebbles represented as a function of heliocentric distance in the case of scenario 1, and compared with those measured in R2’s coma (blue bar) at 0.1, 0.5, and 1 Myr of the PSN evolution. Both ratios measured in R2 are simultaneously matched by our model at a distance of 7 au and 1 Myr.
In the text

Fig. 8 N2/CO (top panel) and CO/H2O (bottom panel) abundance ratios in pebbles represented as a function of heliocentric distance in the case of scenario 2, and compared with those measured in R2’s coma (blue bar) at 0.1, 0.5 and 1 Myr of the PSN evolution. Only the N2/CO ratio is reproduced in R2 at 8 au and 1 Myr. The CO/H2O ratio is however approached by our model at the same location and epoch of the PSN evolution.

In the text

	Fig. 9 Time evolution of the elemental abundances (relatives to the pro-tosolar values) at a heliocentric distance of 4 au in the cases of scenario 1 (top panel) and scenario 2 (bottom panel). The blue area corresponds to the range covered by the elemental abundances derived from spacecraft measurements (see text).
In the text

	Fig. 10 Radial profiles of the C/N (top panel) and C/S (bottom panel) ratios calculated in pebbles at different epochs of the PSN evolution in the case of scenario 1. The horizontal dashed line represents the minimum ratio measured in the tropospheres of Uranus and Neptune. The orange area encompasses the current locations of the ice giants in the solar system.
In the text

	Fig. 11 Radial profiles of the C/N (top panel) and C/S (bottom panel) ratios calculated in pebbles at different epochs of the PSN evolution in the case of scenario 2. The horizontal dashed line represents the minimum ratio measured in the tropospheres of Uranus and Neptune. The orange area encompasses the current locations of the ice giants in the Solar System.
In the text