The role of ice lines in the formation of Uranus and Neptune

We aim at investigating whether the chemical composition of the outer region of the protosolar nebula can be consistent with current estimates of the elemental abundances in the ice giants. To do so, we use a self-consistent evolutionary disc and transport model to investigate the time and radial distributions of H2O, CO, CO2, CH3OH, CH4, N2 and H2S, i.e. the main O-, C-, N and S-bearing volatiles in the outer disc. We show that it is impossible to accrete a mixture composed of gas and solids from the disc with a C/H ratio presenting enrichments comparable to the measurements (approx. 70 times protosolar). We also find that the C/N and C/S ratios measured in Uranus and Neptune are compatible with those acquired by building blocks agglomerated from solids condensed in the 10–20 arb. units region of the protosolar nebula. By contrast, the presence of protosolar C/N and C/S ratios in Uranus and Neptune would imply that their building blocks agglomerated from particles condensed at larger heliocentric distances. Our study outlines the importance of measuring the elemental abundances in the ice giant atmospheres, as they can be used to trace the planetary formation location, the origin of their building blocks and/or the chemical and physical conditions of the protosolar nebula. This article is part of a discussion meeting issue ‘Future exploration of ice giant systems’.


Introduction
Uranus and Neptune are the outermost giant planets of our solar system. The apparent size of these two planets in the sky is roughly a factor of 10 smaller than Jupiter and Saturn, making their physical properties much more difficult to characterize. The bulk compositions of Uranus and Neptune are poorly known [1]. Interior models suggest that these planets contain approximately 10-25% of H-He by mass, where the exact number depends on the assumed composition of the heavy elements [2]. Therefore, although there are estimates for the total heavy-element mass, the available data are insufficient to discriminate among different elemental compositions, implying that the water-to-rock ratio remains unknown. Often models assume a large region of the planet presenting high volatiles enrichments [3][4][5][6], as shown by the measurement of the C abundance, found to be enriched approximately 70 times its protosolar value in both Uranus and Neptune (table 1). By contrast, recent measurements suggest that N/H and S/H ratios might be instead subsolar in their envelopes [11,12]. However, the potential presence of hidden reservoirs of nitrogen and sulfur at deeper levels [13,14] would imply that volatiles are not homogeneously mixed in the envelopes. At present, it remains unclear whether the measured C abundance represents the planetary bulk. Often it is assumed that this is the case, which corresponds to fully mixed (convective) planetary interiors. However, updated formation, evolution and structure models of the ice giants suggest that both Uranus and Neptune are unlikely to be fully mixed [1]. Nevertheless, it remains useful to use the measured abundances and link them to possible formation locations.
To explain the apparent discrepancy observed between the C and N abundances in their envelopes, it has been proposed that Uranus and Neptune have grown from N-depleted building blocks at the location of the carbon monoxide ice line [15]. However, this scenario did not consider the possibility that the measured N abundance could be simply a lower limit in Uranus and Neptune, a hypothesis which has been proposed since then [13,16]. The model also overestimated the redistribution of volatiles around their ice lines for two reasons. First, calculations of volatile transport were based on a stationary PSN model while the thermodynamic conditions in discs evolve significantly over the first hundred thousand years. Second, instantaneous condensation of vapours was assumed, instead of using condensation rates, thus leading to an overly efficient production of solids in the vicinity of the ice lines.
In this paper, we investigate whether the chemical composition of the outer region of the PSN is consistent with current estimates of the elemental abundances in the ice giants by overcoming the afore-mentioned issues. We use a self-consistent evolutionary disc and transport model to investigate the time and radial distributions of H 2 O, CO, CO 2 , CH 3 OH, CH 4 , N 2 and H 2 S, i.e. the main O-, C-, N, and S-bearing volatiles in the outer PSN [17]. The calculated abundance profiles in the outer part of the disc are then compared with the abundance data available for the two planets. Section 2 is dedicated to a short description of the existing abundance data of species of interest in the atmospheres of the ice giants. The disc and transport model used is summarized in §3. Section 4 presents the calculated abundance profiles in the outer PSN and the comparisons with the compositions of Uranus and Neptune. Our conclusion is presented in §5.

