© The Authors 2025

1 Introduction

A considerable number of exo-Kuiper belts, also referred to as debris discs, have been observed over the past few decades. First detected via their infrared excess due to the dust emission they produce (Aumann et al. 1984; Eiroa et al. 2010), they can now be imaged at high resolution with e.g. ALMA (Sepulveda et al. 2019). Probably one of the most unexpected recent results was the discovery of gas (mostly CO and usually referred to as late gas) in relatively old (10–700 Myr) exo-Kuiper belts that orbit around main-sequence stars (Moór et al. 2017; MacGregor et al. 2017). To date, almost 30 debris discs containing gas have been observed, among which some have a considerable mass of CO (up to 0.1 M_⊕, Moór et al. 2019). Observations indicate that CO, C, and O gas species are present (Cataldi et al. 2014; Brandeker et al. 2016) and the gas composition is expected to be dominated by those species rather than hydrogen or helium (Hughes et al. 2017).

Models indicate that this late gas is most likely not a remnant of the protoplanetary disc gas, but rather a secondary phenomenon, where gas is released from volatiles contained in planetesimals in these debris belts (Kral et al. 2017, 2019). Hence, a considerable quantity of gas is produced in the belt over an extended period of time (approximately 100 Myr), as opposed to the typical several millions years for protoplanetary discs. It is anticipated that this gas viscously spreads both inwards and outwards, plausibly because of magnetorotational instability or molecular viscosity (Kral & Latter 2016; Cui et al. 2024).

As these systems are mature, they are expected to contain planets that are already formed, especially giants, which are potentially close to the belts. In fact, thanks to direct imaging, such planets have already been observed in debris disc systems (Lagrange et al. 2009, 2019). It is then expected that the planets are within reach of the viscously spreading gas disc and that some of this late gas is accreted onto the planets. The analytical model developed by Kral et al. (2020) suggests that the accretion of this late gas onto planets is highly efficient, given that it occurs over an extended period and that the atmospheres have sufficient time to cool.

The late gas release observed in relatively young extrasolar systems possibly also happened in the young Solar System. However, today’s Kuiper belt (KB) is considerably less massive than the aforementioned debris discs and gas release may be very low. Although a certain amount of CO could still be released in the current KB, the solar wind would expel this low-density gas from the belt, making it difficult to detect (Kral et al. 2021). However, our Solar System is considerably more ancient than those with extrasolar gas detected. It is reasonable to assume that the Kuiper belt would once have been much more massive and thus produce a larger quantity of gas, sufficient to overcome the solar wind and, like in some extrasolar systems, capable of spreading viscously inwards as far as Neptune and Uranus, and even Saturn and Jupiter. Indeed, the current low-mass Kuiper belt is expected to be a remnant of a much larger belt that was dispersed during the early stages of the Solar System’s formation by the influence of Neptune and Uranus, which makes it likely that it resembled the gas-rich belts observed today (Bottke et al. 2023).

One of the most studied models of the evolution of our Solar System designed to reproduce its current architecture is the Nice model (Tsiganis et al. 2005; Morbidelli et al. 2005; Gomes et al. 2005). It was progressively enriched over time to solve fine-tuning issues (Levison et al. 2011), timing issues (Nesvorný & Morbidelli 2012), and get a more precise representation of the primordial disc (Griveaud et al. 2024). The Nice model suggests that the system was in a more compact configuration than today, with planets within 16 au. The planets were then pushed onto their current orbits after an instability between Jupiter and Saturn via resonant interactions between them. However, to get the model working and push Jupiter and Saturn into resonance, the initial belt needs to have a mass between 20 to 50 M_⊕ in place of the ~0.1 M_⊕ of the current Kuiper belt. The primordial massive belt gets depleted when Neptune migrates outwards, which captures some Kuiper belt objects into resonances and can explain the current KB architecture. Moreover, refined collisional models such as Bottke et al. (2023) show that the Nice model can explain the size distributions of Trojans and of craters on the moons Europa and Ganymede, as well as many features of the current Kuiper belt.

There are also alternative models that may explain the broad strides of the evolution of the Solar System. For instance, Liu et al. (2022) argue that as the disc dissipates from the inside out via photoevaporation, some ‘rebound’ effect may happen that could push the outer ice giants outwards and destabilize the system, thus causing outward migrations for Jupiter, Saturn, Uranus, and Neptune. In the end, after the dissipation of the protoplanetary disc, the system is already in an extended configuration close to the current one, unlike in the Nice model. In this case, a belt mass as low as 5 M_⊕ is enough to reproduce the current planet orbits.

If the primordial Kuiper belt was indeed massive, as currently assumed by most models, the gas release rate could have been similar to what is observed in some extrasolar systems, which may have important consequences on Uranus and Neptune atmospheres compsotion. Indeed, the icy giants are the closest to the Kuiper belt, and they are expected to have accreted this late gas, which may significantly alter the composition of their atmospheres, as is shown in this paper.

The metallicity of the atmospheres of Neptune and Uranus is constrained to be highly super-solar, with a C/H ratio of ~50–80 times the protosolar metallicity (Atreya et al. 2020; Guillot et al. 2023). To explain such an abnormally high metallicity, the current models posit that the planets formed in proximity to the CO ice line by pebble accretion (Ali-Dib et al. 2014; Mousis et al. 2024a,b). The pebbles in this area are expected to be enriched in carbon. Pebbles drift inwards because of their interaction with primordial gas and then sublimate when they cross the ice line. Unless otherwise stated, the term ‘ice’ in this paper refers exclusively to CO ice. The sublimated gas viscously spreads inwards and outwards. Then a part of it condensates onto the icy pebbles that have not yet crossed the CO ice line. When Uranus and Neptune are close to the ice line, they accrete pebbles enriched in carbon, which can enhance the C/H ratio in their atmospheres. However, even though those pebbles might pollute the hydrogenhelium atmospheres of Uranus and Neptune, there are still large uncertainties about their final impact (Mousis et al. 2024b).

The atmospheres of Neptune and Uranus are ~1–4 M_⊕, i.e. much lighter than the total masses of the planets (Guillot et al. 2023). Accreted late CO gas is mixed in those relatively light atmospheres, which might significantly contribute to this observed super-solar metallicity, as we explore in further detail in this paper.

After introducing our numerical model in Section 2, we study the carbon enrichment of Uranus and Neptune with our model for different plausible scenarios of the Kuiper belt evolution in Section 3. Finally, we discuss our results in Section 4 before concluding. Briefly, we find that the primordial belts tested produced a significant amount of gas and that the most efficient Solar System formation model that results in an enrichment of the atmospheres of Uranus and Neptune in carbon at the observed level is that derived from the most recent extension of the Nice model (Griveaud et al. 2024).

2 Method

The aim of this paper is to investigate the accretion of gas produced in the Kuiper belt early in the Solar System evolution onto the Solar System’s outer planets (Jupiter, Saturn, Uranus, and Neptune) over the entire lifetime of the late gas disc. Debris discs may persist for several Gyr, like in our Solar System, and the gas in it is thus expected to be present for a couple hundred million years at the minimum (Matrà et al. 2017; Moór et al. 2017).

2.1 Viscous diffusion

2.1.1 The viscous model

During the large time span of the gas in the debris disc phase, it will have time to viscously evolve (Kral et al. 2016; Kral & Latter 2016). However, given current numerical constraints, it is not realistic to run hydrodynamic simulations over such an extended period of time. In this context, we adopted the numerical methodology proposed for debris disc by, e.g. Moór et al. (2019); Marino et al. (2020), and solve a one-dimensional viscous diffusion equation for a fluid made up of a mixture of chemical elements, which orbits the star. We assume that gas is initially released as CO. Its photodissociation caused by both the interstellar radiation field (ISRF) and stellar photons results in the partial conversion of this gas into carbon and oxygen, hence the use multi-species model. This means that each fluid component has its own dynamic influenced by the other species.

In our model, we assume that the viscous parameter α (Lynden-Bell & Pringle 1974) is identical for all species, that the disc is axisymmetric and that the vertical scale height is defined, for a given radius R, by H(R) = c_s/Ω, where $\mathrm{\Omega }=\sqrt{\frac{G{M}_{\ast }}{R^3}}$ is the angular speed and $c_s=\sqrt{\frac{k_{B}T}{\mu m_p}}$ is the sound speed. M_* is the mass of the central star, T is the gas temperature at radius R, μ is the mean molecular weight, G and k_B are respectively the gravitational and the Boltzmann constants, and m_p is the proton mass. Let us also define Σ the total surface gas density and Σ_i the surface density of the different species i such that $$\mathrm{\Sigma }=\sum _i{\mathrm{\Sigma }}_i.$$(1)

Our simulations will account for three species, namely CO, C and O. The radial velocity of each species has two components: the global fluid radial velocity v_r and the diffusive flux v_ri which is the velocity due to the diffusion between species.

In our model, ${\dot{\mathrm{\Sigma }}}_i(R,t)$ is the gas generation/destruction rate for species i. This term is the sum of the gas generation rate in the belt, and of its destruction rate due to photodissociation and planetary accretion, i.e. all the non-diffusive terms. The prescriptions for ${\dot{\mathrm{\Sigma }}}_i(R,t)$, depending on the species, will be described more precisely in Subsection 2.5. The gas production rate ${\dot{\mathrm{\Sigma }}}_{\mathrm{C}\mathrm{O}}$ is obtained via the total gas mass production rate ${\dot{M}}_{\mathrm{C}\mathrm{O}}$ calculated in Subsections 2.2 and 2.3. We assumed a radial power law distribution for the belt of −3/2 (as in e.g. Kral et al. 2019). Therefore, we get ${\dot{\mathrm{\Sigma }}}_{\mathrm{C}\mathrm{O}}(r)=\frac{{\dot{M}}_{\mathrm{C}\mathrm{O}}}{S_{\text{norm }}}{\left(\frac{r}{a_0}\right)}^{-3/2}$ with a₀ and Δa the belt’s semi-major axis and radial extension, respectively, and where $S_{\text{norm }}={\int }_{a_0-\mathrm{\Delta }a/2}^{a_0+\mathrm{\Delta }a/2}2\pi {\left(\frac{r}{a_0}\right)}^{-3/2}rdr=4\pi a_0^2(\sqrt{1+\frac{\mathrm{\Delta }a}{2a_0}}-\sqrt{1-\frac{\mathrm{\Delta }a}{2a_0}})$ is the normalisation factor.

The mass conservation of gas in the disc leads to (in cylindrical coordinates) $$\frac{\mathrm{\partial }{\mathrm{\Sigma }}_i}{\mathrm{\partial }t}=-\frac{1}{R}\frac{\mathrm{\partial }}{\mathrm{\partial }R}\left(R{v}_r{\mathrm{\Sigma }}_i\right)-\frac{1}{R}\frac{\mathrm{\partial }}{\mathrm{\partial }R}\left(R{v}_{ri}\mathrm{\Sigma }\right)+{\dot{\mathrm{\Sigma }}}_i(R,t),$$(2)

while the radial velocity can be obtained from the momentum conservation since there is a friction torque defined via the viscous coefficient ν = αc_sH, where α is the viscous parameter (Shakura & Sunyaev 1973) assumed to be constant in our disc $$\mathrm{\Sigma }{v}_r=-\frac{3}{\sqrt{R}}\frac{\mathrm{\partial }}{\mathrm{\partial }R}(\nu \mathrm{\Sigma }\sqrt{R}).$$(3)

We consider α values ranging from 10⁻³ to 10⁻¹ (see Figure 1). Indeed, for debris discs, both turbulent viscosity, such as that caused by the magnetorotational instability, and molecular viscosity result in higher values of α than those observed in protoplanetary discs (Kral & Latter 2016; Cui et al. 2024).

