# Formation of secondary atmospheres on terrestrial planets by late disk accretion

Quentin Kral, Jeanne Davoult, Benjamin Charnay

► To cite this version:

Quentin Kral, Jeanne Davoult, Benjamin Charnay. Formation of secondary atmospheres on terrestrial planets by late disk accretion. Nature Astronomy, 2020, 4 (8), pp.769-775. 10.1038/s41550-020-1050-2. hal-03262625

HAL Id: hal-03262625

https://hal.science/hal-03262625v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Formation of secondary atmospheres on terrestrial planets by late disk accretion

Quentin Kral ID, JEANNE DAVOULT and BENJAMIN CHARNAY

Recently, gas disks have been discovered around main-sequence stars well beyond the usual protoplanetary disk lifetimes (that is, ≥10 Myr), when planets have already formed¹⁻⁴. These gas disks, mainly composed of CO, carbon and oxygen⁵⁻⁷, seem to be ubiquitous³ in systems with planetesimal belts (similar to our Kuiper belt), and can last for hundreds of millions of years⁸. Planets orbiting in these gas disks will accrete⁹,¹⁰ a large quantity of gas that will transform their primordial atmospheres into new secondary atmospheres with compositions similar to that of the parent gas disk. Here we quantify how large a secondary atmosphere can be created for a variety of observed gas disks and for a wide range of planet types. We find that gas accretion in this late phase is very important and an Earth's atmospheric mass of gas is readily accreted on terrestrial planets in very tenuous gas disks. In slightly more massive disks, we show that massive CO atmospheres can be accreted, forming planets with up to sub-Neptune-like pressures. Our results demonstrate that new secondary atmospheres with high metallicities and high C/O ratios will be created in these late gas disks, resetting their primordial compositions inherited from the protoplanetary disk phase, and providing a new birth to planets that lost their atmosphere to photoevaporation or giant impacts. We therefore propose a new paradigm for the formation of atmospheres on low-mass planets, which can be tested with future observations (James Webb Space Telescope (JWST), Extremely Large Telescope (ELT), Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL). We also show that this late accretion would show a very clear signature in sub-Neptunes or cold exo-Jupiters. Finally, we find that accretion creates cavities in late gas disks, which could be used as a new planet detection method, for low-mass planets a few to a few tens of astronomical units from their host stars.

The discovery of large amounts of gas around main-sequence stars is recent, with most detections occurring in the past few years³,¹¹. These late gas disks are observed in systems that have planetesimal belts, which are older than 10 Myr and can last for hundreds of millions of years⁸. It is thought that the observed gas is released from volatile-rich planetesimals when they collide with each other in the system's belts⁶,¹². The gas then viscously evolves¹³, spreading inward and outwards¹⁴,¹⁵. Hence the observed gas is likely secondary (rather than of primordial origin) and this late disk (main-sequence) phase is different from the younger (<10 Myr) protoplanetary disks that are much more massive and hydrogen rich, and in which giant planets form within a few millions of years¹⁶.

These late gas disks are nearly ubiquitous around A-type stars; gas has been detected around more than 70% of systems with bright planetesimal belts³. As for other stellar type stars or lower-mass

systems, the statistics are still based on too small a sample as these gas disks are harder to detect, but gas evolution models⁶ predict that all stars surrounded by planetesimal belts should have gas at a certain level that depends on the mass of solids in the system's planetesimal belt. More than 25% of stars have planetesimal belts massive enough to be detected through their infrared excess¹⁷, and it may be that most stars have belts below current detection limits (for instance, an exact equivalent of the Kuiper belt around a nearby star is not massive enough to be detectable with current facilities). These late gas disks might therefore not be the exception but rather the rule.

At >10 Myr in these late gas disks, planets have already formed. For instance, we now observe massive planets as early as 5 Myr (refs. ¹⁸⁻²⁰) and we know that terrestrial planets such as Mars formed very early (<10 Myr from cosmochemical evidence²¹) and the Earth took slightly longer (10–30 Myr), but most of its mass was acquired within 10 Myr (ref. ²²). The planets embedded in these disks will be able to accrete the disk gas in a similar way as when planets accrete gas in younger protoplanetary disks but for much longer timescales because late gas disks can last for tens of millions of years. Here we estimate how much gas can be accreted in this late phase onto the already-formed planets, and whether this gas can create new secondary atmospheres (with a composition similar to the source gas disk, that is, rich in carbon and oxygen and depleted in hydrogen) that would replace their original atmospheres.

To calculate the amount of gas accretion onto terrestrial or super-Earth planets in these late disks, we used an accretion model first developed for protoplanetary disk environments⁹,¹⁰ that we adapted to work in the late disk phase studied in this paper (Methods). A planet embedded in a gas disk will quickly fill its Hill sphere and whether it is able to accrete mass from the Hill sphere to its atmosphere depends on how fast the planet can cool (radiate away energy), and then contract to eventually accrete some more gas into its Hill sphere¹⁰ and grow an atmosphere. In the protoplanetary disk case, it is predicted that gas accretion should not depend strongly on disk density. However, in late gas disks (much more tenuous than protoplanetary disks), the mass available can be much lower than the mass that can be accreted from cooling/contraction and atmospheric growth is therefore limited by the gas disk mass (Methods).

We run the accretion model starting from an Earth-like planet of mass $1 M_⊕$ at 1 au from its host star and assuming that the planet's atmosphere has no, or very little, mass, because of photoevaporation (that is, it is in the desiccated part of the radius valley)²³ or desiccated after a large impact²⁴ (that should happen frequently in the late stages of planetary formation²⁵). Figure 1 shows the temporal evolution of the gas-to-core ratio (GCR) of the accreting planet for different gas crossing rates from the disk $\dot{M}$ at the planet location
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Formation of massive secondary atmospheres. Temporal evolution of the GCR of an Earth-like (1*M*⊕ at 1 au) planet starting with no atmosphere and orbiting in a late gas disk, for disks with gas crossing rates *Ṁ* at the planet varying from 10⁻⁸ to 10⁻² *M*⊕ Myr⁻¹. The dashed line shows the Earth's GCR at 8.6 × 10⁻⁷.</figcaption></figure>

(a more massive belt provides higher $\dot{M}$; Methods). If the disk is at steady state, $\dot{M}$ is the same as the gas release rate in the planetesimal belt¹⁴. For the most massive belts, we observe that the CO input rate in the belt can reach ~10⁻¹ *M*⊕Myr⁻¹ but for less massive belts releasing less CO per unit time, we expect values that can go below 10⁻⁷ *M*⊕Myr⁻¹ (ref. ⁸) and this is why we explore the effect of different $\dot{M}$ (10⁻⁸, 10⁻⁶, 10⁻⁴ and 10⁻² *M*⊕Myr⁻¹) on the total atmospheric mass accreted by a planet. We assume that most of the accretion happens within 100 Myr (Methods). Mass accretion for different accretion times can be straightforwardly extrapolated from Fig. 1.

We find that a planet starting with a low atmospheric mass will very rapidly accrete gas even at low gas input rates $\dot{M}$. Indeed, within 1 Myr, all simulations with $\dot{M} \ge 10^{-6} M_\oplus \text{Myr}^{-1}$ accreted more than an Earth's atmosphere mass. Given enough time, late accretion transforms a terrestrial planet (1$M_\oplus$ at 1 au) with no atmospheres (or even starting with an Earth-like or Venus atmosphere; Extended Data Fig. 8) into a planet with a massive atmosphere with a GCR up to >10⁻² for $\dot{M} \ge 10^{-4} M_\oplus \text{Myr}^{-1}$. We also show that this late accretion is very efficient on both larger (for example, 5$M_\oplus$ super-Earth) and smaller (for example, Mars-size) planets and for both closer (for example, at 0.1 au) and more distant planets (for example, at 10 au) in Fig. 2 and in Extended Data Fig. 6. We note that in our accretion scenario, there is no risk of runaway because the gas the planet can accrete from is limited and creating a Jupiter from a terrestrial mass core in these low-mass disks is simply not possible (Methods).

