Interpreting the atmospheric composition of exoplanets: sensitivity to planet formation assumptions

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.


INTRODUCTION
The distribution of bulk planetary properties such as mass, radius, and orbital parameters encode 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. 2020a;Schlecker et al. 2021). In addition, the chemical composition of planet atmospheres has long been regarded as a key to unlock 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 is determined by the process of planet formation. Consequently, there have been a number of studies which 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 1 , or just overall metal content (e.g., Öberg et al. 2011;Fortney et al. 2013;Marboeuf et al. 2014;Madhusudhan et al. 2014;Cridland et al. 2016;Mordasini et al. 2016;Madhusudhan et al. 2017;Lothringer et al. 2021;Schneider & Bitsch 2021a;Khorshid et al. 2021).
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 towards C/O 1 (Madhusudhan et al. 2011;Crossfield et al. 2012;Swain et al. 2013;Line et al. 2014;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 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. 2021, and the references therein).
Luckily, the quality of observational constraints on the planetary composition is expected to be rapidly improving. The advent of retrieval methods for mediumand high-resolution observations (e.g., Brogi et al. 2017;Brogi & Line 2019;Gibson et al. 2020), may allow 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, ; Pelletier et al. (2021); 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 isotpologues 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, later in the decade 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;Danielski 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., Line & Parmentier 2016;Feng et al. 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;Taylor et al. 2020;Feng et al. 2020;MacDonald et al. 2020;Pluriel et al. 2020;Lacy & Burrows 2020;Changeat et al. 2021;MacDonald & Lewis 2021;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 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, and 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.

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, 1dimensional disk, the authors calculated where important volatile gases such as H 2 O, CO 2 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 2 O freezes out first when moving outward, followed by CO 2 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. 2 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 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.