Following Charnoz et al. (2019), the interspecies diffusion velocity equation is given by $$v_{ri}=-v\frac{\mathrm{\partial }}{\mathrm{\partial }R}\left(\frac{{\mathrm{\Sigma }}_i}{\mathrm{\Sigma }}\right).$$(4)

Using equations (1), (2), (3), and (4), we can retrieve the usual mono-species viscous evolution equation for Σ defined classically in Lynden-Bell & Pringle (1974) for $\dot{\mathrm{\Sigma }}(R,t)=0$.

Fig. 1 Schematic presenting a description of our set of nine simulations. The heavy belt cases are at the top (compact configuration) and middle (extended configuration), and the light cases are at the bottom. Jupiter, Saturn, Uranus, and Neptune are designated by J, S, U, and N, respectively. The different simulations differ in terms of the locations of planets, accretion efficiency, the viscosity of the gas, and the belt’s mass. For the heave belt configurations, the time when depletion starts is also a parameter. Each parameter has its own symbol and colour as can be seen in the bottom part of the image.

2.1.2 The numerical approach

The set of partial differential equations we need to solve has no general analytical solution, and must be solved numerically. In order to achieve this, the spatial domain is discretized on a grid with n cells via the finite difference method using centred numerical derivatives. Subsequently, we obtain n differential equations of order 2 in time, which we turn into 2n order 1 differential equations in time. To solve these equations, we use an adaptive time step RK45 solver, via the solve_ivp function of the SciPy package (Virtanen et al. 2020). After testing our setup, we chose n = 8000 to circumvent numerical instabilities and ensure global mass conservation.

We use boundary conditions for the surface densities that are similar to those set by Marino et al. (2020). For the outermost cell, we use a power-law extrapolation of the surface density based on the nearest cells. Although we must truncate the disc at an arbitrary outer radius, the disc is expected to continuously decrease in surface density rather than stop as expected at steady state (Kral et al. 2019). We thus assume a power-law extrapolation as a good physical approximate of the outer disc. Since the computational cost increases as the step size decreases, it is impractical to model the innermost part of the disc. Therefore, we also use a power-law extrapolation for the innermost cell. However, this extrapolation may be less reliable than for the outer part of the disc, especially at the beginning of the simulation where the gas surface density varies significantly. Thus, we also add the constraint that the flux in the innermost cell cannot be greater than the flux in the previous cell via the condition νΣ = Const. We verify that our numerical model, using these boundary conditions, converges to the expected analytical solution for a constant mass production rate.

2.2 Gas produced from warming ices in the KB

There are two main processes that can release gas in the Kuiper belt. First, gas can be released via collisions as explained in more detail in Subsection 2.3, and second, gas can be released by sublimation of warming KB objects (KBOs) over time (as proposed in Kral et al. 2021), which we study in this subsection. For our different setups, we find that one process always dominates the other as we will mention more specifically later.

For KBOs in the Solar System, the timescale for the sublimation of CO ices and the timescale for the CO diffusion into the porous structure of KBOs is relatively short in comparison to the thermal diffusion characteristic timescale (Lellouch et al. 2013; Kral et al. 2021). Therefore, KBOs will take time to reach their equilibrium temperature, while the Sun’s heat penetrates slowly inside KBOs in an outside-in fashion. Since the sublimation timescale is much shorter than the thermal diffusion timescale, the model will not depend on the distance to the central star, or rather only indirectly, assuming the belt is close enough to the central star for CO ice to sublimate. Let us consider a KBO with a diameter s that is initially composed of refractory materials and CO ices, with a ratio of CO ice mass to refractory mass f_ice = M_ice/M_refr with typical values of 0.1 (e.g. Mumma & Charnley 2011).

At a given time t, the ice within the radius $r=\sqrt{Kt}$ (where K is the thermal conductivity of the KBO) will remain intact but it will turn into gas above r. We assume a thermal conductivity K of ~10⁻¹⁰ m² s⁻¹ as given by observations of these cold objects (Prialnik et al. 2004; Lellouch et al. 2013). Because all the CO above the radius r sublimates, we can calculate that during a timescale dt, all of CO in the layer between $s/2-\sqrt{Kt}$ and $s/2-\sqrt{K(t+dt)}$ sublimates. Hence, the mass loss rate for a KBO of diameter s as a function of time can be easily estimated (e.g. Kral et al. 2021) $$\begin{array}{rl}\frac{\mathrm{d}{M}_{\mathrm{C}\mathrm{O}}}{\mathrm{d}t}(s)& =2\pi {\rho }_{\text{refr }}{f}_{ice}{K}^{3/2}\sqrt{t}\text{ if }t<\frac{s^2}{4K},\\ & =0\text{ if }t\ge \frac{s^2}{4K},\end{array}$$(5)

where ρ_refr is the density of the refractory material.

As explored in the next subsection, the gas production rate from collisions is usually much higher when the belt is young and massive. Hence, we expect the ice sublimation to only dominate when the disc becomes less massive, which happens after the massive primordial belt depletion (see Section 2.6). We assume that the KBO number density will not evolve significantly during this later phase and we can calculate the total sublimation rate at a time t coming from all KBOs in the belt as $$\frac{\mathrm{d}{M}_{\mathrm{C}\mathrm{O}}}{\mathrm{d}t}={\int }_{max(s_{min},\sqrt{Kt})}^{s_{max}}2\pi {\rho }_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{r}}f{K}^{3/2}\sqrt{t}n(s)\mathrm{d}s,$$(6)

where n(s) is the belt’s size distribution, which will be one of the parameters we explore in this study. ${\dot{M}}_{\text{CO}}$ can then be fed in our model via the ${\dot{\mathrm{\Sigma }}}_{\mathrm{C}\mathrm{O}}(R,t)$ variable.

2.3 Gas produced from collisions

In younger, heavier belts, such as the primordial Kuiper belt we consider here, the effect of collisions onto gas production might not be negligible (Thébault & Augereau 2007; Bonsor et al. 2023). The standard model to predict gas production rates in collisionally evolving debris discs is described in Kral et al. (2017). It assumes that the gas production rate is proportional to the solid mass loss rate. The idea is that collisions produce smaller debris and expose new surfaces that can sublimate and release gas. We use this approach to compute the gas release rate due to collisions coupled to an analytic model for the estimation of the solid mass loss rate, able to account for several size distribution slopes (Löhne et al. 2008). The latter model can be adapted to compute the solid mass loss rate in the primordial Kuiper belt using primordial slopes for the young Kuiper belt size distribution from Bottke et al. (2023). We use formula 38 of Löhne et al. (2008) to estimate the belt’s mass as a function of time and then derive it numerically to obtain the solid mass loss rate ${\dot{M}}_{\text{belt}}$. Similar to Kral et al. (2017), we then estimate the gas generation rate to be ${\dot{M}}_{\mathrm{C}\mathrm{O}}=-f_{\text{ice }}{\dot{M}}_{\text{belt }}$, where we assume that all ice contained in KBOs will eventually sublimate as the planetesimals are turned into dust. This is a reasonable assumption since the area of exposed KBOs increases significantly during collisions (since it produces smaller chunks with higher cross sections given assumed size distributions) and the sublimation proceeds at the surface, not deep inside the KBOs.

2.4 Influence of the planet

In order to model three-dimensional hydrothermodynamic effects, such as planetary accretion, within our one-dimensional viscous diffusion code, it is necessary to parameterize the problem. To do so, let us define the accretion efficiency $$f_{\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{r}}=\frac{{\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{r}}}{{\varphi }_{in}},$$(7)

where ${\dot{M}}_{\text{accr}}$ is the planet’s mass accretion rate, and ϕ_in denotes the incoming mass flux. It should be noted that the incoming mass flux is not necessarily the total mass flux at the planet’s semi-major axis. Indeed, the gas disc has a considerable vertical scale height, whereas the planet has a gravitational sphere of influence of radius R_accr, beyond which the gas is not accreted. Accounting for that, the incoming mass flux to a planet of semimajor axis a can be expressed using the minimum function as follows (Kral et al. 2020) $${\varphi }_{in}(t)=min(1,\frac{R_{\text{accr }}}{H})\mathrm{\Sigma }(a,t)v_r(a,t).$$(8)

We estimate the planet’s accretion radius, R_accr, such that R_accr = min(R_Hill, R_B), where $R_{\text{Hill }}=a{\left(\frac{M_{\text{planet }}}{3M_{\ast }}\right)}^{1/3}$ is the Hill radius beyond which the central star’s gravity will be stronger than the planet’s gravity, and $R_B=\frac{G{M}_{planet}}{v_r^2(a,t)}$ the Bondi radius, which is the planet’s escape radius for a particle moving at radial velocity v_r. Note that we do not use the gas sound speed in the disc as the characteristic speed, since the planet’s atmosphere is expected to be decoupled from the disc in low-density environments similar to the late stages of giant planet formation in protoplanetary discs (Mordasini et al. 2012). We note that in our upcoming simulations, R_Hill is always smaller than R_B.

The gas discs under consideration have typically a lower density than protoplanetary discs and in this case, Kral et al. (2020) predict that the accretion efficiency is very high and may be close to 100% (as will be discussed further in Section 4). Consequently, the main limiting factor of the gas accretion rate onto a planet would be determined by hydrodynamic processes in addition to the radius of influence R_accr. As demonstrated by previous studies in protoplanetary disc environments, it is challenging to accurately determine the amount of gas flowing around the planet that gets accreted (Lubow & D’Angelo 2006; Mordasini et al. 2012; Ormel et al. 2015). In addition, numerical hydrodynamic simulations have not yet been run in the debris disc phase and more refined effects may have to be accounted for. Based on previous numerical estimates (Lubow & D’Angelo 2006; Mordasini et al. 2012; Tanigawa et al. 2012; Ormel et al. 2015; Lambrechts et al. 2019), we assume an accretion efficiency due to hydrodynamics processes of f_accr = 0.5. However, to account for uncertainties in this parameter, based on those hydrodynamic studies, we also test two alternative values of 0.1 and 0.8.

Numerically, we use sink cells to account for gas accretion onto planets. Therefore, at each time step, we compute the amount of gas mass that should be accreted onto the different planets and remove it from the gas phase similar to simulations of Kral et al. (2024) for water released in the young asteroid belt.

2.5 Photodissociation

Our model follows the evolution of carbon monoxide outgassed in the belt as well as its photodissociation products C and O. We assume that photodissociation is mainly produced by the interstellar radiation field (ISRF). This is because the Kuiper belt is located at a considerable distance from the Sun and that the disc’s radial extent, due to radial spreading, is significantly larger than its vertical scale height. Consequently, the majority of the disc is optically thick in the radial direction and thus mostly opaque to solar radiation in the UV (e.g. Kral et al. 2016).

The theoretical lifetime of CO, t₀, in the case of impinging photons from the ISRF is estimated to be 120 years (Visser et al. 2009). However, CO can be shielded both by itself (Visser et al. 2009) and by neutral carbon (Kral et al. 2019). We interpolate the CO self-shielding function Θ(Σ_CO) used in our model (when no self-shielding is present, it tends to 1) from Visser et al. (2009). For neutral carbon, we estimate a critical surface density Σ_crit above which the optical depth to incoming vertical radiation would be superior to one. Following Kral et al. (2019); Marino et al. (2020), we expect carbon shielding to have an exponential dependence on neutral carbon surface density, with a critical density Σ_crit ~ 10⁻⁷ M_⊕ au⁻². In this model, we assume that the carbon ionization fraction is small, similar to Marino et al. (2020). This assumption is less reliable for the discs with low masses (Kral et al. 2017). However, even in the extreme case where all carbon is ionized, because the gas disc is of low mass, shielding will not be important and the ionized carbon gas disc will behave exactly as the neutral carbon since it has the same molecular mass. We find that the interspecies diffusion is negligible in our simulations leading to no differences between neutral and ionized carbon diffusion.