In Fig. 2, we show the effect of gas accretion on a planet's pressure and bulk density over time (Methods) and confirm that a variety of pressures from Earth-like (1 bar) to sub-Neptune-like (>10⁴ bar) can be reached on low-mass planets formed in these disks (Fig. 2, top). Even though pressures as high as 10⁵ bar (Neptune-like) can be reached on these planets, their densities never reach values as low as that of Neptune because CO atmospheres have a much higher mean molecular weight compared with hydrogen-dominated

atmospheres. We note that in the new accretion scenario proposed here, CO accretion does not depend much on the core mass (at least for planets with $R_H/H > 1$, where $R_H$ is the Hill radius and H is the scale height of the gas disk; Methods) so that the density is higher for higher-mass cores. We thus predict statistically that rocky planets (with $R_H/H > 1$ and excluding hydrogen-dominated atmospheres) with larger cores will have higher densities, which is in contrast with current models of planet formation and may be in line with current observations not seeing any clear correlation between the core mass and GCR but rather a large spread in GCRs²⁶,²⁷. We also note that as accretion is very efficient, we expect the outermost planets to accrete most disk material before it has time to spread further in, leading to decreasing densities for increasing semi-major axes (for a given core mass but note that for planets with $R_H/H < 1$ accretion will be slightly less efficient; Methods). There is a competing but less important effect, where closer-in planets have higher temperatures and thus lower densities for a given accreted mass (Fig. 2, bottom). For instance, our scenario may explain systems such as TRAPPIST-1 with decreasing densities (5.64 ± 0.4, 4.50 ± 0.20, 4.18 ± 0.19, 3.96 ± 0.6)²⁸ for the four outermost planets (e, f, g and h, respectively), where the low densities of TRAPPIST-1f, g and h could be explained by massive CO atmospheres of ~10⁵ bar.

Another strong prediction for this type of accretion is that planet atmospheres formed with this mechanism should have a high mean molecular weight and be mostly made of carbon and oxygen (rather than hydrogen) with a C/O ratio close to 1 (as CO is expected to be the main volatile released, but see Methods for detail as other molecules may also be released). Current Hubble Space Telescope (HST) observations of super-Earths revealed flat transit spectra, interpreted as the presence of atmospheres with high mean molecular weights and clouds²⁹⁻³¹. JWST and ARIEL will have the power (Methods) to take spectra for many more terrestrial planets (for example, the temperate TRAPPIST-1 planets³²) and super-Earths and find out whether it is common for these planets to have such high metallicities (and high C/O ratio) to test whether this new phase of late accretion really happens widely.

We find that our late accretion scenario is much more efficient and favourable to accrete volatiles on a terrestrial planet than delivery from impacts, that is, even if a late heavy bombardment-like event happens several hundreds of millions of years after the gas disk dissipated, the atmosphere would still be dominated by volatiles accreted by late gas accretion (Supplementary Information). The impact and late gas accretion scenarios could be distinguished based on the final composition of the observed planets, in particular looking at their C/O ratio (Supplementary Information). We also find that our scenario of late accretion works for a wide range of planets going from Mars-like to super-Earths and for close-in as well as distant planets from their stars. It also works if the planet is not initially devoid of atmosphere and has an Earth-like or a Venus atmosphere initially as it would replace the bulk of these atmospheres with gas coming from late disk accretion, even for cases where the disk is not very massive (Methods). Finally, if the gas accretion happens for only 10 Myr or for a period longer than 100 Myr, we find that the quantity of gas accreted by planets always remains considerable, leading to planetary atmospheres with masses at least that of the Earth's atmosphere for disks more massive than typical very-low-mass disks such as the Kuiper belt (Methods).

We also show that for existing (initially hydrogen rich) sub-Nep-tunes or more massive planets, the accretion is also going to affect their atmospheres once they mix with this new gas (Supplementary Information). In Fig. 3 (top), we show that the metallicity in sub-Neptunes could reach >1,000 times the solar metallicity, down to a factor 1.25 for more massive Jupiter-like planets. There are now some direct measurements of atmospheric metallicities in Neptunes and sub-Neptunes³³,³⁴. Some studies find near-solar metallicities (for example, GJ 3470 b)³³, while some others find super-solar
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Pressure and density evolution of an initially desiccated planet embedded in a late gas disk. Top: pressure versus GCR of an Earth-like (1*M<sub>⊕</sub>* at 1 au) planet starting with no atmosphere for different gas input rates, evolving up to 100 Myr. Bottom: bulk density versus *Ṁ* for different planet masses and semi-major axes, plotted at 100 Myr up to a maximum GCR of 0.1.</figcaption></figure>

metallicities (for example, HAT-P-11 b or HAT-P-26 b or K2-18 b)34–37 and more measurements will be welcome to test our scenario. The C/O ratio in these giant planets may also become close to 1 for a great variety of atmospheric masses and *Ṁ*, and even increase by 10% in a Jupiter-like planet for large values of *Ṁ* (Fig. 3, bottom). In Supplementary Fig. 2, we plot the metallicity and C/O ratio after 10 Myr of evolution to show that these effects can occur early in the planet's history. Measurements of C/O ratios are still scarce, but first results show super-solar (for example, HR 8799 c)38 as well as sub-solar values (for example, β Pictoris b)39 and more measurements, especially for smaller sub-Neptune planets, will help to test our scenario.

<figure><img src="image_5.png" /><figcaption>**Fig. 3** | Signature of late gas accretion on giant planets. Temporal variation of metallicity (top) and C/O ratio (bottom) as accretion proceeds for 100 Myr from an initially hydrogen-rich primordial atmosphere for different gas input rates and different initial atmospheric masses (sub-Neptune up to Jupiter).</figcaption></figure>

We also make predictions for detecting ongoing accretion onto young Jupiter-like planets or more massive brown dwarfs in direct imaging (Supplementary Information). When the gas accretes onto the outer envelope of the giant planet, it accumulates and diffuses inwards over time⁴⁰. We show that this accumulation (Supplementary Fig. 3) will be detectable on spectra observed with JWST Near Infrared Camera or Near Infrared Spectrograph as well as instruments on ELT in the M-band around 4.5–5.0 µm (Supplementary Fig. 4) for the coldest giant planets or brown dwarfs (<800 K), which would be a clear signature of this accretion. Only a few spectra of planets cooler than 800 K (for example, GJ 504 b¹¹) have been obtained so far in direct imaging but none targeted the required CO bands⁴²,⁴³.

With our accretion model, we find that accretion onto planets from late gas disks is very efficient, and for most configurations, a large fraction of the incoming gas is accreted rather than passed on further in (Methods). This means that these gas disks will often be very depleted inwards of a planet and one could infer the presence of a planet from the gas distribution. As the gas extends in the inner region of planetary systems, this new planet detection method could allow us to indirectly detect low-mass planets at a few astronomical units or farther from their host star. Using the high spatial resolution of the Atacama Large Millimeter/submillimeter Array (ALMA),
<figure><img src="image_6.png" /><figcaption>**Fig. 4 |** Cavity created by a low-mass planet in a late gas disk. ALMA synthesized [C] image (at 492.16 GHz, band 8) of a late disk whose cavity is clearly resolved and carved by a planet at 0.2 arcsec, that is, 10 au at 50 pc (with 5 h on source and a beam of 0.12 arcsec).</figcaption></figure>