Abundances of heavy elements in Uranus and Neptune's atmospheres
The composition of the deep atmospheres of Uranus and Neptune is shrouded in mystery since most heavy constituents condense at pressures deeper than can readily be probed remotely. Water is expected to condense at very deep pressures ( 100 bar) and ammonia and hydrogen sulfide are expected to combine together to form a cloud of either NH 4 SH or (NH 4 ) 2 S at a pressure of 40-50 bar, leaving the more abundant gas (i.e. either NH 3 or H 2 S) to condense at lower pressures (less than 10 bar). Finally, methane is predicted to condense at a pressure of 1-2 bar (depending on abundance). For Jupiter and Saturn, the abundance of NH 3 [11,12] and assuming S/N > 4.5-5 × solar.
spectrum of Uranus and Neptune showed it to be much brighter than expected, suggesting that the abundance of NH 3 (which is a strong microwave absorber) was depleted [18][19][20], leading to the suggestion that the abundance of H 2 S exceeded that of NH 3 at 40-50 bar, and that H 2 S should be present alone at lower pressures and condense at about three bars. This presence of H 2 S gas absorption features above the clouds of both Uranus and Neptune has recently been detected in Gemini/NIFS observations [11,12], which support this hypothesis. One explanation for the apparently low abundance of ammonia in the atmospheres of Uranus and Neptune is that it is partially dissolved in an aqueous ammonia cloud. However, this is not predicted to be able to absorb sufficient ammonia and it has been suggested that ammonia (and hydrogen sulfide) may instead partially dissolve in a water 'ocean' or even in an 'ionic ocean' at depth (e.g. [13]), modifying the apparent deep abundance of nitrogen and sulfur in the atmospheres of Uranus and Neptune. Alternatively, if we believe the observed abundances of H 2 S and NH 3 to be representative of the bulk composition of these worlds, then it suggests that Uranus and Neptune may have formed in a colder part of the solar nebula than Jupiter and Saturn, in a region where more S than N was accreted into the planets [11,12]. Finally, a potential clue to the deep abundances of Uranus and Neptune comes from the observations of carbon monoxide, which resides mostly in the stratosphere, but which may perhaps have also been detected in the upper troposphere. On Uranus, CO estimates suggest abundances of 7.1-9.0 ppb in the stratosphere [21], and a 3-σ upper limit of 2.1 ppb in the troposphere for pressures 0.1-0.2 bar [22]. The abundance of CO in Neptune's atmosphere is even found to be much higher with [23] determining a step-type profile with 1-2 ppm in the stratosphere and 0-0.3 ppm in the troposphere at pressures >1 bar. The presence of CO in the stratosphere can be explained through the impact of comets or via interplanetary dust particles, but models that can predict such high abundances in the troposphere need to have atmospheric compositions that are heavily enriched in O/H by at least 280 times solar [21,23,24]. However, such a large enrichment of O/H is not compatible with D/H measurements, which suggest more modest O/H enrichments of approximately 50-150 [25] if Neptune's internal water was sourced from protoplanetary ices with D/H comparable to present day comets. These considerations point toward how challenging it is to derive the O/H ratio from CO [13,26]. Given those difficulties, the present work only takes into account the measurements of C, N and S in the envelopes of Uranus and Neptune (table 1).

