# Interpreting the Atmospheric Composition of Exoplanets: Sensitivity to Planet Formation Assumptions

Paul Mollière, Tamara Molyarova, Bertram Bitsch, Thomas Henning, Aaron Schneider, Laura Kreidberg, Christian Eistrup, Remo Burn, Evert Nasedkin, Dmitry Semenov, et al.

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

## To cite this version:

Paul Mollière, Tamara Molyarova, Bertram Bitsch, Thomas Henning, Aaron Schneider, et al. Interpreting the Atmospheric Composition of Exoplanets: Sensitivity to Planet Formation Assumptions. The Astrophysical Journal, 2022, 934, 10.3847/1538-4357/ac6a56. insu-03874922

HAL Id: insu-03874922

https://insu.hal.science/insu-03874922

Submitted on 28 Nov 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

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

# Interpreting the Atmospheric Composition of Exoplanets: Sensitivity to Planet Formation Assumptions

Paul Mollière<fnref n="1" />, Tamara Molyarova<fnref n="2" />, Bertram Bitsch<fnref n="1" />, Thomas Henning<fnref n="1" />, Aaron Schneider<fnref n="3" />,<fnref n="4" />, Laura Kreidberg<fnref n="1" />, Christian Eistrup<fnref n="1" />, Remo Burn<fnref n="1" />, Evert Nasedkin<fnref n="1" />, Dmitry Semenov<fnref n="1" />,<fnref n="5" />, Christoph Mordasini<fnref n="6" />, Martin Schlecker<fnref n="1" />, Kamber R. Schwarz<fnref n="1" />, Sylvestre Lacour<fnref n="7" />,<fnref n="8" />,<fnref n="9" />, Mathias Nowak<fnref n="10" />,<fnref n="11" />, and Matthäus Schulik<fnref n="12" />

*Received 2021 December 10; revised 2022 April 22; accepted 2022 April 24; published 2022 July 26*

## Abstract

Constraining planet formation based on the atmospheric composition of exoplanets is a fundamental goal of the exoplanet community. Existing studies commonly try to constrain atmospheric abundances, or to analyze what abundance patterns a given description of planet formation predicts. However, there is also a pressing need to develop methodologies that investigate how to transform atmospheric compositions into planetary formation inferences. In this study we summarize the complexities and uncertainties of state-of-the-art planet formation models and how they influence planetary atmospheric compositions. We introduce a methodology that explores the effect of different formation model assumptions when interpreting atmospheric compositions. We apply this framework to the directly imaged planet HR 8799e. Based on its atmospheric composition, this planet may have migrated significantly during its formation. We show that including the chemical evolution of the protoplanetary disk leads to a reduced need for migration. Moreover, we find that pebble accretion can reproduce the planet's composition, but some of our tested setups lead to too low atmospheric metallicities, even when considering that evaporating pebbles may enrich the disk gas. We conclude that the definitive inversion from atmospheric abundances to planet formation for a given planet may be challenging, but a qualitative understanding of the effects of different formation models is possible, opening up pathways for new investigations.

*Unified Astronomy Thesaurus concepts: Exoplanet atmospheric composition (2021); Exoplanet formation (492)*

## 1. Introduction

The distribution of bulk planetary properties such as mass, radius, and orbital parameters encodes critical information that constrains planet formation models (see, e.g., Ida & Lin 2004; Alibert et al. 2005; Mordasini et al. 2009; Hasegawa & Pudritz 2011; Lambrechts & Johansen 2012; Bitsch et al. 2015; Nayakshin & Fletcher 2015; Cridland et al. 2016; Emsenhuber et al. 2021; Schlecker et al. 2021). In addition, the chemical composition of planet atmospheres has long been regarded as a key to unlocking the process of planet formation (e.g., Gautier & Owen 1989; Owen & Encrenaz 2003) and is the explicit goal of many atmospheric characterization studies (see, e.g., Madhusudhan 2019, for a recent review).

This is because the chemical abundances of planetary atmospheres are highly complementary to bulk planetary parameters: they likely relate to the composition of the planetary building blocks in the protoplanetary disk, be it planetesimals, pebbles, or

gas. The composition of the building blocks is determined by disk processes, while their relative importance and accretion location for a given planet are determined by the process of planet formation. Consequently, there have been a number of studies that investigate how planet formation may set the composition of an exoplanet, focusing on the planetary carbon-to-oxygen number ratio (C/O), nitrogen content, content in refractory material,<fnref n="13" /> or just overall metal content (e.g., Öberg et al. 2011; Fortney et al. 2013; Madhusudhan et al. 2014, 2017; Marboeuf et al. 2014; Cridland et al. 2016; Mordasini et al. 2016; Khorshid et al. 2021; Lothringer et al. 2021; Schneider & Bitsch 2021a).

Combining observatories such as the Hubble Space Telescope (HST) and Spitzer Space Telescope allowed for a first look at atmospheric C/O values, albeit with large uncertainties (e.g., Line et al. 2014; Benneke 2015; Brewer et al. 2017), or leading to controversial findings, as in the case of WASP-12b, which was claimed to be either carbon or oxygen rich, with finally a firm water detection in transit pointing toward C/O ≲ 1 (Madhusudhan et al. 2011; Crossfield et al. 2012; Swain et al. 2013; Line et al. 2014;

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

<footnote n="1">¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2">²Institute of Astronomy, Russian Academy of Sciences, 48 Pyatnitskaya St., Moscow, 119017, Russia</footnote>

<footnote n="3">³Centre for ExoLife Sciences, Niels Bohr Institute, Øster Voldgade 5, DK-1350 Copenhagen, Denmark</footnote>

<footnote n="4">⁴Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001, Leuven, Belgium</footnote>

<footnote n="5">⁵Department of Chemistry, Ludwig-Maximilians-Universität, Butenandtstraße 5-13, D-81377 Munich, Germany</footnote>

<footnote n="6">⁶Center for Space and Habitability, Universität Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="7">⁷ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="8">⁸Max Planck Institute for extraterrestrial Physics, Giessenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="9">⁹European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching, Germany</footnote>

<footnote n="10">¹⁰Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="11">¹¹Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="12">¹²Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK</footnote>

<footnote n="13">¹³ Refractories generally encompass all those chemical species with a high condensation temperature, such that they are found in the solid phase of the protoplanetary disk, except for at the smallest orbital separations.</footnote>
Stevenson et al. 2014; Benneke 2015; Kreidberg et al. 2015). Studying the bulk atmospheric enrichment of exoplanets, mostly based on water detections in HST WFC3 spectra, has also been attempted, but the community would clearly benefit from data with higher signal-to-noise ratio and larger spectral coverage to improve abundance constraints, break degeneracies with clouds, and probe additional atmospheric absorbers (e.g., Kreidberg et al. 2014; Fisher & Heng 2018; Wakeford et al. 2018; Pinhas et al. 2019; Welbanks et al. 2019). We note here (and discuss later) that a connection between atmosphere and bulk planet composition is far from trivial (also see, e.g., the recent discussion in Helled et al. 2022, and references therein).

Luckily, the quality of observational constraints on the planetary composition is expected to be rapidly improving. The advent of retrieval methods for medium- and high-resolution observations (e.g., Brogi et al. 2017; Brogi & Line 2019; Gibson et al. 2020) may allow us to constrain the atmospheric volatile and refractory content (Lothringer et al. 2021) from the ground and to trace even isotopologues (Mollière & Snellen 2019). For example, Gandhi et al. (2019), Pelletier et al. (2021), and Line et al. (2021) used high-resolution retrievals to constrain planetary C/O values, while the medium-resolution retrievals of Zhang et al. (2021) indeed revealed isotopologues for the first time. Moreover, recent observations of the GRAVITY instrument at the Very Large Telescope Interferometer (VLTI) have led to some of the most precise constraints on C/O for planetary-mass objects to date (Gravity Collaboration et al. 2020; Mollière et al. 2020). Most importantly, the recent launch of the James Webb Space Telescope (JWST) and the launch later in the decade of ARIEL are expected to lead to excellent constraints on the C/O ratio, especially for transiting planets, and may probe the nitrogen and refractory content of cool planets (e.g., Greene et al. 2016; Wang et al. 2017; Daniëlski et al. 2018; Tinetti et al. 2018). With these next-generation telescopes the focus will likely shift from observational uncertainties to uncertainties in the models for atmospheric characterization (e.g., Feng et al. 2016; Line & Parmentier 2016; Blecic et al. 2017). The development of new characterization techniques is therefore necessary for interpreting future observations. This work has already begun (e.g., Caldas et al. 2019; Feng et al. 2020; Lacy & Burrows 2020; MacDonald et al. 2020; Pluriel et al. 2020; Taylor et al. 2020; Changeat et al. 2021; MacDonald & Lewis 2022; Nixon & Madhusudhan 2022).

Now that the atmospheric abundance constraints may become more precise than ever before, it is timely to revisit the justification stated for many observational campaigns: How can a planet's formation history *actually* be constrained, given its atmospheric abundances, and how well? What are the actual formation quantities that are constrainable? What are the major obstacles that would need to be overcome in case this is not possible? And, lastly, if such an inversion process is challenging for a single planet, could the distribution of abundance patterns be used to constrain some aspects of planet formation?

Our study aims at addressing some of the questions stated above. Specifically, we discuss planet formation and its complexities in the context of the inversion challenge (Section 2). In Section 3, we present a methodology that may prove useful for assessing the consequences of a given formation model choice, where we use a nested sampling method to constrain formation parameters, given the atmospheric composition, for different model assumptions. We show example

<figure><img src="image_5.png" /><figcaption>**Figure 1.** C/O distribution of the solid and gas component in a smooth, static, circumstellar disk assuming equilibrium condensation. The locations of the H₂O, CO₂, and CO ice lines are indicated in the plot.</figcaption></figure>

applications, namely, how chemical disk evolution, or pebble drift, evaporation, and accretion, may affect the inferred formation and migration history of the planet HR 8799e. Our method can be used to qualitatively understand differences between planet formation implementations. In Section 4 we summarize which molecular and atomic species can and will be probed by atmospheric observations, as well as how these may serve to broadly inform the process of planet formation. We end with a short discussion and summary of our study in Section 5.

## 2. The Complexity of the Planet Formation Problem

The idea of using planetary composition to constrain planet formation gained traction in the field of exoplanets with the seminal paper by Öberg et al. (2011). Here the authors propose that the C/O value derived from a planet's composition could be used to constrain where in a protoplanetary disk it formed. The general idea is outlined in Figure 1, very similar to the original Figure 1 in Öberg et al. (2011). Assuming a smooth, static, one-dimensional disk, the authors calculated where important volatile gases such as H₂O, CO₂, and CO (sorted by decreasing condensation temperature) freeze out, if present. Because for the temperature gradient in a protoplanetary disk it holds that $dT/dr < 0$, where *r* is the distance from the star, H₂O freezes out first when moving outward, followed by CO₂ and CO. This directly affects the C/O values in the gas and solid phases because water, for example, removes oxygen from the gas phase, when condensing.<fnref n="14" />

The idea, then, is that if the planetary C/O and overall metallicity (here: C and O content) are known, it is possible to determine where in the disk a planet has formed. This method of determining the process of planet formation has since been cited in virtually every study that aims at constraining the

<footnote n="14">¹⁴ We note that we assumed that 33% of all O and 38% of all C is contained in the refractory solids in the example setup shown in Figure 1, which leads to a high solid-phase C/O inside the H₂O ice line.</footnote>
atmospheric composition of an exoplanet. This also has to do with the comparative ease with which the atmospheric C/O value may be constrained, as we discuss in Section 4. The general idea presented in Öberg et al. (2011) is powerful, but it is undeniably so that planet formation is more complicated than assumed in their study. In the following we give a summary of processes that may have to be taken into account, and inverted, when trying to connect planet composition to formation in practice.

## 2.1. Disk Elemental Composition and Structure

Constraining a planet's formation location based on its elemental abundance ratios (e.g., C/O) requires that the elemental composition of the protoplanetary disk is known. A good starting point may be to assume that the protoplanetary disk has an elemental composition that is identical to that of the host star (planet formation may deplete stellar photospheres in metals with respect to the disk, however; see Chambers 2010; Bitsch et al. 2018a; Booth & Owen 2020). It is therefore crucial to have knowledge about the host star's abundances that is as complete as possible, in excess of just [Fe/H], or to at least use existing scaling relations to approximate the stellar metal content for elements other than iron (e.g., Bitsch & Battistini 2020). The disk composition and assumed mass also set an upper limit on the amount of metals that a planet may accrete during its formation (e.g., Baraffe et al. 2008). Alternatively, the retrieved metal content of an exoplanet may also be used to place a lower limit on the disk metal content, or even total mass, analogous to the concept of the minimum-mass solar nebula (e.g., Hayashi 1981).

Moreover, the disk's physical and thermal structure is important to set the radial fractionation of elements into different molecular species, in solid and gaseous form, that can be accreted by a forming planet. The disk structure depends on the assumed (dust) opacities and therefore the dust evolution (Schmitt et al. 1997; Birnstiel et al. 2016; Savvidou et al. 2020). Moreover, disks are viscously (Lynden-Bell & Pringle 1974) evolving (or due to photoevaporation and disk winds; see Clarke et al. 2001; Bai et al. 2016; Suzuki et al. 2016; Chambers 2019), prone to different instabilities (e.g., Flock et al. 2017; Klahr et al. 2018), and will be affected by the presence of (especially massive) planets that may induce spiral density perturbations, lead to the formation of vortices, or open gaps (e.g., Lin & Papaloizou 1986; Crida et al. 2006; Lobo Gomes et al. 2015; Pinilla et al. 2015; Binkert et al. 2021).

Another important effect determining the disk gas composition is the evaporation of pebbles inside of ice lines, which may significantly increase the local volatile content of gas in the disk (Piso et al. 2015; Booth & Ilee 2019; Schneider & Bitsch 2021a). This could lead to planets being enriched in volatiles much more than expected from pure gas accretion in the classical Öberg et al. (2011) setup. It is important to mention that the dynamics of pebbles is likewise determined by the disk structure. The structure sets the pebbles' growth rates, their Stokes parameters (i.e., drift speeds), and may trap pebbles into local pressure maxima (e.g., Paardekooper & Mellema 2006; Ormel & Klahr 2010; Birnstiel et al. 2012; Lambrechts et al. 2014).

## 2.2. Disk Chemistry

The chemical composition of the protoplanetary disk (e.g., Henning & Semenov 2013) is of central importance for

determining the composition of planetary building blocks, that is, the disk's gas and solid phases. In Öberg et al. (2011) and the more recent Öberg & Wordsworth (2019) study, a simplified and static disk chemical model is assumed. In practice, the disk's chemical composition will evolve because both the disk gas and the volatile ices on grain surfaces will undergo chemical processing (e.g., Eistrup et al. 2016, 2018; Molyarova et al. 2017). This means that chemical reactions between atoms and molecules in both gas and ice may alter which molecular species are the dominant carriers of elements such as C, O, and N over time. This is an important effect and is expected to alter the inferred formation history of a planet, as has been pointed out by Eistrup et al. (2016). Examples are the conversion of CO into CO₂ ice over time, or the conversion of N₂ gas into NH₃ ice (Semenov & Wiebe 2011; Molyarova et al. 2017; Eistrup et al. 2018). Another example of the importance of disk chemistry on planet formation are the processes that lead to the observed carbon depletion in the inner solar system (e.g., Mordasini et al. 2016; Cridland et al. 2019), which may be attributable to the irreversible chemical destruction of carbon grains within a disk's so-called soot line, or connected to chondrule formation (e.g., Kress et al. 2010; Gail & Trieloff 2017; van't Hoff et al. 2020; Li et al. 2021).

The disk chemistry itself is sensitive to many processes, such as stellar evolution (Miley et al. 2021), or the cosmic-ray ionization rate (Eistrup et al. 2016; Schwarz et al. 2019). Moreover, whether or not the initial composition of the disk matter is molecular, that is, “inherited” from the composition of the natal molecular cloud, or elemental (“reset” scenario) can strongly influence the C/O values (Eistrup et al. 2016). The disk's physical structure may also have an impact on its chemical evolution. As an example, we highlight the effect of the self-shadowing of the disk, allowing for nominally too cool compositions to occur closer to the star than otherwise expected (Ohno & Ueda 2021).

## 2.3. Planet Formation

The idea that planet formation may be constrained through planet composition, as presented in Öberg et al. (2011), conceptually boils down to comparing the planetary C/O value and total metal enrichment to the C/O of the disks' solid and gaseous phases as a function of orbital distance from the star. However, planet formation involves and connects many complex processes. This means that a forming planet cannot accrete arbitrary amounts of gas or solids at (or from) arbitrary positions in the disk.