it would be possible to pinpoint the planet location. In Fig. 4, we show a synthetic ALMA image of the carbon emission of a late gas disk with a terrestrial planet at 10 au from its host star located 50 pc away from Earth (Methods). The cavity is well resolved and the ALMA sensitivity is high enough to detect the signal at the 0.12" resolution. The carbon emission does not seem to extend all the way to the star for the few carbon gas disks¹⁵,⁴⁴ observed with ALMA so far but the resolution needs to be pushed further with ALMA to make any strong conclusions.

## Methods

**Input rate of gas and lifetime of gas disks.** The input rate of gas from the disk at the planet location is given by *M* in our study. This is the quantity of gas per unit time integrated over the whole disk scale height that is transferred radially inwards. These late gas disks are expected to viscously evolve¹⁴,⁴⁵ (maybe owing to the magnetorotational instability¹³) and there will be a transfer of most of the gas mass inwards (and angular momentum outwards). When steady state is reached, *M* becomes equal to the gas mass released in the planetesimal belt, which flows constantly inwards over time and accretes onto a planet or the central star. In this case, we can relate *M* to the surface density of the gas¹⁶ $Σ$ as $M = 3πvΣ$, where $v = αcH$ is the disk viscosity that can be parameterized with an $α$ value¹⁷, and $c_*$ and $H = c_*/Ω$ are the sound speed and disk scale height, respectively, with $Ω$ the Keplerian frequency. The value of $α$ has been estimated in a few studies by comparing the carbon quantity to what is expected from the measured gas input rate at the planetesimal belt location and it can vary from $10^{-4}$ to 0.1 (refs. ¹⁴,¹⁵,⁴⁸). Using population synthesis of these gas disks, it seems that all observations so far are consistent with $α$ being close to 0.1 (ref. ⁴⁹). An analytic study of the magnetorotational instability in these debris disks¹¹ shows that it could be very active owing to the high ionization fraction in these disks (compared with protoplanetary disks) and lead to $α$ values close to 0.1, or indeed smaller if non-ideal effects such as ambipolar diffusion are at play.

How fast the disk spreads viscously is set by $α$ such that the viscous timescale $t_v$ equals $r^2Ω/(αc_*^2)$, with $r$ the location of the planetesimal belt from which gas is released¹⁶. In Extended Data Fig. 1, we plot $t_v$ as a function of $α$ for different radial locations $r$ of the belt and gas temperatures $T$. The viscous timescales are typically

≥1 Myr (even when using the highest $α$ values), meaning that the spreading is often slow but steady state should be reached within the first tens of millions of years of evolution.

The release rate of gas in late gas disks is expected to decrease with time because planetesimals get destroyed over time and less gas can be released¹⁹⁻²⁵. The *M* we use in this study is an average of the real decreasing input rate over 100 Myr. We suppose that most of the accretion happens when the system is still young (as these belts are mostly observed in systems younger than 100 Myr (ref. ³)) before its planetesimal belt loses too much mass²⁵, that is, we do not follow the evolution beyond 100 Myr in our model, but it is possible that accretion still happens after 100 Myr at a lower rate. We can quantify the gas release rate as a function of time assuming that gas is released together with dust along the collisional cascade in Solar System comet proportions⁶. It is well known that mass loss rates in collisional cascades decrease with time and a simplified model of the evolution of the mass $M_s$ and mass loss rate $\dot{M}_s$ of solids with time is given by²³,⁵⁴

$$M_s(t) = \frac{M_{\text{init}}}{1 + t/t_{\text{col}}} \qquad (1)$$

$$\dot{M}_s(t) = \frac{M_s(t)^2}{M_{\text{init}} t_{\text{col}}} \qquad (2)$$

$$t_{\text{col}} = \frac{1.4 \times 10^{-9} r^{13/3} \frac{dr}{r} D_c Q_D^{*6/5}}{e^{5/3} M_*^{4/3} M_{\text{init}}} \qquad (3)$$

with $t_{col}$ is the initial collisional timescale of largest planetesimals, where we took typical values²² of largest planetesimal size $D_*$ = 10 km, belt width $dr = 0.5r$, mean eccentricity $e = 0.1$, planetesimal strength $Q_D^* = 330$ J kg⁻¹ and stellar mass $M_*$ = $1 M_\odot$ to compute Extended Data Fig. 2. We let the disk mean radius $r$ and the initial mass of solids $M_{\text{init}}$ be free parameters. We choose to vary $r$ from 50 to 100 au, which is typical of extra-solar Kuiper belt distances²⁵ and the initial total mass of solids $M_{\text{init}}$ between a Kuiper belt-like mass of 0.1 and 100 $M_\odot$, which is the maximum mass of solids available in protoplanetary progenitors to debris disks, inferred from submillimetre surveys²⁶. We then plot the temporal evolution of the gas release rate $\dot{M}$ in Extended Data Fig. 2 assuming that the gas release rate is 10% of the dust release rate²⁶. We can see that *M* remains approximately constant over 100 Myr of evolution except for the most massive belts ($M_{\text{init}} = 100 M_\odot$) that are closer-in (50 au) and have a much faster collisional evolution. Our approximation of taking $\dot{M}$ constant over 100 Myr works in most cases. However, for the most massive belts, our value of *M* should be understood as being the mean gas input rate integrated over time and one can trace back the *M* value at a given time $t$ using Extended Data Fig. 2 or the set of three equations above.

**Model.** The accretion model we used was first developed to explain the formation of super-Earths in protoplanetary disks⁹. A good understanding of the way accretion works was obtained through analytical model fitting to numerical simulations¹⁰. It is shown that the accretion rate hardly depends on the gas density from which it accretes but rather on the planet's ability to cool (or radiate away its energy). The more a planet can cool, the more it contracts, emptying parts of its Hill sphere, which get refilled very rapidly and the same process can happen again on a Kelvin-Helmholtz time. We also note that this accretion model is valid even in low gas density environments where planetary atmospheres are optically thin to incident starlight²⁷. Therefore, a good estimate for the accretion timescale is given by the cooling timescale equal to

$$t_{\text{cool}} = |E|/L_{\text{cool}} \qquad (4)$$

where $E$ is the atmosphere energy and $L_{\text{cool}}$ its luminosity¹⁰. From that timescale, we calculate the gas-to-core ratio GCR = $M_{\text{gas}}/M_{\text{core}}$ as a function of time $t$

$$\text{GCR}(t) = \frac{-bt + \sqrt{b^2t^2 - 4abt}}{2a} \qquad (5)$$

where $a = 3f_E k_B^{1+\frac{1}{\gamma-1}}(4\pi)^{\frac{1}{3}(\frac{3}{\gamma-1})} \rho_B^{\frac{1}{3}(\frac{3}{\gamma-1})} \kappa_{\text{rcb}}$ and $b = -4\pi^2 G^{1+1/(\gamma-1)} \sigma T_{\text{rcb}}^{3-1/(\gamma-1)} (\mu_{\text{rcb}} m_{\text{H}})^{1+1/(\gamma-1)} \nabla_{\text{ad}}^{1+1/(\gamma-1)} M_{\text{core}}^{2(\frac{1}{\gamma}-1)} 3^{\frac{1}{3}(\frac{3}{\gamma-1})}$.