Volatile distribution model
The volatile transport and distribution model used in our work is derived from the approach described in [27], to which the reader is referred for details. In our model, the accretion rate onto the star declines over time, and the disc develops a transition radius between inward and outward gas flows. The location of this transition radius is moving outward with time. Our model is governed by the following differential equation [27]: which describes the time evolution of a viscous accretion disc of surface density Σ g and viscosity ν, assuming hydrostatic equilibrium in the z 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Ṁ: where Q = d ln(Σ g ν)/d ln(r). The first equation is a mass conservation law, and the second one is a diffusion equation. The mass accretion rate can be expressed in terms of the gas velocity field v g asṀ = −2π rΣ g v g .
The viscosity ν is computed in the framework of the α-formalism [28]: where Ω K = GM /r 3 is the keplerian frequency with G the gravitational constant, and c s is the isothermal sound speed given by In this expression, R is the ideal gas constant and μ g is the gas mean molar mass. In equation (3.4), α is a non-dimensional parameter measuring the turbulence strength, which also determines the efficiency of viscous heating, hence the temperature of the disc. The midplane temperature T of the disc is computed via the addition of all heating sources, giving the expression [29]: (3.6) The first term corresponds to the viscous heating [30], where σ sb is the Stefan-Boltzmann constant, and τ R and τ P are the Rosseland and Planck mean optical depths, respectively. For dust grains, we assume τ P = 2.4τ R [30]. τ R is expressed as [29]: where κ R is the Rosseland mean opacity, computed as a sequence of power laws of the form κ R = κ 0 ρ a T b . Parameters κ 0 , a and b are fitted to experimental data in different opacity regimes [31] and ρ denotes the gas density at the midplane. The second term corresponds to the irradiation of the disc by the central star of radius R and surface temperature T . It considers both direct irradiation at the midplane level and irradiation at the surface at a scale height H = c s /Ω K . The last term accounts for background radiation of temperature T amb = 10 K. At each time step, Σ g is evolved with respect to equation (3.2), using the forward Euler method. Then equation (3.6) is solved iteratively with equations (3.4), (3.5) and (3.7) to produce the new thermodynamic properties of the disc. The new velocity field is then computed following equation (3.3). The evolution of the disc is started with an initial profile given by Σ g ν ∝ exp(−r 2−p ), with p = 3/2 for an early disc [27]. In our computations, α = 5 × 10 −3 and the disc mass is set equal to 0.1 M . Most of the disc mass (99%) is encapsulated within approximately 200 arb. units and the initial mass accretion rate onto the Sun is set to 10 −7.6 M yr −1 [32]. All quantities are defined on 1000 logarithmically scaled bins. Derivatives are calculated using central differences, except for boundaries where we assume ∂Ṁ/∂r = 0. To test the stability of the numerical scheme, we compute the disc's mass by integrating Σ g over its radius and the mass leaving each boundary. This quantity remains equal to the disc's initial mass within numerical precision.   The size of dust particles used in our model is determined by a two-population algorithm written from the approach depicted in [33]. This algorithm computes the representative size of particles through the estimate of the limiting Stokes number in various dynamical regimes. In our model, dust is initially present in the form of particles of sizes a 0 = 10 −7 m, and grow through mutual collisions. This growth is limited by the maximum sizes imposed by fragmentation or by the drift velocity of the grains (see [34] for details). The dust surface density is the sum over all surface densities of available solids at given time and location, assuming a protosolar ice-to-rock ratio (approx. 2.57) [35], and a bulk density of 1 g cm −3 for rocks and 0.3 g cm −3 for ices.
We follow the approaches of [33,36,37] for the dynamics of trace species in terms of motion and thermodynamics, respectively. Our two-populations algorithm provides the size of particles, as well as the associated Stokes number and diffusion coefficient. These are inputs to the onedimensional radial advection-diffusion depicting the motion of vapours and dust [33,36]. We assume that all trace species are entirely independent in our simulations and that the disc is uniformly filled with H 2 O, CO, CO 2 , CH 3 OH, CH 4 , N 2 and H 2 S. These molecules are considered to be the dominant volatile species in the PSN, assuming protosolar abundances for O, C, N and S [35]. Half S is assumed to be in H 2 S form with the other half forming refractory sulfide components [38], and all C forms CO, CO 2 , CH 3 OH and CH 4 , with the remaining O going into H 2 O. We have set CO : CO 2 : CH 3 OH : CH 4 = 10 : 30 : 1.67 : 1 in the gas phase of the disc. The CO : CO 2 ratio comes from ROSINA observations of the coma of comet 67P/Churyumov-Gerasimenko from which estimates of the production rates of both CO and CO 2 [39]. The CO : CH 3 OH : CH 4 ratio is consistent with the production rates measured in the southern hemisphere of 67P/C-G in October 2014 by the ROSINA instrument [40]. No chemistry is assumed to happen between the trace species. Sublimation of grains occurs during their inward drift when partial pressures of trace species become lower than the corresponding vapour pressures. Once released, vapours diffuse both inward and outward. Because of the outward diffusion, vapours can recondense back in solid form following the rates defined by Drążkowska & Alibert [37], and condensation occurs either until thermodynamic equilibrium is reached or until no more gas is available to condense. The motion of dust and vapour is computed by integrating the onedimensional radial advection-diffusion equation derived from [33,36], and detailed in [34]. The vapour pressures of trace species are taken from [41]. Figure 1 represents the radial profiles of temperature, gas surface density, dust surface density, and absolute value of mass accretion rate as a function of time in our adopted PSN model. Figure 2 shows the evolution of the condensation radii of the considered species as a function of time in our PSN model.   vapours. Peaks of vapours also appear along the ice lines, as a result of the sublimation of drifting ices. With time, the O abundance in vapour form increases in the inner disc, then decreases until it reaches a plateau close to the protosolar value at 1 Myr. By contrast, the O abundance in solid form decreases with time and heliocentric distance beyond the peak formed after the H 2 O ice line in the outer disc. In this case, the O abundance drastically varies with heliocentric distance but always remains subsolar. The vapour and solid carbon enrichment profiles correspond to the sum of CO, CO 2 , CH 3 OH and CH 4 contributions. These profiles present characteristics close to those found in the vapour and solid O enrichment profiles, and evolve in a similar manner.