For example, a limiting factor for planets forming via the core accretion paradigm (e.g., Mizuno 1980; Pollack et al. 1996), specifically when accreting pebbles, is the pebble isolation mass. Pebble accretion, which is the accretion of roughly centimeter-sized solids by the forming planet (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012), can only dominate the solid accretion process until the growing planet reaches this isolation mass $M_{iso}$ (e.g., Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018b), after which pebble accretion stops. This is because the planet induces the formation of a pressure bump in the disk, exterior to its orbit, which traps the inward-drifting pebbles. The isolation mass places an upper limit on the refractory content of a planet. A refractory content higher than allowed by this concept of $M_{iso}$ could point to the importance of accreting planetesimals (e.g., Mordasini et al. 2016; Brügger et al. 2020), unless the planet
formed very close to its star, within the refractories' ice lines (Schneider & Bitsch 2021b).

For planets growing in situ, a planetesimal isolation mass was established by Lissauer & Stewart (1993). It is caused by the planet depleting the local reservoir of planetesimals within the zone of its gravitational influence (the so-called “feeding zone”). In contrast to the pebble isolation mass, the planetesimal feeding zone increases with planet mass. The different ways of how a planet's refractory content and total mass scale in planetesimal and pebble accretion models might thus allow us to put limits on the contributions of the different solid reservoirs for a given planet. The planetesimal isolation mass can also be overcome via (giant) impacts or if a protoplanet migrates into regions of the disk still containing planetesimals (Alibert et al. 2005).