Finally, the CO photodissociation timescale in the disc can be estimated as $$t_{\mathrm{p}\mathrm{h}}=t_0\frac{\exp \left[\frac{{\mathrm{\Sigma }}_{\mathrm{C}}}{{\mathrm{\Sigma }}_{\text{crit }}}\right]}{\mathrm{\Theta }\left({\mathrm{\Sigma }}_{\mathrm{C}\mathrm{O}}\right)},$$(9)

where Σ_C and Σ_CO are the surface density of neutral carbon and that of CO, respectively. The photodissociation contribution to

${\dot{\mathrm{\Sigma }}}_{\mathrm{i}}$ is then given for the different species by $${\dot{\mathrm{\Sigma }}}_{\mathrm{C}\mathrm{O}}=-\frac{{\mathrm{\Sigma }}_{\mathrm{C}\mathrm{O}}}{t_{\mathrm{p}\mathrm{h}}},$$(10) $${\dot{\mathrm{\Sigma }}}_{\mathrm{C}}=\frac{{\mu }_{\mathrm{C}}}{{\mu }_{\mathrm{C}\mathrm{O}}}\frac{{\mathrm{\Sigma }}_{\mathrm{C}\mathrm{O}}}{t_{\mathrm{p}\mathrm{h}}},$$(11) $${\dot{\mathrm{\Sigma }}}_{\mathrm{O}}=\frac{{\mu }_{\mathrm{O}}}{{\mu }_{\mathrm{C}\mathrm{O}}}\frac{{\mathrm{\Sigma }}_{\mathrm{C}\mathrm{O}}}{t_{\mathrm{p}\mathrm{h}}},$$(12)

where μ_CO, μ_C, and μ_O are the molecular weight of CO, C, and O, respectively.

2.6 Kuiper belt modelling and evolution

In our simulations, we consider three different configurations of the young Kuiper belt, each corresponding to a distinct scenario for the belt’s early history (see Fig. 1). The first set of simulations (HKB-fid, 2 and 3 simulations, see Table 2) takes as a basis the Nice model, i.e. we start with a heavy primordial belt and planets closer to the Sun than their current positions. The second set (HKB-4, 5) takes its initial setup from the ‘Rebound’ scenario of Liu et al. (2022), which also assumes a relatively massive initial belt but with planets closer to their current locations than in the Nice model. Last, we consider a set of simulations (LKB-fid, 2, 3, 4, see Table 1) assuming that the KB was born “light”, with a configuration similar to the current KB and planets at their current locations. This light-KB scenario can be considered as a lower limit in terms of gas production within the belt and its effect on Uranus and Neptune atmospheres.

2.6.1 The Nice model scenario

In the Nice model scenario, the primordial belt is massive (≳10 M_⊕) and evolves over a time t_ev of 10–100 million years before being depleted by Neptune’s outer migration. Recent results suggest that this value is closer to 10 Myr, which is what we assume in our model (Bottke et al. 2023). We thus consider that the mass depletion of the belt starts after a time t_ev and is modelled by an exponential decay with an e-folding time t_fold. For the Nice model, we expect t_fold ~ 10 Myr (Bottke et al. 2023), which is the value we take. We assume that the mass tends to zero at the end of the simulations, which is unrealistic but does not affect the final results, dominated by the early phase.

The planets and belt positions are also similar to those assumed in the most recent iterations of the Nice model (e.g. Griveaud et al. 2024, and references therein). The primordial belt is taken to be located at 27 au from the Sun and to be 6 au wide (Bottke et al. 2023). Neptune, Uranus, Saturn, and Jupiter are in a compact configuration, at 16, 14, 8, and 5 au, respectively. Because of gravitational interactions with the primordial KB, the outer planets move outwards and Jupiter and Saturn cross their 2:1 resonance, which excites their eccentricities and triggers an instability, which eventually depletes the KB. We do not model the migration of planets before the instability and subsequent depletion but keep the belt and planet positions fixed throughout. This is because the positions of planets are not an important parameter and do not alter our final results as discussed later. However, because of the migration of planets, planetesimals in the young belt become excited dynamically. To account for that, we assume a mean eccentricity of 0.25 in our collisional model, slightly higher than the 0.05–0.1 typically assumed (e.g. Thébault & Augereau 2007).

The size distribution of the primordial KB is assumed to be different from that of the current KB. We use state-of-the-art models of the early KB for our simulations, which assume that the size distribution can be separated into three regimes: for bodies between 10⁻⁴ m and 100 km in diameter, the size distribution has a slope of −2.1, between 100 and 300 km the slope is −6, and for bodies up to 4000 km the slope is −3.5 (Bottke et al. 2023). We assume that the spatial distribution scales as R^−3/2, independent of time (similar to Bottke et al. 2023). The HKB-fid, 2 and 3 simulations only differ in the assumed accretion efficiency of gas onto planets (0.5 and 0.1) and the initial belt mass (20 and 50 M_⊕ are tested). In these simulations, the gas production mechanism is dominated by collisions (see Sec. 2.3).

2.6.2 Alternative heavy initial belt scenarios

We run another type of simulations to account for scenarios in which the primordial belt starts heavy but the planets, especially Neptune, start farther away, i.e. closer to their current positions. In this case, models do not need as heavy a belt as in the Nice model, and an initial belt under 10 M_⊕ can reproduce the Solar System evolution (e.g. Liu et al. 2022). We assume a typical belt mass of 5 M_⊕ for those scenarios.

The ‘Rebound’ scenario proposed recently in Liu et al. (2022) for forming the Solar System can be modelled via those simulations. In this scenario, the protoplanetary disc is photo-evaporated from the inside out from a certain radius, creating a gap in the young protoplanetary disc. When the outer edge of the gap reaches Saturn, it pushes the giant outwards, which triggers an orbital compression of planets and finally an instability (Liu et al. 2022). In this scenario, Neptune migrates over a smaller distance than in the Nice model. In our simulations taking as a basis this scenario, which shall be referred to as the “extended” configuration, Neptune’s location is 23 au (instead of 16 au in the Nice model runs). For all the other 3 planets, we consider the same locations (5,8, and 14 au) as in the compact Nice-model simulations of Section 2.6.1. As for the compact simulations, the migration of the planets is not dynamically modelled and the planets remain at fixed positions throughout the runs, but we will show that the planets’ locations only weakly affect their level of carbon enrichment. As for the eccentricity, belt position and width, planetesimal size distributions and surface density, they are taken identical as in the compact Nice-model simulations.

We run 2 simulations for this extended configuration, HKB-4 and HKB-5, which differ in their t_ev and t_fold values. The HKB-4 simulation uses the same values as that of the Nice model (10 Myr for both parameters) but for the HKB-5 simulation, we test values of 100 Myr for both t_ev and t_fold. This is because the survival timescale of the primordial KB in those scenarios is not very well constrained and needs to be explored. In these simulations, the gas production mechanism is also dominated by collisions (see Sec. 2.3).

2.6.3 Light early belt scenarios

Some less common scenarios for the formation of the Solar System assume that the KB was born light with a mass and size distribution very similar to the belt we observe today (Shannon et al. 2016). We run 4 simulations (LKB-fid, 2, 3, and 4) based on this scenario, for which we assume an initial belt mass of 0.1 M_⊕ and a size distribution corresponding to the current one. Using Morbidelli et al. (2021) we take a slope of −4 for bodies with radii between 10⁻⁴ and 30 m, of −3 for bodies smaller than 100 km and, −8 for bodies smaller than 4000 km. For these light belt simulations, we assume the KB’s surface density follows a smooth Gaussian distribution centred at R₀ = 44 au with a total width of 8 au. Similar to the heavy primordial belt scenarios, we do not model any dynamical effect so that the locations of the planets remain fixed throughout our simulations and are assigned to their current values of 30, 19.2, 9.5, and 5.2 au for Neptune, Uranus, Saturn, and Jupiter, respectively. We do not introduce a depletion term for the belt, whose mass only decreases because of the CO ice it loses (see Sect. 2.2).

The LKB-fid, 2, 3, and 4 simulations differ in their respective values of α (10⁻³, 10⁻²), and accretion efficiency (0.1, 0.5, 0.8) as presented in Table 1. In these simulations, the gas production mechanism is dominated by the heating of the ice in the planetesimals rather than by collisions (see Sec.2.2). We note that to be realistic, we could take the start of those simulations as the end state of the heavy-belt scenarios after the belt depletion. However, our final results would not be affected (see discussion) and we decided to optimize CPU time by not coupling those two simulations and only use the LKB models for the scenario where the KB was born light.

2.7 Estimation of the planet atmospheric $\frac{\mathrm{C}}{\mathrm{H}}$ ratios in our model

Observations give us estimates of the atmospheric $\frac{C}{H}$ ratios for Uranus, Neptune, Jupiter and Saturn. However, our numerical model provides estimates of the accreted carbon mass and we will now describe how we turn it into a modelled $\frac{C}{H}$ ratio prediction. We consider that carbon mix into the atmospheres of the giants, which consist of a mixture of hydrogen and helium (Guillot et al. 2023). From the accreted carbon mass $M_{\text{accr }}^C$ and the atmospheric mass M_atm, we can then calculate the $\frac{C}{H}$ abundance ratios such that $$\frac{C}{H}=\frac{{\mu }_H}{{\mu }_C}\frac{M_{\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{r}}^C}{M_{\mathrm{a}\mathrm{t}\mathrm{m}}^H},$$(13)

where μ_H = m_p and μ_C = 12 m_p, with m_p the proton mass, and the atmospheric hydrogen mass $M_{\text{atm}}^H=f_H{M}_{\text{atm}}$, with f_H the hydrogen mass fraction in the atmosphere. We assumed f_H ~ 0.746 (Asplund et al. 2021).

The atmospheres of Neptune and Uranus are light compared to the total masses of these planets. We assumed atmospheric masses of 1.25–3.5 M_⊕ and 1.6–4.15 M_⊕ for Uranus and Neptune, respectively (Guillot et al. 2023). We use the error bars on those atmospheric masses to calculate the error bars of the modelled $\frac{C}{H}$ ratios.

For giant planets, the total atmospheric masses are much higher. We expect the late gas to be accreted early in the life of the giant planets (before 1 Gyr), i.e. before the hydrogen-helium phase separation that only happens once the planets have cooled down (Howard et al. 2024). It is therefore not expected that the current H/He transitions at depths corresponding to atmospheric masses of ~30M_⊕ for both Jupiter and Saturn (Markham & Guillot 2024; Guillot et al. 2023) could have trapped carbon in this convective upper layer. One should rather consider the total atmospheric mass. The only uncertainty then comes from the extent of the diluted solid core. We consider a typical 60 M_⊕ diluted solid core for both Jupiter and Saturn (Howard et al. 2023; Mankovich & Fuller 2021; Guillot et al. 2023). Some static models of Jupiter favour a higher dilute mass of ~50% of the total mass (Militzer et al. 2022), leading to a larger but still small enrichment in Jupiter, given its massive envelope. We also use lower limits of 20 and 30 M_⊕ for the solid core masses of Saturn and Jupiter, respectively (Howard et al. 2023; Mankovich & Fuller 2021; Miguel et al. 2022), which we use to derive upper limits of the total atmospheric mass carbon can mix with.

Since the atmospheric carbon enrichment in the planets is usually given compared to the protosolar abundance in the literature, we will use the notation $\left[\frac{C}{H}\right]$ when relative to protosolar abundance such that $\left[\frac{C}{H}\right]\equiv \frac{\frac{C}{H}}{{\frac{C}{H}|}_{\text{proto}}}$, where ${\frac{C}{H}|}_{\text{proto}}$ is the protosun’s $\frac{C}{H}$ ratio equal to (3.33 ± 0.31) × 10⁻⁴ (Guillot et al. 2023).

3 Results