Results
The behaviour of the enrichment profiles of N-, and S-bearing species is simpler in our calculations because they individually only depend on one species. Peaks of vapours and solids quickly form to the left and right of the ice lines, respectively. The two panels share the same features with a moderate enrichment (about two to three times protosolar at most) of the vapours present in the inner disc and a significant depletion of the solids in the outer disc. The figure shows that it is not possible to accrete a mixture composed of gas and solids with a C/H ratio presenting enrichments comparable to the values (approx. 70 times protosolar) measured in the atmospheres of Uranus and Neptune. This suggests that the formation of Uranus and Neptune from disc instability (DI) [42] is rather unlikely, since in this mechanism the formed planets consist of large fractions of H-He, even if clumps formed by DI are expected to accrete solids (pebbles/planetesimals) and could therefore significantly increase their metallicities. However,  since the planets formed by this mechanism are typically gas giants, a post-formation mechanism, such as photoevaporation must occur to explain the current masses and compositions of the ice giants. It is therefore more accepted that the ice giants formed by the concurrent accretion of gas and solids by migrating embryos [43,44], provided the accreted mixtures shared a composition consistent with those of Uranus and Neptune. Subsequent photoevaporation can also be consistent with this scenario. Figure 4 shows the time evolution of the C/N, C/S and S/N ratios radial profiles throughout the PSN. Very high abundance ratios can be obtained in solid forms in the 10-20 arb. units region of the PSN, i.e. the location where the ice lines are at play. By contrast, these ratios, still in solid phase, are those initially postulated in the PSN ((C/N) = (C/N) 3.4, (C/S) = 2 × (C/S) 33.9, (S/N) = 0.5 × (S/N) 0.1) at higher heliocentric distances, irrespective of the epoch considered. While the inner regions of the PSN cannot be considered as a viable formation region of the ice giants' building blocks, it is noticeable that the C/N and C/S ratios, both in vapour phase, are found to be supersolar at these locations. If the building blocks accreted by Uranus and Neptune agglomerated from grains formed in the 10-20 arb. units region of the PSN, they should display C/N and C/S ratios close to those recently measured in their tropospheres (C/N ≥ 175 and C/S ≥ 35; see table 1). On the other hand, if the ice giants' building blocks agglomerated from grains formed beyond these distances, they should display protosolar C/N and C/S ratios, in agreement with some hypotheses formulated regarding the composition of their interiors [11,18].