Disk elemental composition and structure
Constraining a planet's formation location based on its elemental abundance ratios (e.g., C/O) requires that 2 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 H2O iceline. 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 sets 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;Suzuki et al. 2016;Bai 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;Pinilla et al. 2015;Lobo Gomes et al. 2015;Binkert et al. 2021).
Another important effect determining the disk gas composition is the evaporation of pebbles inside of icelines, 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, which sets the pebbles' growth rates, Stokes parameters (i.e., drift speeds), and may trap them into local pressure maxima (e.g., Paardekooper & Mellema 2006;Ormel & Klahr 2010;Birnstiel et al. 2012;Lambrechts et al. 2014).

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;Molyarova et al. 2017;Eistrup et al. 2018). 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 2 ice over time, or the conversion of N 2 gas into NH 3 ice (Molyarova et al. 2017;Eistrup et al. 2018;Semenov & Wiebe 2011). 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 selfshadowing of the disk, allowing for nominally too cool compositions to occur closer to the star than otherwise expected (Ohno & Ueda 2021).

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 cm-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;Bitsch et al. 2018b;Ataiee et al. 2018), 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' icelines (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 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 proto-planet migrates into regions of the disk still containing planetesimals (Alibert et al. 2005).
Moreover, the accretion of gas by the growing planet is a 3-dimensional process (e.g., D'Angelo et al. 2003;Ayliffe & Bate 2009;Szulágyi et al. 2014;Ormel et al. 2015;Schulik et al. 2019Schulik et al. , 2020. This may be especially important as the gas composition and C/O value is thought to be changing above the mid-plane 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 discs viscous resupply of the planetary feeding zone (e.g., Lubow et al. 1999;Lissauer et al. 2009;Ayliffe & Bate 2009;Schulik et al. 2020;Bergez-Casalou 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 icelines, migration may be less of a problem if a forming planet did not migrate across icelines. Interestingly, planets forming in the inner disk may actually be trapped at locations just beyond the water iceline (e.g., Bitsch & Johansen 2016;Cridland et al. 2016;Müller et al. 2021). There also exist other traps not connected to icelines, 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;Izidoro et al. 2021;Emsenhuber et al. 2020a). 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 further out (e.g., Jurić & Tremaine 2008;. 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. 2020b;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 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. (2020a), 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, Neptunelike, 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'.

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);  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 presently investigated (Vazan et al. 2015(Vazan et al. , 2016Moll 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 also rainout processes 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 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 towards 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 which 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. (2021), 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.

Atmospheric evolution
The atmospheric composition of an exoplanet may also evolve due 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 due to thermal or non-thermal 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. Firstly, the atmosphere needs to be of low enough mass to allow for a significant amount to be lost. Secondly, 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 9779 b, 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 years before being removed, assuming impacts of kmsized comets every 20-200 years, 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 τ mix = H 2 P /K zz , where H P is the planetary pressure scale height and K zz the atmospheric eddy diffusion coefficient, and find A K zz value of 10 8 cm 2 /s 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. 2013a;Agúndez et al. 2014). For self-luminous planets the correct value to be chosen is unclear, but K zz = 10 5 cm 2 /s 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 towards 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 τ 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 radiative-convective boundary. By extension, the enrichment of the visible atmosphere may be lower if the impactors deposit their metals below the radiativeconvective boundary: 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 1-d 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: Here,Ṁ 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 RCB is the location of the radiative-convective boundary 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 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 ∆X/X 0 = 3 × 10 −4 to 3 × 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 1-d diffusion approximation, is needed.

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 be either 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 anaylsis setup by starting with the inversion of the formation model used inÖberg et al. (2011);Ö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 the effects of either 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.

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 can be measured with a given uncertainty, we are then interested in the probability distribution of formation parameters ϑ of a formation model M, given this measurement. Using Bayes' theorem, this can be written as Here, P (ϑ|M) is the prior probability of ϑ before considering any data C, while P (C|ϑ, M) is the likelihood for observing C, given that ϑ is true. In practice this  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 represent constraints of the planet formation process. Running inversions for various toy formation models then allows to study the impact of differing model assumptions.
may be written as, for example, where N species is the number of measured atmospheric species, C M i (ϑ) is the formation model prediction of the planet abundance of species i, and ∆C i are the measurement uncertainties. Here we chose 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 (ϑ|C, 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 model evidence ratios allow 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 (2019), assuming their static protoplanetary disk model as described for the young solar nebula, but with the icelines 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 iceline the volatile species is in the gas phase, 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 to map to the bulk composition of the planet: where M P is the total planet mass, and M solid is the mass of the solids (refractory species and volatile ices) accreted by the planet. It holds that M P = M gas + M solid , which we use to determine the amount of gas a planet accreted. The parameters a solid and a gas de- note the orbital distances where the solids and gas were accreted, respectively. For a given value of ϑ our setup then uses the disk model to determine the planet's accretion location with respect to the disk icelines. 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 C.

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;Eistrup et al. 2018;Bosman et al. 2018;Schwarz et al. 2018). If a large fraction of CO gas in a given region of the disk mid-plane 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 2-d axis-symmetric 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 are taken from the disk's chemical evolution calculations. Formally this also allows to add the formation time as an additional parameter to be constrained, but here we chose to initially only study and compare the inversion outcomes for different times during chemical evolution.

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 towards 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 icelines (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, indentical to the treatment of the disk's chemical evolution.

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 VLT, the community is already starting to derive precise C/O values using ground-based instruments. This has to do both with the data, but also the use of state-of-the-art retrieval techniques (Brogi & Line 2019;Mollière et al. 2020;Pelletier et al. 2021;Line 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 λ-Boötis-type star; this means that the abundances measured for its iron-peak elements are subsolar, with values of [Fe/H] = −0.55 ± 0.10 (Sadakane 2006) or [Fe/H] = −0.52 ± 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 ± 0.12, [O/H] = 0.12 ± 0.14, and (C/O)/(C/O) = 0.96 ± 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 volatilerich 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 λ Boo abundances of HR 8799 at face value, however. The iceline 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 with 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-σ 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 as well as relative gas-to-solid accretion fraction sets the O, C, and 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 which 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 was 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)  The inversion process finds a clear preference for the solids to stem from outside the CO iceline, or from within the H 2 O iceline. 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.   (i), but additionally placing a prior on the accreted solid mass, corresponding to a 20-M⊕ 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⊕). An arbitrary offset has been applied to the ellipses of Cases (i) to (iv), for clarity.
the C/O ratio set by the accretion of solids. Accreting solids of stellar C/O is possible at orbital distances outside of the icelines 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 iceline 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 iceline, 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 1-and 2-d 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 2 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 origin of the accreted solids is still outside the CO iceline. The formerly second likely location, inside the H 2 O iceline, 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.

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 yr −1 ), where L and M are the solar luminosity and mass. 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 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 upper right panel of Figure 4, 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 iceline 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 iceline 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 iceline of CO 2 over time. The conversion rate depends on the CO abundance in the ice, which drops rapidly inside the CO iceline. So while the reaction rate increases with temperature, the conversion is more efficient right inside the CO iceline (Bosman et al. 2018). Thus, the process occurs first for larger disk radii, and later for smaller radii, and drives the solid C/O towards the stellar value also inside the static CO iceline. We note the 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 inwards 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. Inreasing 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 iceline. In the inner part of the disk, within the CO 2 iceline, 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 iceline. 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.

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 sub-stellar metallicities and super-stellar C/O values. In contrast, planets with an atmospheric metal enrichment dominated by solid accretion may have have super-stellar metallicities, but sub-stellar C/O ratios (e.g., Öberg et al. 2011;Madhusudhan et al. 2014;Mordasini et al. 2016;Madhusudhan et al. 2017). In the case of pebble drift, however, evaporation of pebbles inside of the CO, CO 2 , and potentially the CH 4 icelines can lead to disk gas that is significantly enriched in these species, allowing for super-stellar metallicities and C/O ratios in the disk's Time (yr) 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 iceline, due to pebble evaporation, followed by its viscous spreading.
gas phase, and therefore in the atmospheres of planets (e.g., Booth et al. 2017;Schneider & Bitsch 2021a,b). For setting up the chemcomp pebble disk model, we used the same initial disk surface density and temperature structure 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 α = 5 × 10 −4 , where α is the usual dimensionless diffusion coefficient, in units of c s H, where c s is the local mid-plane 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 α of the order of 10 −3 − 10 −4 (Pinte et al. 2016;Flaherty et al. 2017Flaherty et al. , 2018. The disk viscosity is a key parameter for the pebble problem, with smaller values of α leading to larger pebbles, thus generally faster inward drift, and longer persistence time scales 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 µ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 iceline 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 iceline. The same effect is visible just outside the the CO 2 and H 2 O icelines, where the solid C/O values reach 0.5 and 0, respectively. Away from the icelines 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 due 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 iceline 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' C/O starts rising again as the CO gas that has evaporated off the pebbles inside the CO icelines 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 2 O iceline; the gas' 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 inwards.
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): like (i), but putting an upper limit of 20 M ⊕ 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): like (ii), but replacing the upper mass limit on the planetary mass from the spectroscopic retrieval (M P < 14 M 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.9 −1.8 M Jup . Case (iv): like (iii), but increasing the pebble isolation mass to 100 M ⊕ . 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 lower panel of Figure 4. In this scenario pebble drift, evaporation, and accretion is 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 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 icelines, the result that significant accretion of solids from outside the CO iceline is likely does not change when compared to the static setup of the disk composition. Just outside the CO iceline 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 iceline 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 iceline, 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 mount 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 super-stellar atmospheric metallicity, leads to a 1-σ upper limit of 570 M ⊕ 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 ⊕ on the accreted solid mass, a more likely scenario. We note that the pebble isolation mass is a function of the disk viscosity α, and that it is very sensitive to the disk aspect ratio (M iso ∝ [H/r] 3 , with H being the disk scale height). The value of 20 M ⊕ is what we derive for the HR 8799 disk model at the location of the CO iceline, using the scaling relations reported in Bitsch et al. (2018b). The compositional fit for Case (ii) is shown in the lower panel of 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 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 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 (cf. 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.  (iii), that is Miso = 20 M , MP = 9.6 +1.9 −1.8 MJup, at solar disk composition (left ellipse) with a setup where the λ Boo-type composition of HT 8799 was assumed for the disk instead (right ellipse). An offset was applied to these ellipses for clarity.
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 lower 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 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 probability distribution of a 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 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 iceline, in the outer parts of the disk. Due to the disk flaring, H/r increases towards the outer disk, and we would find M iso = 50 M ⊕ , 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 iso calculation (α = 0.004, instead of the nominal 0.0005); this results in M iso = 100 M ⊕ . The corresponding enrichment pattern of the planet is shown in the lower panel of Figure 4, rightmost ellipse. Unsurprisingly, note 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. Like expected, a solid values outside the CO iceline are the most likely for this case, with some additional gas accretion from within the CO 2 iceline.
From our investigation it thus becomes evident that pebbles alone, for average M 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 the our Case (ii), where the inversion when imposing M iso = 20 M Jup resulted in an a good fit by constraining the planetary mass to below 3.5 M Jup . Next, if the pebble isolation mass is much higher than our baseline case (e.g., 100 instead of 20 M ⊕ ), which is possible for large disk viscosities and the planet initially forming far out in the disk, pebble enrichment becomes 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 ⊕ 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 λ 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 iceline (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, thus shorter accretion timescales) it may have been less affected by the formation of the outer planets.

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 pile-up at the water iceline. 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 mid-plane 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); . The impact of metallicity gradients could potentially be studied by adding a parameter f mix which 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 super-stellar 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 iceline. 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 icelines). If f mix was added as a free parameter, small f mix values would again have yielded regions outside the CO iceline 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 which 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) which 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 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