Moreover, the accretion of gas by the growing planet is a three-dimensional process (e.g., D'Angelo et al. 2003; Ayliffe & Bate 2009; Szulágyi et al. 2014; Ormel et al. 2015; Schulik et al. 2019, 2020). This may be especially important, as the gas composition and C/O value are thought to be changing above the midplane of the disk (e.g., Molyarova et al. 2017; Cridland et al. 2020a). Similar to the solid accretion processes, the amount of gas a planet can accrete during formation is limited. In contrast to the solids, however, the ultimately limiting factor is the lifetime of the protoplanetary disk. Once a planet enters runaway accretion, it will accrete gas as quickly as can be provided by the viscously evolving disk, until the disk dissipates. More specifically, it has been shown that gas delivery to the planet can be severely limited by gap formation, which in turn is controlled by the disk's viscous resupply of the planetary feeding zone (e.g., Lubow et al. 1999; Ayliffe & Bate 2009; Lissauer et al. 2009; Bergez-Casalou et al. 2020; Schulik et al. 2020). For embedded planets below the gap opening mass, gas accretion may be limited if radiative cooling is counteracted by hot inflowing gas from the ambient disk. In this case even gas that enters the planetary Hill sphere is not accreted (so-called “recycling”; Cimerman et al. 2017).

Another important process of planet formation is migration. The migration history of a planet is critical to its final composition because it determines where in the disk it accretes material. Because the Öberg et al. (2011) approach strives to ultimately constrain formation locations with respect to the disk ice lines, migration may be less of a problem if a forming planet did not migrate across ice lines. Interestingly, planets forming in the inner disk may actually be trapped at locations just beyond the water ice line (e.g., Bitsch & Johansen 2016; Cridland et al. 2016; Müller et al. 2021). There also exist other traps not connected to ice lines, such as the disk's dead-zone inner edge (Bitsch et al. 2014; Cridland et al. 2016). Qualitatively speaking, and when not currently trapped, planets are expected to migrate either by fast type I migration or via slower type II migration, the latter of which ensues once the planet is massive enough to open a gap in the disk (see, e.g., the review by Baruteau et al. 2016). Quantitatively, there is an ongoing debate about the actual magnitude of the torques (and therefore speed of migration), for example, for type II migration (e.g., Dürmann & Kley 2015; Robert et al. 2018).

Further complicating the picture are *N*-body interactions between the forming planets. This process is now regularly included in models of planet formation and population syntheses but comes at an increased numerical cost (e.g., Alibert et al. 2013; Chambers 2016; Lambrechts et al. 2019; Emsenhuber et al. 2021; Izidoro et al. 2021). While increasing

the complexity of describing planet formation, *N*-body interactions may represent an alternative avenue for producing hot Jupiters, which may have been scattered in by planets farther out (e.g., Jurić & Tremaine 2008; Bitsch et al. 2020). *N*-body interactions can also alter the composition of a planet via giant impacts. The amount of solids brought into a giant planet via such impacts might be substantial or even dominant compared to the amount accreted from small bodies like planetesimals or pebbles (e.g., Emsenhuber et al. 2021; Ginzburg & Chiang 2020; Ogihara et al. 2021). If the impactors originally formed in clearly different regions of the disk than the forming planet, this would blur the meaning of a well-defined formation location of a planet. Accounting for the effects of *N*-body interactions when trying to invert planet formation thus seems challenging. Interestingly, it has recently been shown that while *N*-body interactions tend to randomize the process of planet formation, machine-learning techniques such as random forests still allow us to predict the outcome of planet formation quite accurately (Schlecker et al. 2021). In this work the authors show that the initial parameters of the planet formation model described in Emsenhuber et al. (2021), mainly the initial location of the planetary embryo and the dust mass of the disk, may be used to predict which class a forming planet will belong to (super-Earths, Neptune-like, giant planets, etc.). These classes also correspond to certain orbital distances and planetary compositions.

The discussion here focused mostly on planet formation via the core accretion paradigm. Other aspects are of importance if a planet forms via gravitational instability (GI; see Boss 1997). The disk structure (and thus formation environment of the planet) will be quite different in this case. This is because GI planets may form early, when the disk is still massive, in the outer parts of the disk (Boss 2021; Schib et al. 2021). We note that while GI is classically regarded as a way to produce wide-separation gas giant planets, it has also been suggested to allow for the formation of less massive, small-separation planets (Nayakshin 2010), caused by fast inward migration after formation and the associated mass loss, dubbed “tidal downsizing.”

## 2.4. Planetary Bulk—Atmosphere Coupling

When aiming to constrain planet formation based on the results of atmospheric abundance characterizations, one must make assumptions on how the atmospheric composition relates to the bulk composition of a planet. It has been pointed out that planets growing via core accretion may have a layer of heavily metal-enriched gas above their solid cores, due to the evaporation of pebbles and planetesimals that are destroyed when entering the hot planet's proto-atmosphere (e.g., Mordasini et al. 2016; Helled & Stevenson 2017; Brouwers et al. 2018; Brouwers & Ormel 2020). The recent constraints by the Juno spacecraft on the interior of Jupiter are consistent with this assessment, pointing to the existence of a dilute core (Wahl et al. 2017; Debras & Chabrier 2019). We note that the exact distribution of metals in Jupiter' interior is difficult to explain, however, and may involve a formation process stretching over 2 Myr (see Helled et al. 2022, for a recent review). What is more, the planet bulk metallicity constraints for transiting planets derived in Miller & Fortney (2011), Thorngren et al. (2016), and Thorngren & Fortney (2019) tend to be higher than the metallicities reported for planetary atmospheres (albeit with large uncertainties; see Welbanks
et al. 2019), making an enrichment of the interior with respect to the atmosphere a likely scenario, and allowing for a first estimate regarding the efficiency of planetary mixing.

The question that naturally arises from these findings is how representative the inferred atmospheric composition is of the bulk of the planet, even for giant planets, where the gas dominates the mass budget. For this it needs to be understood if and how well the metals can be mixed throughout the planetary envelope. Whether this happens at all is not clear, as a gradient in metallicity (and therefore mean molecular weight) may stabilize the planetary interior against convective motions (e.g., Ledoux 1947; Chabrier & Baraffe 2007; Leconte & Chabrier 2012). How a solid core or compositional gradients tend to mix throughout the planetary envelope for gas giant planets, often taking Jupiter as an example, is currently being investigated (Vazan et al. 2015, 2016; Moll et al. 2017; Vazan et al. 2018; Müller et al. 2020; Ormel et al. 2021). The results reported in these studies do not yet agree, predicting either fully mixed envelopes or ones where a metal gradient (and core) persists in the planet. In addition to imperfect mixing, rainout processes also likely play an important role, potentially depleting both solar system and exoplanet atmospheres in metals (e.g., Spiegel et al. 2009; Wilson & Militzer 2010).

From the discussion above it becomes clear that the metal enrichment inferred from the atmospheric retrievals, which serve as an input for any formation analysis, is only a lower limit for the true planetary metal enrichment. As long as the metals are locked into the invisible interior of the planet homogeneously, and not selectively, this may still allow us to constrain the solids' location of origin in the protoplanetary disk, as long as the atmosphere is still enriched enough for solid accretion to be the dominating factor. In the case where the relative elemental composition in the atmosphere (except for H and He) is different from the deep interior, or when the atmosphere is depleted with respect to the deeper interior to the point where it mimics the atmospheric metallicities expected from pure gas accretion, this poses a problem. Depending on how planet formation ensued, the former case may still trace the origin of the solids that were accreted toward the end of the formation process. The latter case could be resolved by checking the relative atomic abundances in the presumed metal-poor planet. If refractory atomic species are relatively abundant, this may point to a dominant accretion of solids that are mostly locked into the planet's interior. An interesting recent discussion of how to constrain planet formation in the case where the atmospheric and planetary bulk compositions differ can be found in Helled et al. (2022), who come to similar conclusions. We note again that these avenues for analyzing the origin of the metals in a planet's atmosphere will be further complicated if volatile-rich gas from evaporated pebbles was accreted by the planet.

## 2.5. Atmospheric Evolution

The atmospheric composition of an exoplanet may also evolve owing to atmospheric evaporation, or by secular enrichment of infalling comets and asteroids. Evaporation is especially important for close-in low-mass planets (e.g., Jin & Mordasini 2018). The atmosphere may be partially or fully lost owing to thermal or nonthermal processes, where the thermal escape separates into the regimes of Jeans escape or hydrodynamical escape, depending on the local thermal state of the atmosphere (see, e.g., Barman 2018). For the atmosphere

to become relatively enriched in metals by evaporation, two criteria have to be met. First, the atmosphere needs to be of low enough mass to allow for a significant amount to be lost. Second, the atmospheric escape process needs to preferentially retain the heavier atmospheric species, for which the atmosphere needs to be in the Jeans escape regime (e.g., Bourrier & Lecavelier des Etangs 2018). In the hydrodynamic escape regime the heavier metal species would be lost together with hydrogen. We note that this transition is gradual, however, and that there can be mass fractionation in hydrodynamic outflows as well, depending on the magnitude of the total mass flux (e.g., Hu et al. 2015). These authors also report that such outflows can selectively deplete the atmospheres of Neptune- and sub-Neptune-mass planets in hydrogen over multiple Gyr, provided that the initial atmospheric mass is small enough ($<10^{-3}$ of planetary mass). For gas giant planets evaporation may thus be a less relevant process for changing the atmospheric composition. An extreme case that is worth mentioning is the Roche lobe overflow that may affect the closest-in gas giant planets. This process could strip away the upper gaseous envelope, potentially revealing the more metal-enriched layers below. Roche lobe overflow has been discussed as the potential origin of LTT 9779b, a planet in the hot Neptune desert (Jenkins et al. 2020). By extension, significant atmospheric evaporation may lead to similar outcomes if a compositional gradient is present in the atmosphere.

Another possibility of atmospheric evolution is the secular contamination of the planetary atmospheres by infalling comets or asteroids. Here the frequency of cometary impacts and the persistence of the enrichment they cause in the atmosphere need to be estimated. Turrini et al. (2015) find that the additional water a comet may deposit in the visible atmosphere of the hot Jupiter HD 189733b would have to persist for 500–5000 yr before being removed, assuming impacts of kilometer-sized comets every 20–200 yr; otherwise, no significant enrichment is possible. A quantitative assessment of cometary enrichment requires an estimate of the persistence timescale, however. As this appear to be lacking from the literature, we present a simple first-order analysis below. We start by approximating local mixing timescales as $\tau_{mix} = H_p^2/K_{zz}$, where $H_p$ is the planetary pressure scale height and $K_{zz}$ the atmospheric eddy diffusion coefficient, and find

$$ \tau_{\text{mix}} = 0.3 \left( \frac{H_p}{300 \text{ km}} \right)^2 \left( \frac{10^8 \text{ cm}^2 \text{ s}^{-1}}{K_{zz}} \right) \text{yr.} \quad (1) $$

A $K_{zz}$ value of $10^8 \text{ cm}^2 \text{ s}^{-1}$ within the radiative zone is well within the estimates by passive tracers reported from model calculations for HD 189733b and HD 209458b (Parmentier et al. 2013; Agúndez et al. 2014). For self-luminous planets the correct value to be chosen is unclear, but $K_{zz} = 10^5 \text{ cm}^2 \text{ s}^{-1}$ may at least be a useful lower limit (Ackerman & Marley 2001). In the deeper regions of the atmosphere convective overshoot may lead to higher values for $K_{zz}$, smoothly transitioning toward the value expected for fully convective atmospheres as the radiative-convective boundary (RCB) is approached (e.g., Ludwig et al. 2002; Helling et al. 2008).

The $\tau_{mix}$ estimate of less than a year in Equation (1) thus shows that mixing potentially proceeds on faster timescales than those quoted by Turrini et al. (2015), at least locally, meaning that any material added by cometary impacts should be mixed away quickly. This does not preclude a slower,
homogeneous enrichment of the radiative atmosphere by cometary impacts over time. However, the question is how quickly this enrichment will be removed into the bulk interior of the planet by entrainment into overshooting convective blobs at the RCB. By extension, the enrichment of the visible atmosphere may be lower if the impactors deposit their metals below the RCB: for the bulk of Jupiter, for example, the mixing timescale has been estimated to be at most a few years (Debras & Chabrier 2019).

Assuming diffusive mixing in a 1D atmosphere, we derive the following expression for the increase in mass fraction X of a certain species in the planetary atmosphere, due to cometary impacts:

$$ \Delta X = \frac{\dot{M}_g H_p^2 P_{\text{RCB}} - P_i}{4\pi R_p^2 K_{zz} P_{\text{RCB}} P_i} \quad (2) $$

Here $\dot{M}$ is the mass accretion rate of comets, g the planetary gravity, $R_p$ the planetary radius, $H_p$ the atmospheric pressure scale height, and $P_i$ the average location of the destruction of impacting comets in units of pressure. $P_{\text{RCB}}$ is the location of the RCB or, alternatively, the location where the deep mixing becomes fast enough to make the local mass fraction equilibrate to the average of the planetary interior (which is assumed to be well mixed). In any case it must hold that $P_i < P_{\text{RCB}}$ for this expression. A derivation can be found in Appendix A. When assuming pure-water comets and tracking the change in the mass fraction of water, a relative enrichment of $\Delta X/X_0 = 3 \times 10^{-4}$ to $3 \times 10^{-2}$ is found, where $X_0$ is the planet's bulk water mass fraction. See Appendix A for more details and which values to assume for the various quantities. We therefore conclude that secular cometary enrichment of the planetary atmosphere may be unlikely for giant planets, but a better modeling of the process, for example, the lower boundary treatment when using a 1D diffusion approximation, is needed.

### 3. From Planet Composition to Formation Outcomes

From Section 2 it is obvious that a full inversion, from atmospheric composition to planet formation, is still a long way off. However, for qualitatively understanding the ramifications of different planet formation assumptions, it would be useful to possess a tool that compares the inverted outcomes of such models. In this case the effect of a given process may be studied in isolation, allowing the user to get an intuition for the importance of a given assumption. This is in contrast to attempting to invert a full formation model, which may either be too numerically costly or require too many parameters when compared to the limited number of observational characteristics.

In what follows, we will demonstrate such an analysis setup by starting with the inversion of the formation model used in Öberg et al. (2011) and Öberg & Wordsworth (2019) and applying it to the compositional constraints obtained for the directly imaged planet HR 8799e. As a second step, we will introduce either the effects of including chemical evolution of the protoplanetary disk or the effect of pebbles that drift and evaporate in the disk. Comparing the results of these setups for HR 8799e serves to highlight the likely importance of disk chemical evolution for its inferred migration history and studies whether pebble accretion may have been a likely scenario for

this planet. We end this section by suggesting other toy model setups, testing for the influence of various formation model complexities described in Section 2.

#### 3.1. Formation Model Inversion

Formation models produce synthetic planets, or populations thereof, starting from a physical model and set of formation parameters. These can be, for example, the initial disk mass, the disk composition, and the starting position of a planetary embryo in a disk. When attempting to constrain planet formation based on measured planetary compositions, the formation models need to be inverted, because planet composition is an outcome of the formation models. More specifically, if assuming that planetary compositions

$$ C = (\text{O/H}, \text{C/H}, \text{N/H}, \text{Fe/H}, \text{Si/H}, \text{P/H}, \text{S/H},...) \quad (3) $$

can be measured with a given uncertainty, we are then interested in the probability distribution of formation parameters $\vartheta$ of a formation model $\mathcal{M}$, given this measurement. Using Bayes's theorem, this can be written as

$$ P(\vartheta|\mathcal{C}, \mathcal{M}) = \frac{P(\vartheta|\mathcal{M})P(\mathcal{C}|\vartheta, \mathcal{M})}{P(\mathcal{C}|\mathcal{M})} \quad (4) $$

Here $P(\vartheta|\mathcal{M})$ is the prior probability of $\vartheta$ before considering any data $\mathcal{C}$, while $P(\mathcal{C}|\vartheta, \mathcal{M})$ is the likelihood for observing $\mathcal{C}$, given that $\vartheta$ is true. In practice, this may be written as, for example,

$$ \log P(\mathcal{C}|\vartheta, \mathcal{M}) = - \sum_{i=1}^{N_{\text{species}}} [C_i - C_i^{\mathcal{M}}(\vartheta)]^2 / 2\Delta C_i^2 + \text{cst} \\ = -\chi^2/2 + \text{cst}, \quad (5) $$

where $N_{\text{species}}$ is the number of measured atmospheric species, $C_i^{\mathcal{M}}(\vartheta)$ is the formation model prediction of the planet abundance of species i, and $\Delta C_i$ are the measurement uncertainties. Here we choose a simple form of the likelihood for clarity, assuming that the measured abundances of different atmospheric constituents are independent and follow a Gauss distribution. In general, the functional form of the log-likelihood can be arbitrarily complicated. For example, the abundance posterior of an atmospheric retrieval may be used directly, which can be approximated by, say, a Gaussian mixture model. In our application in Section 3 we chose an intermediate step, accounting for the covariance between the atmospheric oxygen and carbon content. We note that such an inversion process does not necessarily only need to consider elemental abundances as input measurements. Any observed property of a planet, such as its orbital parameters, could in principle be included in this analysis, as long as it is predicted by a formation model.

In practice, we will compute samples of the target distribution $P(\vartheta|\mathcal{C}, \mathcal{M})$ by numerically integrating the numerator of the right-hand side of Equation (4) using the so-called nested sampling method (Skilling 2004). In short, nested sampling is a Monte Carlo technique to integrate functions in highly dimensional model parameter spaces. Here it is the space spanned by the formation parameters. When integrating the numerator of Equation (4), the integral value resulting in the so-called model evidence, nested sampling will automatically generate samples of our target distribution. In principle,
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Schematic outlining the mode of operation of a formation model inversion. Exoplanet observations first result in inferred atmospheric compositions. This composition is then used as input for the formation model inversion, which uses nested sampling to compare the inferred composition to the prediction of a toy planet formation model. The resulting posterior distribution of the formation parameters represents constraints of the planet formation process. Running inversions for various toy formation models then allows us to study the impact of differing model assumptions.</figcaption></figure>

model evidence ratios allow us to distinguish between different formation models. However, as long as we cannot invert full state-of-the-art formation models, this may be possible only when considering sets of assumptions that lead to wildly different outcomes, and where one clearly represents a better fit. We use the PyMultiNest (Buchner et al. 2014) package for inverting formation models, which is a Python wrapper of the MultiNest code (Feroz et al. 2019). A schematic illustrating our approach is shown in Figure 2.

### 3.1.1. Inverting the Öberg et al. (2011) Model

As a first toy model we will use arguably one of the simplest formation models conceivable when aiming at studying the usefulness of atmospheric C/O constraints. For this we follow the setup presented in Öberg et al. (2011) and Öberg & Wordsworth (2019), assuming their static protoplanetary disk model as described for the young solar nebula, but with the ice lines and composition adapted to the planet system of interest. For every volatile species we define a constant mass fraction in relation to the total disk mass. Inside of its ice line the volatile species is in the gas phase, whereas outside it is in the solid phase. This leads to the well-known step-like behavior of the C/O values in the solid and gaseous phase of the disk, as shown in Figure 1. More details on our implementation can be found in Appendix B.1.

We then assume that the planet formation process can be fully described by a set of four parameters, which ultimately allow us to map to the bulk composition of the planet:

$$ \vartheta = (M_P, M_{\text{solid}}, a_{\text{solid}}, a_{\text{gas}}), \qquad (6) $$

where $M_P$ is the total planet mass and $M_{\text{solid}}$ is the mass of the solids (refractory species and volatile ices) accreted by the planet. It holds that $M_P = M_{\text{gas}} + M_{\text{solid}}$, which we use to determine the amount of gas a planet accreted. The parameters $a_{\text{solid}}$ and $a_{\text{gas}}$ denote the orbital distances where the solids and gas were accreted, respectively.

For a given value of $\vartheta$ our setup then uses the disk model to determine the planet's accretion location with respect to the disk ice lines. Species in the gas phase will be used to determine the composition of the accreted gas; the analogous is done for the solids. Together with the mass fractionation between solids and gas in the planet, this determines the planet's composition $\bar{C}$.

### 3.1.2. Adding Chemical Time Evolution

The Öberg et al. (2011) disk setup is convenient for conceptually studying the usefulness of C/O. As discussed in Section 2, there are many complicating factors that make a true inversion from C/O to planet formation parameters extremely challenging. In our second toy model setup we single out one of these processes and study how chemical evolution in the protoplanetary disk changes inferences when compared to the static disk model.

As discussed before, chemical reactions may process the initial disk abundances over time, shifting carbon and oxygen atoms to different chemical species. For example, CO, which is very volatile, can be removed from the gas by surface reactions, processing it into less volatile $CO_2$ (Molyarova et al. 2017; Bosman et al. 2018; Eistrup et al. 2018; Schwarz et al. 2018). If a large fraction of CO gas in a given region of the disk midplane undergoes such processing into $CO_2$ ice, then an amount of elemental C and O equivalent to that initially carried in CO will be processed from the gas into the ice in this region. This is one example of how chemical processing in disks can alter the elemental partitioning of, for example, C and O, in the gas and ice. These processes are thus expected to influence inferences on the location of planet formation as first discussed in Eistrup et al. (2018).

To study the effect of disk chemical evolution in formation inversions, we replaced the static disk abundance model with a time-dependent model. For this we calculated the evolution of the disk chemical composition using the ANDES code, which describes a quasi-stationary 2D axisymmetric protoplanetary disk (Akimkin et al. 2013; Molyarova et al. 2017). ANDES solves for the time-dependent chemical composition of the disk with a detailed description of grain surface and gas-phase reactions; see Appendix B.2 for more details. For the initial disk setup we use the abundances given in Table 2, that is, equal to the disk abundances used for our static disk model. Although the chemical model includes other elements, we only consider H-, He-, C-, O-, and N-bearing species in the calculations. The refractories are considered to be chemically inert and condensed at all times. An example for the resulting time evolution of the C/O ratio in the disk is shown as the yellow to dark-blue lines in Figure 3. The process of planet formation is then modeled in the same way as in the static disk case, but with the difference that the composition of the accreted gas and solids is taken from the disk's chemical evolution calculations. Formally this also allows us to add the formation time as an additional parameter to be constrained, but
<figure><img src="image_7.png" /><figcaption>**Figure 3.** C/O values in the protoplanetary disk's gas and solid phases, shown as dashed and solid lines, respectively. C/O ratios shown in red are obtained from the Öberg et al. (2011) disk model (static disk composition). C/O ratios shown in colors changing from yellow to dark blue are obtained from the model including the chemical evolution of the protoplanetary disk.</figcaption></figure>

here we chose to initially only study and compare the inversion outcomes for different times during chemical evolution.

### 3.1.3. *Studying the Effect of Pebble Accretion*

In addition to the two setups described above, we study the effect of pebbles drifting and evaporating in the protoplanetary disk and their accretion, as a third scenario. In general, pebbles will quickly drift toward the central star in an unperturbed disk because of the torque exerted by the head wind of the disk gas. This wind is caused by the radial pressure support of the gas, which makes it orbit the star at sub-Keplerian speeds. Therefore, unless there are local pressure enhancements that trap pebbles (caused by, e.g., a planet), pebbles will drain into the inner parts of the disk quickly, releasing copious amounts of volatiles into the gas phase when crossing their respective ice lines (e.g., Booth et al. 2017; Booth & Ilee 2019; Schneider & Bitsch 2021a). Therefore, in addition to accreting pebbles directly, this process can be crucial for setting the composition of forming planets by gas accretion.

For modeling the effect of pebbles on the formation inversion of HR 8799e, we fed disk compositional structures from the *chemcomp* model (Schneider & Bitsch 2021a) into our inversion framework. In short, *chemcomp* solves for the pebble growth, drift, and evaporation in a viscously evolving disk. In addition, it includes a full planet formation model in the pebble accretion paradigm, handling planet-disk interactions such as gap opening and planet migration. Because *chemcomp* is too slow for inverting it directly, we use its pebble drift and evaporation prescription in an unperturbed disk to obtain a disk compositional structure as a function of time, to which we then apply our inversion framework, identical to the treatment of the disk's chemical evolution.

## 3.2. Application of Toy Model Inversions to HR 8799e

As discussed in Section 1, deriving accurate and precise C/O ratios from current exoplanet observations is challenging, but

this will likely change with JWST. Using high-resolution spectrographs, or the interferometric GRAVITY instrument at the Very Large Telescope, the community is already starting to derive precise C/O values using ground-based instruments. This has to do with the data but also with the use of state-of-the-art retrieval techniques (Brogi & Line 2019; Gandhi et al. 2019; Mollière et al. 2020; Line et al. 2021; Pelletier et al. 2021). Below we will make use of the atmospheric composition derived for HR 8799e from GRAVITY observations (Mollière et al. 2020) and try to constrain how its derived formation history changes when disk chemical time evolution or pebbles are included.

In order to run the formation inversion for HR 8799e, a disk elemental composition needs to be assumed for HR 8799. HR 8799 is a $\lambda$ Boötis-type star; this means that the abundances measured for its iron-peak elements are subsolar, with values of $[Fe/H] = -0.55 \pm 0.10$ (Sadakane 2006) or $[Fe/H] = -0.52 \pm 0.08$ (Wang et al. 2020) having been inferred for iron specifically. A similar depletion is expected for Mg, Si, and other massive iron-peak elements, while the abundances of elements typically found in volatile species (C, N, O) are expected to be close to solar (e.g., Paunzen 2004). Indeed, the latest analysis of Wang et al. (2020) inferred $[C/H] = 0.11 \pm 0.12$, $[O/H] = 0.12 \pm 0.14$, and $(C/O)/(C/O)_\odot = 0.96 \pm 0.19$, which are all consistent with solar but slightly enriched in C and O. According to Wang et al. (2020), the most likely explanation for the observed composition of HR 8799 is recent accretion of volatile-rich material onto the outer layers of the star, for example, from an evaporating hot Jupiter, or of volatile-rich ices scattered into the inner system by the four HR 8799 planets. We will therefore use the composition of the solar system from Öberg & Wordsworth (2019) as our nominal abundance model because it could be unlikely that the star has a bulk elemental composition identical to its observed photospheric values. When relevant, we will also report on how our results change if taking the $\lambda$ Boo abundances of HR 8799 at face value, however. The ice line positions in the HR 8799 disk are set to the values derived from the ANDES chemical evolution model, at $t=0$.

For HR 8799e we make use of the atmospheric retrieval results reported in Mollière et al. (2020). Of relevance for the formation inversion are the derived values for the planetary mass, as well as the atmospheric metallicity and C/O ratio. As the mass is spectroscopically determined, it has large error bars. However, it still results in a constraint on the total amount of solids that the planet incorporated, which can also be estimated from multiplying the atmospheric metallicity by the planetary mass. More specifically, multiplying the planet mass by the inferred atmospheric metallicity results in a lower limit on the metal mass (in solid or gaseous state) that a planet accreted, which may be dominated by solids in cases of high metallicities. We used the actual posterior distributions on planetary gravity and radius from the spectral retrievals to construct the inversion prior for the planetary mass (effectively corresponding to a $1\sigma$ upper limit of $14 M_{Jup}$). We also study the effect of using a tighter mass constraint in the pebble inversion scenario.

Converting the atmospheric C/O ratio for use in the inversion method requires special care. In the spectral retrievals the metallicity is used as a free parameter to scale all elemental abundances except H and He, after which the C/O ratio is set by scaling O. In the formation model the formation location and relative gas-to-solid accretion fraction set the O, C, and
<figure><img src="image_8.png" /><figcaption>**Figure 4.** C/H and O/H distribution of HR 8799e (black ellipses) as obtained from the GRAVITY retrievals reported in Mollière et al. (2020). The three ellipses denote the $1\sigma$, $2\sigma$, and $3\sigma$ uncertainties of the C/H, O/H distributions. Colored ellipses show the distribution of values ($1\sigma$) resulting from running formation inversions as a function of time for different disk scenarios. The slanted dashed lines denote (C/H, O/H) value pairs of constant C/O from 0.1 to 1.2, in steps of 0.1. The solar C/H, O/H value is shown as a blue filled circle. Arrows indicate the direction of increasing atmospheric metallicity and C/O. Top left panel: formation inversion using a static disk model. Top right panel: formation inversion taking into account the chemical evolution of the disk. Bottom panel, from left to right: formation inversions in cases including pebble drift and evaporation. Case (i): replacing the static disk model by the pebble disk model. Case (ii): same as Case (i), but additionally placing a prior on the accreted solid mass, corresponding to a 20 $M_\oplus$ upper limit defined by the pebble isolation mass. Case (iii): same as Case (ii), but including a tighter prior on the mass of HR 8799e, based on a dynamical mass estimate. Case (iv): same as Case (iii), but using a larger value for the pebble isolation mass (100 $M_\oplus$). An arbitrary offset has been applied to the ellipses of Cases (i)–(iv), for clarity.</figcaption></figure>

refractory metal content, from which a C/O ratio can be calculated. Therefore, the C content is no longer strictly coupled to the refractory metal content, in contrast to the spectral retrievals that we use as input for our formation retrievals. This inconsistency has to be kept in mind when we use the C/O and metallicity of the spectral retrieval to obtain atmospheric C/H and O/H values and compare these to the C/H and O/H predicted by the formation model. In general, independently constraining C/H and O/H in atmospheric retrievals is the better avenue for running formation inversion studies. We note that it is also important to account for the amount of oxygen that has been sequestered into atmospheric clouds. Because the C/O constraints from Mollière et al. (2020) include this effect, the atmospheric retrieval results can be used without modification. For C/O constraints from retrievals that constrain absorber abundances independently, corrections would need to be applied.

Moreover, sampling the C/O and metallicity posterior of the spectral retrieval leads to tightly correlated C/H and O/H values, and we take this into account by fitting it with a two-dimensional Gaussian distribution. The distribution's covariance matrix is then used to describe the uncertainties of C/H and O/H during the inversion process. If this were not done, the independent uncertainties in C/H and O/H (as obtained from the diagonals of the covariance matrix) would allow for a spread in C/O values much larger than obtained from the spectral retrieval, rendering a formation inversion meaningless. Taking the spectral retrieval results from Mollière et al. (2020), we used $[\text{Fe/H}] = 0.48 \pm 0.27$ and C/O = $0.60 \pm 0.075$. The

resulting 2D distribution of C/H and O/H is shown in Figure 4.

Lastly, because our adapted Öberg & Wordsworth (2019) setup for the disk abundance led to a slightly subsolar C/O ratio (0.52 instead of 0.55), we scaled the planetary value of 0.6 by 0.52/0.55 prior to generating the planetary C/H, O/H distribution as input for the formation inversion process. This is done to conserve the relative distance in C/O, with respect to solar abundances, of the input planetary composition.

### 3.2.1. Static Disk Chemistry

We start the formation inversion for HR 8799e using the disk model à la Öberg et al. (2011) and Öberg & Wordsworth (2019). The magenta ellipse in the top left panel of Figure 4 shows the distribution of sampled C/H, O/H pairs of the formation inversion, indicating a good fit. The most interesting result from the formation inversion is shown in the left panel of Figure 5: here we plot the 1D posterior of $a_{\text{solid}}$, which is the location where the solids that are enriching the planet originated or were accreted.

The inversion process finds a clear preference for the solids to stem from outside the CO ice line, or from within the H$_2$O ice line. This is intuitively easy to understand: because the spectral retrieval resulted in a superstellar atmospheric metallicity and a C/O ratio consistent with stellar, this means within the Öberg et al. (2011) model that the planet must have accreted solids of stellar C/O ratio. The metal content of the gas, being stellar or substellar, is not high enough to offset the C/O ratio set by the accretion of solids. Accreting solids of stellar C/O is possible at orbital distances outside of the ice
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Left panel: posterior distribution of the location where the solids in HR 8799e were accreted or originate when using the static disk composition in the formation model. Because the planet has an increased metallicity and stellar C/O ratio, its enrichment is likely dominated by solids that originate from outside the CO ice line, as the ices then incorporate all major C- and O-bearing species. HR 8799e's current orbital distance is shown as a blue vertical line, where the ice line positions were computed from chemical disk models around a young HR 8799 host star; see description in text. Right panel: formation inversion of HR 8799e including the chemical evolution of the protoplanetary disk. The ice line positions of NH₃ and N₂ have been omitted for clarity. The x-axis has been log-scaled in both panels.</figcaption></figure>

lines of all major carbon- and oxygen-carrying species, the outermost being CO (see Figure 1). For HR 8799e's current orbital position, ~15 au (Wang et al. 2018), this could mean that the planet underwent some orbital migration after solid accretion, as the CO ice line for a young A5 host star such as HR 8799 is expected to be around 35 au. To illustrate this, we overplot today's orbital location of HR 8799e in the left panel of Figure 5.

Alternatively, due to a high organic carbon content of the refractories in our toy model, a roughly stellar C/O ratio is also attainable if the planet formed within the water ice line and then migrated or scattered outward to its current orbital position. This formation channel for distant giant planets was proposed, for example, in Marleau et al. (2019). While this is an intriguing result, we stress that it is dependent on the disk compositional model we assume and the formation model used in general. We also note that the high carbon content of the refractories in the inner disk may be unlikely; see our discussion in Section 2. The 1D and 2D projection of the full posterior of the formation inversion is shown in the left panel of Figure 12 and discussed in Appendix C.

We also carried out inversions using the λ Boo composition of the star for the disk. This was done by increasing the oxygen and carbon abundance by 30% and decreasing the iron and silicate content of the refractory material by 70%. The oxygen no longer bound in silicates was added to H₂O, which is the dominant reservoir of oxygen in the protoplanetary disk. CO, the second most abundant oxygen reservoir, should not change because the carbon content of the disk is not changed when applying the depletion of the iron-peak elements. Finally, the oxygen abundance is increased until C/O = 0.54 is reached, which is the value reported in Wang et al. (2020). We find that the most likely location of the origin of the accreted solids is still outside the CO ice line. The formerly second likely

location, inside the H₂O ice line, vanished: the solid C/O there is exclusively set by the refractory species, which have a much higher C/O value of 2.6 now because of the strong silicate depletion. The associated posterior is shown in the right panel of Figure 12, in Appendix C.

### 3.2.2. Chemical Disk Evolution

Next, we analyze how robust the above findings are when adding chemical evolution of the disk composition. We thus ran a formation inversion of HR 8799e using the formation model that included chemical evolution. In practice, this was done by inferring the planet formation parameters using the composition of the ANDES disk chemical models, as a function of time. ANDES computes the abundances as a function of altitude above the midplane. We used the resulting surface densities of the disk to determine the composition of the gas and solids. We assumed a young (1 Myr) host star at HR 8799's current mass and $L = 3.58 L_⊙$ (Yorke & Bodenheimer 2008), with a disk that produces an accretion luminosity of 0.233 $L_⊙$ (corresponding to $10^{-8} M_⊙ \text{ yr}^{-1}$), where $L_⊙$ and $M_⊙$ are the solar luminosity and mass, respectively. For a given inversion at time $t$ we assumed that the disk composition is fixed at the value that the chemical evolution predicts at that time. This is an approximation because it implicitly assumes that planet formation happens over a characteristic timescale $<10^5 \text{ yr}$, which is chosen as the time step between the snapshots of the disk composition. Because the goal of the present exercise is to study the zeroth-order effect that chemical evolution may have, we deem this approximation acceptable. Future applications could incorporate the time of formation as another free parameter, also assuming (or trying to infer) the duration of the planet formation process.

The results of the inversion including disk chemical evolution are shown in the top right panel of Figure 4,
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Left panel: C/O in the disk's gas and solid phase when incorporating the effect of pebbles drifting to the inner disk, and their evaporation at the ice lines. Right panel: time-dependent 1D posterior of $a_{\text{solid}}$ for the formation inversion of HR 8799e in the pebble drift scenario, Case (i). The ice line positions of NH$_3$ and N$_2$ have been omitted for clarity, and the x-axis has been log-scaled.</figcaption></figure>

indicating a good fit of the atmospheric composition. The right panel of Figure 5 shows the posteriors on the location where the planet accreted its solids (or, alternatively, where these solids originated in the disk) as a function of time. For reference, the ice line positions of the static disk model are indicated as well.

To understand the results of the inversions with chemical evolution, it is useful to reconsider the underlying C/O distribution in the disk as a function of time, shown in Figure 3. For $t > 0$ it is seen that the disk gas C/O value outside the CO$_2$ ice line decays, while the C/O of the solid component increases. This is because CO is converted into CO$_2$ ice on the surfaces of dust grains outside the ice line of CO$_2$ over time. The conversion rate depends on the CO abundance in the ice, which drops rapidly inside the CO ice line. So while the reaction rate increases with temperature, the conversion is more efficient right inside the CO ice line (Bosman et al. 2018). Thus, the process occurs first for larger disk radii and later for smaller radii and drives the solid C/O toward the stellar value also inside the static CO ice line. We note that the ANDES model also included the formation of CH$_3$OH ice.

Because HR 8799e is found to have a C/O ratio similar to the stellar one, this means that the region for its most likely formation (or the region of origin of its accreted solids) expands inward over time to include smaller disk radii. This effect is clearly visible in the right panel of Figure 5. We therefore confirm the findings presented in Mollière et al. (2020), where it was argued that processing CO gas into CO$_2$ ice may have a significant effect on the formation location of HR 8799e. As stated in Mollière et al. (2020), this also has consequences for how strongly HR 8799e may have migrated to reach its present-day orbit. If chemical evolution was significant in HR 8799e's natal protoplanetary disk, and if the exoplanet formed late enough, it may have migrated much less (or not at all) than in cases where it formed early.

In general, our findings emphasize the importance of disk chemical evolution for planet formation that has been reported in Eistrup et al. (2018). It also shows that any analyses that try

to infer planet formation based on atmospheric compositions should compare the relevant chemical timescales to the timescales of planet formation.

Similar to the static disk case, assuming λ Boo-type elemental abundances for the chemical evolution is not expected to change the results significantly. Increasing the carbon and oxygen abundance by 30% is within the modeling uncertainties of the disk chemistry, and the additional oxygen going into H$_2$O due to the silicate depletion is irrelevant to the evolution of the CO ice line. In the inner part of the disk, within the CO$_2$ ice line, water ice is slowly destroyed to form CO$_2$ ice, which raises the solid C/O in the inner disk over time, similar to the CO condensation within the static CO ice line. The amount of available CO$_2$ that can be formed is independent of the H$_2$O ice fraction to first order, and the timescale over which this happens is set by the cosmic-ray ionization rate, so it is independent of the water concentration.

### 3.2.3. Pebble Drift and Evaporation

In this section we model the effect of pebble drift and evaporation on the formation inversion of HR 8799e. This process can be crucial for setting the composition of forming planets. When neglecting pebble drift, planets whose atmospheric metal content is set by gas accretion are generally expected to have substellar metallicities and superstellar C/O values. In contrast, planets with an atmospheric metal enrichment dominated by solid accretion may have superstellar metallicities but substellar C/O ratios (e.g., Öberg et al. 2011; Madhusudhan et al. 2014, 2017; Mordasini et al. 2016). In the case of pebble drift, however, evaporation of pebbles inside of the CO, CO$_2$, and potentially the CH$_4$ ice lines can lead to disk gas that is significantly enriched in these species, allowing for superstellar metallicities and C/O ratios in the disk's gas phase and therefore in the atmospheres of planets (e.g., Booth et al. 2017; Schneider & Bitsch 2021a, 2021b).

For setting up the chemcomp pebble disk model, we used the same initial disk surface density and temperature structure
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Evolution of the CO concentration of the disk gas, normalized by the initial CO concentration, in the pebble drift and evaporation scenario. The earliest times are characterized by a spike in CO close to the CO ice line, due to pebble evaporation, followed by its viscous spreading.</figcaption></figure>

as for the disk’s chemical evolution case. Likewise, the initial disk composition was fixed to the one described in Table 2. For the disk viscosity we chose an intermediate value of $\alpha = 5 \times 10^{-4}$, where $\alpha$ is the usual dimensionless diffusion coefficient, in units of $c_sH$, where $c_s$ is the local midplane sound speed and $H$ the disk’s pressure scale height (Shakura & Sunyaev 1973). This value is consistent with observational data on turbulence in protoplanetary disks, suggesting $\alpha$ of the order of $10^{-3}-10^{-4}$ (Pinte et al. 2016; Flaherty et al. 2017, 2018). The disk viscosity is a key parameter for the pebble problem, with smaller values of $\alpha$ leading to larger pebbles, thus generally faster inward drift, and longer persistence timescales of the gas locally enriched by pebble evaporation. All solid material is considered to be in the form of pebbles, with an initial particle size of $1 \mu$m, which then evolve by growth and drift (e.g., Birnstiel et al. 2012).

The disk’s resulting C/O values in the solid and gas phase are shown in the left panel of Figure 6. At $t=0$ the disk C/O values reproduce our static setup. At larger times, however, the effect of drifting pebbles becomes noticeable very quickly. Pebbles drifting across the CO ice line will start enriching the gas phase in CO (also see Figure 7). Some of this gas diffuses outward again, condensing on the inward-drifting pebbles and increasing the pebble C/O value to unity just outside the CO ice line. The same effect is visible just outside the $CO_2$ and $H_2O$ ice lines, where the solid C/O values reach 0.5 and 0, respectively. Away from the ice lines the C/O of the solids remains largely unchanged, however. At the same time we note that the solid surface density will drop significantly over the simulated time owing to pebble drift, by up to two orders of magnitude, while the gas surface density only drops by less than one order of magnitude. Inside the $CO_2$ ice line the C/O ratio of the gas immediately drops at $t>0$ as $CO_2$ evaporates off the inward-drifting pebbles. At later times the gas’s C/O starts rising again as the CO gas that has evaporated off the pebbles inside the CO ice lines reaches the inner disk regions, due to the disk’s viscous evolution. An analogous evolution

can be observed for the disk gas inside the $H_2O$ ice line; the gas’s C/O value first drops significantly, due to the water evaporating off the pebbles, but rises again at later times as gas enriched in $CO_2$ and CO viscously spreads inward.

For studying the effect of pebble accretion, drift, and evaporation, we investigated four scenarios with our formation inversion setup. Case (i) is simply applying the disk compositional model, as determined by the pebble drift and evaporation framework, in the formation forward model. Case (ii) is like Case (i), but putting an upper limit of $20 M_\oplus$ on the mass that can be accreted as solids, accounting for the concept of the pebble isolation mass (see Section 2.3 and Bitsch et al. 2018b). Case (iii) is like Case (ii), but replacing the upper mass limit on the planetary mass from the spectroscopic retrieval ($M_P < 14 M_{\text{Jup}}$; Mollière et al. 2020) with a tighter prior from the dynamical mass estimate reported in Brandt et al. (2021), that is, $M_P = 9.6_{-1.8}^{+1.9} M_{\text{Jup}}$. Case (iv) is like Case (iii), but increasing the pebble isolation mass to $100 M_\oplus$. The reasoning for testing these different cases will be discussed below, where we summarize the inversion results obtained for the different cases.

The compositional fit for Case (i) is depicted by the leftmost ellipse in the bottom panel of Figure 4. In this scenario pebble drift, evaporation, and accretion are able to reproduce the observed abundance pattern of HR 8799e. Conceptually, Case (i) simply tests whether the results of the static disk model inversion change when introducing pebbles, but it does not yet apply any prior knowledge on how pebbles are accreted, such as the concept of the pebble isolation mass. The reason for the good compositional fit of Case (i) becomes evident when studying the right panel of Figure 6, which shows the resulting posterior for the most likely accretion location of the solids for HR 8799e. Because the solid C/O values do not change significantly, except for just outside the ice lines, the result that significant accretion of solids from outside the CO ice line is likely does not change when compared to the static setup of the disk composition. Just outside the CO ice line the probability goes down, however, because CO gas recondensing on the pebbles drives up the pebbles’ C/O to values larger than the planetary one. What is noticeable is that the region inside the CO ice line at $t>0$ is somewhat more likely when compared to $t=0$. This is because the disk gas, enriched by CO from evaporating pebbles, and with C/O=1, is of high enough metallicity to somewhat offset the C/O value of solids accreted inside the CO ice line, which is too low when compared to the planet. The enrichment of the disk gas in CO over time is shown in Figure 7. We note that the likelihood for accreting a significant amount of pebbles decreases over time because pebbles will drain to the inner parts of the disk. We neglect this effect here.

In Case (i) the upper limit on the planetary mass from the spectroscopic retrieval, together with a superstellar atmospheric metallicity, leads to a $1\sigma$ upper limit of $570 M_\oplus$ on the accreted pebble mass. Such a high value is inconsistent with the concept of the pebble isolation mass. Therefore, we deem Case (ii), where we set an upper limit of $20 M_\oplus$ on the accreted solid mass, a more likely scenario. We note that the pebble isolation mass is a function of the disk viscosity $\alpha$ and that it is very sensitive to the disk aspect ratio ($M_{\text{iso}} \propto [H/r]^3$, with $H$ being the disk scale height). The value of $20 M_\oplus$ is what we derive for the HR 8799 disk model at the location of the CO ice line, using the scaling relations reported in Bitsch et al. (2018b). The compositional fit for Case (ii) is shown in the bottom panel of
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Time-dependent posterior of $a_{\text{gas}}$, the location where gas was accreted by the forming planet HR 8799e. The ice line positions of NH$_3$ and N$_2$ have been omitted for clarity, and the x-axis has been log-scaled. Left panel: pebble drift scenario when including a pebble isolation mass prior of 20 $M_\oplus$ for the accreted solids (Case ii). Right panel: same as the left panel, but for Case (iii), which is additionally including a dynamical mass prior on the mass of HR 8799e.</figcaption></figure>

Figure 4 (second ellipse from the left). Also in this scenario pebble drift is able to reproduce the observed abundance pattern of HR 8799e but leads to a generally somewhat lower planetary metal enrichment. We note that these results assume that all accreted pebbles are visible in the atmosphere, which is equivalent to full core dissolution and mixing. Moreover, in order to allow pebble enrichment to have a noticeable effect on the planet composition, the inversion constrains the planetary mass to <3.5 $M_{\text{Jup}}$. What is more, as a result of the prior limit on the accreted solid mass and the planetary-mass prior, the inversion deems scenarios more likely where the composition of the accreted gas has more impact than in Case (i). The resulting probability distribution on the locations $a_{\text{gas}}$ where the planet accreted its gas is shown in the left panel of Figure 8. The most likely locations and times for the gas accretion correspond to the situation where the gas enriched by the evaporating pebbles reaches approximately C/O values of 0.6, corresponding to the planet's atmosphere (see left panel of Figure 6). In this scenario the most likely formation accretion location is inside of HR 8799e's current orbital position, which would require some outward migration if taken at face value.

In Case (iii) we study the effect of enforcing an upper limit on solid accretion, due to the pebble isolation mass, and a tighter constraint on the planetary mass. The mass prior stems from a dynamical analysis based on the orbital characterization of the HR 8799 system and accelerations from the Gaia-Hipparcos catalog (Brandt et al. 2021). The compositional fit for Case (iii) is shown in the bottom panel of Figure 4 (second ellipse from the right). In this case the inversion struggles to reproduce the observed enrichment pattern of HR 8799e; while it fits the atmospheric C/O ratio well, the planetary enrichment is generally too low, but it improves at later times. This is explained from the fact that the high mass prior assumed for HR 8799e, together with the low pebble isolation mass, does not allow for the pebbles to play a significant role in the planet enrichment, while especially at early times the disk gas is not enriched enough by gas that has evaporated off the inward-drifting pebbles. This situation is thus alleviated at later times,

when the disk gas enrichment increases, but it is never enough to fully reproduce the planetary metal enrichment. The right panel of Figure 8 shows the probability distribution of $a_{\text{gas}}$. Because only gas accretion is able to affect the planetary composition noticeably in Case (iii), it is essentially a higher-contrast version of the $a_{\text{gas}}$ distribution of Case (ii), shown in the left panel.

Case (iv) essentially studies the case when the planet started forming very far outside the CO ice line, in the outer parts of the disk. Due to the disk flaring, $H/r$ increases toward the outer disk, and we would find $M_{\text{iso}} = 50 M_\oplus$, corresponding to $H/r = 0.07$, at 200 au. Because we are interested in an upper limit on what pebble accretion could contribute, we also assume that the disk viscosity is very high for the $M_{\text{iso}}$ calculation ($\alpha = 0.004$, instead of the nominal 0.0005); this results in $M_{\text{iso}} = 100 M_\oplus$. The corresponding enrichment pattern of the planet is shown in the bottom panel of Figure 4, rightmost ellipse. Unsurprisingly, not only the planet C/O but also its metal enrichment are better fit now, when compared to Case (iii), leading to a good fit overall. As expected, $a_{\text{solid}}$ values outside the CO ice line are the most likely for this case, with some additional gas accretion from within the $CO_2$ ice line.

From our investigation it thus becomes evident that pebbles alone, for average $M_{\text{iso}}$ values, may not be sufficient to fully explain the observed abundance pattern of HR 8799e, even when making the assumption that all pebbles accreted onto the planet (likely onto the forming planetary core) mix into the visible atmosphere. This conclusion hinges on at least three assumptions. First, if the planetary mass was actually lower than reported in Brandt et al. (2021), enriching the planet by the accreted solid pebbles becomes easier. This is shown by our Case (ii), where the inversion when imposing $M_{\text{iso}} = 20 M_{\text{Jup}}$ resulted in a good fit by constraining the planetary mass to below 3.5 $M_{\text{Jup}}$. Next, if the pebble isolation mass is much higher than our baseline case (e.g., 100 instead of 20 $M_\oplus$), which is possible for large disk viscosities and the planet initially forming far out in the disk, pebble enrichment becomes
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Like Figure 4, but comparing the nominal pebble inversion Case (iii), that is, $M_{iso} = 20 M_{\odot}$, $M_p = 9.6_{-1.8}^{+1.9} M_{\oplus}$, at solar disk composition (left ellipse), with a setup where the $\lambda$ Boo-type composition of HT 8799 was assumed for the disk instead (right ellipse). An offset was applied to these ellipses for clarity.</figcaption></figure>

a likely scenario for explaining the abundance pattern of HR 8799 again. Lastly, if the composition of the disk is different from our baseline case, even $M_{iso} = 20 M_{\oplus}$ with the dynamical mass prior of HR 8799e (Case iii) becomes a likely scenario again. This is seen in Figure 9, where we show what happens when running Case (iii) again, but assuming the $\lambda$ Boo-type composition of HR 8799 for the disk composition. Due to the carbon and oxygen content being ~30% higher in this case, the accretion of gas that is enriched by evaporated pebbles leads to a better agreement with the total atmospheric metal enrichment. We note that all of these conclusions are based on the spectroscopic retrieval result for HR 8799e, and a slightly lower retrieved metallicity would make pebble accretion more likely again. Due to the large uncertainties on the atmospheric metallicity, even the pebble scenario with the worst fit (Case iii) is only about one standard deviation away from the mean composition derived in the spectroscopic retrievals.

Lastly, it should be kept in mind that other likely important effects connected to pebbles were not studied here. For example, we neglected the effects that the outer HR 8799 planets may have had on the pebble flux that reaches the inner disk and therefore HR 8799e's position. Outer planets may prevent pebbles from drifting inward and evaporating at the CO ice line (e.g., Bitsch et al. 2021; Schneider & Bitsch 2021a). It is unclear to what degree this effect is important here because the giant planets may have formed late enough that some pebble drift may already have taken place in the disk before shutting off the pebble flux. In addition, if HR 8799e, the innermost planet, formed first (high surface densities and orbital periods, and thus shorter accretion timescales), it may have been less affected by the formation of the outer planets.

### 3.3. Suggested Toy Models to Study Other Formation Aspects

Above it was studied how inferences drawn from a simple formation model change if chemical evolution of the protoplanetary disk, or the drift, evaporation and accretion of pebbles

is included. As discussed in Section 2, planet formation is the combination of quite a number of key processes. A concurrent formation inversion with all the ingredients appears both numerically and conceptually unworkable, at the moment. It will still be instructive, however, to add certain aspects of the planet formation problem to such inversion calculations, to study their influence in isolation, or to assess the magnitude of their importance for atmospheric compositions. In Table 1 we list the way in which many of the aspects mentioned in Section 2 may be studied via inversion of the formation process.

To give an example, it would be straightforward to feed disk compositional models that include the disk's self-shadowing into the inversion framework. This process has been suggested by Ohno & Ueda (2021), where the shadowing is caused by a dust pileup at the water ice line. Depending on the grain properties and densities, such a scenario may allow for very volatile species such as CO, $N_2$, and even noble gases such as Ar to condense at distances from the star that are nominally too hot. To study such an effect, various midplane disk and abundance structures, for differing dust density contrasts, could be explored.

Another setup that would be instructive is to further investigate the effect of incomplete mixing between the deep interior and planetary atmosphere for gas giant planets. As discussed in Section 2, the metal enrichment inferred from atmospheric characterization studies is only a lower limit for the true planetary metal enrichment. Where available, an upper limit could be placed based on the analyses of planetary bulk metallicities, as obtained in Thorngren et al. (2016) and Thorngren & Fortney (2019). The impact of metallicity gradients could potentially be studied by adding a parameter $f_{mix}$ that describes whether the metals accreted during formation fully mix ($f_{mix} = 1$) into the atmosphere or not ($f_{mix} = 0$). As long as a planetary atmosphere is of superstellar metallicity, $f_{mix}$ will simply be inversely correlated with the accreted solid mass (if pebble evaporation is neglected). Once a planetary atmosphere is of stellar or substellar metallicity, $f_{mix}$ may also correlate with the formation location of a planet, depending on the disk's abundance structure. An example for this can be constructed by considering our inversion results for HR 8799e in the static disk picture. Because the atmospheric metallicity is high and the planet has a stellar C/O value, the inferred atmospheric C/O ratio could only be reproduced by accreting solids from outside the CO ice line. If the planet's atmospheric metallicity was stellar, it could instead have formed at any location in the disk, as long as the location of gas accretion is equal to the location of solid accretion (measured with respect to the ice lines). If $f_{mix}$ was added as a free parameter, small $f_{mix}$ values would again have yielded regions outside the CO ice line as the most likely region of origin for the solids accreted by the planet.

It is also conceivable to construct a three-component model for the formation inversion that separates the planetary mass into three reservoirs: solids accreted onto the core, solids accreted and mixed into the gaseous envelope, and the gas itself. Each of these would also be associated with a parameter that described where the corresponding material was accreted. One could then define an $f_{mix}$ for the core (or solids in the deep interior) that would describe the degree to which the deep core dissolves and mixes into the envelope.

Lastly, abundance constraints on the refractory content of a planet, as presented in Lothringer et al. (2021) for WASP-121b, for example, may be used to put an upper limit on the mass a planet accreted through pebbles. As long as the planet did not
<table><caption>**Table 1**
Aspects of Planet Formation and Their Potential Treatment in Toy Formation Model Inversions</caption><thead><tr><th>Aspect</th><th>Potential Tractability in Formation Model Inversions</th></tr></thead><tbody><tr><td><strong>Disk Composition and Structure</strong></td><td></td></tr><tr><td>Unknown disk elemental abundances</td><td>Scale using stellar [Fe/H], try varying composition according to scaling uncertainties.</td></tr><tr><td>Available solid reservoir</td><td>Impose limit based on likely disk mass and dust-to-gas ratio.</td></tr><tr><td>Disk (thermal) structure</td><td>Feed in disk structures from dedicated disk models.<br>Explore whether parameterizing 3D effects in 1D model is possible.<br>Changes in disk structure will affect, e.g., ice line positions, as in Ohno & Ueda (2021).</td></tr><tr><td>Planetary back-reaction on disk</td><td>Use simplified gap opening criteria to limit gas accretion, compare to disk lifetimes.<br>Apply pebble isolation mass (limit refractory reservoir accessible to planet).</td></tr><tr><td>Include pebble drift and evaporation at ice lines</td><td>Increase gas metallicity inside of ice lines as a function of time; also see Section 3.2.</td></tr><tr><td><strong>Disk Chemistry</strong></td><td></td></tr><tr><td>Chemical evolution of disk</td><td>Run formation inversion with chemical composition as a function of time.<br>Also see Section 3.2.</td></tr><tr><td>Inherited or “reset” disk abundances</td><td>Explore impact of differing assumptions on disk abundances for the inversion process.</td></tr><tr><td>Cosmic-ray ionization and stellar irradiation</td><td>Use best guesses for retrievals, otherwise explore different values.</td></tr><tr><td>Refractory carbon depletion in inner disk</td><td>Explore via on/off switch.</td></tr><tr><td><strong>Planet Formation</strong></td><td></td></tr><tr><td>Pebble and planetesimal accretion</td><td>Compare inferred solid (refractory) mass of planet with isolation masses.<br>Constrain upper limit on accreted M<sub>pebbles</sub>, lower limit on accreted M<sub>planetesimals</sub>.</td></tr><tr><td>3D planet accretion</td><td>Test impact of parameterizations, for example, vertically averaged abundances for gas accretion.</td></tr><tr><td>Planet migration</td><td>Test inferring multiple formation locations?<br>Add priors on formation location: planet traps?<br>Add priors enforcing inward migration (e.g., a<sub>gas</sub> &le; a<sub>solid</sub>)?</td></tr><tr><td>Leveraging full complexity of formation models</td><td>Explore use of machine-learning techniques.<br>E.g., random forest predictors as demonstrated in Schlecker et al. (2021).</td></tr><tr><td>Planet formation by gravitational instability</td><td>Likely treatable, but requires changes.<br>E.g., steady state viscous disk model &rarr; infall disk model</td></tr><tr><td><strong>Planet Bulk-Atmosphere Coupling</strong></td><td></td></tr><tr><td>Metallicity gradient inside planet</td><td>Use multicomponent model, infer mixing efficiency f<sub>mix</sub> &isin; [0, 1]<br>to reveal correlations with other formation parameters.</td></tr><tr><td><strong>Atmospheric Evolution</strong></td><td></td></tr><tr><td>Evaporation</td><td>Can be important for lower-mass planets or gas planets with a metallicity gradient.<br>Use inverted evaporation models to reveal correlation with formation parameters?</td></tr><tr><td>Infall of comets/asteroids</td><td>Better quantitative modeling needed.<br>Potentially not important for gas giant planets.</td></tr></tbody></table>

form very close around its star, where also refractories may enter the gas phase, refractories can only be incorporated into the planet by solid accretion. If the inferred amount of accreted refractories is higher than expected from the concept of the pebble isolation mass (see Section 2), an upper limit on the amount of accreted pebbles, as well as a lower limit on the amount of accreted planetesimals (or other impactors such as smaller planets; e.g., Ginzburg & Chiang 2020), may be constrainable. Similar constraints may be obtained from the cases where the volatile enrichment of a planet is higher even than what pebble evaporation in a disk may provide: any additional volatile mass must then be accreted in the form of volatile ices.

## 4. Future Observatories and a Census of Atmospheric Compositions

In the previous sections we discussed that inverting atmospheric compositions to reveal the detailed formation history of a planet is hardly at the moment: the process of planet formation is too complex, with too many unknowns, and likely

too numerically costly to invert. However, the JWST, the class of future ground-based Extremely Large Telescopes (ELTs), and later ARIEL will record high-quality spectra for hundreds of planets. In this section we summarize the compositional constraints that can be extracted from such atmospheric measurements and how the resulting atmospheric enrichment patterns for the planetary population may allow us to constrain planet formation in a broader sense.

### 4.1. C/O

The importance of the planetary C/O ratio for informing planet formation has been discussed in Sections 1 and 2. In these sections we also discuss the complications that may make the picture likely more complex than suggested by the foundational study by Öberg et al. (2011). For example, while dominant solid enrichment is generally expected to lead to stellar or substellar C/O ratios and superstellar planetary metallicities, pebble drift and evaporation may lead to super-stellar enrichments of the gas at C/O values both smaller or larger than stellar. Interestingly, however, superstellar C/O
Atmospheric occurrence of formation tracers in H – He – dominated atmospheres

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

**Figure 10.** Potential atmospheric visibility of various absorbers in planetary atmospheres. Every species or group of species shown here is known to be spectrally active. We searched the literature for the average atmospheric temperatures where these species are visible. Alternatively, we used the equilibrium chemistry code described in Mollière et al. (2017) and checked for which temperature range the species is present in the atmosphere. Our standard assumption was solar metallicity and abundance ratios and a pressure of 0.1 bar, whereas dissociation and ionization values were obtained from assuming pressures from 0.1 to 0.001 bars. We assumed either solar C/O (=0.55) or C/O = 1.1. The temperatures given therefore should only serve as rough guidance and do not necessarily correspond to a planet’s effective temperature. We note that chemical transitions also depend on the metallicity and the pressure at the planetary photosphere (therefore effectively also on the planetary gravity *g*). Moreover, many of these species can be affected by disequilibrium chemistry (see, e.g., Fortney et al. 2020) or be cold-trapped into condensates (e.g., Spiegel et al. 2009; Parmentier et al. 2016). The chemical behavior of the species listed here is described in Section 4 and Appendix D for the refractories.

values and enrichments are difficult to obtain without considering pebble evaporation (also see Section 3.2), so a large population of planets with such abundance characteristics may indicate a dominant role of pebbles for setting planetary abundances. Similarly, a large enough overall metal enrichment of a planet, especially if formed in the outer disk, may be difficult to explain from pebble evaporation, even for low disk viscosities, which would point more toward planetesimal accretion playing an important role.

In general, C/O is also popular because it determines the relative abundances of the spectrally active C- and O-bearing molecules in the exoplanet atmospheres such as H₂O, CH₄, CO, CO₂, HCN, and C₂H₂. C/O therefore regulates the spectral appearance of a planet in the near- to mid-infrared (e.g., Fortney et al. 2005; Seager et al. 2005; Madhusudhan 2012; Moses et al. 2013; Mollière et al. 2015; Molaverdikhani et al. 2019a; Goyal et al. 2020; Hobbs et al. 2021b).

For reference, Figure 10 shows under which atmospheric conditions the absorbing species that trace the C/O ratio in gas-dominated planets may be visible. Also see, for example, Lodders & Fegley (2002) for a detailed description of the atmospheric chemistry. For temperatures below about 1000 K the atmosphere will be rich in H₂O and CH₄; for higher temperatures these species will be converted into CO until either C or O runs out, depending on the C/O ratio. For high temperatures and C/O ≥ 1, CH₄ will thus be visible; for high

temperatures and C/O ≲ 1, H₂O will be visible. For further increasing temperatures and C/O ≥ 1, CH₄ is replaced by increasing amounts of C₂H₂ and HCN (see, e.g., Madhusudhan 2012; Mollière et al. 2015). We note that the chemical transitions mentioned here also depend on the local atmospheric pressure (Mollière et al. 2015; Molaverdikhani et al. 2019a). CO can still be visible in cool atmospheres, especially of self-luminous brown dwarfs and planets, because atmospheric mixing may transport CO-rich gas from the deep (hotter) atmosphere to the photosphere (e.g., Zahnle & Marley 2014; Miles et al. 2020). Whether such disequilibrium abundances are expected for irradiated (often transiting) planets is less clear because the insolation leads to more isothermal atmospheres. For planets that are still strongly cooling (with a high internal temperature) or heated by processes such as eccentricity dampening, CH₄ may be strongly suppressed (Fortney et al. 2020). As mentioned, ground-based high-contrast or high-resolution observations have started to obtain the first useful constraints on C/O. The state of the art will greatly improve once JWST and later ARIEL allow for a larger census of planetary compositions (also see our discussion in Section 1).

## 4.2. N/O, N/C

The importance of atmospheric nitrogen-bearing species such as NH₃ or HCN for constraining exoplanet formation has
been recognized recently, especially for planets forming in the outer solar and extrasolar disks. The reason for this is that nitrogen, predominantly in the form of N₂ in protoplanetary disks, is extremely volatile. Planets forming at increasingly larger distances, when dominated by solid metal enrichment, will therefore exhibit increasingly lower N/O or N/C ratios, and vice versa if dominated by gas metal enrichment. This is because several ice lines of C- and O-bearing species are crossed toward larger orbital radii, while N₂ stays in the gas phase (Turrini et al. 2021). If the planet forms at wide enough orbital distances, eventually N₂ will freeze out as well, leading to an enhanced atmospheric nitrogen content, which will scale similarly with atmospheric metallicity as the abundances of C- and O-bearing species. The high nitrogen content of Jupiter has therefore led to the interpretation that Jupiter formed in the outer regions of the solar system, beyond the location of the N₂ ice line at ~30 au, which is also consistent with the planet's elevated abundance of noble gases (e.g., Owen & Encrenaz 2003; Bitsch et al. 2015; Bosman et al. 2019; Öberg & Wordsworth 2019; Cridland et al. 2020b). Similar to the discussion of C/O, the situation is likely more complicated also for the nitrogen enrichment. Both disk self-shadowing (see Ohno & Ueda 2021, and our discussion in Section 3.3) and pebble drift and evaporation (Schneider & Bitsch 2021b) are likely complicating factors. We also note that a planet that forms late within a protoplanetary disk's lifetime may be less sensitive to N₂, as cosmic-ray ionization may process N₂ to NH₃ ice over Myr timescales, such that the importance of NH₃ and its ice line increases over time, with the ice line of NH₃ being much closer to the star than that of N₂ (Semenov & Wiebe 2011).

In exoplanets the only spectrally active nitrogen-bearing species of relevance are NH₃ and HCN. N₂, which is the dominating nitrogen bearer at larger temperatures, has negligible opacity in the near- and mid-infrared. NH₃, on the other hand, should be detectable in the mid-IR using JWST in exoplanet atmospheres (e.g., Daniëlski et al. 2018). Moreover, evidence for NH₃ has been seen in high-resolution studies (Giacobbe et al. 2021; Sánchez-López et al. 2022). For C/O ~1, NH₃ is only abundant up to ~500 K (Lodders & Fegley 2002). HCN, on the other hand, will be visible for temperatures of 1500 K or larger, if C/O > 1 (e.g., Madhusudhan 2012; Mollière et al. 2015). We indicate these detectability ranges in Figure 10. Chemical disequilibrium may (or may not) allow for NH₃ or HCN to be visible at intermediate temperatures (500 K < T < 1500 K) in irradiated planets as well (see, e.g., MacDonald & Madhusudhan 2017; Fortney et al. 2020; Hobbs et al. 2021a, and references therein). For self-luminous planets disequilibrium chemistry may play less of a role for N, as isoabundance lines are parallel to atmospheric pressure-temperature profiles (Zahnle & Marley 2014). As before, all chemical transition temperatures also depend on the atmospheric pressure.

We note that the chemical behavior of N-, C-, and O-bearing species described here mostly hinges on chemical equilibrium or simple atmospheric disequilibrium treatments, also considering the planetary atmospheres to be one-dimensional and of mostly scaled solar abundances (except for the C/O ratio). The recent and intriguing results of Giacobbe et al. (2021), who detected H₂O, CO, HCN, C₂H₂, NH₃, and CH₄ in the atmosphere of HD 209458b (with an equilibrium temperature of ~1500 K), are a reminder that atmospheric chemistry may be much more complex than discussed above. An important

effect is the horizontal advection of chemical abundances predicted from coupling chemical models to the output of 3D general circulation models (e.g., Agúndez et al. 2014; Baeyens et al. 2021), which could also be connected to condensate rainout (Sánchez-López et al. 2022). Photochemistry is also important, especially in the upper atmospheric layers (e.g., Kopparapu et al. 2012; Venot et al. 2012; Molaverdikhani et al. 2019b). Telescopes such as JWST, current high-resolution spectrographs, and ultimately instruments mounted on ELT-class telescopes will allow us to investigate these effects more thoroughly.

### 4.3. R/O

Measuring the refractory content of an atmosphere could provide unique insight into a planet's formation history. As has been argued recently by Lothringer et al. (2021), measuring the refractory-to-oxygen ratio R/O of a planet constrains the importance of metal enrichment by rocky accretion relative to icy or gaseous accretion. Here R stands for any element that traces the refractory content of the planet (Fe, Na, K, Si, Mg, Ti, ...), or an average of such elements. As argued further, this may allow the placement of constraints on whether the planet (or its solid building blocks) migrated significantly during formation. We argue that R/O may potentially even be useful to constrain the relative importance of pebble and planetesimal accretion, in the core accretion paradigm; see also Section 3.3 and the discussion in Schneider & Bitsch (2021b).

In Figure 10 we indicate the temperature ranges over which various refractory-tracing atmospheric absorbers are visible. We refer the reader to Appendix D for a discussion of the chemistry of the refractory-tracing absorbers. Lothringer et al. (2021) put emphasis on ultrahot Jupiters, for which various refractory elements exist as molecules (metal oxides or hydrides), atoms, or ions in the gas phase. We also note that species such as H₂S and PH₃ may be useful refractory tracers at intermediate atmospheric temperatures (Wang et al. 2017; Öberg & Wordsworth 2019). While H₂S and PH₃ are volatile species, the dominant carrier of P and S atoms in a protoplanetary disk appears to be refractory species (Öberg & Wordsworth 2019). Moreover, measuring the abundances of Na and K in planetary atmospheres may be worthwhile to trace of the refractory content (Welbanks et al. 2019).

Refractory cloud species may affect planetary spectra by muting molecular features and reddening the spectral energy distribution. Silicate particles like MgSiO₃ and Mg₂SiO₄ are especially interesting, as they may lead to visible absorption features around 10 µm (e.g., Cushing et al. 2006; Wakeford & Sing 2015). Due to the complex microphysical problem of cloud formation (e.g., Rossow 1978; Powell et al. 2018; Woitke et al. 2020), measuring a refractory abundance from observed cloud absorption may prove difficult, however. Moreover, clouds may complicate measuring and interpreting the abundances of gas refractory species due to cold trapping by condensate rainout (e.g., Spiegel et al. 2009; Parmentier et al. 2016). An interesting alternative to silicate clouds could be searching for the absorption of gaseous SiO at ~7 µm with JWST's MIRI instrument. SiO is promising because it is the most abundant Si-bearing gas species after the silicates evaporate (Visscher et al. 2010) and is more stable than H₂O against dissociation (by about 500 K). This should allow detection of this species in ultrahot Jupiters.
Other metal oxides such as TiO, VO, AlO, and CaO have features in the optical and near-infrared (e.g., Sharp & Burrows 2007; Gandhi & Madhusudhan 2019; Lothringer et al. 2020). Similarly, metal hydrides may be useful refractory tracers, at similar temperatures to the metal oxides. Species such as FeH, CaH, MgH, NaH, CrH, and TiH all have absorption features in the optical and near-infrared (Sharp & Burrows 2007; Gandhi & Madhusudhan 2019). Metal atoms are visible in the atmosphere once the refractory clouds are no longer present (e.g., Mg, Fe), or once the dominant molecular species (such as SiO for Si) have been dissociated. Mg, Fe, Ca, Cr, Ni, V, Na, and maybe Co have been detected in the ultrahot atmospheres of KELT-9b and WASP-121b in the optical (Hoeijmakers et al. 2019, 2020). Finally, metal ions become visible in the hottest atmospheres as soon as the atoms have been ionized. This has led to the detection of Fe+, Ti+, Cr+, Sc+, Y+, and maybe Sr+ in the hottest known exoplanet KELT-9b in the optical (Hoeijmakers et al. 2019).

## 5. Discussion and Summary

Inferring the formation history of a planet, based on its atmospheric composition, is one of the most cited goals of the atmospheric characterization community. In our work we take a look at what obstacles need to be overcome to make such an inversion feasible.

Summarizing the complex and interconnected processes that govern planet formation (see Section 2), we conclude that actually inverting planet formation in this way is still a long way off, if even possible at all. Current formation models are likely too complex (too many free parameters), too uncertain (which processes to consider, which assumptions to make for them), and too numerically costly (N-body interactions, dust evolution, hydrodynamical evolution, disk chemical evolution, etc.). Many of these problems may actually be alleviated in the coming years or decades, but the degree to which such a full formation inversion will ever become possible is difficult to assess, at the moment. As an interesting avenue for inverting full, state-of-the art formation models, we want to highlight the recent work by Schlecker et al. (2021), where a random forest technique was used to predict planetary formation outcomes based on formation model input parameters. We will have to wait to see how far this method can be used to predict planetary abundances.

Apart from this conclusion, we also introduce a method that allows us to study and compare the qualitative impact of different assumptions made in the modeling of planet formation; see Section 3. Assuming some measured planetary compositions as observations, we use nested sampling to invert simplified formation models, constraining their corresponding formation parameters. Due to the challenges mentioned above, such invertible formation models cannot be complex enough to yield reliable results on a given planet's formation process. However, they may allow us to study the importance of various formation aspects in isolation. As an example, we show how the deduced formation history of the directly imaged planet HR 8799e changes if the composition of the protoplanetary disk in which it forms is allowed to evolve chemically. We find that chemical evolution may significantly affect the migration history inferred for this planet; the planet may have migrated much less if chemical evolution is taken into account. What is more, we show that the drift, evaporation, and accretion of pebbles are able to reproduce the planetary C/O value, but whether it can reproduce the inferred high atmospheric

metallicity depends on the assumptions made for the disk viscosity, pebble isolation mass, and disk composition. We end this section by suggesting a number of other formation processes that could be studied in a similar way, for example, metallicity gradients and ineffective mixing of the planetary interior.

While the detailed inversion of planet formation may still be in the far future, it is clear that the atmospheric abundance constraints obtained with new and upcoming instruments will be crucial to inform planet formation models in a broader sense. In Section 4 we summarize under which atmospheric conditions various spectrally active atmospheric species that trace the C/O value (H₂O, CO, CH₄, CO₂, C₂H₂, HCN), nitrogen content (NH₃ and HCN), and refractory content (H₂S, PH₃, alkalis, refractory clouds, metal oxides, hydrides, atoms and ions) may be observable in H/He-dominated atmospheres. Instruments such as GRAVITY, CRIRES+ (or other high-resolution spectrographs), JWST, and facilities further in the future like ARIEL and the ELTs will obtain abundance constraints for many of these species. We discuss how the C/O values derived for the atmospheric composition of exoplanets may allow us to constrain the importance of pebble drift and evaporation, as well as how the refractory content of a planet may constrain the relative contribution of planetesimal and pebble accretion.

Making the connection between atmospheric abundances and formation a reality seems daunting, but the likely transformative nature of observations of many upcoming observational facilities will lead to more precise atmospheric abundance constraints for exoplanets than ever before. The constraints obtained from these observations will require being put into context, to assess what information on planet formation may possibly be gleaned from them. With these data one may begin assessing the degree to which planet formation can indeed be informed by the atmospheric composition of exoplanets.

We would like to thank the anonymous referee for their detailed report, which greatly improved the quality of this manuscript. P.M. and Th.H. acknowledge support from the European Research Council under the European Union's Horizon 2020 research and innovation program under grant agreement No. 832428-Origins. T.M. acknowledges support of the Ministry of Science and Higher Education of the Russian Federation under grant 075-15-2020-780 (N13.1902.21.0039; Sections 2.2 and 3.1.2). B.B. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. A.S. acknowledges funding from the European Union H2020-MSCA-ITN-2019 under grant agreement No. 860470 (CHAMELEON). R.B. acknowledges the support from the Swiss National Science Foundation under grant P2BEP2_195285. C.M. acknowledges the support from the Swiss National Science Foundation under grant BSSG10_155816 "PlanetsInTime." D.S. acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1833: "Building a Habitable Earth" (SE 1962/6-1). Parts of this work have been carried out within the frame of the National Center for Competence in Research PlanetS supported by the SNSF.

## Appendix A
### Atmospheric Evolution

Here we derive Equation (2), which estimates the change in mass fraction of a given species due to the enrichment of a
planet's atmosphere by impacts; see also Section 2.5. We also discuss the case of pure-water comets increasing the atmospheric water content of a Jovian planet.

To begin, we estimate the mixing in the atmosphere by 1D diffusion. We thus write, using the usual 1D diffusion equation for concentrations (e.g., Parmentier et al. 2013),

$$ \rho \frac{\partial X}{\partial t} = \frac{\partial}{\partial z} \left( K_{zz} \rho \frac{\partial X}{\partial z} \right) + W(z, z_i) \frac{\dot{M}}{4\pi R_p^2 \Delta z_i}, \quad (A1) $$

where $\rho$ is the atmospheric density, $X$ the atmospheric water mass fraction, $t$ the time, $z$ the atmospheric altitude, $\dot{M}$ the mass accretion rate of pure-water comets, and $R_p$ the planetary radius. $W(z, z_i)$ is defined as

$$ W(z, z_i) = \Theta(z - z_i + \Delta z_i / 2) - \Theta(z - z_i - \Delta z_i / 2), \quad (A2) $$

with $\Theta$ being the Heaviside step function. This means that we assume that the comets are destroyed in a narrow layer of width $\Delta z_i$, at altitude $z_i$ in the atmosphere. Using the equation of hydrostatic equilibrium ($\partial P / \partial z = -\rho g$), together with the equation of state of an ideal gas ($P = \rho k_B T / \mu$) and that $H_p = k_B T / \mu g$, one finds that one can express Equation (A1) as

$$ \frac{\partial X}{\partial t} = \frac{\partial}{\partial P} \left( \frac{K_{zz}}{H_p^2} P^2 \frac{\partial X}{\partial P} \right) + W(P, P_i) \frac{\dot{M} g}{4\pi R_p^2 \Delta P_i}, \quad (A3) $$

where $P$ is the atmospheric pressure, $g$ the gravity, $k_B$ the Boltzmann constant, $T$ the atmospheric temperature, and $\mu$ the atmospheric mean molecular weight. For a steady-state ansatz and integration over $P$ one obtains that

$$ \frac{\partial X}{\partial P} = - \frac{\dot{M} g H_p^2}{4\pi R_p^2 K_{zz} P^2} \quad (A4) $$

for $P \gg P_i + \Delta P_i / 2$ and $\partial X / \partial P = 0$ for $P \ll P_i - \Delta P_i / 2$ and a linear transition between these two cases for $P \in (P_i - \Delta P_i / 2, P_i + \Delta P_i / 2)$. In the following we assume that $\Delta P_i \ll P_i$. Thus, it holds that the enrichment in units of mass fractions $\Delta X = X - X_0$ at $P_i$ and at lower pressures (higher altitudes) is

$$ \Delta X = \frac{\dot{M} g H_p^2 P_{RCB} - P_i}{4\pi R_p^2 K_{zz} P_{RCB} P_i}, \quad (A5) $$

where $K_{zz}$ was assumed to be constant for simplicity. $P_{RCB}$ denotes the location of the RCB or, more generally, the altitude in the planet below which the planet is well mixed, with $X(P > P_{RCB}) = X_0$ and $P_{RCB} > P_i$.

From Equation (2) we see that $\Delta X$ will become large for small $P_i$ and large $P_{RCB}$. Pinhas et al. (2016) found that water ice comets of 1 km in size will be destroyed by 100 bars if impacting Jupiter at terminal velocity. From our condition that $P_{RCB} > P_i$ this means that we need to consider $P_{RCB} > 100$ bars. For warm Jupiters the maximum depth of the RCB may be as deep as $P_{RCB} = 200$ bars (Thorngren et al. 2019; Sarkis et al. 2021). Then, assuming $X_0 = 10^{-3}$, $M_P = 1 M_{Jup}$, $R_P = 1 R_{Jup}$, $H_P = 200$ km, $K_{zz} = 10^8 \text{ cm}^2 \text{ s}^{-1}$, $P_{RCB} = 200$ bar, and $P_i = 100$ bar results in a relative enrichment of $\Delta X / X_0 = 3 \times 10^{-4}$, if a very high impact rate of $10^5$ comets of 1 km size per year are assumed. This would correspond to 275 impacts per day (and would double the total water content of the planet in $5 \times 10^6$ yr). Therefore, combining reasonable estimates for the parameters describing the atmosphere with a

<figure><img src="image_15.png" /><figcaption>Figure 11. Verification of Equation (A5) by comparing to the numerical integration of Equation (A3) to very long times $t$, to obtain the limit $t \to \infty$. To obtain nonzero $\Delta X / X_0$ values, an unrealistically low value of $P_i = 1$ bar was chosen, and 20 cometary impacts per year with comets of 20 km in size were assumed.</figcaption></figure>

very high cometary impact rate would not really change the water content of the planet's atmosphere (neglecting the change in $X_0$ over $5 \times 10^6$ yr). An obvious caveat of our toy model that comes to mind is the assumption that the planet will mix any pollution away instantaneously at pressures larger than the RCB. The deep convective $K_{zz}$, estimated from mixing-length theory, may be in the range of $K_{zz} = 10^9 \text{ cm}^2 \text{ s}^{-1}$ (see, e.g., Equation (4) of Zahnle & Marley 2014). Extending the integration domain to $2 \times 10^4$ bars and setting $K_{zz} = 10^9 \text{ cm}^2 \text{ s}^{-1}$ for $P > P_{RCB}$ leads to $\Delta X / X_0 = 3 \times 10^{-2}$ when numerically integrating Equation (A4), that is, 100 times higher, but still too low.

Finally, to demonstrate the good agreement of our analytical Equation (A5), we show a comparison to Equation (A3), numerically integrated to very long times $t$, to obtain the limiting case $t \to \infty$, in Figure 11. For this comparison 20 comets per year of 20 km size were assumed, which are destroyed at an unrealistically low $P_i = 1$ bar (to obtain nonnegligible $\Delta X / X_0$ values).

## Appendix B: Toy Formation Models

### B.1. Öberg et al. (2011) Model

Here we describe the disk setup used for inverting the Öberg et al. (2011) model (see Section 3.1.1), which assumes a static disk composition. More specifically, our disk model is based on the description in Öberg & Wordsworth (2019), which assumes a static power law for the disk temperature and density of the young solar nebula. Inside its ice line position a given volatile species is in the gas phase, whereas outside it is in the solid phase. For every species the mass fraction compared to the total disk mass is tabulated, in addition to the mass fractions of the constituent atoms within a volatile species. We also account for refractory material, which we include in our framework by
<table><caption>**Table 2**
Composition of the Protoplanetary Disk Model and Assumed Ice line Positions, Adapting the Abundances Given for the Young Solar Nebula in Öberg & Wordsworth (2019) to the Case of HR 8799</caption><thead><tr><th>Disk Species</th><th>Constituent Atoms</th><th>Relative Atom Mass Contribution to Disk Species</th><th>Disk Mass Fraction</th><th>Ice Line (au)</th></tr></thead><tbody><tr><td>Refractories</td><td></td><td></td><td>6.2 × 10<sup>-3</sup></td><td>0</td></tr><tr><td></td><td>Fe</td><td>0.21</td><td></td><td></td></tr><tr><td></td><td>Mg</td><td>0.15</td><td></td><td></td></tr><tr><td></td><td>Si</td><td>0.17</td><td></td><td></td></tr><tr><td></td><td>O</td><td>0.29</td><td></td><td></td></tr><tr><td></td><td>C</td><td>0.13</td><td></td><td></td></tr><tr><td></td><td>P</td><td>8.9 × 10<sup>-4</sup></td><td></td><td></td></tr><tr><td></td><td>S</td><td>0.05</td><td></td><td></td></tr><tr><td>H<sub>2</sub>O</td><td></td><td></td><td>1.6 × 10<sup>-3</sup></td><td>2.3</td></tr><tr><td></td><td>H</td><td>2/18</td><td></td><td></td></tr><tr><td></td><td>O</td><td>16/18</td><td></td><td></td></tr><tr><td>NH<sub>3</sub></td><td></td><td></td><td>8.5 × 10<sup>-5</sup></td><td>6.5</td></tr><tr><td></td><td>N</td><td>14/17</td><td></td><td></td></tr><tr><td></td><td>H</td><td>3/17</td><td></td><td></td></tr><tr><td>CO<sub>2</sub></td><td></td><td></td><td>1.3 × 10<sup>-3</sup></td><td>8.0</td></tr><tr><td></td><td>C</td><td>12/44</td><td></td><td></td></tr><tr><td></td><td>O</td><td>32/44</td><td></td><td></td></tr><tr><td>CO</td><td></td><td></td><td>2.3 × 10<sup>-3</sup></td><td>35.0</td></tr><tr><td></td><td>C</td><td>12/28</td><td></td><td></td></tr><tr><td></td><td>O</td><td>16/28</td><td></td><td></td></tr><tr><td>N<sub>2</sub></td><td></td><td></td><td>6 × 10<sup>-4</sup></td><td>73.0</td></tr><tr><td></td><td>N</td><td>1</td><td></td><td></td></tr><tr><td>Background</td><td></td><td></td><td>0.987</td><td>1000</td></tr><tr><td></td><td>H</td><td>0.75</td><td></td><td></td></tr><tr><td></td><td>He</td><td>0.25</td><td></td><td></td></tr></tbody></table>

**Note.** The ice lines of the refractory and background species were set to arbitrarily small/large values to ensure that they stay condensed/gaseous within the simulation domain.

setting its ice line position to zero. The background gases H₂ and He are included by setting their ice line to 1000 au. This ensures that the refractory and background species stay condensed/gaseous within the simulation domain. Table 2 lists the mass fractions, ice line positions, and atomic composition of all considered disk species and their constituent atoms. We note that the ice line positions given here are those expected for the disk around HR 8799, which have been obtained from our ANDES disk model; see Section 3.2.

The mass fractions were obtained using the provisional protosolar nebula composition from Öberg & Wordsworth (2019) with some modifications; see below. The abundance of a species *i* is given as number fractions *x*<sub>*i*</sub> in Öberg & Wordsworth (2019), compared to the number of hydrogen atoms *n*<sub>H</sub>. This was converted into mass fractions *m*<sub>*i*</sub> relative to the total disk mass using

$$m_i = \frac{\mu_i x_i}{1.4}, \qquad (B1)$$

with $\mu_i$ being the molecular mass of species *i* in atomic mass units. This expression is obtained from the disk's total mass

$M_{\text{tot}}$, from the mass of species $M_i$ of species *i*, and by setting

$$m_i = \frac{M_i}{M_{\text{tot}}} \qquad (B2)$$

$$\approx \frac{\mu_i x_i n_H}{\mu_H n_H + \mu_{\text{He}} n_{\text{He}}} \qquad (B3)$$

$$= \frac{\mu_i x_i}{\mu_H + \mu_{\text{He}} n_{\text{He}} / n_H}, \qquad (B4)$$

where it was assumed that most of the disk mass is contributed by H and He atoms. Setting $\mu_H = 1$ and $\mu_{\text{He}} = 4$ and assuming that $n_{\text{He}}/n_H = 0.1$ (see Table 8 in Lodders 2019, for the recommended protosolar abundances) leads to the relation given in Equation (B1).

The refractory composition model was likewise constructed using the information given in Öberg & Wordsworth (2019). In their model, this results in 30% of all oxygen in the form of refractory silicates, identical to the amount of oxygen in H₂O. We assumed that the silicates consist purely of MgSiO₃, which also conserves the solar Mg/Si abundance ratio, which is close to unity (Asplund et al. 2009; Lodders 2019). The refractory carbon component plus volatile organics account for 50% of all carbon atoms, with a 3:1 ratio between the two. To simplify, we added the carbon of the volatile organics by increasing the CO abundance, taking the required oxygen from the water mass reservoir. Because the fate of organic carbon, especially in the inner part of the disk, is uncertain anyway (see, e.g., Mordasini et al. 2016; Cridland et al. 2019, and references therein), we decided to forgo a more careful treatment of the organic carbon reservoir for our conceptual study here. Iron, sulfur, and phosphorous atoms were assumed to only be present in the refractory phase. The C, Fe, S, and P abundances, relative to H, were taken from Asplund et al. (2009). Our resulting C/O ratio distribution for the disk solid and gas components is shown in Figure 1.

## B.2. Disk Chemical Evolution

Here we describe the ANDES chemistry model used for inverting the formation model including disk chemical evolution (see Section 3.1.2), which includes the evolution of the disk's chemical composition. In ANDES, surface reactions are described by the Langmuir-Hinshelwood mechanism and are not limited to hydrogenation. H and H₂ tunneling through reaction barriers is also included. Any dynamical effects on the distribution of C and O, such as drifting grains, are omitted. The surface density profile is described by a power law with the exponent equal to 1.5. The chemical network is based on ALCHEMIC (Semenov et al. 2010; Semenov & Wiebe 2011), with updated binding energies from Cuppen et al. (2017). It incorporates the effects of extreme ultraviolet irradiation, cosmic rays, and radionuclides as ionization sources. The dust size distribution is described by a power law with $p = -3.5$ between 0.005 and 25 µm, which reflects dust growth in disks compared to the interstellar medium. It is used to calculate the radiation field and dust temperature in the disk's upper layers. An average grain radius of 0.35 µm is used for calculating surface reaction rates. The disk abundances are initialized assuming that all volatile species are in the gas phase, using the same abundances as used for the static disk model; see Table 2.
## Appendix C
### Full Inversion Posterior of HR 8799e

In the left panel of Figure 12 we show the 1D and 2D projections of the full posterior resulting from the formation inversion of HR 8799e, as discussed in Section 3.2, using the static disk composition. The posterior of the planetary mass closely follows the mass prior, which we sampled by using the spectral retrieval results for HR 8799e, namely, $log(g) = 4.0 \pm 0.5$, $R_P = 1.12 \pm 0.09 R_{Jup}$, as reported in Mollière et al. (2020). We neglected the error on $R_P$, assumed that $log(g)$ follows a Gaussian distribution, and converted to mass via $10^{log(g)} R_P^2/G$, where G is the gravitational constant. A flat prior was assumed on the formation/accretion locations. The prior on the accreted planetesimal mass was taken to be flat, ranging from 0 to 1000 $M_\oplus$. The posterior of the accreted solid mass can be explained considering the atmospheric

metallicity that was used as an input to the formation inversion ($[Fe/H] = 0.48 \pm 0.25$) and the total mass of the planet. The increased probability of the planet having accreted solids from outside the CO ice line or inside the H$_2$O ice line (discussed in Section 3.2) is visible. There also exists a less likely solution with lower total metallicity (lower solid mass) where both the solids and the gas were accreted between the CO$_2$ and CO ice lines. This branch of solutions can be explained by studying Figure 1, showing the variation in C/O in the disk gas and solids as a function of distance: between these two ice lines the solids' C/O is substellar and the planetary C/O can be raised to higher values by accreting gas that has an increased C/O ratio. In our current model setup this only works if the planet has a low metallicity; otherwise, the gas enrichment cannot compete with the metal enrichment from the solids. In the right panel we show the corresponding posterior in the case when $\lambda$ Boo-type abundances are assumed for the disk of HR 8799. The solution

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Full posterior of the formation inversion of HR8799e assuming the static disk model à la Öberg et al. (2011). Left panel: nominal disk composition (solar). Right panel: same as the left panel, but assuming a $\lambda$ Boo-type composition for HR 8799.</figcaption></figure>
inside the H₂O ice line is no longer valid for *a*<sub>solid</sub>, due to the high local solid C/O ratio.

## Appendix D
### Refractory Chemistry

Here we give a short description of the chemical behavior of atmospheric species that trace the planetary refractory content, as shown in Figure 10. We outline their behavior as a function of temperature, assuming a pressure of 0.1 bars, whereas dissociation and ionization values were obtained from assuming pressures of 0.1–0.001 bars. We assumed either solar C/O (>= 0.55) or C/O = 1.1. We only roughly determine the transition temperatures, as these may also depend on the atmospheric gravity and metallicity. Moreover, disequilibrium chemistry, internal luminosity, insolation flux, and cold trapping can play important roles (see, e.g., Spiegel et al. 2009; Parmentier et al. 2016; Fortney et al. 2020). The temperatures given here therefore do not necessarily directly translate into planetary effective temperatures. If no reference is given, we use the equilibrium chemistry code described in Mollière et al. (2017) to determine the chemical behavior.

#### D.1. H₂S

H₂S condenses into NH₄SH at ~200 K; the higher-temperature condensates MnS, ZnS, and Na₂S are of minor importance (Lodders 2010). For C/O < 1, H₂S dissociates at ~2000 K, while it moves into species like CS for C/O > 1 at temperatures around 1500–2000 K.

#### D.2. PH₃

PH₃ condenses into H₃PO₄ at 500 K. Its presence in Jupiter's atmosphere at lower temperature indicates a deep quenching point, however, such that PH₃ may still be visible at lower temperatures (see, e.g., Baudino et al. 2017). For temperatures approaching 1000 K, PH₃ is increasingly converted into PH₂.

#### D.3. Na

Na condenses into Na₂S at ~900 K. In principle, alkalis such as Na could also be sequestered into high-temperature condensates such as feldspars. However, this likely does not occur owing to the rainout of silicates (e.g., Line et al. 2017) that deplete Si from the atmosphere, which is needed for feldspar formation. Above 900 K Na is thus in the gas phase, until it gets ionized at around 2500 K.

#### D.4. K

K condenses into KCl at ~900 K. In analogy to Na, sequestration of K into feldspars likely does not occur owing to silicate rainout. Thus, K is in the gas phase from ~900 to ~2000 K, after which it is ionized. Ionization occurs at temperatures roughly 500 K cooler than for Na.

#### D.5. Refractory Clouds

As mentioned above, the refractory cloud species Na₂S and KCl likely form at temperatures below 900 K. Here we focus on the remaining cloud species forming at intermediate to hot temperatures and concentrate on those carrying the largest mass and/or opacity, often using the data given in Wakeford et al. (2017), or our own equilibrium chemistry calculations.

Refractory clouds can exist if the atmospheric temperature is below their respective evaporation temperature. Like all chemical transitions discussed here, this temperature depends on the elemental abundance and local atmospheric pressure. Under our adopted standard conditions silicates such as MgSiO₃ and Mg₂SiO₄ evaporate at temperatures around ~1600 K, while iron clouds are stable until ~1700 K. VO and calcium titanates are stable until 1600–1800 K, respectively. Aluminum-bearing condensates such as Al₂O₃, which are among the most stable ones, evaporate around 1900 K. Among the species listed here, potentially only silicates, Al₂O₃, and KCl may actually form in the visible part of the atmosphere, as these species have low surface energies, leading to high nucleation rates (Gao et al. 2020). The cloud bases will reside deeper inside the planetary atmosphere for lower temperatures, with the cloud particles entering from above, or settling below the photosphere. For brown dwarfs this temperature-dependent removal of silicate clouds is thought to cause the L–T transition, which typically occurs at T<sub>eff</sub> = 1200–1400 K (e.g., Best et al. 2021, and references therein), while for planets and low-gravity brown dwarfs this limiting temperature may be as low as approximately 1000 K (e.g., Marley et al. 2012; Morley et al. 2012; Charnay et al. 2018).

#### D.6. SiO

SiO is an especially interesting molecule for tracing the abundance of the refractory silicates in the atmosphere such as MgSiO₃ or Mg₂SiO₄. As soon as the silicates evaporate (around 1600 K), their constituent atoms move into the gas phase. While atomic Mg is then the preferred gaseous form of Mg, Si will move into SiO (Visscher et al. 2010). For C/O ≥ 1, SiO enters the gas phase at around 1300 K, which is when SiC evaporates. SiO then starts to dissociate around 3500 K for C/O ≤ 1, while moving into SiS around 2000 K for C/O ≥ 1. The local evaporation-dependent temperatures given here could be lower than the observed transition as a function of planetary effective temperature, where high-pressure cloud formation could cold-trap Si into silicates.

#### D.7. Metal Oxides

Similar to SiO, the other metal oxides form as soon as the refractory clouds evaporate. Possible species of interest are TiO, VO, SiO, AlO, and CaO (e.g., Sharp & Burrows 2007; Gandhi & Madhusudhan 2019), with the relevant evaporation temperatures of the clouds ranging from ~1600 to 1900 K at our adopted standard conditions. Again, these are then only expected to be visible in the atmosphere if not cold-trapped into condensates at lower altitudes, that is, higher pressures. Except for SiO (see discussion in the SiO section above), most of these metal oxides are not expected to form in atmospheres with C/O ≥ 1 (Madhusudhan 2012; Gandhi & Madhusudhan 2019). In general, metal oxides will be destroyed by dissociation at high enough temperatures, with TiO and VO dissociating at temperatures similar to water (around 3000 K). As stated above, SiO is a bit more stable, dissociating from temperatures higher by about 500 K.

#### D.8. Metal Hydrides

Similar to metal oxides, metal hydrides such as FeH, CaH, MgH, NaH, CrH, and TiH may form as soon as the refractory
clouds have been evaporated, thus at local atmospheric temperatures of around 1600–1900 K and cooler temperatures for NaH, as Na₂S evaporates at ~900 K already. Of course, the cold trapping statement from above holds here as well. The hydrides will be destroyed by thermal dissociation at high enough temperatures, for example, MgH and FeH dissociate around 3000 K or so (Lothringer et al. 2018). For these two species it should also be noted that the main gas-phase carrier is atomic Mg and Fe in (ultra)hot Jupiter atmospheres and that MgH and FeH are less abundant by about a factor of 10⁴ (Visscher et al. 2010).

## D.9. Metal Atoms

Metal atoms can exist in the gas phase as soon as the sequestering refractory condensates evaporate (modulo cold trapping). As mentioned above, gaseous Fe and Mg are the main gas carriers of these elements once the silicates and iron condensates are gone. Si takes over as the main gas species only after SiO is dissociated. Fe, Mg, Ti, Ca, and Ni all are ionized between 3500 and 4000 K or so, with Al being ionized at somewhat lower temperatures (Kitzmann et al. 2018; Lothringer et al. 2018; Hoeijmakers et al. 2019).

## D.10. Metal Ions

Metal ions become visible in the atmosphere as soon as the atoms have been ionized; see immediately above.

## ORCID iDs

Paul Mollière https://orcid.org/0000-0003-4096-7067

Tamara Molyarova https://orcid.org/0000-0003-0448-6354

Bertram Bitsch https://orcid.org/0000-0002-8868-7649

Thomas Henning https://orcid.org/0000-0002-1493-300X

Aaron Schneider https://orcid.org/0000-0002-1448-0303

Laura Kreidberg https://orcid.org/0000-0003-0514-1147

Christian Eistrup https://orcid.org/0000-0002-8743-1318

Remo Burn https://orcid.org/0000-0002-9020-7309

Evert Nasedkin https://orcid.org/0000-0002-9792-3121

Dmitry Semenov https://orcid.org/0000-0002-3913-7114

Christoph Mordasini https://orcid.org/0000-0002-1013-2811

Martin Schlecker https://orcid.org/0000-0001-8355-2107

Kamber R. Schwarz https://orcid.org/0000-0002-6429-9457

Sylvestre Lacour https://orcid.org/0000-0002-6948-0263

## References

Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872

Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A, 564, A73

Akimkin, V., Zhukovska, S., Wiebe, D., et al. 2013, ApJ, 766, 8

Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109

Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Ataiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110

Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 393, 49

Baeyens, R., Decin, L., Carone, L., et al. 2021, MNRAS, 505, 5603

Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152

Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315

Barman, T. S. 2018, in Planetary Evaporation Through Evolution, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer), 29

Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, SSRv, 205, 77