Conclusion
In this work, assuming the carbon abundances determined in the atmospheres of Uranus and Neptune are representative of their bulk composition (the case of a fully mixed interior, which is likely to be oversimplified), we have shown that it is impossible to accrete a mixture composed of gas and solids from the PSN with a C/H ratio presenting enrichments comparable to the measurements (approx. 70 times protosolar). By contrast, we found that the C/N and C/S ratios in Uranus and Neptune are compatible with those acquired by grains condensed in the vicinities of N 2 and CO ice lines in the PSN. If the measurements of the C/N and C/S ratios in the ice giants are representative of their bulk compositions, then they could be explained by the formation of their building blocks from grains and pebbles condensed in this region of the PSN. On the other hand, the presence of protosolar C/N and C/S ratios in Uranus and Neptune would imply that their building blocks agglomerated from particles condensed at higher heliocentric distances.
Variations of the viscosity parameter α in the 10 −4 -10 −2 range, which is commonly adopted for PSN models, show that the results are qualitatively close to those presented here and that the conclusions remain unchanged. The lowest value α = 10 −4 corresponds to the assumption of an MRI dead zone in our simulation box [45][46][47]. Our results are quite different from those derived by Ali-Dib et al. [15] who found very high abundances of solid CO in the vicinity of its iceline. As mentioned above, the reasons for these differences are twofold: the PSN model used by Ali-Dib et al. [15] is stationary and assumes instantaneous condensation of vapours, instead of using a time-dependent model with prescriptions for condensation rates, as we did here.
Our study demonstrates the importance of measuring the elemental abundances in the ice giant atmospheres, as they can be used to trace the planetary formation location and/or the chemical and physical conditions of the PSN. Future facilities such as the James Webb Space Telescope should provide better determinations of the elemental abundances in Uranus and Neptune in the near future [48]. A question that arises is whether the atmospheric elemental abundances indeed represent the bulk. It is possible that the measured composition is affected by a recent impact of comets/asteroids, as it has been invoked for the detection of CO. In addition, as both Uranus and Neptune are unlikely to be fully mixed, it is possible that the ratios between the different elements change with depth due to chemistry. This argues in favour of the in situ measurement by an entry probe of the noble gases in the envelopes of Uranus and Neptune since these species should be insensitive to any form of chemistry and are not expected to condense at depth. Even if the absolute abundances of noble gases should be affected by the compositional gradient with depth due, for example, to primordial composition gradients and/or layered convection, their abundance ratios (He/Ne/Ar/Kr/Xe) are expected to remain constant.
Although it has been suggested that Uranus and Neptune could be formed by disc instability [42], this is rather unlikely given that the planetary composition is dominated by the presence of heavy elements. In order to have the ice giants formed by this mechanism a substantial amount of solids must have been captured [49], followed by significant loss of H-He gas via photoevaporation [50]. Core accretion remains the preferred model for the formation of Uranus and Neptune in a way consistent with their high metallicities, despite the clear challenges in terms of formation timescale and expected composition [1]. Our study shows that if the heavy elements delivered to the ice giants at their formation location are mostly composed of H 2 O, CO, N 2 and H 2 S, it is very difficult to explain the measured C, N and S atmospheric abundances via classical structure models. This emphasizes the possibility that the classical 3-layer structure model is not consistent with the composition of the ice giants. It is also possible that the inferred atmospheric composition does not represent the bulk composition and/or that the planets are not fully mixed. Alternatively, the metallicity of the PSN could have been much higher than the protosolar value assumed here and/or that there was an accumulation of solids near the formation locations of the planets. However, the present study shows that the solids enhancement at the locations of the ice lines is not as efficient as was initially anticipated [15]. Finally, this work emphasizes the link between the early stages of the formation of planets in the disc phase, their evolution (e.g. mixing) and internal structure.