The term $f_E$ equals $G(4\pi\rho_B/3)^{1/3}$, and $\rho_B$ is the bulk density of the core, $M_{\text{core}}$ is the core mass, $\nabla_{\text{ad}} = (\gamma-1)/\gamma$ is the adiabatic gradient, $T_{\text{rcb}}$, $\mu_{\text{rcb}}$ and $\kappa_{\text{rcb}}$ are the temperature, mean molecular weight and atmospheric opacity at the radiative/convective boundary (rcb). Finally, $m_{\text{H}}$ is the proton mass, $k_B$ the Boltzmann constant, $G$ the gravitational constant and $\sigma$ the Stefan-Boltzmann constant. This new analytic expression given by equation (5) is otherwise general and do not make the assumption that GCR ≪ 1 as in previous work¹⁰. Unless otherwise stated, we fix all values as given in the reference paper¹⁰. To be realistic, we use real opacities derived for highly metal-enriched (100 times solar abundances) super-Earths⁵⁸ and in this case GCR($t$) must be obtained numerically (because now opacities vary with GCR). We fix $\gamma = 1.4$ and the mean molecular weight as being 28, which corresponds to an atmosphere accreting a majority of CO as would be the case
in our scenario¹⁵. We note that in non-shielded disks (that is, disks not massive enough for carbon production to be sufficiently high that neutral carbon shields CO from photodissociating¹⁶), μ is closer to 14, that is, dominated by carbon and oxygen¹⁴⁻¹⁶ rather than CO but our final results are not so much affected by the value of μ (except for the very high *M* values) as demonstrated further in Extended Data Fig. 3. There could also be some other molecules released in small quantities in these belts, such as CN, but only CO is detected so far.

We find that GCR(*t*) can be simplified in the limit where $-a/(bt) \gg 1$ or $-a/(bt) \ll 1$. For $-a/(bt) \gg 1$, we find GCR = $\sqrt{-bt/a}$ and when $-a/(bt) \ll 1$, GCR = $-bt/a$. The latter regime has not been studied so far as it only appears at large *t* and large μ, which is typical of late gas disks but not of protoplanetary disks.

The maximum rate at which a planet can accrete an atmosphere is given by *M*<sub>gas</sub>, which we obtain by deriving *M*<sub>gas</sub> = GCR(*t*)*M*<sub>core</sub>. We find

$$ \dot{M}_{\text{gas}} = \frac{M_{\text{core}}}{2a} \times \left[ -b + \frac{1}{2} \times (b^2 t^2 - 4abt)^{-1/2} \times (2b^2 t - 4ab) \right] \quad (6) $$

where *a* and *b* are defined above. In equation (6), we have assumed a constant opacity, but in reality the opacity varies with time as the GCR increases. In our code, we compute this derivative (M̃<sub>gas</sub>) numerically to take this complexity into account. We also derive M̃<sub>gas</sub> analytically for the case where κ depends on *t*. The opacity can be parameterized³⁸ as κ ∝ κ₀ρ<sub>rcb</sub>τ<sup>β</sup><sub>rcb</sub>Z<sup>δ</sup>, where the time dependence comes in the density at the rcb ρ<sub>rcb</sub> which is proportional to GCR¹⁰. We find that there are also two regimes depending on the value of $-a/(bt)$. If $-a/(bt) \gg 1$ then M̃<sub>gas</sub> = *M*<sub>core</sub>/(2 + 3α)√$-b/(at)$ and when $-a/(bt) \ll 1$ then M̃<sub>gas</sub> = $-M*\text{core}/(1 + α)(b/a)$. For the case where $-a/(bt) \gg 1$, we thus find

$$ M_{\text{gas}} \propto \left[ T_{\text{rcb}}^{3-\beta/(1+\alpha)/(\gamma-1)} (\mu_{\text{rcb}} \nabla_{\text{ad}})^{1+(1+\alpha)/(\gamma-1)} \rho_{\text{b}}^{-4/3-\alpha+1/3(1+\alpha)/(\gamma-1)} \times M_{\text{core}}^{2/3(1+\alpha)/(\gamma-1)+4/3+2\alpha} Z^{-\delta} \kappa_0^{-1} t^{1/(2+\alpha)} \right] \quad (7) $$

$$ \dot{M}_{\text{gas}} \propto \left[ T_{\text{rcb}}^{3-\beta/(1+\alpha)/(\gamma-1)} (\mu_{\text{rcb}} \nabla_{\text{ad}})^{1+(1+\alpha)/(\gamma-1)} \rho_{\text{b}}^{-4/3-\alpha+1/3(1+\alpha)/(\gamma-1)} \times M_{\text{core}}^{2/3(1+\alpha)/(\gamma-1)+4/3+2\alpha} Z^{-\delta} \kappa_0^{-1} t^{-1-2\alpha} \right]^{1/(2+\alpha)} \quad (8) $$

and for the case $-a/(bt) \ll 1$, we find

$$ M_{\text{gas}} \propto \left[ T_{\text{rcb}}^{3-\beta/(1+\alpha)/(\gamma-1)} (\mu_{\text{rcb}} \nabla_{\text{ad}})^{1+(1+\alpha)/(\gamma-1)} \rho_{\text{b}}^{-4/3-\alpha+1/3(1+\alpha)/(\gamma-1)} \times M_{\text{core}}^{2/3(1+\alpha)/(\gamma-1)+1/3+\alpha} Z^{-\delta} \kappa_0^{-1} t^{1/(1+\alpha)} \right] \quad (9) $$

$$ \dot{M}_{\text{gas}} \propto \left[ T_{\text{rcb}}^{3-\beta/(1+\alpha)/(\gamma-1)} (\mu_{\text{rcb}} \nabla_{\text{ad}})^{1+(1+\alpha)/(\gamma-1)} \rho_{\text{b}}^{-4/3-\alpha+1/3(1+\alpha)/(\gamma-1)} \times M_{\text{core}}^{2/3(1+\alpha)/(\gamma-1)+1/3+\alpha} Z^{-\delta} \kappa_0^{-1} t^{-\alpha} \right]^{1/(1+\alpha)} \quad (10) $$

Taking an opacity with α=0.6, β=2.2 and δ=1 valid for dust free atmospheres beyond 1 au (ref. ¹⁰) and that γ=1.4 (for CO atmospheres), we find for the case $-a/(bt) \gg 1$ that

$$ M_{\text{gas}} \propto T_{\text{rcb}}^{-1} \mu_{\text{rcb}}^{1.6} \rho_{\text{b}}^{-0.19} M_{\text{core}}^{1.6} Z^{-0.3} \kappa_{0}^{-0.6} t^{0.3} \quad (9) $$

$$ \dot{M}_{\text{gas}} \propto T_{\text{rcb}}^{-1} \mu_{\text{rcb}}^{1.6} \rho_{\text{b}}^{-0.19} M_{\text{core}}^{1.6} Z^{-0.3} \kappa_{0}^{-0.6} t^{-0.7} \quad (10) $$

and for the case $-a/(bt) \ll 1$, we find

$$ M_{\text{gas}} \propto T_{\text{rcb}}^{-2} \mu_{\text{rcb}}^{3.1} \rho_{\text{b}}^{-0.38} M_{\text{core}}^{2.3} Z^{-0.6} \kappa_{0}^{-0.6} t^{0.6} \quad (11) $$

$$ \dot{M}_{\text{gas}} \propto T_{\text{rcb}}^{-2} \mu_{\text{rcb}}^{3.1} \rho_{\text{b}}^{-0.38} M_{\text{core}}^{2.3} Z^{-0.6} \kappa_{0}^{-0.6} t^{-0.4} \quad (12) $$