Aspect
Potential tractability in formation model inversions Disk composition and structure Unknown disk elemental abundances Scale using stellar [Fe/H], try varying composition according to scaling uncertainties. Available solid reservoir Impose limit based on likely disk mass and dust-to-gas ratio. Disk (thermal) structure Feed in disk structures from dedicated disk models. Explore if parameterizing 3-d effects in 1-d model is possible. Changes in disk structure will affect, e.g., iceline positions, as in Ohno & Ueda (2021). Planetary back-reaction on disk Use simplified gap opening criteria to limit gas accretion, compare to disk lifetimes. Apply pebble isolation mass (limit refractory reservoir accessible to planet).

Include pebble drift & evaporation at icelines
Increase gas metallicity inside of icelines as function of time, also see Sect. 3.2.

Disk chemistry
Chemical evolution of disk Run formation inversion with chemical composition as function of time. Also see Section 3.2. Inherited or 'reset' disk abundances Explore impact of differing assumptions on disk abundances for the inversion process. Cosmic ray ionization and stellar irradiation Use best guesses for retrievals, otherwise explore different values.

Refratory carbon depletion in inner disk
Explore via on/off switch.

Planet formation
Pebble and planetesimal accretion Compare inferred solid (refractory) mass of planet with isolation masses. Constrain upper limit on accreted M pebbles , lower limit on accreted M planetesimals .