Baudino, J.-L., Mollière, P., Venot, O., et al. 2017, ApJ, 850, 150

Benneke, B. 2015, arXiv:1504.07655

Bergez-Casalou, C., Bitsch, B., Pierens, A., Crida, A., & Raymond, S. N. 2020, A&A, 643, A133

Best, W. M. J., Liu, M. C., Magnier, E. A., & Dupuy, T. J. 2021, AJ, 161, 42

Binkert, F., Szulágyi, J., & Birnstiel, T. 2021, arXiv:2103.10177

Birnstiel, T., Fang, M., & Johansen, A. 2016, SSRv, 205, 41

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Bitsch, B., & Battistini, C. 2020, A&A, 633, A10

Bitsch, B., Forsberg, R., Liu, F., & Johansen, A. 2018a, MNRAS, 479, 3690

Bitsch, B., & Johansen, A. 2016, A&A, 590, A101

Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112

Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018b, A&A, 612, A30

Bitsch, B., Morbidelli, A., Lega, E., Kretke, K., & Crida, A. 2014, A&A, 570, A75

Bitsch, B., Raymond, S. N., Buchhave, L. A., et al. 2021, A&A, 649, L5

Bitsch, B., Trifonov, T., & Izidoro, A. 2020, A&A, 643, A66

Blecic, J., Dobbs-Dixon, I., & Greene, T. 2017, ApJ, 848, 127

Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994