Here, we explore how a primordial KB may have enhanced the metallicity of the giants of our Solar System. In particular, the metallicities of Uranus and Neptune are highly super-solar with $\left[\frac{C}{H}\right]$ values around 50–80 times the protosolar abundance. However, some part of this enriched metallicity may come from enrichment during planet formation in the protoplanetary disc, due to pebble or planetesimal accretion phases (Guillot & Hueso 2006; Mousis et al. 2024b). Therefore, to give an idea of the extra metallicity needed to fully explain the current C/H values of Uranus and Neptune, we start this result section by computing an estimate of the potential early enrichment (i.e. not due to late gas accretion) using the S/H ratios in Uranus and Neptune atmospheres as a good tracer of this early enrichment. Indeed, S is not expected to be released as gas (unlike CO) and it may rather be accreted together with pebbles or planetesimals when forming the planets. The S/H ratio and the corresponding C/H coming from early accretion are therefore not expected to vary after the formation of the planets, contrary to the global C/H ratio that could increase due to late gas accretion. We can then deduce the amount of late gas needed to explain observations while accounting for some potential early formation enrichment (see Section 3.1). In this section, we will also try the same calculation based on the D/H rather than the S/H to explore if we can extract more information about this early enrichment. In our simulations, we explore whether we can explain the extra potential enrichment as coming from gas released in the young primordial Kuiper belt and then spreading towards the giants with some of it captured by the different planets. Because the S/H and D/H values are not necessarily reliable tracers of early accretion (as error bars are quite large, see later), we also assess whether the observed $\left[\frac{C}{H}\right]$ values can be explained from late gas without considering any primordial enrichment. As a matter of fact, we will see that the answer is positive.

We estimate the gas production rate originating from slowly warming young KBOs above the sublimation temperature of CO to be significantly lower than that originating from collisions. Hence, to get a lower estimate of the atmospheric carbon enrichment expected for Uranus and Neptune (as well as for Jupiter and Saturn), we will continue this result section by first showing the results for the light Kuiper belt scenario (Section 3.2), i.e. a belt very similar to the current KB. We will finish this section with Subsection 3.3, which is the heart of this paper showing the main results concerning the plausible important enrichment of Uranus and Neptune by late gas when assuming a massive primordial KB.

3.1 Estimation of the C enrichment coming from early planetary formation in the protoplanetary disc

We start this results section by estimating the carbon enrichment from early accretion of planetesimals or pebbles into the protoplanetary disc using two different observables: a) the S/H ratio, and b) the D/H ratio. It is necessary to do this before presenting the results of our simulations because we want to know how much carbon enrichment is required from late gas compared with the amount that could have been supplied earlier during planet formation. The following estimations rely on the assumptions that during late gas accretion, only CO and its photodissociation products (C and O) are accreted, whereas many other elements are accreted with carbon during the planets’ formation. During planet formation, solids pollute the atmosphere with heavy elements such as carbon and sulfur, as well as deuterium. The late gas accreted later is produced by the sublimation of CO ices and we do not consider any water sublimation. We do not expect any dust accretion since Neptune is outside the dust disc, and the drag between dust and gas is small in debris discs because the mass of the secondary gas disc is much smaller than that of a protoplanetary disc. Therefore, we do not expect dust to be able to diffuse inwards with gas up to the planet locations (Takeuchi & Artymowicz 2001; Olofsson et al. 2022).

3.1.1 Using the S/H ratio to estimate carbon enrichment from planetary formation

In our model, the accreted material consists only of carbon and oxygen, since gas is produced by sublimation of CO ices. However, during their formation in the protoplanetary phase, Uranus and Neptune were also enriched in carbon as well as other heavy elements such as phosphorus (P) or sulphur (S) because of infalling planetesimals or pebbles (e.g. Guillot & Hueso 2006; Mousis et al. 2024b). This is a key difference in being able to determine the amount of carbon that may have come from early planetary formation rather than late gas. Therefore, we can use the amount of S observed in Uranus and Neptune and extract the corresponding amount of carbon associated with S at that early time. By comparing this with the current carbon enrichment, we can then access the amount of extra carbon needed beyond the protoplanetary disc phase to explain the observations, which is the value that our simulations should try to reproduce rather than the currently observed C/H ratio. This is the essence of the equations presented below.

Estimates of the S quantity are based on measurements of H₂S molecules detected at radio wavelengths deep in the atmospheres of these planets (Molter et al. 2021; Tollefson et al. 2021). Estimates of the carbon quantity are based on CH₄ abundances at ~1 bar (Sromovsky et al. 2019; Karkoschka & Tomasko 2011). CH₄, H₂S, and NH₃ are the only major equilibrium species detected in Uranus and Neptune. The NH₃ abundance is very low and probably not a reliable guide because it is removed by meteorological processes linked to water condensation (“mushballs”, see Guillot et al. 2020).

For Uranus, observations lead to a $\left[\frac{C}{H}\right]$ ratio between 44 and 74, and a $\left[\frac{S}{H}\right]$ ratio between 29 and 47. For Neptune, the $\left[\frac{C}{H}\right]$ ratio is measured to be between 55 and 92, and the $\left[\frac{S}{H}\right]$ ratio between 41 and 68 considering the case of wet adiabats (Tollefson et al. 2021; Guillot et al. 2023). The $\left[\frac{S}{H}\right]$ ratios are smaller than the $\left[\frac{C}{H}\right]$ ratios or wet adiabats, which could be a sign of late gas accretion onto Uranus and Neptune. This is not the case for dry adiabats, but since this model is less favoured in the analysis of Tollefson et al. (2021), we only use the case of wet adiabats to estimate $\left[\frac{C}{H}\right]$ and $\left[\frac{S}{H}\right]$ in this paper. From these values, an estimate of carbon enrichment due to the early formation can be obtained and we can then work out an estimate of the atmospheric carbon enrichment needed from late gas when comparing to observed C/H ratios, as we now show.

Let us take an arbitrary volume V of atmosphere. This volume contains a number n_s of sulphur atoms, n_C of carbon and n_H of hydrogen. We assume that all the hydrogen atoms come from the protosolar nebula gas, neglecting the hydrogen enrichment due to the accretion of solid materials in the protoplanetary phase. Sulphur is assumed to come from the protosolar nebula and from early accreted materials during the planetary formation in the protoplanetary disc such that $n_S^{atm}=n_S^{proto}+n_S^{ppd}$. For carbon, there is an extra term from the accretion of late gas produced in the KB such that $n_C^{atm}=n_C^{proto}+n_C^{ppd}+n_C^{KB}$. Since we assume the number of hydrogen atoms to be constant, we obtain ${\frac{S}{H}|}_{atm}=\frac{n_S^{atm}}{n_H}$ and ${\frac{C}{H}|}_{atm}=\frac{n_C^{atm}}{n_H}$ ratios or equivalently ${\frac{S}{H}|}_{atm}={\frac{S}{H}|}_{proto}+{\frac{S}{H}|}_{ppd},{\frac{C}{H}|}_{atm}={\frac{C}{H}|}_{proto}+{\frac{C}{H}|}_{ppd}+{\frac{C}{H}|}_{KB}$. Finally, when rewriting those terms relatively to the protosolar $\frac{S}{H}$ and $\frac{C}{H}$ ratios, we find $$\frac{{\frac{S}{H}|}_{atm}}{{\frac{S}{H}|}_{\text{proto }}}=\frac{{\frac{S}{H}|}_{\text{ppd }}+{\frac{S}{H}|}_{\text{proto }}}{{\frac{S}{H}|}_{\text{proto }}},$$(14) $$\frac{{\frac{C}{H}|}_{atm}}{{\frac{C}{H}|}_{\text{proto }}}=\frac{{\frac{C}{H}|}_{\mathrm{p}\mathrm{p}\mathrm{d}}+{\frac{C}{H}|}_{\mathrm{K}\mathrm{B}}+{\frac{C}{H}|}_{\text{proto }}}{{\frac{C}{H}|}_{\text{proto }}}.$$(15)

To get an estimate of ${\frac{C}{H}|}_{\mathrm{K}\mathrm{B}}$ (which is our final goal), we use the values of ${\frac{S}{C}|}_{\text{ppd}}$, the $\frac{S}{C}$ ratio of accreted material during planet formation in the protoplanetary disc, and of ${\frac{S}{C}|}_{\text{proto}}$, the protosolar $\frac{S}{C}$ ratio. We find $$\frac{{\frac{C}{H}|}_{\mathrm{K}\mathrm{B}}}{{\frac{C}{H}|}_{\text{proto }}}=\frac{{\frac{C}{H}|}_{atm}}{{\frac{C}{H}|}_{\text{proto }}}-\frac{{\frac{S}{C}|}_{\text{proto }}}{{\frac{S}{C}|}_{\text{ppd }}}(\frac{{\frac{S}{H}|}_{atm}}{{\frac{S}{H}|}_{\text{proto }}}-1)-1.$$(16)

While ${\frac{S}{C}|}_{\text{proto}}$ is measured to be 4.6 × 10⁻² (Guillot et al. 2023), we do not have precise estimates of ${\frac{S}{C}|}_{\mathrm{p}\mathrm{p}\mathrm{d}}$. Based on models of planet formation we can estimate that its value may be close to that of the protosolar abundance (e.g. Mousis et al. 2024b). Therefore, we assume $\frac{{\frac{S}{C}|}_{\text{proto }}}{{\frac{S}{C}|}_{\text{ppd }}}\sim 1$ to compute our predictions on early enrichment.

Finally, using Eq. (16), we find that ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^U={23.68}_{-26.74}^{20.04}$ for Uranus, and ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^N={21.13}_{-33.99}^{29.84}$ for Neptune. For Jupiter and Saturn, we obtain ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^J=0.64\pm 1.54$, and ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^S=1.33\pm 3.41$, respectively. The results are summarized in Figure 2.

Fig. 2 Estimations of the contribution of late gas to the atmospheric [C/H] ratios for Jupiter (J), Saturn (S), Uranus (U), and Neptune (N) calculated from the observed [C/H] and [S/H] (see Section 3.1.1). The crosses show the mean values and the bars show the uncertainties. We note that the minimal values go to zero.

3.1.2 Using the D/H ratio to estimate carbon enrichment from planetary formation

The D/H ratios of Uranus and Neptune might also be used to get another estimation of the gas-rich Kuiper belt contribution to the atmospheric carbon enrichment for Uranus and Neptune. In fact, the D/H values for Uranus and Neptune are intermediate between the protosolar and cometary values and the difference with the protosolar abundances may come from the first materials accreted (planetesimals, pebbles) in the atmosphere of these planets during the protoplanetary phase. The calculation presented here below is very similar to the S/H ratio described above. In this case, the late gas only contributes carbon and oxygen, while the early planetary formation in the protoplanetary disc contributes all the other elements, such as deuterium. From the abundance of deuterium accreted at the start of planetary formation, we can estimate the corresponding amount of carbon accreted at the same time and quantify the additional carbon accretion required once the protoplanetary disc phase is over. This additional carbon accretion is smaller than the C/H currently observed, and it is this value that must be compared with the results of our simulations.

Consider an arbitrary atmospheric volume V_atm. This volume contains $N_H^f$ atoms of hydrogen. Of these $N_H^f$ atoms, $N_H^{\text{proto}}$ come from the primordial gas, and $N_H^{\mathrm{p}\mathrm{p}\mathrm{d}}$ from the accretion of solid bodies during the protoplanetary disc phase. A fraction of this hydrogen is deuterium, but the ratio $\frac{D}{H}$ is expected to be different for the hydrogen of the primordial nebula ${\frac{D}{H}|}_{\text{proto}}$, the solid bodies ${\frac{D}{H}|}_{\text{solids}}$, and the final atmospheres ${\frac{D}{H}|}_{\text{atm}}$. Knowing their respective values, we can obtain an estimate of $N_H^{\text{ppd}}$ by calculating the number of deuterium atoms as follows $${\frac{D}{H}|}_{\text{atm }}{N}_H^f={\frac{D}{H}|}_{\text{solids }}{N}_H^{\text{ppd }}+{\frac{D}{H}|}_{\text{proto }}{N}_H^{\text{proto }},$$(17)