In late gas disks, the total gas mass is much smaller than in protoplanetary disks and it may happen that the gas available per unit time is smaller than M̃<sub>gas</sub>. Therefore, in our code, at each time step, we compare the gas crossing rate M̃ to M̃<sub>gas</sub>. For most cases studied in this paper (except the Mars-mass planet case; Extended Data Fig. 6), we find that M̃ is indeed lower than M̃<sub>gas</sub> and the accretion is limited by the quantity available rather than by the planet cooling. For the cases where M̃ < M̃<sub>gas</sub>, the mass that accumulates is lower than the theoretical mass M<sub>gas</sub> given by equation (5) so that when computing M̃<sub>gas</sub> one should take the theoretical accretion rate for the mass that actually accumulated rather than the theoretical mass, which we do in our code. For the regime where $-a/(bt) \gg 1$ this means that at a given time *t*, one should take the accretion rate at time *t* = *t*(*M*<sub>atm</sub>/M<sub>gas</sub>)<sup>2+α</sup>, which gives an accretion rate which is higher by a factor M<sub>gas</sub>/(M̃<sub>atm</sub>). We show the theoretical M̃<sub>gas</sub> in Extended Data Fig. 4, where we plot the temporal evolution of M̃<sub>gas</sub> (taking into account that the opacity varies with time) for different values of planet semi-major axes *a*<sub>pl</sub>, atmospheric mean molecular weight μ and core mass *M*<sub>core</sub>, corresponding to the typical range of values that we are interested in in our study.

We also compare the Hill radius R<sub>H</sub> (equal to *a*<sub>pl</sub>[M<sub>pl</sub>/(3*M*<sub>⋆</sub>)]<sup>1/3</sup>) of the planet (of semi-major axis *a*<sub>pl</sub> and mass *M*<sub>pl</sub>) with the scale height *H* of the disk. For very small Earth-like planets (Extended Data Fig. 5), it can happen that R<sub>H</sub> < *H* and in this case, some gas crossing at a rate M̃ cannot be accreted by the planet. We then recalculate a new M̃ that can be accreted by only considering the M̃ that crosses the planet's Hill sphere rather than the whole scale height. To calculate the quantity of gas that cannot be accreted, we assume a sphere of radius *H* on top of the sphere of

radius R<sub>H</sub>, and take out the parts of the large *H* sphere with a height greater than R<sub>H</sub>. This gives a new M̃<sub>new</sub> value that is 3/2(R<sub>H</sub>/H) - 1/2(R<sub>H</sub>/H)<sup>3</sup> of the full M̃. Note that for cases that reach GCR values approaching 1, R<sub>H</sub> can become larger and R<sub>H</sub>/H becomes greater, hence enhancing accretion. Therefore, we find that M̃<sub>new</sub> scales roughly with R<sub>H</sub>/H rather than its square as would be expected from Bondi accretion or dividing the volume of two tori of radii R<sub>H</sub> and *H*, respectively. This is because the orbital timescale of the planet is much faster than the viscous drift timescale of the gas flow and the new M̃ will be approximately equal to the ratio of cylindrical collisional cross-section of the planet 2πa<sub>pl</sub>R<sub>H</sub> and the cylindrical flow cross-section 2πa<sub>pl</sub>*H*, which equals R<sub>H</sub>/H.

Furthermore, we compare our typical accretion rate M̃ with a Bondi-like accretion that may be relevant in the regime where R<sub>H</sub> < *H*. Bondi-like accretion happens at a rate M̃ = 4πR<sub>H</sub><sup>2</sup>ΣΩ, where Σ is the local secondary disk gas surface density (with Σ = M̃/(3πν)) at steady state) and Ω is the Keplerian frequency. We find that M̃<sub>B</sub>/M̃ = 4/(3α)(R<sub>H</sub>/H)<sup>2</sup> or M̃<sub>B</sub>/M̃<sub>new</sub> ≈ (1/α)(R<sub>H</sub>/H), with α typically between 10<sup>-4</sup> and 0.1 as explained earlier¹³,¹⁴. It means that in most cases as R<sub>H</sub>/H > 0.1, the Bondi accretion rate is higher than the rate at which gas is delivered, and one is still limited by M̃. For instance, for a 1 M<sub>☉</sub> planet at 1 au (with R<sub>H</sub>/H ≈ 0.9, as shown in Extended Data Fig. 5), we obtain a new input rate that is ~0.98 times of the full M̃. This value is therefore usually very close to 1 and does not make important differences even when taking extreme cases such as a Mars-mass (0.1 M<sub>☉</sub>) planet at Mars distance (1.5 au) or a 0.5 M<sub>☉</sub> planet at 10 au, for which we find R<sub>H</sub>/H ≈ 0.37 for both cases, and the new input rate M̃<sub>new</sub> is 0.53 of the full M̃, therefore lowering the final masses accreted by a factor of 2 on such planets. We show the results for these two planets in Extended Data Fig. 6. We end up with GCR values very similar to the case in Fig. 1 (1 M<sub>☉</sub>) for the case 0.5 M<sub>☉</sub> because although the final gas mass is divided by two, the core mass is also twice smaller, hence the GCR is similar. Owing to the much lower core mass of the Mars-mass case, the final GCR is five times higher than in the 1 M<sub>☉</sub> case except for the highest input rates where the theoretical cooling accretion rate becomes smaller than M̃ and the accreted mass becomes smaller, explaining the shallower slope of GCR at large *t*. For cases with R<sub>H</sub>/H > 1, for example, planets with more massive cores, the results are similar to the case shown in Fig. 1 after a correction of the GCR to take account of the *M*<sub>core</sub> value (that is, the GCR goes down by a factor ∝ *M*<sub>core</sub>) but we note that in this case the hydrodynamics of the flow can be complicated owing to shock waves developing near the planet/disk interface³⁹.

To verify that the Hill radius is the radius of interest, which sets the length scale of accretion as assumed here, we also compute the Bondi radius (equal to 2GM<sub>pl</sub>/c<sup>2</sup>) for a large variety of planet masses and planet semi-major axes. Extended Data Fig. 7 shows that indeed the Bondi radius is always greater than the Hill radius in our study, confirming that the Hill radius should be used in previous computations.

To compute the planet's bulk density given its GCR and core mass, we improve on previous work³⁰ and do not assume that the gravity *g* is constant with pressure so that the thickness of the adiabatic convective region of the atmosphere is given by (when R<sub>core</sub> > c<sub>p</sub>/g<sub>batm</sub>(T<sub>batm</sub> - T<sub>atm</sub>))

$$ d_{\text{atm}} = \frac{c_p}{g_{\text{batm}}} (T_{\text{batm}} - T_{\text{atm}}) \frac{1}{1 - c_p / (R_{\text{core}} g_{\text{batm}}) (T_{\text{batm}} - T_{\text{atm}})} \quad (11) $$