Booth, R. A., & Ilee, J. D. 2019, MNRAS, 487, 3998

Booth, R. A., & Owen, J. E. 2020, MNRAS, 493, 5079

Bosman, A. D., Cridland, A. J., & Miguel, Y. 2019, A&A, 632, L11

Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182

Boss, A. P. 1997, Sci, 276, 1836

Boss, A. P. 2021, ApJ, 911, 146

Bourrier, V., & Lecavelier des Etangs, A. 2018, in Characterizing Evaporating Atmospheres of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer), 148

Brandt, G. M., Brandt, T. D., Dupuy, T. J., Li, Y., & Michalik, D. 2021, AJ, 161, 179

Brewer, J. M., Fischer, D. A., & Madhusudhan, N. 2017, AJ, 153, 83

Brogi, M., Line, M., Bean, J., Désert, J. M., & Schwarz, H. 2017, ApJ, 839, L2

Brogi, M., & Line, M. R. 2019, AJ, 157, 114

Brouwers, M. G., & Ormel, C. W. 2020, A&A, 634, A15

Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65

Brügger, N., Burn, R., Coleman, G. A. L., Alibert, Y., & Benz, W. 2020, A&A, 640, A21

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Caldas, A., Leconte, J., Selsis, F., et al. 2019, A&A, 623, A161