where we assume $\frac{N_H^f}{N_H^{\text{proto}}}\sim 1$ because the atmospheres of Uranus and Neptune are mainly composed of primordial gas despite their high metallicities (i.e. the extra amount of hydrogen from solid materials is negligible). This assumption will help simplify our equations further to reach the final contribution of early accretion to carbon enrichment.

The solids are composed of water and refractory materials with a mass ratio of $\chi =\frac{M_{H_{2}O}}{M_{\text{refr }}}\sim 1$ (Marschall et al. 2025), with M_refr the mass of refractory materials and M_H2O the water ice mass.$N_H^{\text{ppd}}$ can be separated into two distinct contributions that are $N_H^{\mathrm{p}\mathrm{p}\mathrm{d},\mathrm{w}}$ due to water ice (i.e. leading to no carbon enrichment), and $N_H^{\text{ppd,refr}}$ due to refractory materials, which produces carbon enrichment. Let $r_{ref}^w=\frac{N_H^{\text{ppd,w}}}{N_H^{\text{ppd,refr}}}$, then $N_H^{\text{ppd}}=N_H^{\text{ppd,refr}}(1+r_{ref}^w)$. Assuming that the $\frac{C}{H}$ ratio in refractory materials is ${\frac{C}{H}|}_{\text{refr }}\sim 1$ (Isnard et al. 2019; Hänni et al. 2022), we can finally get an estimate of the atmospheric carbon enrichment during the protoplanetary disc phase using that $$N_C^{\mathrm{p}\mathrm{p}\mathrm{d}}={\frac{C}{H}|}_{\text{refr}}{N}_H^{\mathrm{p}\mathrm{p}\mathrm{d},\text{refr}},$$(18)

which can be used to estimate the associated $\frac{C}{H}$ ratio $${\frac{C}{H}|}_{\mathrm{p}\mathrm{p}\mathrm{d}}\sim \frac{{\frac{C}{H}|}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{r}}}{1+r_{ref}^w}\frac{{\frac{D}{H}|}_{\mathrm{a}\mathrm{t}\mathrm{m}}-{\frac{D}{H}|}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{o}}}{{\frac{D}{H}|}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}\mathrm{s}}}.$$(19)

The measurements of the D/H ratios of Uranus and Neptune lead to ${\frac{D}{H}|}_{\text{atm }}=(4.4\pm 0.4)\times {10}^{-5}$ and (4.1 ± 0.4) × 10⁻⁵, respectively (Guillot et al. 2023). The protosolar D/H ratio is equal to ${\frac{D}{H}|}_{\text{proto }}=(1.67\pm 0.25)\times {10}^{-5}$ (Asplund et al. 2021). However, the mean D/H in accreted solids is not very well constrained with comet-like objects having D/H ratios between 10⁻⁴ and 10⁻³ (Altwegg et al. 2015). Let us use the direct measurement of the solid $\frac{D}{H}$ ratios for the comet 67P/Churyumov–Gerasimenko as a reasonable assumption. In-situ measurements lead to ${\frac{D}{H}|}_{\text{solids}}$ equal to 1.57 ± 0.54 × 10⁻³ (Paquette et al. 2021).

Now, let us estimate the value of $r_{ref}^w=\frac{N_H^{\text{ppd,w}}}{N_H^{\text{pop,refr}}}$ to finally obtain the C/H contribution from solid accretion. Let ${\mathrm{N}}_H^{H_{2}O}$ and $N_H^{\text{ref}}$ be the number of hydrogen atoms for a molecule of water and a molecule of refractory, μ_H2O, μ_ref the average molecular masses of water and refractory, respectively, and ρ_H2O, ρ_ref their respective densities. We can then write $$r_{ref}^w=\frac{N_H^{H_{2}O}\frac{{\rho }_{{\mathrm{H}}_{2}O}}{{\mu }_{{\mathrm{H}}_{2}O}}}{N_H^{\text{ref }}\frac{{\rho }_{\text{ref }}}{{\mu }_{\text{ref }}}}=\frac{N_H^{H_{2}O}}{N_H^{\text{ref }}}\frac{{\mu }_{\text{ref }}}{{\mu }_{H_{2}O}}{\chi }^{-1},$$(20)

with $N_H^{H_{2}O}=2$ and μ_H2O = 18 proton mass. We estimate $N_H^{\text{ref}}$ and μ_ref from the average solid composition C₁H_1.56O_0.134N_0.046S_0.017 measured by ROSINA for comet the 67P/C-G (i.e. μ_ref = 16.62 m_p and $N_H^{\text{ref}}=1.56$, Hänni et al. 2022).

Finally, from Eq. (19), we can compute the carbon enrichment from accreted solids in the protoplanetary disc phase to be ${\left[\frac{C}{H}\right]}_{\mathrm{p}\mathrm{p}\mathrm{d}}\sim 18$ for both planets, which is smaller than the observed values for Uranus and Neptune closer to 50–80 times the protosolar abundance (Guillot et al. 2023). From the difference, we can estimate the potential late gas extra carbon enrichment for Uranus to be ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^U\sim 40\pm 15$ and ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^N\sim 55\pm 18$ for Neptune, which is slightly higher than the values derived with the S/H method, though of the same order of magnitude. We note that in this case, we find that ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^U>25$ and ${\left[\frac{C}{H}\right]}_{\mathrm{K}\mathrm{B}}^N>37$ assuming the largest value in our range of ${\frac{D}{H}|}_{\text{solids}}$ values comprised between 10⁻⁴ and 10⁻³. This may indicate that late gas is indeed necessary if the extent of average values of ${\frac{D}{H}|}_{\text{solids}}$ is indeed within that range.

Fig. 3 CO mass production rate in the belt as a function of time for all LKB simulations.

3.2 The light (present-day) Kuiper Belt (LKB)

In the case of debris discs, it is anticipated that the viscous parameter will be somewhat elevated in comparison to that observed in protoplanetary discs, with a likely range of values between α = 10⁻³ to 10⁻¹ due to higher ionisation fraction in the former (Kral & Latter 2016; Cui et al. 2024). The fiducial simulation (LKB-fid) has an α value of 10⁻³, however, a simulation of a belt with a higher α value was also run to check the effect of this parameter (see Table 1). In the LKB-fid simulation, we fixed the accretion efficiency to be 0.5; however, we also tested an efficiency of 0.1 and 0.8 (see Table 1).

The gas production rate in the belt is independent of α and f_accr but depends on the initial belt mass, which we take equal to its current value of ~0.1 M_⊕ for the present-day KB (Vitense et al. 2012). Hence, the quantity of CO gas released in the belt is equal for all LKB’s simulations and never exceeds 6 × 10⁻⁴ M_⊕/Myr as seen in Figure 3. We find that CO is not shielded from photodissociation for the α values that were explored in our set of LKB simulations (e.g. see Figure 4 for the LKB-fid simulation where CO remains co-located with the KB). As a consequence CO gets photodissociated before being able to spread far outside its birthplace (the KB) and the gas accreted by planets is composed of neutral carbon and atomic oxygen.

In Fig. 3, we can see that there is a peak of gas production rate initially at t ~ 500 yr. This is because all the small grains lose their ices very rapidly. However, most of the ice mass is in the large solids and the majority of it is released over large timescales as can be seen in both Figs. 4 and 5. The rate of gas production throughout most of the simulations is around 10⁻⁶ M_⊕/Myr between 10–100 Myr and decreases slowly over time. This is because heat takes time to diffuse in the largest bodies (the thermal conductivity K value is much lower than in asteroids for example), where most of the mass is located given the size distribution assumed. Hence, there is still a significant amount of ice in the belt after 1 Gyr (see Figure 4).

Once produced, the gas viscously evolves and when the steady state is reached after roughly 10 Myr (i.e. when gas had time to spread both inwards and outwards), we find that there are no significant differences in the gas mass accreted by the different giant planets as a function of time (see Figure 5). This might appear counter-intuitive because Jupiter and Saturn being located farther out from the KB than the icy giants, one may expect that they receive less gas as they may get starved by Uranus and Neptune. However, given their higher masses and corresponding Hill radii, they also have greater incoming mass fluxes (see Eq. (8)). Hence, the greater degree of carbon enrichment due to late gas accretion for the atmospheres of Uranus and Neptune observed in our simulations (Fig. A.2) over that of Jupiter and Saturn (Fig. A.3) is due to their smaller atmospheric masses (in which carbon can mix) rather than more mass accreted.

Changing α only affects the time to reach steady state. However, the final accreted masses remain similar for all values of α tested (see Figure A.1). It is expected because α does not influence the amount of gas mass released from the belt nor its accretion efficiency onto planets.

We do not observe a significant neutral carbon enrichment/depletion due to interspecies diffusion (see Fig. 4). The C/O ratio of the accreted gas is then always equal to 1 in our simulations, since we assume that only CO ices sublimate. We note that it is unclear whether some CO₂ ices may also sublimate in the KB in addition to CO, which would then decrease the gas disc C/O ratio (De Prá et al. 2024).

In Figure A.1, we can observe that there is no direct proportionality between the accreted mass and the accretion efficiency f_accr. The accreted mass for the icy planets in the LKB-fid simulation (with f_accr = 0.5) is about 20 times higher (i.e. much greater than five times) than the accreted mass in LKB-3 (with f_accr = 0.1). As can be seen in the same figure, the icy planets demonstrate a greater capacity to accrete mass compared to giants when the efficiency of the accretion process is high. This is because, in the case of a higher accretion efficiency, the gas giants are deprived of gas that is accreted onto the icy planets, which results in a lower accreted mass than in the fiducial simulation.

Finally, regardless of the circumstances, we find that for the LKB simulations the accreted masses on both Uranus and Neptune turn into a small C/H value when compared to the observed values. For instance, for our LKB-fid simulation, we can anticipate a maximum of one percent of the protosolar C/H ratio in the atmospheres of icy planets due to accretion of late gas from a light Kuiper Belt (see Figure A.2). The results obtained for the other simulations are analogous, given the comparable accretion rates (See Table A.1).

Fig. 4 Surface density of CO (in variations of blue) and neutral carbon (in variations of orange) at different times for the LKB-fid simulation. The colour goes from darker to brighter for increasing time, as indicated by the colour bar. The grey vertical bars represent the planets’ accretion zones (from left to right: Jupiter, Saturn, Uranus, and Neptune).

Fig. 5 Neutral carbon accreted mass for each planet as a function of time for the LKB-fid simulation.

3.3 The heavy (primordial) Kuiper belt (HKB)

The primordial Kuiper belt is expected to have been much more massive than the current KB with a much higher gas production rate than that predicted in the previous subsection. To evaluate the efficiency of gas accretion onto the giants from gas released in the young primordial KB, a series of simulations were conducted. The predictions for the accretion of gas onto the icy giants will vary depending on the assumed primordial KB model (as described in Section 2.6), with the most important parameters being the initial mass of the KB (5, 20, 50 Earth masses), its lifetime (10 to 100 Myr), and the system’s architecture (compact or extended). The different setups of our heavy KB (HKB) simulations are listed in Table 2. The final $\left[\frac{C}{H}\right]$ for all simulations are summarized in Table 3.

3.3.1 The compact architecture

The compact architecture represents a scenario akin to that of the Nice model (Tsiganis et al. 2005; Levison et al. 2011; Griveaud et al. 2024). The planets begin in a compact configuration, with a substantially massive belt of 20–50 M_⊕. The outer planets migrate outwards, resulting in the depletion of the primordial belt. This is the most comprehensive model to date of the evolution of the young Solar System, capable of reproducing both the current configuration of planets and the present-day Kuiper belt. Accordingly, the characteristics of our HKB-fid, 2, and 3, simulations have been selected to align with this architecture (see Section 2.6).