where *P*<sub>atm</sub>=0.1 bar typical of the pressure at which the atmosphere transitions from adiabatic to isothermal<sup>51</sup>, T<sub>batm</sub> = *T*<sub>atm</sub>(P<sub>batm</sub>/P<sub>atm</sub>)<sup>*κ*</sup>, where κ=2/(2+n), with *n* the number of degrees of freedom equal to 5 for diatomic gases, with *T*<sub>atm</sub> = *T*<sub>*</sub>/2√R<sub>*</sub>/(2α) (*T*<sub>*</sub> and *R*<sub>*</sub> are the temperature and radius of the host star), and we calculate *P*<sub>batm</sub> by integrating *dm* = 4π(R<sub>core</sub> + *z*)<sup>2</sup>D*P*/*g*(*P*) from *P*<sub>atm</sub> up to a value *P* for which the atmosphere mass is *M*<sub>atm</sub> = GCRM<sub>core</sub>. In the integral, we use that *g*(*z*) = g<sub>batm</sub>/(1 + *z*/R<sub>core</sub>)<sup>2</sup> and *z* = *T*<sub>atm</sub>c<sub>p</sub>/g<sub>batm</sub>((P<sub>batm</sub>/P<sub>atm</sub>)<sup>*κ*</sup> - (P/P<sub>atm</sub>)<sup>*κ*</sup>)/((1 - c<sub>p</sub>*T*<sub>atm</sub>/(R<sub>core</sub>g<sub>batm</sub>))(P<sub>batm</sub>/P<sub>atm</sub>)<sup>*κ*</sup> - (P/P<sub>atm</sub>)<sup>*κ*</sup>), where c<sub>p</sub> is the heat capacity, *T*<sub>atm</sub>, the atmosphere surface temperature, g<sub>batm</sub> and *T*<sub>batm</sub> the gravity and the temperature at the bottom of the atmosphere, respectively. Since *z* depends on *P*<sub>batm</sub>, we calculate *P*<sub>batm</sub> by iteration, starting from an initial value equal to *M*<sub>atm</sub>g<sub>batm</sub>/(4πR<sub>core</sub><sup>2</sup>). We then add the isothermal part of the atmospheric thickness from 100 to 20 mbar: *R*<sup>*</sup>*T*<sub>atm</sub>/g(100 mbar)μ)ln(100/20) (where *R*<sup>*</sup> is the universal gas constant) assuming that the transit radius is typically observed at around 20 mbar (ref. ⁶⁰).

**Accretion onto a pre-existing Earth-like atmosphere.** In Extended Data Fig. 8, we show what happens when an atmosphere starts growing from a pre-existing Earth-mass atmosphere rather than an empty atmosphere as in Figs. 1 and 2. The results are the same, with all the lines being just shifted up, starting at the pre-existing level.

This example shows that the final atmosphere will be dominated by volatiles delivered from late gas accretion for all cases with M̃ > 10<sup>-8</sup> M<sub>☉</sub> Myr<sup>-1</sup>, that is, in all cases with belts more massive than the Kuiper belt (Extended Data Fig. 2), which is the least massive belt we know of so far. One could start with an even more massive atmosphere and just shift the lines up by the pre-existing mass or pressure to predict the final atmospheric masses and pressures. For instance, starting with a Venus-like atmosphere would still up with an atmosphere being dominated by late accretion after 100 Myr for M̃ > 10<sup>-6</sup> M<sub>☉</sub> Myr<sup>-1</sup>.
**Formation of cavities owing to planets.** Our accretion models for both Earth-like and giant planets (Supplementary Information) suggest that accretion onto planets is very efficient and in most cases, the rate of gas that can be accreted is higher than what is available (that is *M*). For this reason, when the gas disk spreads inwards and crosses a planet's orbit, it may not be able to spread further in as all of the inflowing gas is accreted onto the planet that is being crossed. However, this depends on $R_{11}/H$. As shown in Extended Data Fig. 5, for small or distant planets, $R_{11}$ can be smaller than the disk scale height *H* and a certain fraction of the gas will flow inwards. We note that for atmospheres with low *μ*, high *γ* or accreting for a very long time, the theoretical cooling accretion rate may become smaller than *M̄* and gas would flow inwards. We also note that the gas model used for low-mass planets is one dimensional and thus assumes an axisymmetric gas flow but the gas flow geometry around the planet may be more complicated in three dimensions⁶², which may lead to gas flowing inwards anyway⁶³. However, it is still unclear how much these complications affect the overall one-dimensional accretion rates⁶⁴, but recent three-dimensional simulations find gas accretion rates that are comparable to one-dimensional derivations⁶⁵. In the end, we expect accretion to be efficient and a gap in density should be seen in the gas distribution after crossing a planet, which would pinpoint the accreting planet location.

**Gas distribution with ALMA to infer the planet position.** In Fig. 4, we used carbon observations as a good tracer for these cavities instead of CO. This is for two main reasons. First, carbon emission in band 8 seems to be a better tracer of this gas than CO in either band 6 or 7 according to models⁶⁶ and to the first few observations of neutral carbon in these disks²,¹⁵,⁴⁴. Second, carbon is always expected to spread in the inner region, given enough time, while CO photodissociates in about 100 yr in unshielded disks and remains colocated with the parent belt of planetesimals, implying that the CO cavity observed is then due to photodissociation rather than accreting planets. Only in the case of massive gas disks can CO be shielded by carbon¹⁵ and may have time to viscously spread further inwards than the planetesimetal belt. For the latter case, for systems where CO had time to spread in the inner region where planets are located, then CO cavities could also be used to infer planets but carbon is much more general as it is not subject to photodissociation.

We now explain the details of how we produced the synthetic ALMA image shown in Fig. 4. We first created a density profile for a late gas disk at steady state that scales as $r^{-1}$ inwards¹⁵. The centre of the belt is at 50 au where the gas temperature is 20 K (scaling as $r^{-0.5}$). An input rate of $10^{-3} M_{\oplus} Myr^{-1}$ of CO has been chosen. We put the planet at 10 au and therefore impose that the gas density drops to zero within 10 au.

We then used the radiative transfer code RADMC-3D⁶⁶ to obtain the emission of the carbon fine structure line at 492.16 GHz assuming that the inclination is 30° and the position angle is 45°. We then used the CASA software to simulate an ALMA observation at a resolution of 0.12" in band 8 for 5 hours on source (we used the C43-5 configuration together with the C43-2 to get more extended emission). Finally, we created a moment-0 image of the final cube and obtained Fig. 4. We see that the gas disk is well detected and the cavity within 10 au is well resolved.

**Data availability**

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

**Code availability**

The particular scripts used for the analysis are made in python and available on reasonable request from the corresponding author.

Received: 8 August 2019; Accepted: 18 February 2020;

Published online: 6 April 2020

**References**

1. Dent, W. R. F. et al. Molecular gas clumps from the destruction of icy bodies in the β Pictoris debris disk. *Science* **343**, 1490–1492 (2014).

2. Lieman-Sifry, J. et al. Debris disks in the scorpius-centaurus OB association resolved by ALMA. *Astrophys. J.* **828**, 25 (2016).

3. Moór, A. et al. Molecular gas in debris disks around young A-type stars. *Astrophys. J.* **849**, 123 (2017).

4. Matrà, L., Óberg, K. I., Wilner, D. J., Olofsson, J. & Bayo, A. On the ubiquity and stellar luminosity dependence of exometary CO gas: detection around M dwarf TWA 7. *Astron. J.* **157**, 117 (2019).

5. Cataldi, G. et al. Herschel/HIFI observations of ionised carbon in the β Pictoris debris disk. *Astron. Astrophys.* **563**, A66 (2014).

6. Kral, Q., Matrà, L., Wyatt, M. C. & Kennedy, G. M. Predictions for the secondary CO, C and O gas content of debris discs from the destruction of volatile-rich planetesimals. *Mon. Not. R. Astron. Soc.* **469**, 521–550 (2017).

7. Higuchi, A. E. et al. Detection of submillimeter-wave [C1] emission in gaseous debris disks of 49 Ceti and β Pictoris. *Astrophys. J.* **839**, L14 (2017).

8. Matrà, L. et al. Detection of exometary CO within the 440 Myr old Fomalhaut belt: a similar CO+CO₂ ice abundance in exocomets and Solar System comets. *Astrophys. J.* **842**, 9 (2017).

9. Lee, E. J., Chiang, E. & Ormel, C. W. Make super-Earths, not Jupiters: accreting nebular gas onto solid cores at 0.1 au and beyond. *Astrophys. J.* **797**, 95 (2014).

10. Lee, E. J. & Chiang, E. To cool is to accrete: analytic scalings for nebular accretion of planetary atmospheres. *Astrophys. J.* **811**, 41 (2015).