Chabrier, G., & Baraffe, I. 2007, ApJ, 661, L81

Chambers, J. 2019, ApJ, 879, 98

Chambers, J. E. 2010, ApJ, 724, 92

Chambers, J. E. 2016, ApJ, 825, 63

Changeat, Q., Al-Refaie, A. F., Edwards, B., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 913, 73

Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172

Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Crida, A., Morbidelli, A., & Masset, F. 2006, Icar, 181, 587

Cridland, A. J., Bosman, A. D., & van Dishoeck, E. F. 2020a, A&A, 635, A68

Cridland, A. J., Eistrup, C., & van Dishoeck, E. F. 2019, A&A, 627, A127

Cridland, A. J., Pudritz, R. E., & Alessi, M. 2016, MNRAS, 461, 3274

Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2020b, A&A, 642, A229

Crossfield, I. J. M., Barman, T., Hansen, B. M. S., Tanaka, I., & Kodama, T. 2012, ApJ, 760, 140

Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, SSRv, 212, 1

Cushing, M. C., Roellig, T. L., Marley, M. S., et al. 2006, ApJ, 648, 614

D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540

Danielski, C., Baudino, J.-L., Lagage, P.-O., et al. 2018, AJ, 156, 276

Debras, F., & Chabrier, G. 2019, ApJ, 872, 100