3-d planet accretion
Test impact of parameterizations, for example vertically averaged abundances for gas accretion Planet migration Allow to fit for multiple formation locations? Add priors on formation location: planet traps? Add priors enforcing inward migration (e.g., agas ≤ a solid )? Leveraging full complexity of formation models Explore use of machine learning techniques.
E.g., random forest predictors as demonstrated in Schlecker et al. (2021). Planet formation by gravitational instability Likely treatable, but requires changes. E.g., steady state viscous disk model → infall disk model Planet bulk -atmosphere coupling Metallicity gradient inside planet Use multi-component model, infer mixing efficiency f mix ∈ [0, 1] to reveal correlations with other formation parameters.

Atmospheric evolution Evaporation
Can be important for lower mass planets or gas planets with a metallicity gradient. Use inverted evaporation models to reveal correlation with formation parameters? Infall of comets / asteroids Better quantitative modeling needed.
Potentially not important for gas giant planets. Table 1. Aspects of planet formation and their potential treatment toy formation model inversions.
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.

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 James Webb Space Telescope, 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 summa-rize the compositional constraints that can be extracted from such atmospheric measurements and how the resulting atmospheric enrichment patterns for the planetary population may allow to constrain planet formation in a broader sense.

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.  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 bar. We either assumed 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.
drift and evaporation may lead to super-stellar enrichments of the gas at C/O values both smaller or larger than stellar. Interestingly, however, super-stellar C/O values and enrichments are difficult to obtain without considering pebble evaporation (also see Section 3.2), so a large population of planets which 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 towards 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 2 O, CH 4 , CO, CO 2 , HCN, C 2 H 2 . C/O therefore regulates the spectral appearance of a planet in the near-to mid-infrared (e.g., Seager et al. 2005;Fortney 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 2 O and CH 4 , 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 4 will thus be visible; for high temperatures and C/O 1, H 2 O will be visible. For further increasing temperatures and C/O 1, CH 4 is replaced by increasing amounts of C 2 H 2 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 which are still strongly cooling (with a high internal temperature) or heated by processes such as eccentricity dampening, CH 4 may be strongly suppressed . As mentined, ground-based highcontrast or high-resolution observations have started to obtain the first useful constraints on C/O. The state-ofthe-art will greatly improve once JWST and later Ariel will allow for a larger census of planetary compositions (also see our discussion in section 1).

N/O, N/C
The importance of atmospheric nitrogen-bearing species such as NH 3 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 2 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 icelines of C and O-bearing species are crossed towards larger orbital radii, while N 2 stays in the gas phase (Turrini et al. 2021). If the planet forms at wide enough orbital distances, eventually N 2 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 2 iceline 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;Öberg & Wordsworth 2019;Bosman et al. 2019;Cridland et al. 2020b). Similar to the 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) or 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 2 , as cosmic ray ionization may process N 2 to NH 3 ice over Myr timescales, such that the importance of NH 3 and its iceline increases over time, with the iceline of NH 3 being much closer to the star than the one of N 2 (Semenov & Wiebe 2011).
In exoplanets the only spectrally active nitrogen bearing species of relevance are NH 3 and HCN. N 2 , which is the dominating nitrogen bearer at larger temperatures, has negligible opacity in the near-and mid-infrared. NH 3 , on the other hand, should be detectable in the mid-IR using JWST in exoplanet atmospheres (e.g., Danielski et al. 2018). Moreover, evidence for NH 3 has been seen in high-resolution studies (Giacobbe et al. 2021;Sánchez-López et al. 2022). For C/O 1, NH 3 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 3 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 the references therein). For self-luminous planets disequilibrium chemistry may play less of a role for N, as iso-abundance 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 2 O, CO, HCN, C 2 H 2 , NH 3 , and CH 4 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 3-d general circulation models (e.g., Agúndez et al. 2014;Baeyens et al. 2021), which could also be connected to condensate rain-out (Sánchez-López et al. 2022). Also photochemistry is important, especially in the upper atmospheric layers (e.g., Venot et al. 2012;Kopparapu 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.

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, also see 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 2 S and PH 3 may be useful refractory tracers at intermediate atmospheric temperatures (Wang et al. 2017;Öberg & Wordsworth 2019). While H 2 S and PH 3 are volatile species, the dominant carrier of P and S atoms in a protoplanetary disk appear to be refractory species (Öberg & Wordsworth 2019). Moreover, measuring the abundances of Na and K in planetary atmospheres may be worthwhile tracers 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 3 and Mg 2 SiO 4 are especially interesting, as they may lead to visible absorption features around 10 micron (e.g., Cushing et al. 2006;Wakeford & Sing 2015). Due to the complex micro-physical 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 micron with JWST's MIRI instrument. SiO is promising because it is the most abundant Si-bearing gas species after the silicates evaporate (Viss-cher et al. 2010), and is more stable than H 2 O against dissociation (by about 500 K). This should allow detecting this species in ultra hot Jupiters.
Other metal oxides such as TiO, VO, AlO, CaO have features in the optical and near-infrared (e.g., Sharp & Burrows 2007;Lothringer et al. 2020). Similarly, metal hydrides may be useful refractory tracers, at similar temperatures as the metal oxides. Species such as FeH, CaH, MgH, NaH, CrH, TiH all have absorption features in the optical and near-infrared (Sharp & Burrows 2007;. 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 ultra hot atmospheres of KELT-9b and WASP-121b in the optical (Hoeijmakers et al. 2019(Hoeijmakers et al. , 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).

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. It will have to be seen in how far this method can be used to predict planetary abundances.
Apart from this conclusion, we also introduce a method that allows 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 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 is 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 the 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 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 2 O, CO, CH 4 , CO 2 , C 2 H 2 , HCN), nitrogen content (NH 3 and HCN), and refractory content (H 2 S, PH 3 , 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 con-strain the importance of pebble drift and evaporation, and 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  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, also see 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 1-d diffusion. We thus write, using the usual 1-d diffusion equation for concentrations (e.g., Parmentier et al. 2013b): where ρ is the atmospheric density, X the atmospheric water mass fraction, t the time, z the atmospheric altitude,Ṁ the mass accretion rate of pure-water comets, and R P the planetary radius. W (z, z i ) is defined as with Θ being the Heaviside step function. This means that we assume that the comets are destroyed in a narrow layer of width ∆z i , at altitude z i in the atmosphere. Using the equation of hydrostatic equilibrium (∂P/∂z = −ρg), together with the equation of state of an ideal gas (P = ρk B T /µ) and that H P = k B T /µg one finds that once can express Equation A1 as where P is the atmospheric pressure, g the gravity, k B the Boltzmann constant, T the atmospheric temperature and µ the atmospheric mean molecular weight. For a steady state ansatz and integration over P one obtains that for P ≥ P i + ∆P i /2, and ∂X/∂P = 0 for P ≤ P i − ∆P i /2, and a linear transition between these two cases for P ∈ (P i − ∆P i /2, P i + ∆P i /2). In the following we assume that ∆P i P i . Thus it holds that the enrichment in units of mass fractions ∆X = X − X 0 at P i and at lower pressures (higher altitudes) is where K zz was assumed to be constant for simplicity. P RCB denotes the location of the radiative-convective boundary 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 ∆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 bar if impacting Jupiter at terminal velocity. From our condition that P RCB > P i this means that we need to consider P RCB > 100 bar. For warm Jupiters the maximum depth of the RCB may be as deep as P RCB = 200 bar 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 cm 2 s −1 , P RCB = 200 bar and P i = 100 bar results in a relative enrichment of ∆X/X 0 = 3 × 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 × 10 6 yr). Therefore, combining reasonable estimates for the parameters describing the atmosphere with a 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 × 10 6 years). 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 radiative convective boundary. The deep convective K zz , estimated from mixing length theory, may be in the range of K zz = 10 9 cm 2 s −1 (see, e.g., Equation 4 of Zahnle & Marley 2014). Extending the integration domain to 2 × 10 4 bar and setting K zz = 10 9 cm 2 s −1 for P > P RCB leads to ∆X/X 0 = 3 × 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 → ∞, 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 non-negligible ∆X/X 0 values). 20-km comets, 20 impacts/yr, P i = 1 bar Eq. E.3, numerical integration, t → ∞ Eq. E.5, analytic  (2019), which assumes a static power-law for the disk temperature and density of the young solar nebula. Inside its iceline position a given volatile species is in the gas phase, 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 setting its iceline position to zero. The background gases H 2 and He are included by setting their iceline to 1000 au. This ensures that the refractory and background species stay condensed/gaseous within the simulation domain. Table 2 lists the mass fractions, iceline positions, and atomic composition of all considered disk species and their constituent atoms. We note that the iceline 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 proto-solar nebula composition fromÖberg & Wordsworth (2019) with some modifications, see below. The abundance of a species i is given as number fractions x i inÖberg & Wordsworth (2019), compared to the number of hydrogen atoms n H . This was converted into mass fractions m i relative to the total disk mass using with µ i being the molecular mass of species i in atomic mass units. This expression is obtained from the disk's total mass M tot , the mass of species M i of species i, and setting where it was assumed that most of the disk mass is contributed by H and He atoms. Setting µ H = 1, µ He = 4 as well as assuming that n He /n H = 0.1 (see Table 8 in Lodders 2019, for the recommended proto-solar abundances) leads to the relation given in Equation B6. 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 2 O. We assumed that the silicates consist purely of MgSiO 3 , 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 the references within) we decided to forego a more careful treatment of the organic carbon reservoir for our conceptual study here. Iron, sulphur 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 2 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 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 XUV 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 ISM. 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.