11. Kóspál, Á. et al. ALMA observations of the molecular gas in the debris disk of the 30 Myr old star HD 21997. *Astrophys. J.* **776**, 77 (2013).

12. Zuckerman, B. & Song, I. A 40 Myr old gaseous circumstellar disk at 49 Ceti: massive CO-rich comet clouds at young A-type stars. *Astrophys. J.* **758**, 77 (2012).

13. Kral, Q. & Latter, H. The magnetorotational instability in debris-disc gas. *Mon. Not. R. Astron. Soc.* **461**, 1614–1620 (2016).

14. Kral, Q. et al. A self-consistent model for the evolution of the gas produced in the debris disc of β Pictoris. *Mon. Not. R. Astron. Soc.* **461**, 845–858 (2016).

15. Kral, Q. et al. Around HD 131835: reinterpreting young debris discs with protoplanetary disc levels of CO gas as shielded secondary discs. *Mon. Not. R. Astron. Soc.* **489**, 3670–3691 (2019).

16. Testi, L. et al. Hunting for planets in the HL tau disk. *Astrophys. J.* **812**, L38 (2015).

17. Thureau, N. D. et al. An unbiased study of debris discs around A-type stars with Herschel. *Mon. Not. R. Astron. Soc.* **445**, 2558–2573 (2014).

18. Keppler, M. et al. Discovery of a planetary-mass companion within the gap of the transition disk around PDS 70. *Astron. Astrophys.* **617**, A44 (2018).

19. Müller, A. et al. Orbital and atmospheric characterization of the planet within the gap of the PDS 70 transition disk. *Astron. Astrophys.* **617**, L2 (2018).

20. Haffert, S. Y. et al. Two accreting protoplanets around the young star PDS 70. *Nat. Astron.* **3**, 749–754 (2019).

21. Nimmo, F. & Kleine, T. How rapidly did mars accrete? Uncertainties in the Hf–W timing of core formation. *Icarus* **191**, 497–504 (2007).

22. Jacobsen, S. B. The Hf–W isotopic system and the origin of the Earth and Moon. *Annu. Rev. Earth Planet. Sci.* **33**, 531–570 (2005).

23. Owen, J. E. Atmospheric escape and the evolution of close-in exoplanets. *Annu. Rev. Earth Planet. Sci.* **47**, 67–90 (2019).

24. Yalinewich, A. & Schlichting, H. Atmospheric mass-loss from high-velocity giant impacts. *Mon. Not. R. Astron. Soc.* **486**, 2780–2789 (2019).

25. Quintana, E. V., Barclay, T., Borucki, W. J., Rowe, J. F. & Chambers, J. E. The frequency of giant impacts on Earth-like worlds. *Astrophys. J.* **821**, 126 (2016).

26. Lopez, E. D. & Fortney, J. J. Understanding the mass–radius relation for sub-Neptunes: radius as a proxy for composition. *Astrophys. J.* **792**, 1 (2014).

27. Lee, E. J. & Chiang, E. Breeding super-Earths and birthing super-puffs in transitional disks. *Astrophys. J.* **817**, 90 (2016).

28. Grimm, S. L. et al. The nature of the TRAPPIST-1 exoplanets. *Astron. Astrophys.* **613**, A68 (2018).

29. Kreidberg, L. et al. Clouds in the atmosphere of the super-Earth exoplanet GJ1214b. *Nature* **505**, 69–72 (2014).

30. Charnay, B., Meadows, V., Misra, A., Leconte, J. & Arney, G. 3D modeling of GJ1214b's atmosphere: formation of inhomogeneous high clouds and observational implications. *Astrophys. J.* **813**, L1 (2015).

31. Morley, C. V. et al. Thermal emission and reflected light spectra of super Earths with flat transmission spectra. *Astrophys. J.* **815**, 110 (2015).

32. Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T. & Fortney, J. J. Observing the atmospheres of known temperate Earth-sized planets with JWST. *Astrophys. J.* **850**, 121 (2017).

33. Benneke, B. et al. A sub-Neptune exoplanet with a low-metallicity methane-depleted atmosphere and Mie-scattering clouds. *Nat. Astron.* **3**, 813–821 (2019).

34. Fraine, J. et al. Water vapour absorption in the clear atmosphere of a neptune-sized exoplanet. *Nature* **513**, 526–529 (2014).

35. Wakeford, H. R. et al. HAT-P-26b: a neptune-mass exoplanet with a well-constrained heavy element abundance. *Science* **356**, 628–631 (2017).

36. Benneke, B. et al. Water vapor and clouds on the habitable-zone sub-Neptune exoplanet K2-18b. *Astrophys. J.* **887**, L14 (2019).

37. Tsiaras, A., Waldmann, I. P., Tinetti, G., Tennyson, J. & Yurchenko, S. N. Water vapour in the atmosphere of the habitable-zone eight-Earth-mass planet K2-18 b. *Nat. Astron.* **3**, 1086–1091 (2019).

38. Konopacky, Q. M., Barman, T. S., Macintosh, B. A. & Marois, C. Detection of carbon monoxide and water absorption lines in an exoplanet atmosphere. *Science* **339**, 1398–1401 (2013).

39. GRAVITY Collaboration et al. Peering into the formation history of β Pictoris b with VLTI/GRAVITY long baseline interferometry. *Astron. Astrophys.* **633**, A110 (2020).

40. Charnay, B. et al. A self-consistent cloud model for brown dwarfs and young giant exoplanets: comparison with photometric and spectroscopic observations. *Astrophys. J.* **854**, 172 (2018).

41. Janson, M. et al. Direct imaging detection of methane in the atmosphere of GJ 504 b. *Astrophys. J.* **778**, L4 (2013).
42. Macintosh, B. et al. Discovery and spectroscopy of the young jovian planet 51 Eri b with the gEmini Planet Imager. *Science* **350**, 64–67 (2015).

43. Skemer, A. J. et al. The LEECH exoplanet imaging survey: characterization of the coldest directly imaged exoplanet, GJ 504 b, and evidence for superstellar metallicity. *Astrophys. J.* **817**, 166 (2016).

44. Cataldi, G. et al. ALMA resolves C1 emission from the β Pictoris debris disk. *Astrophys. J.* **861**, 72 (2018).

45. Xie, J.-W., Brandeker, A. & Wu, Y. On the unusual gas composition in the β Pictoris debris disk. *Astrophys. J.* **762**, 114 (2013).

46. Lynden-Bell, D. & Pringle, J. E. The evolution of viscous discs and the origin of the nebular variables. *Mon. Not. R. Astron. Soc.* **168**, 603–637 (1974).

47. Shakura, N. I. & Sunyaev, R. A. Reprint of 1973A&A....24. 337S. Black holes in binary systems. Observational appearance. *Astron. Astrophys.* **500**, 33 (1973).

48. Moór, A. et al. New millimeter CO observations of the gas-rich debris disks 49 Cet and HD 32297. *Astrophys. J.* **884**, 108 (2019).

49. Marino, S. et al. Population synthesis of exometary gas around A stars. *Mon. Not. R. Astron. Soc.* (in the press).

50. Thébault, P. & Augereau, J.-C. Collisional processes and size distribution in spatially extended debris discs. *Astron. Astrophys.* **472**, 169 (2007).

51. Kral, Q., Thébault, P. & Charnoz, S. LIDT-DD: a new self-consistent debris disc model that includes radiation pressure and couples dynamical and collisional evolution. *Astron. Astrophys.* **558**, A121 (2013).

52. Wyatt, M. C. Evolution of debris disks. *Annu. Rev. Astron. Astrophys.* **46**, 339–383 (2008).