In our fiducial simulation, probably the most realistic, we have selected a relatively substantial primordial belt of 50 Earth masses in agreement with recent work showing that assuming a more realistic low viscosity for the protoplanetary disc leads to a more massive primordial KB (Griveaud et al. 2024). In Figure 6 (bottom), we can see that the HKB-fid simulation leads to a total CO mass released greater than 0.1 M_⊕ in roughly 10 Myr. On the same figure (top), we see that the gas production rate reaches values greater than 0.1 M_⊕/Myr, i.e. orders of magnitude higher than for the LKB simulations. For the HKB-3 simulation with a belt mass of 20 Earth masses, the total CO mass released and gas production rates are lower but of the same order of magnitude.

In Figure 7, we see that for the HKB-fid, 2, 3 simulations, CO gets shielded from the impinging interstellar radiation. Hence, CO has time to spread viscously within the system and it can reach the giants in molecular form. Therefore, the accreted gas onto the icy planets is mostly made up of CO rather than atomic carbon and oxygen as was the case for the LKB simulations explored in the previous section. However, CO is not stable in the upper layers of the atmospheres of Uranus and Neptune, and it probably reacts with the primordial hydrogen in the atmospheres to create CH₄ and H₂O (Atreya et al. 2020). The extra carbon brought via late gas could then be seen when observing the C/H ratio in the atmospheres. As discussed later, we note that the O/H ratio is not very informative because water is not only confined to the atmospheres and the interpretation of the O/H ratio is more complex.

For the HKB simulations, the α value was fixed to 10⁻³. We do not vary it as it does not have an important impact on the final predictions of the captured gas mass by planets as previously explored for the LKB. However, there is still plenty of gas at the onset of Neptune’s depletion of the primordial belt (second phase of the simulation, see Figure 7). We do not model the exact planet’s dynamic crossing the disc during its outward migration, and the accretion rates obtained after this time might then be over or underestimated. However, to account for that, we explore the final masses we obtain for different belt lifetimes going from 10 to 100 Myr, which can lead to a factor 2 difference in total mass captured onto planets.

Figure 8 shows the CO mass captured onto the 4 giants for our HKB-fid, 2, 3 simulations. Note that the HKB-fid and HKB-2 only differ by their f_accr values of 0.5 and 0.1, respectively, and the HKB-3 tests a less massive belt of 20 M_⊕. Figure 8 shows that the CO mass accreted onto icy giants is greater than 10⁻² M_⊕, which may be significant when injected into the atmospheres of Neptune and Uranus. Indeed, assuming that the accreted CO only mixes in the atmospheres of Uranus (1.25–3.5 M_⊕, Guillot et al. 2023), and Neptune (1.6–4.15 M_⊕), we can easily compute the corresponding [C/H] ratio.

Figure 9 illustrates that after 10 Myr of late gas accretion, the [C/H] ratio of Uranus reaches 5 to 13, while for Neptune it reaches 9.5 to 25. These are the values only accounting for late gas accretion and assuming that no enhanced primordial metallicity was present. We note that if the gas disc survives 100 Myr instead, the resulting [C/H] ratio is 13.5 to 33 for Uranus and 20 to 58 for Neptune. We remind the reader that observations show that [C/H] is of 44–74 and 55–92 for Uranus and Neptune (Atreya et al. 2020), i.e. our fiducial model may not be far from potentially explain it all by itself.

However, some primordial enhancement of C/H may be present due to the planet formation process (e.g. pebble accretion) and a smaller C/H may be needed to account for observations as explained in the beginning of the result section. To account for that, in Figure 9 we plot two horizontal lines showing the maximum extra carbon enrichment needed to explain current observations with the primordial enrichment calculated using measurements of S/H. Our calculations show that Uranus should accrete on average an extra [C/H] of 24 and Neptune 21. However, given uncertainties, we can only reliably derive upper values of extra [C/H] of 51 for Neptune and 44 for Uranus (see Section 4 for more details), which are plotted on Figure 9. We conclude that the fiducial simulation is able to explain both observed [C/H] values of Uranus and Neptune, with most of the metallicity enrichment coming from accretion of late gas.

The HKB-fid simulation represents the most favourable case for the accretion of carbon-rich gas onto icy planets, characterized by the most massive belt and the highest accretion efficiency. The HKB-2 simulation was run with an accretion efficiency five times smaller than the HKB fiducial value to test the impact of f_accr on our results. It leads to a C/H ratio lower than that of the fiducial simulation, by a factor of between 1.4 and 9.2 for Uranus and 1.5 and 11.2 for Neptune, after 100 million years.

The HKB-3 simulation was run with a less massive belt because even though recent Nice model simulations (Griveaud et al. 2024) suggest it may have reached 50 M_⊕, the collisional models of the belt by Bottke et al. (2023) indicate the possibility of a slightly lighter belt of about 20 M_⊕. We find results very similar to the HKB-2 simulation as seen in Figure 10. The C/H values reached by both simulations are still within the right range to explain plausible extra carbon-enrichment from late gas in Uranus and Neptune, but although significant, it would become a minor source compared to primordial enrichment.

Furthermore, our model also allows us to assess the quantity of carbon-rich gas that can accumulate on the gas giants Jupiter and Saturn. Despite the similarity in absolute accreted mass, the atmospheric mass in which the accreted gas will mix is significantly higher. We account for the effect of dilution of the solid core for Jupiter and Saturn that can change the atmospheric mass into which carbon can mix with accretion. We consider lower limits for the atmospheric mass using a diffuse core of 60 M_⊕ for both Jupiter and Saturn, and upper limits of 20–30 M_⊕ for the solid cores of Saturn and Jupiter, respectively.

The contribution of late gas released in the primordial Kuiper belt to the [C/H] ratio of Jupiter and Saturn is, at most, of 0.16 and 2.0, respectively (see Figure 11). If we consider a dilute core of 50% of the total mass for Jupiter the [C/H] is at most 0.26. However, the observed [C/H] ratios are also much lower for Jupiter and Saturn (Guillot et al. 2023, respectively 3.56 ± 0.86 and 7.50 ± 0.32), hence this contribution might also be non-negligible at least for Saturn. Indeed, the upper limits for the contribution of late gas (calculated using the S/H ratios as a tracer of accretion from planet formation) are 2.19 and 4.7 for Jupiter and Saturn, respectively.

We conclude that our compact configuration simulations (mimicking the Nice model) provide very compelling results and show that late gas could be the major contributor to the observed C/H values in both Uranus and Neptune in the best-case scenario, or at least lead to non-negligible enhancement in the worst-case scenario, which is still needed to be accounted for. We will discuss the consequences in a broader context in Section 4.

Fig. 6 Top: CO mass production rate in the belt for the HKB simulations. Bottom: Total mass of CO gas released in the belt for the HKB simulations. The initial mass of 5 M⊕ (HKB-4, 5) is for the extended configuration, whereas the initial masses of 20 (HKB-3) and 50 M⊕ (HKB-fid, 2) are for the compact configuration.

Fig. 7 Surface density of CO (in variations of blue) and neutral carbon (in variations of orange) at different times for the HKB-fid simulation. The colour goes from darker to brighter for increasing time as indicated by the colour bar. The grey vertical bars indicate the planets’ accretion zones (from left to right: Jupiter, Saturn, Uranus, and Neptune).

Fig. 8 CO accreted mass for each planet as a function of time for the compact configuration of HKB simulations (HKB-fid left, HKB-2 middle, HKB-3 right).

Fig. 9 [C/H] for the HKB-fid simulation as a function of time. The filled areas correspond to predictions for Uranus (purple) and Neptune (blue). The uncertainties are shown via the extension of the filled areas and are due to uncertainties on the respective atmospheric masses. The solid and dashed black lines are the maximum estimation of the gas disc contribution obtained from observations for Uranus and Neptune (see Section 3.1.1).

Fig. 10 [C/H] for the simulation HKB-2 (top) and HKB-3 (bottom) as a function of time.

3.3.2 The extended architecture

In the compact architecture context, the planets and KB need to start in a compact configuration with a massive planetary belt to end up with the present-day Solar System configuration at the end of the system’s evolution. However, some mechanisms allow the planets to migrate earlier, directly in the protoplanetary disc phase as, for instance, the ‘Rebound’ model by Liu et al. (2022). In this case, once the protoplanetary disc fully dissipates, the system has its planets and belt at locations very similar to the present-day Solar System, i.e. more extended than in the compact configuration. In this context, the belt mass can be smaller (e.g. 5 M_⊕) while still allowing the current Solar System architecture to be obtained (Liu et al. 2022).

We ran two simulations mimicking the extended architecture approach with a 5 M_⊕ primordial belt and varying belt lifetimes (10 Myr for HKB-4 and 100 Myr for HKB-5, see Table 2). The latter is important because it is not yet clear what is the exact dynamics of belts in the extended configuration and it may take longer to deplete that in the compact configuration (Liu et al. 2022).

Figure 12 shows the C/H predictions for late gas capture on Uranus and Neptune for the HKB-4 and HKB-5 simulations. As expected, because of the lower initial belt mass, we end up with smaller C/H ratios than the HKB-fid simulation. In this case, irrespective of the lifetime of the belt, the contribution of late gas to the atmospheric carbon enrichment for Uranus and Neptune is small compared to observed ratios though not negligible.

We conclude that the extended configuration (mimicking other models than the Nice model) does not lead to a significant contribution to the metallicity of Uranus and Neptune. We note, however, that if it is confirmed by future observations that extra metallicity is needed to explain the C/H in both ice giants, our work could be used as an argument in disfavour of those extended models, thus reinforcing the Nice model.

Fig. 11 [C/H] for the simulation HKB-fid as a function of time for Jupiter and Saturn. The uncertainties are shown by the extent of the filled-in areas and are due to the uncertainties over the respective atmospheric mass in which the CO can mix (i.e. between a diluted 60 M⊕ diffuse core for both planets and a 20, 30 M⊕ solid core for Saturn and Jupiter, respectively).

4 Discussion

4.1 The O/H ratio

Together with carbon, some oxygen will be accreted onto the planets, either in the form of CO, or directly as atomic oxygen. We then naively expect the atmospheres of Uranus and Neptune to be enriched in oxygen compared to solar abundances. However, water condenses deep into the atmosphere (e.g. Cano Amoros et al. 2024), and we only have an indirect constraint on the O/H abundance through the detection of CO, and then using models of chemical reaction rates converting CO into CH₄ and accounting for upward transport of this CO (Lellouch et al. 2005; Cavalié et al. 2014; Venot et al. 2020). Therefore, no reliable O/H measurement can be used to estimate the amount of solids accreted during planet formation or the contribution from late gas, and we do not explore this aspect further though it may be different for exoplanets with different temperatures as we discuss later in Subsection 4.6.

Fig. 12 [C/H] for the simulation HKB-4 (top) and HKB-5 (bottom) as a function of time.

4.2 The primordial belt model to explain the metallicity of Uranus and Neptune

4.2.1 Gas viscosity

The α viscosity in debris discs is expected to be rather high compared to protoplanetary discs with values that can reach 10⁻³−0.1 both because of the magnetorotational instability (Kral & Latter 2016) or other instabilities, or even via the molecular viscosity which may not be negligible in debris disc conditions (Cui et al. 2024). It is thus expected that the gas will spread from the KB to Neptune in t_vis = ΔR²/ν ~ 0.3 to 30 Myr for α = 0.1 to 10⁻³. As one can see the viscous timescale decreases linearly with α and the only difference in our simulations will thus be the time to reach planets but the accreted mass will remain the same as shown by the LKB-2 simulation and explained in further detail in Kral et al. (2024). Hence, in the long term, we do not expect α to play a major role for the final accreted mass of gas onto planets.