C. FULL INVERSION POSTERIOR OF HR 8799E
In the left panel of Figure 12 we show the 1-and 2-d 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 ± 0.5, R P = 1.12 ± 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 2 P /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 ⊕ . 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 ± 0.25), and the total mass of the planet. The increased probability of the planet having accreted solids from outside the CO iceline or inside the H 2 O iceline (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 within between the H 2 O and CO icelines. 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: within these two icelines the solids' C/O is sub-stellar and the planetary C/O can be raised to higher values by accreting gas which 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 λ Boo-type abundances are assumed for the disk of HR 8799. The solution inside the H 2 O iceline is no longer valid for a solid , due to the high local solid C/O ratio.

D. REFRACTORY CHEMISTRY
Here we give a short description of the chemical behaviour 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 bar, whereas dissociation and ionization values were obtained from assuming pressures of 0.1 to 0.001 bar. We either assumed 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;Fortney et al. 2020;Parmentier et al. 2016). The temperatures given here do therefore 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.
H 2 S H 2 S condenses into NH 4 SH at ∼200 K, the higher temperature condensates MnS, ZnS, Na 2 S are of minor importance (Lodders 2010). For C/O < 1, H 2 S dissociates at ∼ 2000 K, while it moves into species like CS for C/O > 1 at temperatures around 1500-2000 K.
PH 3 PH 3 condenses into H 3 PO 4 at 500 K. Its presence in Jupiters atmosphere at lower temperature indicates a deep quenching point however, such that PH 3 may still be visible at lower temperatures (see, e.g., Baudino et al. 2017). For temperatures approaching 1000 K, PH 3 is increasingly converted into PH 2 .