53. Dominik, C. & Decin, G. Age dependence of the Vega phenomenon: theory. *Astrophys. J.* **598**, 626 (2003).

54. Wyatt, M. C. et al. Steady state evolution of debris disks around a stars. *Astrophys. J.* **663**, 365 (2007).

55. Matrà, L. et al. An empirical planetesimal belt radius-stellar luminosity relation. *Astrophys. J.* **859**, 72 (2018).

56. Williams, J. P. & Cieza, L. A. Protoplanetary disks and their evolution. *Annu. Rev. Astron. Astrophys.* **49**, 67–117 (2011).

57. Lee, E. J., Chiang, E. & Ferguson, J. W. Optically thin core accretion: how planets get their gas in nearly gas-free discs. *Mon. Not. R. Astron. Soc.* **476**, 2199 (2018).

58. Freedman, R. S. et al. Gaseous mean opacities for giant planet and ultracool dwarf atmospheres over a range of metallicities and temperatures. *Astrophys. J. Suppl. Ser.* **214**, 25 (2014).

59. Tanigawa, T., Ohtsuki, K. & Machida, M. N. Distribution of accreting gas and angular momentum onto circumplanetary disks. *Astrophys. J.* **747**, 47 (2012).

60. Dorn, C., Mosegaard, K., Grimm, S. L. & Alibert, Y. Interior characterization in multiplanetary systems: TRAPPIST-1. *Astrophys. J.* **865**, 20 (2018).

61. Robinson, T. D. & Catling, D. C. Common 0.1 bar tropopause in thick atmospheres set by pressure-dependent infrared transparency. *Nat. Geosci.* **7**, 12–15 (2014).

62. Ormel, C. W., Shi, J.-M. & Kuiper, R. Hydrodynamics of embedded planets' first atmospheres—II. A rapid recycling of atmospheric gas. *Mon. Not. R. Astron. Soc.* **447**, 3512–3525 (2015).

63. Béthune, W. & Rafikov, R. R. Envelopes of embedded super-Earths—II. Three-dimensional isothermal simulations. *Mon. Not. R. Astron. Soc.* **488**, 2365–2379 (2019).

64. Fung, J., Artymowicz, P. & Wu, Y. The 3D flow field around an embedded planet. *Astrophys. J.* **811**, 101 (2015).

65. D'Angelo, G. & Bodenheimer, P. Three-dimensional radiation-hydrodynamics calculations of the envelopes of young planets embedded in protoplanetary disks. *Astrophys. J.* **778**, 77 (2013).

66. Dullemond, C. P. et al. RADMC-3D: a multi-purpose radiative transfer tool. *Astrophysics Source Code Library* ascl:1202.015 (2012).

## Acknowledgements

We thank G. Rosotti, P. Thebault and A. Shannon for discussions. Q.K. deidcates this paper to Michaël.

## Author contributions

Q.K. led the work, proposed the original idea, wrote the manuscript and produced Figs. 3 and 4. J.D. coded the model and produced Figs. 1 and 2. B.C. provided atmospheric parameters to input into the model, expertise in atmosphere observations and produced the synthetic spectra shown in the paper. All authors contributed to the interpretation of the results and commented on the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at https://doi.org/10.1038/s41550-020-1050-2.

*Supplementary information* is available for this paper at https://doi.org/10.1038/s41550-020-1050-2.

*Correspondence and requests for materials* should be addressed to Q.K.

*Reprints and permissions information* is available at www.nature.com/reprints.

*Publisher's note* Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2020
<figure><img src="image_7.png" /><figcaption>**Extended Data Fig. 1** | Typical viscous evolution timescales. We plot the viscous timescale $t_v$ as a function of the viscous parameter $\alpha$ for different belt locations (50, 100 and 150 au) and different gas temperatures (10, 30 and 100 K).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Extended Data Fig. 2 | Evolution of $\dot{M}$ with time.** We plot the temporal evolution of $\dot{M}$ for different belt locations (50 and 100 au) and different initial belt masses (0.1, 1, 10, and 100 $M_\oplus$).</figcaption></figure>
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 3** | GCR for different values of mean molecular weight $\mu$. Temporal evolution of the gas-to-core ratio (GCR) with varying $\mu$. We note that the GCR grows more slowly than expected for the cases $\mu = 7$ and 14 when $\dot{M} = 10^{-2} M_{\oplus}/\text{Myr}$, which is because the theoretical cooling accretion rate becomes smaller than $10^{-2} M_{\oplus}/\text{Myr}$ for lower values of $\mu$ (see Extended Data Fig. 4). We also see that for lower values of $\mu$, the accretion is less efficient from the start because the gas disk scaleheight is higher and less gas is accreted (see Extended Data Fig. 5).
<figure><img src="image_10.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 4** | Potential accretion rate on a planet Vs. available accretion rates. We plot the temporal evolution of the potential theoretical accretion rate on a planet $\dot{M}_{\text{gas}}$ (numerical derivative with $\kappa$ varying with time) for different values of planet semi-major axes $a_{pl}$, atmospheric mean molecular weight $\mu$ and core mass $M_{\text{core}}$. The fiducial model is $a_{pl} = 1$ au, $\mu = 14$, and $M_{\text{core}} = 1 M_{\oplus}$. We overplot horizontal lines with different input rate values of our parameter $\dot{M}$, including the case of $10^{-2} M_{\oplus}/\text{Myr}$ over 100 Myr to verify whether in the cases studied in this paper the theoretical accretion rate is higher than our $\dot{M}$ parameter (which is always the case for $\mu = 28$, large $M_{\text{core}}$ or distant planets). We note that for the $\mu = 28$ case, the green lines become less steep at large $t$. This is because as $t$ increases, one reaches the second regime for which $-a/(bt) \ll 1$, where $M_{\text{gas}}$ scales as $t^{-0.4}$ instead of $t^{-0.7}$ in the other regime (see Eqs 9 and 10). We also note that for the case at 0.1 au (for which $T = 1000$ K), the opacity varies more slowly with $T$ for high enough densities and $\beta$ becomes smaller$^{10}$, hence leading to a higher $\dot{M}_{\text{gas}}$.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 5** | Is a planet accreting all gas flowing through the disk? We plot the Hill sphere radius-to-scaleheight ratio Vs. the distance to host star. For the lowest-mass planets, the Hill sphere radius can be smaller than the disk scaleheight and gas can flow inwards rather than being accreted.
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 6** | GCR and pressures for a Mars-mass and a distant planet. Temporal evolution of the gas-to-core ratio (left) and pressure (right) for a Mars-mass (0.1 M<sub>⊕</sub> at 1.5 au) planet and a distant planet (10 au) with a core mass of 0.5 M<sub>⊕</sub> up to a GCR of 0.5. In the pressure plot (right), the thick solid line is for the Mars-like planet case and the thinner line is for the core of mass 0.5 M<sub>⊕</sub> at 10 au. We note that for the Mars-mass case, GCR grows more slowly than expected when M&#775; = 10<sup>-2</sup> M<sub>⊕</sub>/Myr, which is because the theoretical cooling accretion rate becomes smaller than 10<sup>-2</sup> M<sub>⊕</sub>/Myr in this case (see Extended Data Fig. 4).
<figure><img src="image_16.png" /><figcaption>Extended Data Fig. 7 | Hill Vs. Bondi radii. Hill and Bondi radii Vs. planet mass for different planet semi-major axes.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 8** | GCR and pressures for a planet starting with an Earth atmospheric mass. Temporal evolution of the gas-to-core ratio (left) and pressure (right) starting from a pre-existing atmosphere with an Earth atmospheric mass.