We note that one difference it could make is on the composition of the gas disc. Indeed, if α is small, the gas disc will spread more slowly and become denser (Kral & Latter 2016; Kral et al. 2019) and some shielding may appear (e.g. from carbon or CO). However, it is not very relevant for our study because the final quantity of C and O accreted into the atmospheres of Neptune and Uranus will not change whether it is captured as CO or C+O. This gas mass will then be turned into CH₄ and water into the atmospheres of Uranus and Neptune as explained in Subsection 3.3.

4.2.2 Effects of the belt’s mass

The gas production rate we use is based on the solid mass loss rate in the KB, which we estimate using Eq. (38) in Löhne et al. (2008). This collisional model accounts for varying slopes in the size distribution with e.g., a primordial slope (for bodies that had not had time to collide significantly yet), and then two more slopes to account for the gravity and strength regimes, respectively. We note that even though simplified one-slope collisional models show that the mass loss rate scales as the square of the belt mass, it is indeed more complicated in our case. Numerical calculations show for instance that the HKB-3 belt of 20 M_⊕ leads to less than 2 / 5 of the gas mass produced by the HKBfid belt of 50 M_⊕, due to this non-linearity in our model. As the belt mass increases, the collision timescale decreases as it typically scales with the inverse of total belt mass (Wyatt 2008). Increasing the belt mass leads to more collisions, faster, and thus more gas produced when the system is younger. This explains that models with the less massive belts initially (e.g. HKB-5) take longer to release their gas content (e.g. see Fig. 12).

It is therefore essential to pay attention to the exact size distributions and belt mass given how they may change the final gas production rates.

The outer planets might dynamically heat the belt, especially given the potential outward migration of Neptune (Kenyon & Bromley 2004; Nesvorný 2018). Therefore, to account for this effect, we assume a mean eccentricity of 0.25 for the planetesimals in the belt. The rate of mass production in such a belt in collisional evolution is typically proportional to e^5/3 (Wyatt 2008). Thus, a colder belt with a typical eccentricity of 0.1 will produce approximately 5 times less gas than assumed here. A belt with very high dynamical heating and a typical eccentricity of 0.5 will produce approximately 3 times more gas than in our fiducial model. We note that this collisional model can explain most observations (Kral et al. 2017, 2019) but it would be even more satisfying to develop a physical model able to predict directly the amount of gas released from physical mechanisms. However, more data are needed and this is work left for the next decade, which will allow us to refine our predictions.

4.2.3 Effects of the accretion efficiency

The accretion efficiency of a planet capturing gas is susceptible to two effects: (1) thermal cooling of the atmosphere (Piso & Youdin 2014; Piso et al. 2015; Lee & Chiang 2015; Kral et al. 2020), and (2) hydrodynamical effects (Lubow & D’Angelo 2006; Mordasini et al. 2012; Tanigawa et al. 2012; Ormel et al. 2015; Lambrechts et al. 2019). The analytical model developed by Kral et al. (2020) shows that in a debris disc environment cooling is very efficient. If we apply their Eq. (5)¹ for Uranus and Neptune, we obtain a maximum theoretical capture rate onto the planets that is higher than accretion rates obtained from the simulations ran in this paper. This means that theoretically, all incoming gas can be captured if only cooling is accounted for as shown in Figure 13 for the HKB-fid simulation, which leads to the highest accretion rates. We note that the theoretical prediction is analytical and does not account for subtle effects such as the decoupling between the disc and the planet’s atmosphere that may happen at low gas densities such as those found in debris disc environments (Mordasini et al. 2012). Implementing numerical simulations of cooling in debris disc conditions is not trivial and is left for future work.

As stated earlier, the other physical mechanism that may impact the theoretical accretion rate is the complex hydrodynamical flow. Previous numerical hydrodynamical simulations in protoplanetary disc conditions show that the accretion efficiency (f_accr) accounting for hydrodynamic effects is at least 0.1 and most likely higher than 0.5. For instance, 2D hydro simulations by Lubow & D’Angelo (2006) provide f_accr values between 0.38 and 0.88. We note that we derive their accretion efficiency from the ratio ${\dot{M}}_i/{\dot{M}}_p$ they compute. Those results are confirmed by 3D hydrodynamical simulations (e.g. Schulik et al. 2019; Li et al. 2023; Bailey et al. 2024).

Gas recycling may be important in protoplanetary discs (Ormel et al. 2015; Moldenhauer et al. 2022) but it is tightly linked to the cooling timescale, which is very small in our case and this may not operate. Dedicated hydrodynamic simulations in debris disc conditions would provide further insights concerning the right range of values to use for f_accr even though it is most likely of the same order of magnitude as in protoplanetary discs.

We note that gap opening may also impact accretion rates. It is expected that Uranus and Neptune are not massive enough to open a gap (Bergez-Casalou et al. 2020) but Saturn and Jupiter would. In this case, accretion rates would be slightly lowered as shown in Bergez-Casalou et al. (2020), which should be accounted for in future hydrodynamical simulations. Another effect, which should be less important would be to consider the migration of Neptune through the gas disc as it accretes. Given that the migration timescale is small compared to the accretion timescale, we do not expect a large change in our conclusions.

For all those reasons, we consider the simulation HKB-2 with an accretion efficiency of 0.1 to most likely provide a lower bound on the atmospheric carbon enrichment for the icy planets. In that worst-case scenario, we find [C/H] values of 9 and 13 (at most) at the end of the simulation for Uranus and Neptune, respectively. We note that after 10 Myr the carbon enrichment is already non-negligible because accretion is much faster at the beginning of the simulation. We consider that it reinforces our conclusion that late gas may have a significant impact on the metallicity of giant ices and may turn out to be important in extrasolar systems as well as discussed in Subsection 4.6.

Fig. 13 Accretions rates for the HKB-fid simulations compared to the maximum capture rate onto planets due to atmospheric cooling.

4.3 Using our results to constrain the formation model of the Solar System

Given the uncertainties on the contribution of early enrichment during planet formation on the observed C/H ratios of Uranus and Neptune, it is not yet possible to give firm conclusions concerning the need for late gas enrichment and hence infer the initial belt mass needed. However, assuming the central values of the derived extra C/H as more likely, we see that an extra [C/H] of ~20 would be needed for both Uranus and Neptune. Using our model, we see that we need belts with masses greater than 20–30 M_⊕ to be able to reproduce this order of magnitude. If confirmed, it would be an indication that the heavy belt needed in the Nice model to reproduce the current Solar System is needed.

Another constraint that may be relaxed is where the ice giants formed relative to the CO ice line. Indeed, current models postulate that the two planets formed at the CO ice line within the protosolar nebula to be able to explain the high C/H ratio (Irwin et al. 2019; Sromovsky et al. 2019). However, if the carbon contribution comes from late gas, this constraint can be lifted, which could be interesting to consider in future models of planet formation.

4.4 On the results concerning Saturn and Jupiter

One of our results is that the carbon enrichment for Saturn is higher than that for Jupiter (see Figure 11). This difference between the two planets is greater for higher accretion efficiencies. This is expected because for a higher accretion efficiency, there is less gas arriving at Jupiter, thus reducing its accretion rate. There is also a dilution effect due to the highest atmospheric mass for Jupiter. In any case, we find that the predicted atmospheric [C/H] ratios due to late gas accretion in the HKB-fid simulation are 2.0 and 0.16 for Saturn and Jupiter, respectively, in agreement with the potential maximum extra C/H needed from late gas accretion of 4.7 and 2.2 (using the S/H ratios as presented in Section 3.1.1), respectively. Thus, we find that observations of giant planets are compatible with late gas capture, but given uncertainties on the ${\left[\frac{C}{H}\right]}_{ppd}$ coming from accreted solids materials during the protoplanetary phase, it is hard to use this information to constrain the disc properties such as its mass or lifetime, similar to Uranus and Neptune as explained in Section 4.3.

4.5 Late gas accretion onto the inner terrestrial planets

Up until now, we have focused on the accretion of the four giant planets in the Solar System. However, one notices that in our simulations the gas disc reaches the inner terrestrial planets (see for example Figure 7 for the HKB-fid simulation). Even though we do not include those inner planets directly, we can post-process the amount of CO arriving at the planets’ Hill spheres and how much gets accreted using the same technique as described for giant planets. Only accounting for the Earth (i.e. neglecting Mars), we derive a maximum carbon mass that could be accreted onto it in the HKB-fid simulation, that is ~10⁻³ M_⊕ in 10 Myr and ~10⁻² M_⊕ in 100 Myr.

However, we do not expect such a high carbon mass accreted for the Earth. In fact, our simulations do not model any gap formation but Jupiter and Saturn are much larger than Uranus and Neptune and they could both open gaps (Bergez-Casalou et al. 2020). Accounting for this effect would reduce the mass flux crossing both Jupiter and Saturn’s orbits, thus reducing the potential flux reaching Earth’s orbit. We note that the accretion efficiency could also be higher than the assumed 0.5, starving the inner planets, similar to what happens to Jupiter and Saturn in the simulation LKB-4. We also note that in most simulations, gas is accreted onto planets very early, roughly within 10 Myr of protoplanetary disc dissipation where Earth has not yet accreted its final mass (Kleine & Walker 2017; Clement et al. 2018) and is still very hot, thus reducing the accretion efficiency (Piso & Youdin 2014; Piso et al. 2015; Lee & Chiang 2015). Another effect that may be of importance for the Earth is that of giant impacts (Schlichting et al. 2015) including the moon-forming impact (Kleine & Walker 2017) that may remove part of the accreted CO after it has been accreted. Due to those uncertainties, we decided not to include terrestrial planets in our model and do not attempt to conclude further on late CO gas accretion to terrestrial planets. However, we note that in extrasolar systems with no massive giants, late gas may easily reach those inner planets and alter their atmospheres, which could be an important factor to account for to obtain a proper prediction of their final atmospheric makeup.

4.6 Application to extrasolar systems

This paper focuses on the Solar System, however, similar conditions may happen in extrasolar systems. Indeed, it is now accepted that exo-Kuiper belts releasing CO gas may be the norm rather than the exception (Moór et al. 2017). Released gas will also spread inwards and outwards of the belt and deliver CO or C+O gas to exoplanets in the system. Our mechanism seems universal and would naturally increase exoplanet metallicities in carbon and oxygen.

We note that giant planets of Jupiter-mass or greater can be detected in the outer regions of planetary systems using the direct imaging technique (e.g. Vigan et al. 2021). Some debris disc systems have giants detected, such as the archetypal β Pic system with its two planets (Lagrange et al. 2009, 2019). For instance, it is interesting to note that some long-period planets detected, such as GJ 504 b or AF Lep b, have super-stellar metallicities, i.e. they appear to be enriched in CH₄. Those planets are more massive than Jupiter and we expect our scenario to only provide little carbon compared to their global masses, unless their primordial exo-KBs were significantly more massive than in our Solar System, or the atmosphere in which carbon can mix are relatively light. However, the effect on sub-Jupiter planets is expected to be much more significant and one can expect that late gas may contribute to the observed C/H in their atmospheres. Using our simulations, we can extrapolate the extra C/H ratio for an exoplanet as a function of its atmospheric mass (see Figure 14).