NA
Na condenses into Na 2 S at ∼ 900 K. In principle alkalis such as Na could also be sequestered into high-temperature condensates such as feldspars. However, this does likely not occur due to the rainout of silicates (e.g., Line et al. 2017) which 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.
K K condenses into KCl at ∼ 900 K. In analogy to Na, sequestration of K into feldspars does likely not occur due to silicate rainout. Thus, K is in the gas phase from ∼900 K to ∼2000 K, after which it is ionized. Ionization occurs at temperatures roughly 500 K cooler than for Na.

REFRACTORY CLOUDS
As mentioned above, the refractory cloud species Na 2 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 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 3 and Mg 2 SiO 4 evaporate at temperatures around ∼1600 K, while iron clouds are stable until ∼1700 K. VO and calcium titanates are stable until 1600 to 1800 K, respectively. Aluminum-bearing condensates such as Al 2 O 3 , which are among the most stable ones, evaporate around 1900 K. Among the species listed here, potentially only silicates, Al 2 O 3 , 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 dwarf this temperature-dependent removal of silicate clouds is thought to cause the L-T transition, which typically occurs at T eff = 1200 to 1400 K (e.g., Best et al. 2021, and the references therein), while for planets and low-gravity brown dwarfs this limiting temperature may be as low as approximately 1000 K (e.g., Morley et al. 2012;Marley et al. 2012;Charnay et al. 2018).

SIO
SiO is an especially interesting molecule for tracing the abundance of the refractory silicates in the atmosphere such as MgSiO 3 or Mg 2 SiO 4 . 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.

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, CaO (e.g., Sharp & Burrows 2007;, with the relevant evaporation temperatures of the clouds ranging from ∼1600-1900 K at our adopted standard condidtions. 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;. 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.

METAL HYDRIDES
Similar to metal oxides, metal hydrides such as FeH, CaH, MgH, NaH, CrH, TiH may form as soon as the refractory clouds have been evaporated, thus at local atmospheric temperatures of around 1600 to 1900 K and cooler temperatures for NaH, as Na 2 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 10 4 (Visscher et al. 2010).

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, Ni all are ionized between 3500 to 4000 K or so, with Al being ionized at somewhat lower temperatures (Lothringer et al. 2018;Hoeijmakers et al. 2019;Kitzmann et al. 2018).

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