Dürmann, C., & Kley, W. 2015, A&A, 574, A52

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14

Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A69

Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A70

Feng, Y. K., Line, M. R., & Fortney, J. J. 2020, AJ, 160, 137

Feng, Y. K., Line, M. R., Fortney, J. J., et al. 2016, ApJ, 829, 52

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, OJAp, 2, 10

Fisher, C., & Heng, K. 2018, MNRAS, 481, 4698

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131

Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69

Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80
Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, AJ, 160, 288

Gail, H.-P., & Trieloff, M. 2017, A&A, 606, A16

Gandhi, S., & Madhusudhan, N. 2019, MNRAS, 485, 5817

Gandhi, S., Madhusudhan, N., Hawker, G., & Piette, A. 2019, AJ, 158, 228

Gao, P., Thorngren, D. P., Lee, G. K. H., et al. 2020, NatAs, 4, 951

Gautier, D., & Owen, T. 1989, in The Composition of Outer Planet Atmospheres, ed. S. K. Atreya, J. B. Pollack, & M. S. Matthews (Tucson, AZ: Univ. Arizona Press), 487

Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Natur, 592, 205

Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215

Ginzburg, S., & Chiang, E. 2020, MNRAS, 498, 680

Goyal, J. M., Mayne, N., Drummond, S., et al. 2020, MNRAS, 498, 4680