For exoplanets, the question of the magnitude at which the different elements will exhibit an enrichment or not thus depends on the mass of the fluid envelope and that of the outer planetesimal belt as shown in Fig. 14 but also on the atmospheric temperature and entropy, something that is essentially linked to the planet’s distance to the star. As we have seen the evaporation of CO gas from an outer belt should lead to a global increase of the deep atmospheric C and O abundances, but whether these are detectable depends on the atmospheric temperature itself. A strong difference between exoplanets and Solar System giant planets is that for the latter we can probe their deep interior with measurements at microwave wavelengths that probe to tens or hundreds of bars (de Pater et al. 1991; Bolton et al. 2017), something impossible for exoplanets. For planets too hot for water clouds to form, both C and O should be enriched compared to other species such as N and S. For cooler exoplanets, we expect that the oxygen signature should be suppressed at photospheric temperatures lower than about 250 K because of water condensation. For exoplanets with hydrogen atmospheres, the next major condensate is NH₄SH, which, mostly removes S from the upper atmosphere, due to its lower abundance for a solar-composition gas. This should occur around 150 K. For still lower atmospheric temperatures, close to Jupiter’s effective temperature (124 K), NH₃ condenses and also removes N from the atmosphere. At this temperature and lower, interaction between water ice crystals and ammonia vapor leads to the formation of “mushballs” which can transport N even deeper (Guillot et al. 2020). When the photospheric temperature has decreased between Saturn and Uranus’s values, i.e. between 95 and 60 K, we expect that the transport by mushballs should decrease the abundance of N in the upper atmosphere below that of S (Guillot 2021). This effectively implies that the NH₄SH condensation is not removing all of the S atoms, possibly leaving a possibility to observe it in the upper atmosphere. However, the condensation of H₂S, expected for planets with Uranus’ effective temperature and below (60 K) (Irwin et al. 2018) again removes S from the upper atmosphere, with H₂S which remains observable at depths for ice giants in the Solar System using microwave wavelengths but would not be observable in an exoplanet. As can be seen, measuring a specific enrichment in carbon in exoplanets will require dedicated studies and a large sample.

As a conclusion, it appears that the simplest sample to study would be exoplanets with photospheric temperatures higher than 250 K. For these it is be expected that C, O, N and S would be observable in the upper atmosphere. In fact, observations of H₂S have already been made by the JWST for the super-Earth L 98-59 d (Banerjee et al. 2024). Furthermore, HCN can also be observed with current observatories (e.g. the detection of HCN in HR 8799 c using VLTI/Gravity, Nasedkin et al. 2024). With such observations of C, O, N, and S available, and likely to increase in number over the next few years, one can expect to be able to quantify the amount of carbon accreted in the debris disc stage for a given planet. This is the way forward to test late gas accretion and decide on whether it could be the dominant mechanism that increases metallicity in planets.

The direct detection of less massive exoplanets at large distances from their stars will also grow thanks to the JWST sensitivity (using NIRCam and MIRI). Indeed, it may be possible to go down to Saturn-mass exoplanets or even lower for the closest planets to Earth reaching close to Neptune’s mass (Beichman et al. 2019). When the first results arrive for those lower-mass planets, it will be interesting to see whether there is indeed a trend in metallicity and whether it may be linked to late gas accretion. However, we note that the limits in terms of detecting lines (e.g. water, methane, CO) for directly imaged planets with masses less than Jupiter are not yet clearly defined. Only future work or observations with the JWST will help us to draw conclusions (e.g. Mâlin et al. 2023). The ELT is predicted to be able to image and characterize Saturn to Neptune-like exoplanets, notably with METIS, a first generation instrument (Brandl et al. 2021). It is also expected that second generation instruments such as PCS will be able to go down to exo-Earths to look for biosignatures and detecting Neptune-like planets may become routine (Kasper et al. 2021).

On a longer timeframe, the Habitable Worlds Observatory (HWO, National Academies of Sciences 2023) and the Large Interferometer for Exoplanets (LIFE, Quanz et al. 2022) will be able to characterize those planets and further test the contribution of late gas as well as decipher any trend that may be present in terms of metallicity with the different characteristics of the systems.

Fig. 14 Estimation of the extra C/H due to late gas accretion onto an exoplanet as a function of its atmospheric mass for a belt of respectively 50 M⊕ and 5 M⊕ (purple and blue) based on the HKB-fid and HKB-4 simulations. The left axis shows the atmospheric C/H ratio relative to the protosolar abundance, and the right axis shows the absolute atmospheric C/H ratio. The red line represents the protosolar C/H ratio. The blue shaded area represents Neptune’s atmospheric mass.

5 Conclusion

It is now well known that gas (mostly CO) is released in exo-Kuiper belts around other stars. The gas originates from CO ices that sublimate over time from exo-Kuiper belt objects. It is likely that gas was also released in the young Kuiper belt in our Solar System. In this paper, we investigate the consequences of such gas release in the Solar System’s youth on the atmospheres of Uranus and Neptune. Gas released at the Kuiper belt would have spread inwards towards the giants, which would have captured some of it. We have developed a numerical model capable of calculating the final masses of gas accreted onto Jupiter, Saturn, Uranus, and Neptune, assuming that the Solar System formed in a similar way to that proposed by the Nice model, but also for other scenarios that may explain the formation of the Solar System.

Assuming a massive primordial Kuiper belt similar to what is expected from state-of-the-art models, we find that up to ~0.1 M_⊕ of CO gas can be accreted onto each giant planet. The gas mixes within the upper atmospheres of giants and we find that it can totally change their primordial metallicities, in terms of observed C/H ratios. The effect is most impressive for Uranus and Neptune, as our model can explain [C/H] values on the order of 20 and 50, respectively, as coming from the capture of late gas. It means that our model may explain by itself the observed values of [C/H] on the ice giants.

However, some of the observed [C/H] may originate from initial planet formation (e.g. accretion of planetesimals or pebbles). Using the S/H ratio as a tracer of this early accretion, we find that the extra [C/H] that may be needed to explain current observations is on the order of 20 for both Uranus and Neptune, which aligns very nicely with the results of our late gas capture model. We find that such a contribution can be provided with a primordial Kuiper belt more massive than 20–30 M_⊕, similar to that proposed in the Nice model. We also find that constraints on the [C/H] in Jupiter and Saturn are compatible with such Nice-model scenarios. We conclude that late gas may play an important role in shaping the atmospheres of Uranus and Neptune, as its contribution to the observed [C/H] may be dominant over early accretion during planet formation and is far from negligible even in worst-case scenarios.

Furthermore, although this article has been devoted to the study of the implications of late gas in the Solar System, the late gas mechanism we study is universal. Indeed, exo-Kuiper belts are present around more than 1/4 of stars and most of them could have possessed a gaseous component in their youth. Hence, exoplanets in planetary systems suffer gas accretion (CO or C+O), which inevitably alters their [C/H] ratios. We suggest some ways forward to study the effect of late gas using currently accessible data (e.g. CH₄, HCN, and H₂S observations in sub-Neptunes or super-Earths) of sub-Jupiter planets close to their host star. Future observations of distant sub-Jupiter planets using direct imaging are challenging but may provide even more evidence of this late gas accretion. The tip of the iceberg is offered by the JWST, and then more refined imaging from the ground will be possible with the ELT, or with space missions like HWO or LIFE.

Acknowledgements

We thanks Philippe Thébault and the anonymous referee for interesting discussions on this paper.

Appendix A Extra figures and table

Fig. A.1 Neutral carbon accreted mass for each planet as a function of time for the simulations LKB-fid, 2, 3 and 4 (from left to right and top to bottom). The only differences between those simulations are the viscous parameter α and the accretion efficiency faccr (see Table 1).

Fig. A.2 [C/H] for the LKB-fid simulation as a function of time for the ice giants. The uncertainties represented by the filled areas are due to the uncertainties on the atmospheric masses of Neptune and Uranus.

Fig. A.3 C/H to protosolar ratio for the LKB-fid simulation as a function of time for the giant planets. The uncertainties represented by the filled areas are due to the uncertainties on the atmospheric masses of Jupiter and Saturn that may mix with carbon.

References

All Tables

All Figures

Fig. 1 Schematic presenting a description of our set of nine simulations. The heavy belt cases are at the top (compact configuration) and middle (extended configuration), and the light cases are at the bottom. Jupiter, Saturn, Uranus, and Neptune are designated by J, S, U, and N, respectively. The different simulations differ in terms of the locations of planets, accretion efficiency, the viscosity of the gas, and the belt’s mass. For the heave belt configurations, the time when depletion starts is also a parameter. Each parameter has its own symbol and colour as can be seen in the bottom part of the image.

In the text

	Fig. 2 Estimations of the contribution of late gas to the atmospheric [C/H] ratios for Jupiter (J), Saturn (S), Uranus (U), and Neptune (N) calculated from the observed [C/H] and [S/H] (see Section 3.1.1). The crosses show the mean values and the bars show the uncertainties. We note that the minimal values go to zero.
In the text

	Fig. 3 CO mass production rate in the belt as a function of time for all LKB simulations.
In the text

	Fig. 4 Surface density of CO (in variations of blue) and neutral carbon (in variations of orange) at different times for the LKB-fid simulation. The colour goes from darker to brighter for increasing time, as indicated by the colour bar. The grey vertical bars represent the planets’ accretion zones (from left to right: Jupiter, Saturn, Uranus, and Neptune).
In the text

	Fig. 5 Neutral carbon accreted mass for each planet as a function of time for the LKB-fid simulation.
In the text

	Fig. 6 Top: CO mass production rate in the belt for the HKB simulations. Bottom: Total mass of CO gas released in the belt for the HKB simulations. The initial mass of 5 M⊕ (HKB-4, 5) is for the extended configuration, whereas the initial masses of 20 (HKB-3) and 50 M⊕ (HKB-fid, 2) are for the compact configuration.
In the text

	Fig. 7 Surface density of CO (in variations of blue) and neutral carbon (in variations of orange) at different times for the HKB-fid simulation. The colour goes from darker to brighter for increasing time as indicated by the colour bar. The grey vertical bars indicate the planets’ accretion zones (from left to right: Jupiter, Saturn, Uranus, and Neptune).
In the text

	Fig. 8 CO accreted mass for each planet as a function of time for the compact configuration of HKB simulations (HKB-fid left, HKB-2 middle, HKB-3 right).
In the text

Fig. 9 [C/H] for the HKB-fid simulation as a function of time. The filled areas correspond to predictions for Uranus (purple) and Neptune (blue). The uncertainties are shown via the extension of the filled areas and are due to uncertainties on the respective atmospheric masses. The solid and dashed black lines are the maximum estimation of the gas disc contribution obtained from observations for Uranus and Neptune (see Section 3.1.1).

In the text

	Fig. 10 [C/H] for the simulation HKB-2 (top) and HKB-3 (bottom) as a function of time.
In the text

	Fig. 11 [C/H] for the simulation HKB-fid as a function of time for Jupiter and Saturn. The uncertainties are shown by the extent of the filled-in areas and are due to the uncertainties over the respective atmospheric mass in which the CO can mix (i.e. between a diluted 60 M⊕ diffuse core for both planets and a 20, 30 M⊕ solid core for Saturn and Jupiter, respectively).
In the text

	Fig. 12 [C/H] for the simulation HKB-4 (top) and HKB-5 (bottom) as a function of time.
In the text

	Fig. 13 Accretions rates for the HKB-fid simulations compared to the maximum capture rate onto planets due to atmospheric cooling.
In the text

Fig. 14 Estimation of the extra C/H due to late gas accretion onto an exoplanet as a function of its atmospheric mass for a belt of respectively 50 M⊕ and 5 M⊕ (purple and blue) based on the HKB-fid and HKB-4 simulations. The left axis shows the atmospheric C/H ratio relative to the protosolar abundance, and the right axis shows the absolute atmospheric C/H ratio. The red line represents the protosolar C/H ratio. The blue shaded area represents Neptune’s atmospheric mass.

In the text

	Fig. A.1 Neutral carbon accreted mass for each planet as a function of time for the simulations LKB-fid, 2, 3 and 4 (from left to right and top to bottom). The only differences between those simulations are the viscous parameter α and the accretion efficiency faccr (see Table 1).
In the text

	Fig. A.2 [C/H] for the LKB-fid simulation as a function of time for the ice giants. The uncertainties represented by the filled areas are due to the uncertainties on the atmospheric masses of Neptune and Uranus.
In the text

	Fig. A.3 C/H to protosolar ratio for the LKB-fid simulation as a function of time for the giant planets. The uncertainties represented by the filled areas are due to the uncertainties on the atmospheric masses of Jupiter and Saturn that may mix with carbon.
In the text