Gravity Collaboration, Nowak, M., Lacour, S., et al. 2020, A&A, 633, A110

Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17

Hasegawa, Y., & Pudritz, R. E. 2011, MNRAS, 417, 1236

Hayashi, C. 1981, PThPS, 70, 35

Helled, R., & Stevenson, D. 2017, ApJL, 840, L4

Helled, R., Stevenson, D. J., Lunine, J. I., et al. 2022, Icar, 378, 114937

Helled, R., Werner, S., Dorn, C., et al. 2022, ExA, 53, 323

Helling, C., Ackerman, A., Allard, F., et al. 2008, MNRAS, 391, 1854

Henning, T., & Semenov, D. 2013, ChRv, 113, 9016

Hobbs, R., Rimmer, P. B., Shorttle, O., & Madhusudhan, N. 2021a, MNRAS, 506, 3186

Hobbs, R., Shorttle, O., & Madhusudhan, N. 2021b, arXiv:2112.04930

Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165

Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020, A&A, 641, A123

Hu, R., Seager, S., & Yung, Y. L. 2015, ApJ, 807, 8

Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388

Izidoro, A., Bitsch, B., Raymond, S. N., et al. 2021, A&A, 650, A152

Jenkins, J. S., Díaz, M. R., Kurtovic, N. T., et al. 2020, NatAs, 4, 1148

Jin, S., & Mordasini, C. 2018, ApJ, 853, 163

Jurić, M., & Tremaine, S. 2008, ApJ, 686, 603

Khorshid, N., Min, M., Désert, J. M., Woitke, P., & Dominik, C. 2021, arXiv:2111.00279

Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183

Klahr, H., Pfeil, T., & Schreiber, A. 2018, in Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetesimal Formation, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer), 138

Kopparapu, R. K., Kasting, J. F., & Zahnle, K. J. 2012, ApJ, 745, 77

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, ApJL, 793, L27

Kreidberg, L., Line, M. R., Bean, J. L., et al. 2015, ApJ, 814, 66

Kress, M. E., Tiélens, A. G., & Frenklach, M. 2010, AdSpR, 46, 44

Lacy, B. I., & Burrows, A. 2020, ApJ, 905, 131

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35

Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83

Leconte, J., & Chabrier, G. 2012, A&A, 540, A20

Ledoux, P. 1947, ApJ, 105, 305

Li, J., Bergin, E. A., Blake, G. A., Ciesla, F. J., & Hirschmann, M. M. 2021, SciA, 7, eabd3632

Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846

Line, M. R., Brogi, M., Bean, J. L., et al. 2021, Natur, 598, 580

Line, M. R., Knutson, H., Wolf, A. S., & Yung, Y. L. 2014, ApJ, 783, 70

Line, M. R., Marley, M. S., Liu, M. C., et al. 2017, ApJ, 848, 83

Line, M. R., & Parmentier, V. 2016, ApJ, 820, 78

Lissauer, J. J., Hubickyj, O., D'Angelo, G., & Bodenheimer, P. 2009, Icar, 199, 338

Lissauer, J. J., & Stewart, G. R. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine (Tucson, AZ: Univ. Arizona Press), 1061

Lobo Gomes, A., Klahr, H., Uribe, A. L., Pinilla, P., & Surville, C. 2015, ApJ, 810, 94

Lodders, K. 2010, Exoplanet Chemistry (New York: Wiley), 157

Lodders, K. 2019, arXiv:1912.00844

Lodders, K., & Fegley, B. 2002, Icar, 155, 393

Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27

Lothringer, J. D., Fu, G., Sing, D. K., & Barman, T. S. 2020, ApJL, 898, L14

Lothringer, J. D., Rustamkulov, Z., Sing, D. K., et al. 2021, ApJ, 914, 12

Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001

Ludwig, H.-G., Allard, F., & Hauschildt, P. H. 2002, A&A, 395, 99

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

MacDonald, R. J., Goyal, J. M., & Lewis, N. K. 2020, ApJL, 893, L43

MacDonald, R. J., & Lewis, N. K. 2022, ApJ, 929, 20

MacDonald, R. J., & Madhusudhan, N. 2017, ApJL, 850, L15

Madhusudhan, N. 2012, ApJ, 758, 36

Madhusudhan, N. 2019, ARA&A, 57, 617

Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJL, 794, L12

Madhusudhan, N., Bitsch, B., Johansen, A., & Eriksson, L. 2017, MNRAS, 469, 4102

Madhusudhan, N., Harrington, J., Stevenson, K. B., et al. 2011, Natur, 469, 64

Marboeuf, U., Thiabaud, A., Alibert, Y., Cabral, N., & Benz, W. 2014, A&A, 570, A36

Marleau, G.-D., Coleman, G. A. L., Leleu, A., & Mordasini, C. 2019, A&A, 624, A20

Marley, M. S., Saumon, D., Cushing, M., et al. 2012, ApJ, 754, 135

Miles, B. E., Skemer, A. J. I., Morley, C. V., et al. 2020, AJ, 160, 63

Miley, J. M., Panić, O., Booth, R. A., et al. 2021, MNRAS, 500, 4658

Miller, N., & Fortney, J. J. 2011, ApJL, 736, L29

Mizuno, H. 1980, PThPh, 64, 544

Molaverdikhani, K., Henning, T., & Mollière, P. 2019a, ApJ, 873, 32

Molaverdikhani, K., Henning, T., & Mollière, P. 2019b, ApJ, 883, 194

Moll, R., Garaud, P., Mankovich, C., & Fortney, J. J. 2017, ApJ, 849, 24

Mollière, P., & Snellen, I. A. G. 2019, A&A, 622, A139

Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131

Mollière, P., van Boekel, R., Bouwman, J., et al. 2017, A&A, 600, A10

Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47

Molyarova, T., Akimkin, V., Semenov, D., et al. 2017, ApJ, 849, 130

Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139

Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41

Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172

Moses, J. I., Madhusudhan, N., Visscher, C., & Freedman, R. S. 2013, ApJ, 763, 25

Müller, J., Savvidou, S., & Bitsch, B. 2021, A&A, 650, A185

Müller, S., Helled, R., & Cumming, A. 2020, A&A, 638, A121

Nayakshin, S. 2010, MNRAS, 408, L36

Nayakshin, S., & Fletcher, M. 2015, MNRAS, 452, 1654

Nixon, M. C., & Madhusudhan, N. 2022, arXiv:2201.03532

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194

Ogihara, M., Hori, Y., Kunitomo, M., & Kurosaki, K. 2021, A&A, 648, L1

Ohno, K., & Ueda, T. 2021, A&A, 651, L2

Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43

Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512

Ormel, C. W., Vazan, A., & Brouwers, M. G. 2016, A&A, 647, A175

Owen, T., & Encrenaz, T. 2003, SSRv, 106, 121

Paardekooper, S. J., & Mellema, G. 2006, A&A, 453, 1129

Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91

Paunzen, E. 2004, in IAU Symp. 224, The A-Star Puzzle, ed. J. Zverko et al. (Cambridge: Cambridge Univ. Press), 443

Pelletier, S., Benneke, B., Darveau-Bernier, A., et al. 2021, AJ, 162, 73

Pinhas, A., Madhusudhan, N., & Clarke, C. 2016, MNRAS, 463, 4516

Pinhas, A., Madhusudhan, N., Gandhi, S., & MacDonald, R. 2019, MNRAS, 482, 1485

Pinilla, P., de Juan Ovelar, M., Ataiee, S., et al. 2015, A&A, 573, A9

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Piso, A.-M. A., Öberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109

Pluriel, W., Zingales, T., Leconte, J., & Parmentier, V. 2020, A&A, 636, A66

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Powell, D., Zhang, X., Gao, P., & Parmentier, V. 2018, ApJ, 860, 18

Robert, C. M. T., Crida, A., Lega, E., Méheut, H., & Morbidelli, A. 2018, A&A, 617, A98

Rossow, W. B. 1978, Icar, 36, 1

Sadakane, K. 2006, PASJ, 58, 1023

Sánchez-López, A., Landman, R., Mollière, P., et al. 2022, A&A, 661, A78

Sarkis, P., Mordasini, C., Henning, T., Marleau, G. D., & Mollière, P. 2021, A&A, 645, A79

Savvidou, S., Bitsch, B., & Lambrechts, M. 2020, A&A, 640, A63

Schib, O., Mordasini, C., Wenger, N., Marleau, G. D., & Helled, R. 2021, A&A, 645, A43

Schlecker, M., Pham, D., Burn, R., et al. 2021, A&A, 656, A73

Schmitt, W., Henning, T., & Mucha, R. 1997, A&A, 325, 569

Schneider, A. D., & Bitsch, B. 2021a, A&A, 654, A71

Schneider, A. D., & Bitsch, B. 2021b, A&A, 654, A72

Schulik, M., Johansen, A., Bitsch, B., & Lega, E. 2019, A&A, 632, A118

Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2019, ApJ, 877, 131
Seager, S., Richardson, L. J., Hansen, B. M. S., et al. 2005, ApJ, 632, 1122

Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42

Semenov, D., & Wiebe, D. 2011, ApJS, 196, 25

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Sharp, C. M., & Burrows, A. 2007, ApJS, 168, 140

Skilling, J. 2004, in AIP Conf. Ser. 735, Bayesian Inference and Maximum
Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, &
U. V. Toussaint (Melville, NY: AIP), 395

Spiegel, D. S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487

Stevenson, K. B., Bean, J. L., Madhusudhan, N., & Harrington, J. 2014, ApJ,
791, 36

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016,
A&A, 596, A74

Swain, M., Deroo, P., Tinetti, G., et al. 2013, Icar, 225, 432

Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65

Taylor, J., Parmentier, V., Irwin, P. G. J., et al. 2020, MNRAS, 493, 4342

Thorngren, D., & Fortney, J. J. 2019, ApJL, 874, L31

Thorngren, D., Gao, P., & Fortney, J. J. 2019, ApJL, 884, L6

Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016,
ApJ, 831, 64

Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, ExA, 46, 135

Turrini, D., Nelson, R. P., & Barbieri, M. 2015, ExA, 40, 501

Turrini, D., Schisano, E., Fonte, S., et al. 2021, ApJ, 909, 40

van't Hoff, M. L. R., Bergin, E. A., Jørgensen, J. K., & Blake, G. A. 2020,
ApJL, 897, L38

Vazan, A., Helled, R., & Guillot, T. 2018, A&A, 610, L14

Vazan, A., Helled, R., Kovetz, A., & Podolak, M. 2015, ApJ, 803, 32

Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118

Venot, O., Hébrard, E., Agúndez, M., et al. 2012, A&A, 546, A43

Visscher, C., Lodders, K., & Fegley, B., Jr. 2010, ApJ, 716, 1060

Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, GeoRL, 44, 4649

Wakeford, H. R., & Sing, D. K. 2015, A&A, 573, A122

Wakeford, H. R., Sing, D. K., Deming, D., et al. 2018, AJ, 155, 29

Wakeford, H. R., Visscher, C., Lewis, N. K., et al. 2017, MNRAS, 464,
4247

Wang, D., Miguel, Y., & Lunine, J. 2017, ApJ, 850, 199

Wang, J., Wang, J. J., Ma, B., et al. 2020, AJ, 160, 150

Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192

Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, ApJL, 887, L20

Wilson, H. F., & Militzer, B. 2010, PhRvL, 104, 121101

Woitke, P., Helling, C., & Gunn, O. 2020, A&A, 634, A23

Yorke, H. W., & Bodenheimer, P. 2008, in ASP Conf. Ser. 387, Massive Star
Formation: Observations Confront Theory, ed. H. Beuther, H. Linz, &
T. Henning (San Francisco, CA: ASP), 189

Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41

Zhang, Y., Snellen, I., Bohn, A. J., et al. 2021, Natur, 595, 370