With four companions at separations from 16 to 71 au HR 8799 is a unique target for direct imaging presenting an opportunity for a comparative study of exoplanets with a shared formation history. Combining new VLTIGRAVITY observations obtained within the ExoGRAVITY program with archival data we performed a systematic atmospheric characterisation across all four planets. We explored different levels of model flexibility to understand the temperature structure chemistry and clouds of each planet using both petitRADTRANS atmospheric retrievals and fits to selfconsistent radiativeconvective equilibrium models. Using Bayesian model averaging to combine multiple retrievals a total of 89 across all four planets we find that the HR 8799 planets are highly enriched in metals with MH 1 and have stellar to superstellar atmospheric CO ratios. The CO ratio increases with increasing separation from 0.55SUB0.10SUBSUP0.12SUP for d to 0.78SUB0.04SUBSUP0.03SUP for b with the exception of the innermost planet which has a CO ratio of 0.87  0.03. Such high metallicities are unexpected for these massive planets and challenge planetformation models. By retrieving a quench pressure and using a disequilibrium chemistry model we derive vertical mixing strengths compatible with predictions for highmetallicity selfluminous atmospheres. Bayesian evidence comparisons strongly favour the presence of HCN in HR 8799 c and e as well as CHSUB4SUB in HR 8799 c with detections at gt 5 confidence. All of the planets are cloudy with no evidence of patchiness. The clouds of c d and e are best fit by silicate clouds lying above a deep iron cloud layer while the clouds of the cooler HR 8799 b are more likely composed of NaSUB2SUBS. With welldefined atmospheric properties future exploration of this system is well positioned to unveil further details of these planets extending our understanding of the composition structure and formation history of these siblings.
['radiative transfer', 'instrumentation: interferometers', 'methods: observational', 'planets and satellites: atmospheres', 'planets and satellites: composition', 'Astrophysics - Earth and Planetary Astrophysics']
Fourofakind Comprehensive atmospheric characterisation of the HR 8799 planets with VLTIGRAVITY
{'Four-of-a-kind? Comprehensive atmospheric characterisation of the HR 8799 planets with VLTI/GRAVITY': 'E. Nasedkin 1, ⋆ , P. Mollière 1 , S. Lacour 2 , 3 , M. Nowak 4 , L. Kreidberg 1 , T. Stolker 5 , J. J. Wang 6 , W. O. Balmer 7 , 8 , J. Kammerer 3 , J. Shangguan 9 , R. Abuter 3 , A. Amorim 10 , 11 , R. Asensio-Torres 1 , M. Benisty 12 , J.-P. Berger 12 , H. Beust 12 , S. Blunt 6 , A. Boccaletti 2 , M. Bonnefoy 12 , H. Bonnet 3 , M. S. Bordoni 9 , G. Bourdarot 9 , W. Brandner 1 , F. Cantalloube 13 , P. Caselli 9 , B. Charnay 2 , G. Chauvin 14 , A. Chavez 6 , E. Choquet 13 , V. Christiaens 15 , Y. Clénet 2 , V. Coudé du Foresto 2 , A. Cridland 5 , R. Davies 9 , R. Dembet 2 , J. Dexter 16 , A. Drescher 9 , G. Duvert 12 , A. Eckart 17 , 18 , F. Eisenhauer 9 , N. M. Förster Schreiber 9 , P. Garcia 11 , 19 , R. Garcia Lopez 20 , 1 , E. Gendron 2 , R. Genzel 9 , 21 , S. Gillessen 9 , J. H. Girard 8 , S. Grant 9 , X. Haubois 22 , G. Heißel 23 , 2 , Th. Henning 1 , S. Hinkley 25 , S. Hippler 1 , M. Houllé 14 , Z. Hubert 12 , L. Jocou 12 , M. Keppler 24 , P. Kervella 2 , N. T. Kurtovic 9 , A.-M. Lagrange 12 , 2 , V. Lapeyrère 2 , J.-B. Le Bouquin 12 , D. Lutz 9 , A.-L. Maire 12 , F. Mang 9 , G.-D. Marleau 26 , 27 , 28 , 1 , A. Mérand 3 , J. D. Monnier 29 , C. Mordasini 28 , T. Ott 9 , G. P. P. L. Otten 30 , C. Paladini 22 , T. Paumard 2 , K. Perraut 12 , G. Perrin 2 , O. Pfuhl 3 , N. Pourré 12 , L. Pueyo 8 , D. C. Ribeiro 9 , E. Rickman 31 , J.B. Ru ffi o 32 , Z. Rustamkulov 33 , T. Shimizu 9 , D. Sing 7 , 33 , J. Stadler 34 , 35 , O. Straub 35 , C. Straubmeier 17 , E. Sturm 9 , L. J. Tacconi 9 , E. F. van Dishoeck 5 , 9 , A. Vigan 13 , F. Vincent 2 , S. D. von Fellenberg 18 , F. Widmann 9 , T. O. Winterhalder 3 , J. Woillez 3 , ¸S. Yazıcı 9 , and the GRAVITY Collaboration  \n(A ffi liations can be found after the references)  \nReceived 2024-01-24; accepted 2024-04-25', 'ABSTRACT': 'With four companions at separations from 16 to 71 au, HR 8799 is a unique target for direct imaging, presenting an opportunity for a comparative study of exoplanets with a shared formation history. Combining new VLTI / GRAVITY observations obtained within the ExoGRAVITY program with archival data, we performed a systematic atmospheric characterisation across all four planets. We explored di ff erent levels of model flexibility to understand the temperature structure, chemistry, and clouds of each planet using both petitRADTRANS atmospheric retrievals and fits to selfconsistent radiative-convective equilibrium models. Using Bayesian model averaging to combine multiple retrievals (a total of 89 across all four planets), we find that the HR 8799 planets are highly enriched in metals, with [M / H] ≳ 1, and have stellar to superstellar atmospheric C / O ratios. The C / O ratio increases with increasing separation from 0 . 55 + 0 . 12 -0 . 10 for d to 0 . 78 + 0 . 03 -0 . 04 for b, with the exception of the innermost planet, which has a C / O ratio of 0 . 87 ± 0 . 03. Such high metallicities are unexpected for these massive planets, and challenge planet-formation models. By retrieving a quench pressure and using a disequilibrium chemistry model, we derive vertical mixing strengths compatible with predictions for high-metallicity, self-luminous atmospheres. Bayesian evidence comparisons strongly favour the presence of HCN in HR 8799 c and e, as well as CH4 in HR 8799 c, with detections at > 5 σ confidence. All of the planets are cloudy, with no evidence of patchiness. The clouds of c, d, and e are best fit by silicate clouds lying above a deep iron cloud layer, while the clouds of the cooler HR 8799 b are more likely composed of Na2S. With well-defined atmospheric properties, future exploration of this system is well positioned to unveil further details of these planets, extending our understanding of the composition, structure, and formation history of these siblings.  \nKey words. Methods: observational - planets and satellites: atmospheres - instrumentation: interferometers - radiative transfer', '1. Introduction': "Directly imaged exoplanets provide an ideal laboratory for understanding the formation and evolution of planetary systems. These young systems provide unique insight into widely separated companions ( ≳ 10 au): by directly measuring their emission spectra, we can peer into regions of their atmospheres inaccessible through other techniques. While spectroscopically similar to their brown dwarf cousins, these young, low-surfacegravity exoplanets display unique spectral shapes and colours (Faherty 2018), indicative of di ff erences in their chemistry, clouds, and formation history (Marley et al. 2010, 2012; Charnay et al. 2018).  \nThe HR 8799 system is a benchmark target for directly imaged exoplanets, containing four planets (Marois et al. 2008,  \n2010), an inner debris disc (Boccaletti et al. 2023), and an outer Kuiper-belt-like disc (Su et al. 2009). This is among the beststudied systems of exoplanets, with a wide range of photometric and spectroscopic data. The spectroscopic data cover the nearinfrared region (1-4 µ m) at varying spectral resolution, while the photometric data extend out to 15 µ m with the recent addition of JWST / MIRI observations (Boccaletti et al. 2023). Most of these studies, together with extensive modelling work, have tried to answer the following main questions:  \n- 1. How did the HR 8799 system form?\n- 2. What are the dynamics of the system? Is it stable, and how do the planets and disc interact?\n- 3. What are the atmospheres of each planet made of, and how have they evolved through time?  \nIn the present work, we attempt to directly answer question (3), which has implications for question (1). Using Bayesian atmo-  \nspheric retrievals (e.g. Madhusudhan 2019) as well as fits to 1D self-consistent models, we infer the atmospheric properties of each of the four planets. To date, the only comprehensive retrieval study of all four of the HR 8799 planets was by Lavie et al. (2017). New, high-precision K-band spectra obtained with the VLTI / GRAVITY as part of the ExoGRAVITY program (Gravity Collaboration et al. 2019; Lacour et al. 2020), together with updated atmospheric models and opacity databases, provide motivation and the means to perform a systematic reanalysis of this system.  \nWith e ff ective temperatures in the range of 1000-1400 K, the HR 8799 planets sit in the middle of the L / T transition (Kirkpatrick et al. 1999). This spectral transition is marked by changes in chemistry between L- and T-type objects, from CO-dominated carbon chemistry in the hotter objects to methane chemistry as the temperature falls below ∼ 1300 K. While this transition is well established for brown dwarfs, detections of CH4 in exoplanets remain elusive: there have been tentative detections in the atmosphere of HR 8799 b (Barman et al. 2015; Ru ffi o et al. 2021), but the only convincing detections have come from JWST observations of cool, transiting exoplanets (Bell et al. 2023; Madhusudhan et al. 2023) and the coldest directly imaged exoplanets, such as 51 Eridani b (Brown-Sevilla et al. 2023; Whiteford et al. 2023). This 'missing methane' is thought to be driven by convective upwelling in the atmospheres, dredging material from deeper, hotter regions of the atmosphere where CO is more favoured by equilibrium chemistry (Fegley & Lodders 1996). Precise constraints on the abundance of both CO and CH4 would allow better constraints on this vertical mixing, which is typically parameterised by the vertical di ff usion coe ffi cient K zz.  \nThe sharp change in colour in the L / T transition is thought to be caused by the sinking of silicate clouds through the atmosphere, as the cloud base shifts deeper as the e ff ective temperature decreases (Burrows & Sharp 1999). Once the cloud base sinks below the photosphere, the impact of the cloud opacity is increasingly removed from the spectrum, causing the blueward shift characteristic of T dwarfs as the e ff ective temperature falls below 1300 K. Following the mid-infrared observations with Spitzer / IRS (Cushing et al. 2006, 2008), Suárez & Metchev (2022) identified a trend in the silicate absorption feature at 9 µ mas a function of temperature. The strength of this absorption feature was found to correlate positively with the nearinfrared colour for L dwarfs, which is often used as a proxy for cloudiness. The HR 8799 planets lie comfortably below the temperature at which silicate clouds are expected to occur entirely below the photosphere, yet their red colour and near-infrared spectral shape are thought to be clear hallmarks of thick silicate cloud coverage (Mollière et al. 2020). However, Line et al. (2015) and Suárez & Metchev (2023) find that these clouds are not only sensitive to temperature, but also to surface gravity, which plays a role in determining the size and therefore settling speed of the aerosol particles. As young, giant exoplanets still retain significant heat from formation, their atmospheres remain inflated due to low surface gravity, which will in turn result in cloud properties that are unique to this class of object; observations of VHS 1256 b (Miles et al. 2023) remain the only spectroscopic observations of a silicate feature in a directly imaged planet. Burningham et al. (2021) and Vos et al. (2023) use atmospheric retrievals to identify the detailed structure and composition of the clouds, providing for the first time evidence to support the use of particular cloud compositions in these substellar atmospheres.  \nThe mechanism through which four super-Jupiter planets can form in a single system is unclear. The presence of both an inner  \nand an outer debris disc implies that the planets formed within a circumstellar disc; that is to say they did not form like stars. However, it is still unclear whether these objects formed through gravitational instability (GI) (Bodenheimer 1974; Adams et al. 1989) or via core accretion (Pollack et al. 1996). Evolutionary models (Saumon & Marley 2008) suggest that the current temperatures of the planets suggests hot-start boundary conditions for their evolution, which is more typically associated with GI (but also see Mordasini et al. 2017). GI models, such as that of Helled & Bodenheimer (2010), find that the amount of heavy elements accreted by the planets should be small, implying nearly stellar compositions for all four planets. Likewise, current composition estimates suggest that the planets share a C / Oratio with their host star (Hoch et al. 2023), but may be slightly enriched in metals, leading to tension with the predictions of the GI models.  \nIn addition to understanding the formation mechanism, Mollière et al. (2022) present a framework through which we can infer the conditions of the formation environment from measured atmospheric parameters. However, these authors, and many others (e.g. Eistrup et al. 2018; Cridland et al. 2019, 2020; Turrini et al. 2021; Pacetti et al. 2022), demonstrate that this is not a straightforward task. The Öberg et al. (2011) model links the planet C / O ratio to the location of formation relative to snow lines in the disc. This model provides a simplified view through which we can understand the impact of disc conditions on the outcomes of planet formation, but the complex and time-evolving physics and chemistry of discs and forming planets make solving the inverse problem challenging. Nevertheless, the best hope for linking the atmospheric properties back to the protoplanetary disc is to infer robust atmospheric elemental abundances and link these to interior models to determine the bulk planetary composition (Guillot 2005; Fortney et al. 2011), thereby determining what disc conditions could lead to the diversity of planet-formation outcomes.  \nWhile new data and modelling techniques are beneficial, interpreting such model-data comparisons for exoplanet spectra is far from trivial. Multiple techniques must be studied simultaneously to paint a consistent portrait of these worlds. Biases in inferred planet parameters are a common challenge in directimaging analyses: fits to emission spectra often find unphysically small radii that are inconsistent with evolutionary tracks (Marley et al. 2012). Retrievals using free molecular abundances tend to find higher C / Oratios than when disequilibrium chemistry models are considered (Lavie et al. 2017; Wang et al. 2020), possibly due to additional oxygen sequestered in refractory clouds (Fonte et al. 2023). The inferred e ff ective temperatures ( T e ff ) of each planet can vary by hundreds of kelvin, within a region of parameter space where the chemical timescales can vary by orders of magnitude over tens of kelvin (Zahnle & Marley 2014). Complicating matters further are the known discrepancies between spectral measurements (such as between SPHERE and GPI in the Hband; see Mollière et al. 2020), leading to uncertainties in both the shape and overall flux calibration of the spectra that are not reflected in the formal uncertainties. Attempting to address this problem, Nixon et al. (2023) demonstrate the use of Bayesian model averaging (BMA), which can be used to combine the posterior distributions of multiple models, allowing some degree of model uncertainty to be formally incorporated into the inferred parameter uncertainties. Finally, Greco & Brandt (2016) and Nasedkin et al. (2023) demonstrate the importance of properly accounting for the covariance in low-resolution IFS data a thorough treatment of IFS data is necessary to ensure meaningful and unbiased posterior probability distributions.  \n1e 6  \nFig. 1. HR 8799 planets as imaged in the H-band with the Gemini / GPI IFU, originally published in Greenbaum et al. (2018). The IFU cube was processed using KLIP , and the image is mean combined along the spectral axis. HR 8799 b is outside the field of view of GPI.  \n<!-- image -->  \nMany of the questions of chemistry and formation will be addressed through the use of the various instruments aboard JWST. This telescope will open new observational windows, extending out to the mid-infrared, and allow new characterisation methods, such as molecular mapping of the system (Patapis et al. 2022). Simultaneous measurements of CO and CH4 features between 3 and 5 µ m will allow constraints to be placed on the vertical mixing in the atmosphere, and more precise estimates of the C / O and metallicity. Nevertheless, ground-based observations remain crucial: the innermost companion will remain challenging to measure spectroscopically without a coronagraph; across most of its spectra, HR 8799 e is below the 2 × 10 -5 contrast threshold at 300 mas obtained in Ru ffi o et al. (2023).  \nThe present study provides a comprehensive examination of the atmospheres of the HR 8799 companions, making use of new, highS / N observations obtained with VLTI / GRAVITY, together with a combination of retrieval methods and selfconsistent modelling. We present further context and background information on the HR 8799 system in Section 2. The data used in this work are described in Section 3, while the details of the petitRADTRANS ( pRT ) forward model are described in Section 4, with the self-consistent models introduced in Section 4.7. The results of the atmospheric retrievals and selfconsistent grid fits are presented in Section 5. We discuss the limitations of this study, additional work to validate our results, and the implications of our findings in Section 6. The appendices contain details of the data and data analysis (A, B), model validation (C), implementation details (D), and tables of the complete set of retrieval results (E).", '2. The planetary system of HR 8799': 'While HR 8799 is one of the most well studied exoplanetary systems (as seen in Figure 1), there remains significant disagreement in the literature both with respect to the spectroscopic measurements and the inferred planet properties. In general, these super jupiters all host very cloudy atmospheres, with significant impacts of disequilibrium chemistry. The spectra of these planets show characteristics of low surface gravity and have been  \nclassified at the L / T transition, though the variability typical of these objects has not yet been observed in the companions. The composition of the companions is generally found to be moderately enriched compared to the host star, and while H2O and CO are the dominant absorbers, the C / O ratio estimates vary significantly between models. The measurements of the bulk planet properties discussed in this section, together with the results of this work are compiled in Tables 1, 2, 3, and 4 for planets b, c, d, and e respectively.', '2.1. HR 8799 A: The host star': "HR8799 A is an A5V (Gehren 1977; Cannon & Pickering 1993) to F0 + VkA5mA5 (Gray et al. 2003) type star, host to four detected planets and an inner and outer debris disc (Su et al. 2009). It was one of the first identified γ Doradas pulsators (Kaye et al. 1999; Zerbi et al. 1999), and has been classified as a λ Boötis star (Sadakane 2006; Moya et al. 2010) due to its depletion of heavy elements in the atmosphere. Due to this depletion, the [Fe / H] of HR8799 A is measured to be subsolar, with measurements ranging from -0 . 47 ± 0 . 10 (Gray & Kaye 1999; Sadakane 2006) to between -0 . 32 ± 0 . 1 and -0 . 12 ± 0 . 1, depending on the inclination angle (Moya et al. 2010). TESS photometry allowed the measurement of this inclination angle, finding a core rotation period of ≈ 0.7 days, which combined with v sin i and stellar radius measurements would result in a preliminary stellar inclination of ≈ 28° (Sepulveda et al. 2023), and would favour the higher metallicity case presented by Moya et al. (2010). Using high resolution spectroscopy from the LBT / PEPSI and HARPS instruments, Wang et al. (2020) found a very subsolar iron metallicity of -0 . 52 ± 0 . 08, but found the relative carbon (C / H = 0 . 11 ± 0 . 12) and oxygen (O / H = 0 . 12 ± 0 . 14) abundances to be consistent with solar as is characteristic of λ Boötis stars. The derived C / O ratio from their measurements was 0 . 54 + 0 . 12 -0 . 09 . For this work, we use a BT-Nextgen stellar model as fitted in Nasedkin et al. (2023) with best-fit parameters of T e ff = 7200 K, log g = 3.0, and [Fe / H] = 0.0, slightly cooler than the models used in previous studies (Zurlo et al. 2016; Greenbaum et al. 2018). However, this temperature is in line with Sepulveda & Bowler (2022), though based on their dynamical mass estimate of the host star and the radius measurement of Baines et al. (2012) they find a higher surface gravity of 4.28.  \nMost indicators place the HR 8799 system between 25 and 60 Myr in age. Using the debris disc as evidence, Zuckerman & Song (2004) and Rhee et al. (2007) estimate an age of 30 Myr. Zuckerman et al. (2011) identified HR 8799 as 'a likely member of the ∼ 30 Myr old Columba Association', thus providing an age and an association of stars with a shared formation history to which we can compare HR 8799 and its companions. However, using asteroseismology Moya et al. (2010) found that an age of ∼ 1 Gyr or greater is also compatible with their measurements. While we continue to use the standard ∼ 30 Myr age for the system, we acknowledge that there remains some uncertainty in the age and activity of the host star.", '2.2. Photometric studies': "The HR 8799 system has been the subject of extensive photometric characterisation, from the red-optical out to the mid-infrared with JWST / MIRI. The outer three companions were originally detected in Marois et al. (2008), with HR 8799 e following in Marois et al. (2010). Many photometric studies (e.g. Lafrenière et al. 2009; Fukagawa et al. 2009; Metchev et al. 2009; Currie  \nTable 1. Summary of literature and derived planet properties for HR 8799 b.  \nNotes on clouds: both the 'Thick' and 'Slab' clouds are based on Burrows et al. (2006), and are vertically extended throughout the atmosphere above a base pressure, with a decaying mass fraction, though the slab clouds account for a greater range of aerosol opacities. The 'Power law' clouds parameterise the vertical extent and position using a power law, and fix the base pressure to the location at 2300 K. 'AM01' clouds balance the cloud sedimentation and vertical mixing to determine the particle size, and use the sedimentation fraction to determine the vertical extent. 'ER4' is the Exo-Rem4 model from Bonnefoy et al. (2016). 'Mie' clouds do not use physical optical constants, but directly fit for mie scattering parameters. The ' BT-Settl ' cloud model is based on radiation hydrodynamical simulations that solve for the di ff usion and mixing of aerosol particles (Allard et al. 2012).  \nReferences: B11a: Barman et al. (2011); C11: Currie et al. (2011); G11: Galicher et al. (2011); M11: Madhusudhan et al. (2011); M12: Marley et al. (2012); L13: Lee et al. (2013); B15: Barman et al. (2015); B16: Bonnefoy et al. (2016); L17: Lavie et al. (2017); W21: Wang et al. (2021); R21: Ru ffi o et al. (2021).  \n- a Only T int, a model parameter, is reported.\n- b W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7 . 2 ± 0 . 7 M Jup for the inner three planets and 5 . 8 ± 0 . 5 M Jup for HR 8799 b. A radius of 1 . 2 ± 0 . 1 R Jup was used for all planets.\n- c R21 uncertainties were limited by the coarseness and boundaries of their model grid.\n- d Single best retrieval parameters.\n- e Bayesian model averaged parameters from group A ∩ B .\n- f Unweighted average parameters from group A ∩ B .  \net al. 2011; Bergfors et al. 2011; Galicher et al. 2011; Soummer et al. 2011; Skemer et al. 2012; Currie et al. 2012; Esposito et al. 2013; Skemer et al. 2014; Currie et al. 2014; Maire et al. 2015; Rajan et al. 2015; Petit dit de la Roche et al. 2019; Biller et al. 2021; Boccaletti et al. 2023) have identified the companions as L / T transition objects, with near-infrared colours compatible with more extended clouds than L-dwarfs of similar temperatures. This is generally explained as a result of the young age and low surface gravity, where the lower gravity allows the condensate particles to remain aloft above the photosphere at lower temperatures. Even in the earliest studies, disequilibrium chemistry was used as an explanation for the drop in the continuum flux due to CO absorption at 4 µ m (Currie et al. 2011; Janson et al. 2010).  \nIn addition to the four known companions there have been many searches for a fifth companion, interior to HR 8799 e. Thompson et al. (2023) used long time baseline astrometry and deep L' imaging with Keck / NIRC2 to search for this hypothesised companion, finding that an additional companion fits both the astrometry and photometry better than a four planet solution, but does not result in a significant detection. For now, we only examine the four confirmed companions further.  \nBonnefoy et al. (2016) explores the implications of the nearinfrared photometry for all four of the companions, comparing them to spectrally similar field objects from the SpeX PRISM  \nlibrary. Empirically, the HR 8799 planets are much more red in colour than field dwarfs of similar spectral type. They also show that using the dereddening coe ffi cients for corundum (Al2O3), iron (Fe), enstatite (MgSiO3), and forsterite (Mg2SiO4) from Marocco et al. (2014), the colours of the companions more closely match those of field dwarfs. However, they cannot quantitatively distinguish the chemical composition of the clouds, which requires mid-infrared spectroscopic observations of condensate absorption features (Burningham et al. 2021; Miles et al. 2023).  \nMarley et al. (2012) and Bonnefoy et al. (2016) use estimates of the surface gravity and radius from spectroscopic fits to constrain the overall mass and luminosity of the planets, which they in turn compared to planetary evolution models, such as those of Bara ff e et al. (2003) and Saumon & Marley (2008). With selfconsistent, radiative equilibrium models, the planet radius is often di ffi cult to fit, with the radius underestimated by over 30% compared to expectations from the evolutionary models (e.g. Bonnefoy et al. 2016).", '2.2.1. Variability': 'Young brown dwarfs are known to be highly variable (Radigan et al. 2014; Vos et al. 2019, 2023). L / T transition objects display stronger photometric variability - up to 30% - , though  \nTable 2. Summary of literature and derived planet properties for HR 8799 c  \nReferences: C11: Currie et al. (2011); G11: Galicher et al. (2011); M11: Madhusudhan et al. (2011); M12: Marley et al. (2012); L13: Lee et al. (2013); K13: Konopacky et al. (2013); B16: Bonnefoy et al. (2016); L17: Lavie et al. (2017); G18: Greenbaum et al. (2018); W20: Wang et al. (2020); W21: Wang et al. (2021); R21: Ru ffi o et al. (2021); W23:Wang et al. (2023).  \n- a Only T int, a model parameter, is reported.\n- b W20 compared using strong and weak radius priors to enforce physicality.\n- c W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7 . 2 ± 0 . 7 M Jup for the inner three planets and 5 . 8 ± 0 . 5 M Jup for HR 8799 b. A radius of 1 . 2 ± 0 . 1 R Jup was used for all planets.\n- d R21 uncertainties were limited by the coarseness and boundaries of their model grid.\n- e W23 report the temperature at 3.3 bar rather than the e ff ective temperature.\n- f The metallicity of W23 was found by averaging the retrieved C / H and O / H ratios.  \nthis amplitude is rare outside of the transition regime (Radigan 2014). However, as we view the HR 8799 system nearly pole on, it is di ffi cult to see the e ff ects of rotational variation, in addition to the technical challenges of observing variability with high-contrast imaging instruments. Apai et al. (2016) and Biller et al. (2021) have placed upper limits on the photometric variability of the two outermost HR 8799 planets: 10% for b and 25% for c. Wang et al. (2022) used the Subaru / CHARIS instrument to attempt to monitor H-band variability in HR 8799 c and d, placing upper limits of 10% and 30% respectively. The atmospheric turbulence, stellar contamination, and significant postprocessing required to measure the innermost companion has so far prevented measurements of variability for HR 8799 e.', '2.2.2. Orbital dynamics': 'Within the context of directly imaged exoplanets, the HR 8799 companions orbit relatively near to their host star, from a projected separation of 16 au for e out to 71 au for b. Astrometric monitoring has allowed for the precise characterisation of the orbits of the companions, demonstrated in such studies as Wang et al. (2018b); Brandt et al. (2021) and Thompson et al. (2023). Using such orbital fitting techniques, Zurlo et al. (2022) inferred dynamical masses for each of the companions. While they explore a range of models, we use the fit assuming the planets are in a near-resonant 8:4:2:1 configuration and a host star mass of 1.47 M ⊙ . This model produced results typical of the range of models explored; the mass estimates for each companion are: e = 7 . 6 ± 0 . 9 M Jup , d = 9 . 2 ± 0 . 1 M Jup , c = 7 . 7 ± 0 . 7 M Jup , and  \nb = 5 . 8 ± 0 . 4 M Jup . These dynamical mass estimates, as well as those of Brandt et al. (2021), are broadly consistent with mass estimates from evolutionary models, assuming hot start conditions (Marley et al. 2012). Further astrometric analysis of the GRAVITY data will be examined in a forthcoming paper from Chavez et al. (in prep).', '2.3. Spectroscopic characterisation': 'In addition to the multitude of photometric observing campaigns, the spectroscopic characterisation of the HR 8799 planets has traced the development of dedicated exoplanet instrumentation, from long-slit spectrographs (Janson et al. 2010) to high-contrast integral field spectrographs (IFS) (Ingraham et al. 2014; Zurlo et al. 2016) to fibre-fed high resolution spectrometers (Wang et al. 2021). These observations cover a broad swath of wavelength ranges and spectral resolving powers, leading to often conflicting photometric calibration and inferred atmospheric parameters. In particular the H-band spectra as observed with SPHERE (Zurlo et al. 2016), GPI (Greenbaum et al. 2018) and CHARIS (Wang et al. 2022) display di ff erent flux peaks and different H-band shapes. As several atmospheric parameters such as log g and the water abundance are strongly impacted by the shape of this band, they have remained challenging to measure.  \nMany results for individual planets have been presented in the literature. HR 8799 b was first explored in Bowler et al. (2010) with Keck / OSIRIS, where they identify an L5-T2 spectral type, moderate levels of cloudiness and potential impacts of disequilibrium chemistry. Barman et al. (2011) added addi-  \nTable 3. Summary of literature and derived planet properties for HR 8799 d  \nReferences: C11: Currie et al. (2011); G11: Galicher et al. (2011); M11: Madhusudhan et al. (2011); M12: Marley et al. (2012); B16: Bonnefoy et al. (2016); L17: Lavie et al. (2017); G18: Greenbaum et al. (2018); W21: Wang et al. (2021); R21: Ru ffi o et al. (2021).  \n- a Only T int, a model parameter, is reported.\n- b W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7 . 2 ± 0 . 7 M Jup for the inner three planets and  \n.  \n5  \n8  \n±  \n0  \n.  \n5 M  \nJup  \nfor HR 8799 b. A radius of 1  \n.  \n2  \n±  \n0  \n.  \n1  \nR  \nJup  \nwas used for all planets.  \ntional H-band OSIRIS observations, and inferred the low temperature, low surface gravity, and low CH4 abundance of HR 8799 b through the triangular shape of the H-band feature. They also suggest that higher metallicity grids, up to 10 × solar, may be able to better fit the data and provide more plausible radii than their solar metallicity models. This data was augmented with additional wavelength coverage in the K band in Barman et al. (2015), where they claim simultaneous detections of H2O, CO, and tentatively CH4. Oppenheimer et al. (2013) obtained low resolution spectra for all four planets in the Y, J and H bands using the Project 1640 instrumentation suite at the Palomar Hale Telescopes, and thus providing the only additional measurement for HR 8799 b. However these spectra are very low S / N , and are significantly discrepant from subsequent measurements.  \nJanson et al. (2010) were the first to spectroscopically explore HR 8799 c, using the VLT / NACO L-band spectrometer. While they were limited in the available S / N , they still discussed the impact of disequilibrium chemistry on the overall shape of the spectrum, finding that there was strong CO absorption beyond 4 µ m. These L-band observations were later succeeded by LBT / ALES observations (Doelman et al. 2022; Liu et al. 2023), where low resolution spectra at moderate S / N were obtained for the c, d, and e planets. Konopacky et al. (2013) presented the first conclusive evidence of CO and water absorption lines in a directly imaged exoplanet through K-band observations of HR 8799 c using Keck / OSIRIS, measuring the C / O ratio to be slightly above the stellar value at 0 . 65 + 0 . 10 -0 . 05 . The Gemini / GPI instrument provided the first spectra obtained using a coronagraphic instrument in Ingraham et al. (2014), measuring both the c and d planets in the H and K bands. This was followed up with additional post-processing using KLIP (Soummer et al. 2012; Pueyo 2016) in Greenbaum et al. (2018), where HR 8799 e was also detected. All three of the planets were found to best match mid-to-late L-type spectra, with HR 8799 c being most consistent with an L6 dwarf. Consistent with the photometric models, they found HR 8799 c to have a temperature between 1100 K and 1300K, with a log g around 4.0. As with the photometry, the self-consistent models they used to infer the planet  \nproperties struggled to obtain radius estimates consistent with predictions of evolutionary models. Wang et al. (2021) use high resolution spectroscopy to measure the rotation of c, d, and e, finding an upper limit of 14 km / s for c, and measurements of 10 . 1 + 2 . 8 -2 . 7 km / s for d and 15 . 0 + 2 . 3 -2 . 6 km / s for e. Wang et al. (2023) combine several of these datasets and perform a retrieval analysis to constrain the composition of HR 8799 c, finding [C / H] = 0 . 55 + 0 . 36 -0 . 39 , [O / H] = 0 . 47 + 0 . 31 -0 . 32 , and C / O = 0 . 67 + 0 . 12 -0 . 15 . These results depended strongly on the details of the forward model used in the retrieval, and the [C / H] parameter could vary from 0.55 to 0.95, while the [O / H] from 0.47 to 0.80, though they all represent significant enrichment compared to the host star abundances. These results are also significantly discrepant from those of Wang et al. (2020), who found elemental abundance ratios for HR 8799 c of [C / H] = 0 . 16 + 0 . 12 -0 . 13 , [O / H] = 0 . 13 + 0 . 08 -0 . 08 , and C / O = 0 . 58 + 0 . 06 -0 . 06 , though they also found that enforcing strong mass priors led to both the metallicities and C / O ratio being subsolar. Ru ffi o et al. (2019, 2021) and Wang et al. (2021) explore HR 8799 c using moderate and high resolution spectroscopy respectively. Both works characterise the dynamics of the planets, with Ru ffi o et al. (2021) measuring the radial velocities for planets b, c, and d, finding them to be -9 . 1 ± 0 . 4 km / s, -11 . 1 ± 0 . 4 km / s, and -11 . 6 ± 0 . 8 km / s respectively, placing important constraints on the allowed orbits for the planets. They also confirm the presence of water and CO, but are unable to significantly detect CH4, which was consistent with the results of Wang et al. (2018a).  \nThe first reliable spectroscopic measurements of HR 8799 d and e were published by Zurlo et al. (2016). These were obtained using the VLT / SPHERE instrument, and were the first YJH band observations of the inner two planets, and remain the highest quality observations in this band. Together with the modelling work in Bonnefoy et al. (2016), they classify both planets as L6L8 dusty dwarfs, and confirm that only thick cloud models based on the Exo-REM self-consistent modelling code provide reasonable fits to the data, finding e ff ective temperatures of 1200 K, log g in the range of 3.0-4.5, and metallicities of 0.5 for both planets. Compared to previous modelling work of Madhusudhan et al. (2011) and Barman et al. (2011), the Exo-REM models pro-  \nTable 4. Summary of literature and derived planet properties for HR 8799 e  \nB16: Bonnefoy et al. (2016); L17: Lavie et al. (2017); G18: Greenbaum et al. (2018); M20: Mollière et al. (2020); W21:  \nReferences: Wang et al. (2021);  \n- a Only T int, a model parameter, is reported.\n- b W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7 . 2 ± 0 . 7 M Jup for the inner three planets and 5 . 8 ± 0 . 5 M Jup for HR 8799 b. A radius of 1 . 2 ± 0 . 1 R Jup was used for all planets.  \nvided better fits to the data, due to improvements in the opacity databases, cloud treatments, and the inclusion of disequilibrium chemistry. Subsequent SPHERE observations, such as those in Wahhaj et al. (2021) have maintained consistent spectral shapes with these earlier observations. The Gravity Collaboration et al. (2019) performed the first interferometric observations of an exoplanet, measuring the K-band spectrum of HR 8799 e. HR 8799 e was detected as well, and they performed atmospheric analyses on all three planets using the full spectra at 1-2.5 µ m. They found that the spectrum of HR 8799 d has a substantially di ff erent shape than the other two planets, but that all three shared supersolar metallicities and e ff ective temperatures around 1100 K.', '2.4. Retrieval studies': 'Atmospheric retrievals (e.g. Madhusudhan & Seager 2009; Benneke & Seager 2012; Waldmann et al. 2015; Burningham et al. 2017; Mollière et al. 2019) are widely used to solve the inverse atmosphere problem, inferring planet properties such as the thermal structure, chemical composition, and cloud properties from spectroscopic observations. The HR 8799 planets are among the first directly imaged planets to have retrieval methods applied to their spectra. Lee et al. (2013) performed the first retrieval study of HR 8799 b, using the spectrum published of Barman et al. (2011). This pilot study explored various levels of cloudiness, particle sizes, and compositions, finding that the planet is likely cloudy, with relatively large particle sizes (1.5 µ m) and a supersolar metallicity. They note the long-standing degeneracies between the cloud level and the planet radius, making it di ffi cult to distinguish between di ff erent levels of cloudiness in the models. The first systematic characterisation of all four planets was performed in Lavie et al. (2017) using the HELIOS-Retrieval package, with the key goal of constraining the composition of all four planets using the data of Zurlo et al. (2016). After fitting for molecular abundances, they infer the elemental C / H and O / H ratios for each planet, finding oxygen enrichment for b, c, and e, and carbon enrichment for b and c. They find a strongly superstellar C / O ratio for b of 0.9, a stellar value for c, but were unable to constrain the ratio for the inner two planets. While previous works on HR 8799 e were limited due to a lack of high S / N K-band data, Mollière et al. (2020) made use of the GRAVITY spectrum obtained in Gravity Collaboration et al. (2019), to-  \nher with the SPHERE data of Zurlo et al. (2016) and the GPI data of Greenbaum et al. (2018). Using the pRT retrieval framework and a novel temperature profile, they inferred a highly cloudy atmosphere, implementing clouds with multiple scattering. They infer modest enrichement of [M / H] = 0 . 48 + 0 . 25 -0 . 29 and a C / O ratio of 0 . 60 + 0 . 07 -0 . 08 . Finally, Wang et al. (2020) and Wang et al. (2023) both perform pRT retrieval studies of HR 8799 c. The latter study is unique in including high resolution data in the retrieval framework, allowing precise measurements of the elemental abundance ratios, finding modest enrichment of both carbon and oxygen.', '2.5. Self-consistent atmospheric modelling': 'Motivated by the considerable volume of observations, extensive theoretical modelling work has been performed to better understand the physics of the atmospheres of the HR 8799 planets and similar substellar objects. Brown dwarf atmospheres saw extensive 1D modelling e ff orts (e.g. Chabrier et al. 2000; Allard et al. 2001; Burrows et al. 2006; Saumon & Marley 2008), driven largely by the need to trace the evolution of these continuously cooling objects over time. These studies demonstrated the necessity of accounting for silicate clouds in the atmospheres of L / T dwarfs, used to explain the red colour of these objects in the near infrared. Applied specifically to the young, low-gravity companions, Madhusudhan et al. (2011) developed one of the first models specifically for the HR 8799 companions to constrain their mass and age. They identify forsterite and iron as being the important contributors to the clouds, and infer planetary ages between 10 and 150 Myr, consistent with stellar measurements. Marley et al. (2012) provides a deep review of the state of modelling of the atmospheres of the HR 8799 planets, further developing the model of Saumon & Marley (2008). They find masses and ages for the planets consistent with the stellar properties, and that the companions share approximately consistent properties with L / T dwarfs of similar e ff ective temperatures and surface gravities. Using the cloud model of Ackerman & Marley (2001, hereafter AM01), they infer an f sed parameter of 2, implying that the clouds are moderately extended throughout the atmosphere. More recent self consistent models such as petitCODE (Mollière et al. 2015), ATMO (Tremblin et al. 2015; Phillips et al. 2020; Petrus et al. 2023), Exo-REM (Charnay et al.  \n2018), and Sonora (Marley et al. (2021); Karalidi et al. (2021); Morley et al. in prep) have been developed specifically to understand the thermal structure and clouds of directly imaged planets. There remain degeneracies between reddening and damping of spectral features via continuum opacity sources and through reductions in the temperature gradient, hypothesised to be due to diabatic convection (Tremblin et al. 2019).  \nZahnle & Marley (2014) provide an in-depth exploration of the impacts of disequilibrium chemistry on cool, self-luminous atmospheres, providing predictions for the CO, CH4, and NH3 abundances as a function of vertical mixing and e ff ective temperature, identifying the key transition between CH4 and CO dominated chemistry at around 1100 K. Moses et al. (2016) uses a disequilibrium model including photochemistry to predict the chemical composition for a range of surface gravities and e ff ective temperatures, and provides column abundance predictions for HR 8799 b, finding that the CO abundance should dominate over CH4, assuming a solar composition. Soni & Acharyya (2023) extend this to superstellar metallicities and vertical mixing strengths, using the constraints on the CO and CH4 abundances from Barman et al. (2015) to infer a vertical mixing strength of log K zz ∈ [7 , 10] for the 10 × solar metallicity case. To better understand the planet structure, Thorngren et al. (2016) derive a mass-metallicity relationship. As the mass of the object increases, the metallicity tends to decrease, consistent with predictions of core accretion formation, as heavier objects accrete and retain more H2 and He relative to a lower mass object. From their relationship, they predict that a 6 M Jup planet should have a Z pl / Z ⋆ ratio of between 3 to 5 (in a 68% confidence interval).  \nIn addition to the 1D modelling e ff orts, global circulation models (GCMs) of self-luminous, substellar objects, such as those of Showman & Kaspi (2013); Tan & Showman (2021a,b) have been developed. These 3D models allow for the exploration of atmospheric dynamics, longitudinal variations, and time variability. Recent observations are beginning to validate these 3D models: Suárez et al. (2023) finds that brown dwarfs are cloudier when viewing the equator, which is consistent with the cloudiness predictions of rapidly rotating brown dwarfs in Tan & Showman (2021b). Likewise, the prediction of patchy clouds in the photosphere region leading to variability (Showman & Kaspi 2013) seems to match the observations of high variability in lowgravity atmospheres with silicate clouds (Vos et al. 2023).', '2.6. Formation': "With four massive planets on wide orbits, HR 8799 provides a unique system with which to test formation scenarios. In general, these fall under the categories of either gravitational instability models (e.g. Perri & Cameron 1974; Cameron 1978; Adams et al. 1989; Laughlin & Rozyczka 1996; Boss 1997), where the planets form via the direct collapse of the gas into a substellar object, or core accretion (Pollack et al. 1996; Bodenheimer et al. 2000), where a dense core of heavy material grows slowly until it is massive enough to experience runaway accretion and gather an extended hydrogen-helium envelope. GI models tend to produce larger planets on wider orbits with solar compositions, while core accretion scenarios form closer-in planets on more circular orbits, with the possibility of greater metal enrichment. Dodson-Robinson et al. (2009) tested both of these scenarios for HR 8799, finding that while core accretion may better explain the near-orbital resonances of the system, it struggled to form planets on such wide orbits (beyond 30 au), and could not rule out the possibility of direct gravitational collapse. Similarly, Nero & Bjorkman (2009) find that while HR 8799 b may have  \nformed through gravitational instability, it is unlikely that disc fragmentation could have formed the inner three companions.  \nIn addition to constraints from the mass and location of the companions, the present-day planet composition provides insight into the formation and evolution history. The template for this was developed in Öberg et al. (2011), demonstrating how the C / O ratio in the gas and dust varies as a function of position in the disc, which would in turn impact the outcome of the formation process. Eistrup et al. (2018) extended this model to include time evolution, and Mollière et al. (2022) presented a framework to link the measured planet properties to the disc environment in a Bayesian framework, which allows testing the e ff ect of various formation assumptions. However, due to the uncertainty in the atmospheric measurements, combined with the many outstanding questions in formation modelling, this link remains tenuous.  \nThe di ff erent formation scenarios can lead to dramatically di ff erent amounts of energy retained in the planet following the formation process. So-called 'hot-start' models result in young planets retaining the gravitational potential energy as internal heat, to be radiated and cooled over time (Marley et al. 2007; Mordasini et al. 2017). This scenario is typically associated with formation due to gravitational instability. In cold-start scenarios, often tied to core-accretion models, this energy is radiated away by accretion shocks as the gas flows from the circumstellar disc onto the forming planet, resulting in a lower internal energy (Mordasini et al. 2012; Szulágyi & Mordasini 2017). This is a useful, though simplified picture of planet formation. Additional complication comes from the energetics of the accretion shock during core accretion, where di ff erent radiative e ffi ciencies can lead to di ff erent initial entropies of the forming planet (Marleau et al. 2017). These shock-resolving models find typical internal energies that are an order of magnitude higher than in typical cold-start scenarios (Marleau et al. 2019), thus lying somewhere between the hot and cold start scenarios. Over time, all of these scenarios converge to the same cooling rate, though precise mass and luminosity estimates can distinguish between the two scenarios for the first ∼ 100 Myr (Bara ff e et al. 2003; Saumon & Marley 2008). Current measurements of planet masses, temperatures, and radii generally favour hot or warm start models, but can only definitely exclude the coldest initial conditions, such as the cold-start models of Marley et al. (2007). The hot-start models of Bara ff e et al. (2003) led to predictions of 7 M Jup for the inner three planets, and 5 M Jup for HR 8799 b, which are approximately consistent with the current dynamical mass estimates of Zurlo et al. (2022). Using the hot-start model of Saumon & Marley (2008), Marley et al. (2012) finds that the radii of all of the planets should be slightly larger than 1 R jup, and that even assuming very cold initial conditions the planet radii should never fall below 1 R jup, though this claim did not account for significantly nonsolar composition.  \nFinally, HR 8799 is home to both an inner and outer debris disc, imaged with Spitzer (Su et al. 2009), Herschel (Matthews et al. 2014), JWST (Boccaletti et al. 2023), and ALMA in the millimeter (Wilner et al. 2018). The inner debris disc has a temperature of around 150 K and is confined to within 10 au, while the cold outer debris disc is analagous to the Kuiper belt in our own Solar System (Geiler et al. 2018), but at a much wider separation (90-300 au). The structure of the outer disc appears to be sculpted by an additional gravitational component, though it is unclear whether this is due to HR 8799 b or an additional unseen companion (Contro et al. 2015; Faramaz et al. 2021). The inner disc has been detected in thermal emission (Su et al. 2009) and resolved using MIRI coronagraphic imaging (Boccaletti et al. 2023). Modelling e ff orts have placed tentative lim-  \nFig. 2. Flux-calibrated VLTI / GRAVITY K-band spectra for each of the HR8799 planets, normalised to 10 pc. Each spectrum above that of HR 8799 b has an additional 1.5 × 10 -15 W / m 2 / µ m o ff set. The panels on the right show the empirical correlation matrices for each of the four planets. The colour bar is scaled to highlight weak correlations in the GRAVITY data.  \n<!-- image -->  \nits of ∼ 1 M Jup on the allowed mass of companions interior to HR 8799 e (Go'zdziewski & Migaszewski 2018).", '3. Observations': 'While the new GRAVITY spectra represent the best available K-band observations of the HR 8799 system, additional data are required to constrain planetary properties such as surface gravity and C / O ratio. We combine published datasets across a wide wavelength range from a variety of sources in order to present the most complete possible picture of this system. Archival photometric data of the companions are also included in our analysis, the details of which can be found in appendix A. Also included in appendix A is the stellar photometry used in fitting the BT-Nextgen model, with which the companion contrast measurements are flux calibrated. In this section we present a brief overview of the spectroscopic datasets included in the retrieval analysis, with the key observational parameters listed in Table A.3. All of the observational data, together with the complete set of retrievals results is available on Zenodo 1 .', '3.1. GRAVITY data': "In Figure 2 we present new VLTI / GRAVITY observations of HR 8799 e, together with the first interferometric observations of d, c, and b taken as part of the ExoGRAVITY project (Lacour et al. 2020), under ESO program ID 1104.C-0651. GRAVITY is a K-band spectroscopic interferometer that combines light from either the four 8 m Unit Telescopes (UTs) of the VLT, or the 1.8-m Auxiliary Telescopes (Gravity Collaboration et al. 2017). With baselines of up to 134 m, GRAVITY provides unprecedented spatial resolution, allowing for the detection of companions close to their host stars and the measurement of relative astrometry with a precision of few tens of µ as. All observations of the HR8799 system were obtained using the UTs, with the dual-field mode of GRAVITY. The medium resolution mode was used, which o ff ers a resolution of R ∼ 500 over a nominal wavelength range of 2.0 to 2.4 µ m.  \nTwo di ff erent strategies were used for the observations and data-reduction. Observations of HR 8799 e at all dates, except on 2 dates (11 November 2019 and 02 July 2023) were obtained  \nusing the on-axis strategy, in which a 50 / 50 beam-splitter is used to separate the field to between the science and fringe-tracking channels of GRAVITY. In this mode, observations with the science channel pointing at the location of the planet are interleaved with observations obtained with the fibre pointed at the central star. The on-star observations are then used to calibrate both the interferometric phase and amplitudes. This is similar to the observations reported in Nowak et al. (2020). The second strategy is the dual-field / o ff -axis strategy, in which the roof-mirror is used to split the field. The use of the roof-mirror is required to observe planets at larger separation, because the field of view of the beam-splitter does not reach these targets. In this case, the metrology zero point is calibrated using observations of the dedicated calibrator HD 25535, and the interferometric amplitude using an on-axis observation of the central star, typically done at the end of the observation sequence. This strategy is similar to the observation of Sgr A* by GRAVITY Collaboration et al. (2020).  \nThe data-reduction was performed using the tools developed for the ExoGRAVITY large program 2 . The main steps of the reductions are as follows:  \n- 1. All data are first reduced with the GRAVITY pipeline (Lapeyrère et al. 2014), up to the 'astroreduced' data product, which keeps individual DITs separated.\n- 2. For the on-axis observations, the phase reference is extracted from the on-star observations and subtracted from the onplanet observations. For the o ff -axis observations, this phasereference is extracted from the observations of the binarycalibrator HD 25535. In both cases, the amplitude reference is taken using the on-star observations.\n- 3. The stellar light (also called stellar speckle) is subtracted from the reduced data by fitting a fourth-order polynomial in wavelength multiplied by the amplitude reference. The astrometry of the planet is then extracted from the observations.\n- 4. The contrast spectrum is then extracted using a model that also takes into account the residual starlight, and the planet astrometry previously extracted.  \nThis procedure yields a contrast spectrum for each planet, at each observation date. The spectrum extraction, which consists  \nFig. 3. Data of the HR 8799 planets. The OSIRIS b and ALES d datasets have scaling factors of 1.0 and 1.2 applied, respectively. Not shown are the MIRI photometry points from Boccaletti et al. (2023). References: OSIRIS (Barman et al. 2011); SPHERE (Zurlo et al. 2016); GPI (Greenbaum et al. 2018); CHARIS (Wang et al. 2020, 2022); ALES (Doelman et al. 2022).  \n<!-- image -->  \nentirely of linear operations on the complex coherent flux, also propagates the errors reported by the pipeline as covariance matrices. These covariance matrices allow for correlations over the wavelength channels and between the real and imaginary parts of the coherent flux. However, it should be noted that the GRAVITY pipeline does not report such covariances, and so the extraction code starts with fully diagonal covariance matrices.  \nFor each planet, all the available spectra are then combined using a covariance-weighted combination. The final contrast spectrum C = ( c λ 1 , c λ 2 , . . . , c λ n ) T and its associated covariance matrix W arre given by:  \nW =        X t Wt -1        -1 (1) C = W ·        X t Wt -1 · Ct        , (2)  \nwhere Ct and Wt represent the contrast spectrum and its associated covariance matrix on a given observing date t .  \nThe contrast spectra are then converted to fluxes using a model of the stellar flux. For HR 8799, we used a BT-Nextgen model fit to the near infrared photometry, the details of which are more thoroughly discussed in Section 2.1 and are based on Nasedkin et al. (2023).  \nThe faintest companion, HR 8799 b, was detected with a mean S / N of 3.4 per wavelength channel. HR 8799 c was observed with a mean S / N of 27.5 per channel, while HR 8799 d and HR 8799 e had a mean S / N ≈ 20 and S / N ≈ 10 respectively. These observations were taken over a 5 year period. With the 50 microarcsecond astrometric precision of GRAVITY, this will allow the detection of planet-planet orbital perturbances within a few years (Covarrubias et al. 2022), and we leave such analysis to future work.", '3.2. Archival data': 'In addition to the new GRAVITY spectra, we also include archival data covering a broad wavelength range, presented in Fig. 3. Mollière et al. (2020) noted that the SPHERE (Zurlo et al. 2016) and GPI (Greenbaum et al. 2018) data are inconsistent with each other in the H-band. In order to reduce systematic variation and to account for correlations, we rereduce the data with up-to-date pipelines, and reprocesses the datasets optimally as described in Nasedkin et al. (2023) using KLIP (Soummer et al. 2012; Pueyo 2016). However, in order to best extract the planet signal we use KLIP in ADI + SDI mode, in comparison to ADI only mode as described in the previous study. Both the reprocessed SPHERE and GPI spectra can be found in Figs. B.1 and B.2. In total, our dataset includes nearly 400 data points for each planet: Nb = 297 , Nc = 391 , Nd = 387 , Ne = 388.', '3.2.1. SPHERE': 'Two sets of VLT / SPHERE (Beuzit et al. 2008, 2019) data are considered in this study: the first was taken during the commissioning run of the SPHERE instrument on 12 August 2014, and was originally published in Zurlo et al. (2016). This is still the deepest SPHERE observation of HR 8799 covering the full YJH range, but due to the orientation of the field of view does not include HR 8799 c. This dataset was reprocessed as in Nasedkin et al. (2023) using KLIP in ADI + SDI mode, and we extract spectra and covariance matrices for both the e and the d companions. The second SPHERE dataset was published in Flasseur et al. (2020), who processed the dataset using the PACO-ASDI algorithm and were able to extract a spectrum for HR 8799 c in addition to d and e.  \nAdditional SPHERE observations, such as presented in Biller et al. (2021) or Wahhaj et al. (2021) are available. However, in the case of Biller et al. (2021) the observations of the host star used for photometric calibration that were taken before and after the science observations are of insu ffi cient S / N . While we attempted to calibrate the companion spectra using  \nthe satellite spots, this was unreliable. Finally, these observations only cover the Y and J bands, and lack the overlap with the GPI H-band spectrum, which is important for ensuring compatibility across instruments. Therefore we continue with only the datasets of Zurlo et al. (2016) and Flasseur et al. (2020).', '3.2.2. GPI': 'Gemini / GPI (Macintosh et al. 2014) observations of HR8799, originally published in Greenbaum et al. (2018), were taken on 17 November 2013, 18 November 2013, and 19 September 2016 for the K1, K2 and H bands respectively. These were reduced using the standard GPI reduction pipeline (version 1.4.0), and reprocessed with KLIP using the same methods as the SPHERE data. As the new GRAVITY observations supersede the GPI data in the K-band, we only consider the GPI H-band data for this work.', '3.2.3. CHARIS': 'Subaru / CHARIS (Gro ff et al. 2015, 2017) observations of HR 8799 c, d, and e were presented in Wang et al. (2020) and Wang et al. (2022). These observations cover 1.2-2 . 4 µ m range at low resolution. Wang et al. (2022) primarily examined these data for temporal variability, while here we combine the full two nights of observations in order to obtain the highest precision spectrum for each of the three planets. We take the mean spectrum for both nights, and add the errors in quadrature, dividing by the square root of the number of observations (i.e. by √ 2) to obtain a spectrum for each planet.', '3.2.4. ALES': 'Doelman et al. (2022) presented L-band observations of HR 8799 c, d, and e obtained using the LBT / ALES instrument (Skemer et al. 2015). These supersede the VLT / NACO L-band observations of HR 8799 c of Janson et al. (2010), and are the first L-band spectra of HR 8799 d and e. These data also include covariance matrices, estimated using the analytic method of Greco & Brandt (2016).', '3.2.5. OSIRIS': 'Archival Keck / OSIRIS (Larkin et al. 2006) data taken between 2009 and 2010 is included for HR8799b, as published in (Barman et al. 2011). HR 8799 b falls outside the field of view of most high-contrast-imaging IFUs, so OSIRIS is joined only by GRAVITYin measuring the near infrared spectrum of the planet. With an unbinned spectral resolution of R ≈ 4000, and an integration time of 2700s in the H-band and 1800s in the K-band, when binned to a spectral resolution of R ≈ 60 the OSIRIS data achieves a per-channel S / N comparable to or better than that of the GRAVITY observations. As the OSIRIS data were not taken using standard ADI observing modes, we did not attempt any rereduction or reprocessing of the archival data, apart from rescaling the flux and uncertainty by the current Gaia distance estimate of 41 . 2925 pm 0 . 15 pc (Gaia Collaboration 2020).  \nWe also include the K-band spectra of Konopacky et al. (2013) (Figure 2 of that work). As published, this spectrum is not flux calibrated, and so we always fit for a flux-scaling term. For HR 8799 b and c, these K-band spectra allow us to explore the impact of di ff erent measurements on the retrieved atmospheric  \nTable 5. Retrieval prior; temperature profile priors are included in Table 6. N ( µ, σ ), U (low, high).  \nparameters, and to determine if our methods can reproduce the results of earlier work.  \nMore recent observations of the HR 8799 planets using OSIRIS have been explored in Ru ffi o et al. (2021), but these spectra are continuum subtracted, requiring a somewhat di ff erent modelling framework than the rest of the data considered in this work. As such we do not fit these data, but we do use them as an additional check on the robustness of our fits when examining the best-fit models at higher resolution.', '4. Atmospheric modelling': 'The forward models of our atmospheric retrieval setup were computed using pRT version 2.7 (Mollière et al. 2019), a fast, open-source radiative transfer code with which we calculate the emission spectrum of a planetary atmosphere 3 . Our fiducial setup was based on that of Mollière et al. (2020), used to retrieve the atmospheric properties of HR 8799 e, though substantial improvements to the code have been made and are detailed further in Nasedkin et al. (2024). We explore a wide range of model parameterisations, summarising the parameters and prior distributions used in Table 5. As we consider several thermal profile parameterisation, we compare their prior distributions separately in Table 6.  \nTo allow for both a data-driven and physically motivated approach, we retrieved either log g and R pl with uniform priors or R pl and M pl, with Gaussian priors set by the dynamical mass estimates of Zurlo et al. (2022) and broad Gaussian priors centred at 1.1 R Jup, in line with estimates from evolutionary models (Marley et al. 2012).  \nAs the computational cost of a retrieval varied greatly between the planets, it was unfeasible to run every model for every planet. As our primary point of comparison we explored the different temperature profile parameterisations for each planet, and ran both disequilibrium and free chemistry retrievals for each planet. Due to its low computational run time, we ran additional models for HR 8799 e, focusing on di ff erent cloud parameterisations.', '4.1. Thermal structure': 'ets:  \nWe compared a set of four temperature structures in our model comparison in order to distinguish the amount of model flexibility justified by the data and the impact of the temperature structure on other retrieved atmospheric parameters. While the thermal structure of these self-luminous objects is thought to be well-understood from 1D and 3D atmospheric models, this comparison will validate these predictions using an independent, data-driven methodology. At the same time, using the best physical understanding of the thermal structure may help constrain other parameters with greater accuracy and precision; it is necessary to compare both approaches to ensure consistent results. For all di ff erent temperature profiles we computed an e ff ective temperature after the spectrum computation, by integrating F λ over wavelength and applying the Stephan-Boltzmann law. To do this, we integrated a low resolution spectrum from 0.8 to 250 µ m. The lower limit is set by the wavelength coverage of the data; the optical band is unconstrained and leads to unrealistically large uncertainty on the e ff ective temperature. The long wavelength limit is set by the wavelength coverage of the opacity databases.', '4.1.1. Spline profile': 'To allow the data to fully determine the temperature profile of the atmosphere, we used a Piecewise Cubic Hermite Interpolating Polynomial as implemented in the scipy.interpolate. PchipInterpolator function. Following the prescription of Line et al. (2015), we penalised curvature in the temperature profile by adding an additional term to the likelihood function,  \nlog p ( T ) = 1 2 γ N -1 X i = 1 ( Ti + 1 -2 T + Ti -1) 2 -log (2 πγ ) . (3)  \nThis is the additional penalty term, which we found by taking the sum of the discrete second derivative of the temperature profile. An additional hyperparameter, γ , was also included, with an inverse gamma distribution prior. If γ is large, (disfavoured by the prior), then the data truly demands strong curvature in the profile, while if γ is small, the data favours smoother profiles. Following Line et al. (2015), we set the parameters of the prior distribution on γ based on the work of Lang & Brezger (2004); Rahman (2005) and Jullion & Lambert (2007):  \nΓ -1 ( γ ) = β α Γ ( α ) 1 γ ! α + 1 exp -β γ ! , (4)  \nfor fixed α and β parameters given in Table 6. We repeated the retrievals and varied the number of nodes in the profile, which allowed us to use a Bayes factor comparison to determine the allowable level of complexity. This also allowed us to explore how the pressure-temperature profile can compensate for the presence of clouds by reducing the temperature gradient in the photospheric region.', '4.1.2. Guillot profile': 'The Guillot (2010) (G10) profile is a simple analytical model, constructed to estimate the thermal structure of irradiated plan-  \nT 4 Guillot = 3 T 4 int 4 2 3 + τ ! + 3 T 4 irr 4       2 3 + 1 γ √ 3 +       γ √ 3 -1 γ √ 3       e -γτ √ 3       , (5)  \n√  \nwhere T irr = 2 T equ and τ = P × κ IR / g . T equ is the equilibrium temperature of an irradiated body, T int is the intrinsic temperature of the planet, and g is the surface gravity. κ IR is the mean infrared opacity, and γ is the ratio between the optical and infrared opacities. While these parameters can be physically interpreted, we treat them as nuisance parameters that control the shape of the profile, rather than self-consistently linking them to the chemical opacities. As the HR 8799 planets are widely separated, T irr is small, and thus the profile reduces to an Eddington (1930) profile, which corresponds to keeping only the first term on the righthand side of Equation (5).', '4.1.3. Mollière profile': 'Introduced in Mollière et al. (2020) (M20), this is a physically motivated temperature profile, split into three distinct regions in altitude.  \nThe middle level of the atmosphere contains the photosphere. In this region the temperature profile follows an Eddington profile, as in the first term of the Guillot profile in Equation 5. However, for this profile we parameterise the opacity τ as a function of pressure ( P ):  \nτ = δ P α (6)  \nand retrieve parameters of log δ and α , together with T int, as in the G10 profile.  \nThe upper atmosphere is defined as the region above τ = 0 . 1. Above this level, four pressure points are defined, equidistant in log P . The deepest pressure point, at τ = 0 . 1 is fixed to the temperature of the Eddington profile of the middle atmospheric region, while the remaining temperature points are freely retrieved parameters, subject to the constraint that the temperature decreases with altitude (Kitzmann et al. 2020), as inversions are not expected in self-luminous objects. The temperature profile is then interpolated from a cubic spline between the three points. Combined with the Eddington profile parameters, this results in a total of 6 parameters to describe the temperature profile.  \nThe base of the atmosphere is defined as a moist adiabat, up to the radiative-convective boundary. This boundary occurs when the temperature gradient of the Eddington profile is Schwarzchild unstable:  \ndT dr < T P dP dr 1 -1 γ ad ! . (7)  \nThe moist adiabatic gradient is a function of the temperature, pressure, and chemical composition, and as such is interpolated from the disequilibrium chemistry table, discussed further in Section 4.2.2. Once the atmosphere is unstable to convection, the temperature profile is forced onto the moist adiabat.', '4.1.4. Zhang profile': 'Zhang et al. (2023) (Z23) introduced a novel P-T parameterisation, incorporating the results of radiative-convective equilibrium models into the retrievals via careful prior selection. This is  \nTable 6. Priors for temperature profiles.', 'Notes': 'References: B11: Barman et al. (2011); K13: Konopacky et al. (2013); Z16: Zurlo et al. (2016); G18: Greenbaum et al. (2018);  \nG19: Gravity Collaboration et al. (2019); F20: Flasseur et al. (2020); M20: Mollière et al. (2020); W22: Wang et al. (2022); D22: Doelman et al. (2022).  \net al. (2024), we verified that the results of Mollière et al. (2020) could be reproduced. We independently tested updates to the c-k mixing implementation, the adaptive mesh refinement implementation, updated opacity sources for H2O and CO, bug fixes for convergence on multiple scattering in the clouds, the inclusion of photometric data, the inclusion of scaling factors on the SPHERE and GPI datasets, including or excluding the GPI K-band spectra, updates to each of the SPHERE, GPI, and GRAVITYdatasets, di ff erent prior widths and the number of live points used in the retrieval. All of the posterior distributions were fully consistent to within 2 σ with most falling well within 1 σ of the published results, apart from the inclusion of new epochs  \nof GRAVITY data, which led to a significantly higher retrieved metallicity ([M / H] = 1 . 1 ± 0 . 32) and an f sed of 5 ± 2 . 6.  \nWe verified several model assumptions through retrievals that only include the GRAVITY datasets, or the GRAVITY data and photometry. We find that the GRAVITY data alone could not distinguish between clear and cloudy models ( ∆ log 10 Z < 1), while cloudy models were strongly favoured once the broad wavelength coverage of the photometry was included ( ∆ log 10 Z > 10).  \nModels that use the dynamical mass estimates as priors for calculating the surface gravity were marginally favoured over those that freely retrieve log g and R pl, but that the posteriors parameter distributions were generally consistent, with T int and R pl  \nFig. B.2. KLIP and ANDROMEDA extractions from GPI for HR 8799 c, d, and e compared to the spectra published in Greenbaum et al. (2018)  \n<!-- image -->  \nshowing the greatest discrepancy. Using the dynamical mass as a prior and setting a Gaussian prior on the radius led to more reasonable estimates of the radius of HR 8799 e (0 . 97 ± 0 . 04) compared to the free retrieval (0 . 79 ± 0 . 05). However, the composition of the planet and the degree of cloudiness did not vary significantly between the two models.  \nUsing the full dataset for HR 8799 b, we verified that retrievals including scattering clouds are strongly favoured over those without scattering ( ∆ log 10 Z > 10). Without scattering, both the temperature and composition ([M / H] and C / O) are significantly discrepant from retrievals that include scattering clouds.', '4.2. Chemistry': 'Understanding the atmospheric chemistry of the HR8799 planets is one of the key goals of this work. We compared a simplified disequilibrium chemistry model to a free chemistry retrieval with vertically constant abundances. We primarily used opacities from the ExoMol database (Tennyson & Yurchenko 2012; Chubb et al. 2021), and included H2O (Polyansky et al. 2018), CO (Rothman et al. 2010), CH4 (Yurchenko et al. 2017), CO2 (Yurchenko et al. 2020), NH3 (Coles et al. 2019), HCN (Barber et al. 2014), H2S (Azzam et al. 2016), PH3 (Sousa-Silva et al. 2015), FeH (Wende et al. 2010), Na (Allard et al. 2019), K (Allard et al. 2016), SiO (Barton et al. 2013), TiO (McKemmish et al. 2019), and VO (McKemmish et al. 2016).', '4.2.1. Free chemistry': "In the free chemistry retrievals, we assumed a vertically constant mass fraction for each species, and retrieved the log mass fraction abundance (log Xi for each of H2O, CO, CH4, CO2, HCN, H2S, NH3, FeH, Na, and K), subject to the constraint that the sum of the mass fractions is less than one. Due to the lack of spectroscopic data in in the Y and J bands, we did not retrieve FeH, Na or K for HR 8799 b in the free retrievals to reduce the number of free parameters. The hydrogen and helium mass fractions were calculated by using the solar abundances (0.766 for H2, 0.234 for He), and multiplying by one minus the sum of the  \nretrieved mass fractions (1 -Σ i Xi ). The set of molecules included covers the most abundant trace species in the atmosphere, and in an atmosphere with strong vertical mixing in the photosphere the assumption of a vertically constant abundance is reasonable for H2O, CO, and CH4, though other species, such as FeH, have been found to have nonvertically constant abundances (Rowland et al. 2023). To measure the significance of a detection, we performed a 'leave-one-out' retrieval, and calculated the Bayes factor between the complete retrieval and the retrieval when excluding a single chemical species. This comparison was performed using the Zhang temperature profile and clouds condensing at their equilibrium saturation condition. While the detection significant may vary with di ff erent setups, this setup is representative of typical retrievals, and ensures consistent comparisons.  \nTo determine the bulk properties of [M / H] and C / O for the free retrievals, we converted the retrieved mass fraction abundances to volume mixing ratios. [M / H] is defined as the ratio the planetary elemental abundances from the measurements of H2O, CO, CH4, CO2, NH3, and H2S to the solar elemental abundances. The C / O ratio was likewise found from the number ratio of the carbon and oxygen atoms in the same set of molecular species.", '4.2.2. Interpolated (dis)equilibrium': 'The disequilibrium model used a grid of equilibrium chemical abundances interpolated along dimensions of pressure and temperature, as well as [M / H] and C / O, which were freely retrieved parameters. The metallicity parameter scaled all of the elemental abundances, while the C / O scaled only the oxygen abundance. Initial test retrievals used a prior range of [-1.5,1.5], but the high metallicity demanded by the data led to the choice of a prior range of [-0.5,2.5] for the retrievals included in this work. The model of disequilibrium chemistry of the CO, CH4, H2O system was based on transport-induced quenching, resulting in a vertically constant abundance above a given pressure. This (log) quench pressure was one of the retrieved parameters. The equilibrium abundances used to build the grid were computed using easyChem (Mollière et al. 2017), which minimises the Gibbs free energy for the system at a given pressure, temperature, and atomic composition. We included all of the species listed in Section 4.2 as opacity sources, though 95 species are included in equilibrium chemical network used to determine the molecular abundances. For the alkali metals we used the wing profiles of Allard et al. (2016, 2019).', '4.4. Retrieval setup': 'We considered three cloud parameterisations in this work. The first was based of the model of Ackerman & Marley (2001), where cloud particles are lofted into the atmosphere through eddy di ff usion ( K zz), and settle back down with a speed proportional to the parameter f sed. We retrieved each of these parameters independently, together with σ LN, which is the width of the log-normal particle size distribution.  \nWhile these parameters determine the structure of the clouds, we also determined the cloud opacity through the use of cloud optical constants for di ff erent cloud compositions, allowing us to di ff erentiate between compositions, grain structure (amorphous or crystalline), and whether the particle shapes are spherical or based on the distribution of hollow spheres (DHS). In addition to the standard log-normal particle sized distribution used with AM01 clouds, we incorporated the Hansen (1971) particle size distribution, which has been proposed to be a more accurate representation of the particle size than a log-normal distribution. Instead of σ ln, we retrieved the mean e ff ective width b h. The details of how these parameters shape the distribution, together with how they were incorporated into the AM01 model are described in Appendix D.  \nWe tested a range of cloud compositions, including MgSiO3 (both crystalline and amorphous particle shapes), Mg2SiO4, Fe, Al2O3, KCl, and Na2S, as well as several combinations of these compositions. Each cloud composition had a mass fraction abundance at the cloud base that could either be scaled from an equilibrium value or freely retrieved, together with a unique f sed, which allowed for di ff erent vertical extents for di ff erent cloud compositions.  \nNominally, the cloud base occurs at the intersection of the cloud condensation curve and the temperature pressure profile. However, for our second parameterisation we included the cloud base pressure as a freely retrieved parameter, to determine if the clouds form where expected in the atmosphere. We followed the derivation of AM01 and Mollière et al. (2017) to obtain the cloud abundance throughout the atmosphere. The abundance of the cloud species Xi was defined at this base pressure P 0, and decreases with altitude to the power of f sed:  \nXi ( P ) = Xi , 0 P P 0 ! f sed . (10)  \nAt pressures higher than P 0, the cloud is not condensed, and thus Xi = 0.  \nFinally, we also tested a simple grey cloud deck, where the cloud top pressure was freely retrieved, and the cloud acts as a source of opacity at the base of the photosphere.  \nFor each of these models we could retrieve a cloud patchiness fraction, f c. In this setup, we first calculated the usual cloudy spectrum, S cd. We then turned o ff the cloud opacity sources, and calculated another clear atmosphere spectrum, S cl. We then combined the two spectra, weighted by f c:  \nS = f c S cd + (1 -f c) S cl . (11)  \nThis approach divides the atmosphere into only clear and cloudy components. Other approaches, such as those of Vos et al. (2023) or McCarthy et al. (2024) have di ff erent patchiness fractions for di ff erent cloud layers, allowing for di ff erent degrees of cloudiness. Our approach reduces the number of parameters and is simple to implement in a retrieval framework, but future work should explore the patchiness of individual layers of clouds in the atmosphere.  \nWe used the pyMultiNest (Buchner et al. 2014) wrapper of MultiNest (Feroz & Hobson 2008; Feroz et al. 2019) as the basis for our nested sampling routine, as testing showed that it runs significantly faster than the UltraNest sampler, which may provide more accurate estimates of the Bayesian evidence. For all retrievals we used 4000 live points to ensure dense posterior sampling and coverage of the parameter space. We set the sampling e ffi ciency to 0.05, and used constant e ffi ciency mode in order to reduce computation time. Comparisons to retrievals using 4000 live points and a sampling e ffi ciency of 0.8 without constant e ffi ciency showed that this choice does not bias the posterior estimates, and that the importance nested sampling evidence estimate is of su ffi cient precision for model comparison. The (log Z ) evidence tolerance was set to 0.1, ensuring precise estimates of the evidence and ensuring convergence of the retrievals.', '4.5. Retrieval ranking': "Considering the number and range of models run, we must devise a system to systematically evaluate the quality of the retrieval. A true Bayesian approach would be to exclusively use the Bayes factor to evaluate the model fits. However, without well-defined prior odds for each model, we cannot quantitatively account for the prior probability of a given model. For example, based on the current understanding of these objects, the prior probability of a clear atmosphere model should be less than that of a cloudy model, but there is no clear way of assigning an objective probability. Instead we subjectively grouped some models into a 'low odds' category. While these are useful for validating our assumptions about the planets and testing the inclusion of di ff erent datasets, they should not contribute significantly to a final combined parameter estimate. We further sorted the retrievals a posteriori, creating in total three tiers of retrieval results, illustrated in Fig. 4. We focused our overall analysis on a subset of retrievals that are both plausible models, use consistent datasets, and produce physically reasonable results.  \nGroup A : This set of retrievals is defined as those with physically reasonable posterior values. Based on evolutionary models, it is expected that the HR 8799 planets have radii greater  \nFig. 4. Illustration of how retrievals are grouped in this work. Group A retrievals are selected based on having physically plausible posterior distributions. Group B retrievals are models that are subjectively chosen to have a 'reasonable' prior probability. Group C includes the entire set of retrievals run in this work.  \n<!-- image -->  \nTable 7. Self-consistent grid boundaries and step sizes.  \nthan 1 RJup. Even high-mass, cold brown dwarfs over 1 Gyr in age are found to have minimum radii of ∼ 0.88 RJup, with the minimum radius increasing with increasing metallicity (Burrows et al. 2011). Thus we exclude from group A any retrievals with a median retrieved radius less than 0.9 RJup. We additionally enforce that the mass estimate should be broadly consistent with the dynamical mass estimates: the median retrieved mass must be greater than 1 MJup and less than 22 MJup, approximately double the highest dynamical mass estimate of any of the planets (Zurlo et al. 2022). The exact positioning of these cuts does not significantly impact the results.  \nGroup B : this set of retrievals includes those that we consider to have a high prior model probability P ( M ); equivalently we are assigning a model prior probability P ( M ) = 0 to those models that we believe do not describe these atmospheres well. Specifically, we exclude models with a clear atmosphere, those with with poorly parameterised temperature profiles used during validation studies (e.g. retrievals using only 2 nodes to define a spline temperature profile), and those using data inconsistent with our fiducial dataset. Thus while the retrievals using OSIRIS data for HR 8799 c are highly ranked by the Bayes factor, we exclude them from our analysis and from the Bayesian Model Average, as the Bayes factor is only a relevant metric when comparing like datasets. Likewise, a clear atmosphere would require a diabatic temperature profile to explain the reddening of the emission spectra, which we do not include in the retrievals and therefore the clear models are unlikely to be physically meaningful. As the Bayes factors are weighted heavily towards the best retrievals, a weighted posterior distribution e ff ectively reduces to that of Group A.  \nGroup C : the complete set of retrievals included in this work, regardless of prior or posterior likelihood. As the full set of retrievals includes highly unrealistic atmospheric models by design, we do not present combined posterior distributions, but only explore specific comparison retrievals used to validate different model assumptions. Ultimately we found that all of the retrievals fall into group A, universally finding reasonable estimates for the planet masses and radii.  \nThe best set of retrievals is the intersection of groups A and B(indicated by A ∩ B when used to refer to a particular retrieval), which are retrievals that have physically plausible posterior values, and whose model we believe is a reasonable representation of the atmosphere. Tables E.1 to E.4 list the complete set of retrieval results, classifying the individual retrievals by group and sorting by the Bayes factor. We turn to Kass & Raftery (1995)  \nfor an interpretation of the Bayes factor in terms of frequentist statistical significance. Thus a ∆ log 10 Z > 1 is considered substantial evidence, and ∆ log 10 Z > 2 is considered strong evidence, equivalent to > 5 σ significance. Table 2 of Benneke & Seager (2013) present a similar, albeit slightly more conservative interpretation of the Bayes factor, with a similar threshold of log 10 Z = 2 . 1 for 'strong' evidence in favour of one hypothesis over another, equivalent to 3.6 σ significance.", '4.6. Bayesian model averaging': 'We used the techniques of Bayesian Model Averaging (BMA) in order to combine estimates of a single parameter over a range of models, following the review of Fragoso et al. (2018). This has recently been applied to exoplanet spectroscopy in Nixon et al. (2023), demonstrating that these methods provide more realistic posterior uncertainties. They highlight that to naively use BMA, the use of multiple duplicate models must be avoided to avoid the repeated contribution of that model to the average. As we do not have any identical models in our retrieval suite, BMA remains a valid approach.  \nConsider Bayes theorem for the i th model Mi for data D , with parameters θ i :  \nP ( θ i | D , Mi ) = P ( D | θ i , Mi ) P ( θ i | Mi ) P ( D | Mi ) . (12)  \nWe are interested in obtain a joint posterior probability distribution P ( θ | D , M ) for the subset of parameters θ that are shared between the set of models. From each model we require posterior probability distribution P ( θ i | D , Mi ), the likelihood P ( D | θ i , Mi ), and the evidence P ( D | Mi ) ≡ Z = R P ( D | θ i , Mi ) P ( θ i | Mi ) d θ i . We then assume a prior probability distribution over the full set of models under consideration, and therefore each model has an associated prior probability P ( Mi ). The choice of this prior probability should reflect the prior knowledge of the system under consideration. For example, the prior probability of a clear atmosphere model should be lower than that of a cloudy atmosphere model for the HR 8799 planets. However, quantifying this degree of certainty is highly subjective. We choose instead to use an uninformative prior distribution across N models for each planet:  \nP ( Mi ) = 1 N . (13)  \nThis allows the data to determine which models should be favoured based on the evidence.  \nConsidering all models in the range 1 to N , the posterior model probabilities given the data are  \nP ( Mi | D ) = P ( D | Mi ) P ( Mi ) P N j = 1 P ( D | Mj ) P ( Mj ) . (14)  \nThe marginal posterior distribution for a single parameter θ present in all of the models is thus  \nP ( θ | D ) = N X j = 1 P ( θ | D , Mj ) P ( Mj | D ) . (15)  \nThis combined posterior distribution folds in both the uncertainty from the data and prior distributions, but from the model uncertainty as well, providing a more robust estimate of the overall uncertainty on the inferred parameter.', '4.7. Self-consistent forward modelling': "In order to ensure that the retrieval results are robust and insensitive to the details of pRT , we fit each of the planet's spectra using several grids of 1D self-consistent models: ATMO (Phillips et al. 2020; Petrus et al. 2023), Sonora Bobcat, Cholla ,and Diamondback (Marley et al. 2021; Karalidi et al. 2021), Morley et al. (in prep), Exo-REM (Charnay et al. 2018), and petitCODE (Mollière et al. 2015, 2017). These models represent the current state-of-the-art in both cloudy and cloud-free self-consistent 1D models. The boundaries and intervals of each of these grids is presented in Table 7.", '4.7.1. ATMO': 'We used an up-to-date grid of ATMO of models from Petrus et al. (2023), which in turn is based on prior versions from Tremblin et al. (2015) and Phillips et al. (2020). ATMO is a clear atmosphere model, based on the hypothesis that diabatic convection (Tremblin et al. 2016, 2017), not clouds, are responsible for the reddening of the near-infrared spectra of directly imaged exoplanets and brown dwarfs. This convection is instigated by disequilibrium chemical processes that reduce the temperature gradient, thus reddening the atmosphere. In ATMO, this is parameterised through an e ff ective adiabatic index γ ad, which modifies the temperature gradient. The inclusion of this parameter meant that ATMO is the only clear atmosphere grid that produced a reasonable fit to the spectra of the HR 8799 companions.', '4.7.2. Exo-REM': 'The Exoplanet Radiative-convective Equilibrium Model (ExoREM, Baudino et al. (2015); Baudino et al. (2017)) is a selfconsistent model used to study directly imaged exoplanets and brown dwarfs (Charnay et al. 2018), but has also been extended to lower mass transiting planets (Blain et al. 2021). It implements a cloud microphysics model by combining AM01 with the timescale approach of Rossow (1978), which allows it to reproduce the L-T brown dwarf spectral sequence as a function of e ff ective temperature. Bonnefoy et al. (2016) used this grid to explore the atmospheres of the HR 8799 companions, finding atmospheres mildly enriched in metals ([M / H] = 0.5) and well constrained e ff ective temperatures. However, they developed a set of custom grids that implement detailed cloud properties to model the atmospheres, which are likely more suited to the HR 8799 planets than the more general publicly available grid.', '4.7.3. Sonora': 'The newly developed suite of Sonora models are designed to model the spectra and evolution of substellar atmospheres, covering the L-T-Y spectral sequence (Marley et al. 2021). Sonora comes in several flavours, implementing equilibrium chemistry in Sonora Bobcat (Marley et al. 2021), disequilibrium chemistry in Sonora Cholla (Karalidi et al. 2021), and cloudy atmospheres in Sonora Diamondback (Morley et al. 2024). Like Exo-REM and ATMO, the Sonora models are a 1D, radiative-convective equilibrium model that couples hydrostatic and thermochemical equilibrium temperature structure with a radiative transfer scheme to compute the atmospheric emission spectrum.  \nSonora Bobcat and Cholla did not fit the HR 8799 spectra at all, validating the necessity of cloudy (or similar) atmospheric models. Thus we continued only with the cloudy Sonora Diamondback models in order to interpret the HR 8799 atmospheres. While similar to Exo-REM in implementing clouds, Diamondback currently fixes the C / Oratio to a solar value of 0.458 (Lodders et al. 2009), preventing the measurement of this parameter, and potentially leading to biases in the remaining parameters.', '4.7.4. petitCODE': 'petitCODE is a radiative-convective and chemical equilibrium code used to compute the structures and spectra of exoplanet atmospheres (Mollière et al. 2015, 2017). We used the cool-cloudy and hot-cloudy grids computed for Stolker et al. (2020), spanning temperatures from 500-850 ( cool-cloudy ) and 1000-2000 K ( hot-cloudy ). The code setups are based on the work presented in Samland et al. (2017); Linder et al. (2019). Both grids implement the cloud model described in Ackerman & Marley (2001). While the cool grid only assumes Na2S and KCl clouds, and the hot grid adds Mg2SiO4 and Fe clouds.', '4.7.5. Grid fits': 'We performed Bayesian fits using species to interpolate the grids (Stolker et al. 2020), and MultiNest to sample the parameter space. 400 live points were used for these fits, with uniform priors on all parameters covering the grid ranges as described in Table 7, and an additional Gaussian prior on the planet mass. We fit for covariance width and strength for all IFS datasets following the method of Wang et al. (2020), as the empirical covariance matrices cannot be incorporated in species , other than for GRAVITY data. Fitting for the covariance parameters was universally favoured by the Bayes factor, and thus we only present the full fits. Consistent with the expectations of Greco & Brandt (2016) and Nasedkin et al. (2023), including these parameters also tended to broaden the posterior distributions, though posterior widths remain far narrower than the variation between the models.  \nIn addition to the Bayesian fits, we performed a simple χ 2 minimisation over each grid to avoid potential issues with interpolating the spectra along the di ff erent parameter axes. Using this framework, we identify the single best-fit spectra, as presented in Table 9.', '5. Results': 'Based on both the atmospheric retrievals and the self consistent grid fits, the HR 8799 planets share enriched atmospheres,  \nFig. 5. Bayesian-averaged posterior parameter distributions for each of the HR 8799 planets based on the group A ∩ B set of retrievals. In faint lines beneath the total posterior distribution are the individual contributions from di ff erent retrievals. The coloured dashed and dotted lines indicate the median and ± 34 . 1% confidence regions respectively. The vertical grey lines indicate typical parameter values (e.g. solar metallicity and C / O) and serve as a visual reference for each parameter. For the planet mass, they indicate the dynamical mass estimate from Zurlo et al. (2022).  \n<!-- image -->  \nwith stellar-to-superstellar C / O ratios. The properties of each atmosphere for subsets of the ensemble of retrievals are summarised in tables 1-4 for planets b-e, respectively. These tables also contain ranges of plausible parameter values for each planet based on the aggregate of the self-consistent models, while Fig. 5 shows the distributions for a subset of the key parameters. These estimates are synthesised from the results of the Bayesian fits and the χ 2 minimisation, rejecting solutions with unphysical masses and radii. The full results of the grid fits are found in Table 8 for the Bayesian fits and Table 9 for the single-best χ 2 fits. In these tables, an index is assigned to each retrieval, with the format planet.group.index , which serves as the retrieval identifier throughout the text.  \nNearly 100 retrievals were performed for this analysis across the four companions in order to derive robust constraints on primary planetary properties. The aggregate results of the Bayesian model average of group A ∩ B retrievals are presented in Fig. 5, with the best fit models for the same sample of retrievals compared to the data in Fig. 6. Overall we find consistently good fits for all four planets, with best-fit χ 2 /ν < 2 for each planet. The self-consistent grid-fits incorporate additional physical processes and require fewer free parameters than the retrievals, making them less flexible. The additional processes, such as radiative-convective equilibrium, chemistry, and cloud physics act as narrow priors on the bulk atmospheric parameters. While the lack of flexibility leads to worse fits when compared to retrievals, the fits of multiple self-consistent grids still result in mutually consistent results. Notably, petitCODE gives consistent parameter estimates even though it is strongly disfavoured by the goodness-of-fit metrics.  \nBoth the pRT retrievals and the self-consistent grid fits show that all of the atmospheres are strongly enriched in metals, finding median [M / H] ≳ 1.0 for each planet. This is driven by  \nthe carbon and oxygen abundances as measured through H2O and CO, and particularly through the CO absorption feature at 2.3 µ m. The C / Oratios decrease with decreasing separation from 0 . 78 ± 0 . 04 in HR 8799 b to 0 . 60 ± 0 . 05 for HR 8799 d, with HR 8799 e breaking the trend with carbon rich C / O ratio of 0 . 87 ± 0 . 02. In addition to water and CO, HCN is confidently detected in HR 8799 c and e, while CH4 is detected in HR 8799 c.  \nWe find that cloudy atmospheres are universally favoured over clear atmospheres. Extended silicate clouds, together with an iron cloud deck are the preferred models for the three inner planets, with the cooler HR 8799 b showing evidence for Na2S clouds.  \nThe retrievals found masses, radii, and surface gravities consistent with evolutionary models. By construction, these form Group A of our retrievals. In general, retrieving the planet radius is challenging: the radius of directly imaged exoplanets is commonly underestimated by both retrievals and grid fits (e.g. Bonnefoy et al. 2016). These unphysical solutions are found in several of the grid fits, while nearly every retrieval finds plausible values of the mass and radius. The small radius is often compensated for by adjusting the temperature, metallicity, or cloud properties. Nevertheless, in all cases the most favoured retrieval produced masses consistent with dynamical mass measurements and radii consistent with evolutionary models. In the group A ∩ B retrievals that used dynamical mass priors, the retrieved mass estimate is within 1 σ of the dynamical mass for b, d, and e, with the posterior width consistent with the prior width which. For HR 8799 c the retrieved mass is moderately higher than the dynamical mass estimate at 8 . 5 ± 0 . 5 M Jup compared to the dynamical mass of 7 . 7 ± 0 . 7 M Jup . However, if the dynamical mass is not used as a constraint and log g is freely retrieved, the resulting mass estimate is found to be much larger than the dynamical  \nFig. 6. Best-fit temperature profiles and spectra for the group A ∩ B retrievals. From top to bottom are HR 8799 b, c, d, and e.  \n<!-- image -->  \n07  \nmass estimate, highlighting the importance of including additional constraints in a retrieval framework. The estimate of the e ff ective temperature of each planet is consistent with previous findings, with e, c, and d sharing temperatures around 1100 K, and b being cooler at around 950 K. While we do not perform a comparison to spectral templates due to the inferred high metallicity, low surface gravity, and potential complications from the viewing angle of the planets, we find that the spectral types found in Bonnefoy et al. (2016) remain a good description of the four planets. HR 8799 d and e are similar to late-L type brown dwarfs, consistent with their inferred e ff ective temperature. HR 8799 c is more likely fit by an early-T spectral type as it is a few kelvin  \ncooler than e and d and is beginning to show spectral features from CH4. At 950 K HR 8799 b is solidly in the T dwarf regime.  \nIn this section, we present the measured properties for each of the four planets individually. Following that, in Section 5.5 we present a detailed discussion of the results and challenges of comparing di ff erent thermal structures, chemistry and cloud parameterisations, and data inclusion.', '5.1. HR 8799 b': 'For HR 8799 b we included the GRAVITY and OSIRIS spectra, together with the full set of photometry, allowing the OSIRIS  \nFig. 7. [M / H] posterior distributions for all retrievals in group A ∩ B . From left to right are the distributions for b, c, d, and e. The vertical line indicates solar metallicity. Model keys are as in tables E.1-E.4, and are sorted by the Bayes factor from top to bottom.  \n<!-- image -->  \ndata to float as the published data are not flux calibrated. For the grid fits we included additional parameters to describe the covariance of the OSIRIS data.  \nHR 8799 b is the coldest and lowest mass planet in the HR 8799 system. The best estimate of these parameters via Bayesian averaging of group A ∩ B retrievals finds an e ff ective temperature of 942 + 12 -16 K and a mass of 6 . 0 + 0 . 4 -0 . 3 MJ, driven by the use of the dynamical mass estimate of Zurlo et al. (2022) as a prior. The radius is slightly inflated compared to jupiter, with ( Rpl = 1 . 10 + 0 . 03 -0 . 03 ); combining the mass and radius estimates leads to log g = 4 . 10 + 0 . 03 -0 . 04 . This is consistent with the Bayesian grid fits, though the single-best χ 2 fits found both lower (log g = 3 . 5) and higher (log g = 4.5) solutions. The Bayesian averaged results (for group A ∩ B ) are driven by a single retrieval, with a ∆ log 10 Z of 4 relative to the next best retrieval. This single best retrieval  \nuses the Z23 temperature profile and free chemistry, finding a metallicity of 0 . 96 + 0 . 08 -0 . 08 and a C / O of 0 . 78 + 0 . 03 -0 . 04 .  \nThe grid fits find temperatures between 850 K -1100 K, with the single best fit, found using Exo-REM, finding T e ff = 850 K. Using the Bayesian fit, Exo-REM is again the most favoured model by the Bayes factor, and finds T e ff = 862 + 21 -6 K and a radius of 1 . 19 + 0 . 03 -0 . 03 R Jup, consistent with expectations from evolutionary models (e.g. Marley et al. 2012). The ATMO and Sonora Diamondback models finds a somewhat higher temperature and a smaller radius, but are disfavoured by the Bayes factor. petitCODE is divided into cool-cloudy and hot-cloudy grids, with the cool grid extending up to 850 K, and the hot grid beginning at 1000 K. As the e ff ective temperature of HR 8799 b likely falls between these grids, it poorly fit the data, although  \nFig. 8. C / O posterior distributions for all retrievals in group A ∩ B . From left to right are the distributions for b, c, d, and e. The vertical line indicates stellar C / O. Model keys are as in tables E.1-E.4, and are sorted by the Bayes factor.  \n<!-- image -->  \nthe remaining parameters of log g and f sed are compatible with the other self-consistent models.  \nFigure 7 shows the variation of the metallicity across the different models considered. We see that strong enrichment solution is almost always found, particularly by models preferred by the Bayes factor. The degree of enrichment does not systematically vary between disequilibrium and free chemistry retrievals. All of the self-consistent models favour high metallicity solutions, reaching the upper bounds of the grid in all cases.  \nSimilarly to the metallicity, the C / O ratio is shown for the di ff erent retrievals in Fig. 8 and T e ff in 9. The C / O ratio is generally constrained to between 0.6 and 0.8. However, for free chemistry retrievals, the inferred C / O ratio only indicates the gas-phase composition. Additional oxygen is sequestered in the silicate clouds: accounting for this sequestration would result in  \na lower C / Oratio. Among the grid fits, shown in Fig. 10, the C / O ratio shows more variation, ranging from the lower bound of the ATMO grid at 0.3 to 0.55 from Exo-REM. The best fit models from Exo-REM are consistent with the stellar value of 0.54, and additional data covering the near infrared water features is likely necessary to improve these constraints. From Fig. 11 we find that both the free retrieval and disequilibrium retrievals display similar trends in the retrieved chemical abundances, finding that nearly 10% of the atmosphere is CO by mass, while water has a lower abundance of around 1% by mass. The best-fit free retrieval finds systematically lower abundances for both of these species compared to the best-fit disequilibrium retrieval. H2S and CH4 are found to be the next most abundant species in both the disequilibrium and free chemistry retrievals. However, even though their abundances are constrained by the pos-  \nFig. 9. Te ff posterior distributions for all retrievals in groups A ∩ B . From left to right are the distributions for b, c, d, and e. The vertical line indicates 1000 K. Model keys are as in tables E.1-E.4, and are sorted by the Bayes factor.  \n<!-- image -->  \nterior distribution, there is no evidence for their detection when comparing between the full free chemistry retrieval (b.AB.1) and retrievals without these species (b.A.0 and b.A.2). The b.AB.1 retrieval uses the Z23 temperature profile, free chemistry, and clouds condensing at their equilibrium position. This same setup is used for b.A.0 and b.A.2, apart from the exclusion of CH4 and H2S respectively. While the inferred CH4 abundance of log X CH4 = -5 . 0 ± 0 . 4 is compatible with Barman et al. (2015), there is no statistical evidence to support the detection.  \nIn order to compare the cloud composition for HR 8799 b, we considered a set of retrievals using the same temperature structure (M20) and cloud parameterisation, and vary the cloud optical properties and condensation curve. We found that Na2S clouds are preferred over silicate, iron, and KCl clouds ( ∆ log 10 Z ≥ 2). Patchy clouds were also explored, but no evi-  \ne was found for patchiness, regardless of cloud composition. At the temperature of HR 8799 b silicate clouds are expected to occur below the photosphere, with Na2S or KCl clouds becoming the primary aerosol opacity source. This result is likely driven by the condensation temperature of Na2S rather than the optical properties; unlike silicate clouds which have strong absorption features in the mid-infrared, crystalline Na2S is featureless out to 15 µ m, apart from a characteristic scattering slope.', '5.2. HR 8799 c': 'In our standard retrieval setup, we included data from SPHERE, GPI, CHARIS, GRAVITY, and ALES for HR 8799 c. We omitted the OSIRIS data, as it overlaps nearly completely with the GRAVITY data, is not flux calibrated, and requires either bin-  \nFig. 10. Grid fits from Exo-REM, Sonora Diamondback, ATMO, and petitCODE. ALES data for HR 8799 d are scaled by the overall best-fit scaling parameter (Exo-REM, 1.18). Fits are the single best-fit χ 2 model from the grid.  \n<!-- image -->  \nning to lower resolution to be fit with the c-k opacity tables, or the use of the higher-resolution line-by-line opacities to fit the full resolution data, dramatically increasing computation time. Nevertheless, we performed several retrievals incorporating the OSIRIS data rather than the GRAVITY data to determine how this choice impacts the retrieved chemistry and clouds, and to determine if we can reproduce the findings of Konopacky et al. (2013). Overall, HR 8799 c proved challenging to fit: many retrieval setups either required far more model computations before convergence than the other three planets, or failed to converge entirely. The grid-fits for HR 8799 c also displayed the greatest variation between models.  \nLike HR 8799 b, the group A ∩ B retrievals of HR 8799 c are dominated by a single retrieval, with ∆ log 10 Z = 2 compared to the next best retrieval. This retrieval used the Z23 temperature profile and free chemistry, and requires high-altitude silicate clouds (log P base = -3 . 4 ± 1 . 8 bar). The inferred e ff ective temperature of 1159 + 11 -12 K is consistent with the range of temperatures found by the self-consistent grids, which spans from 1100 K to 1200 K, with the Bayesian fits averaging around ∼ 1200 K. The retrieved mass (8 . 5 + 0 . 4 -0 . 5 M Jup) is slightly higher than the dynami-  \ns estimate, though radius (1 . 10 ± 0 . 01 R Jup) is compatible with evolutionary models, and from these we derive a log g of 4 . 26 + 0 . 02 -0 . 03 . The grid-fits also found plausible radii, favouring values slightly larger than 1 R Jup . All of the self-consistent models found temperatures of 1100 -1200 K, and log g between 3.5 and 4.5. From the single-best χ 2 fits, the ATMO model found a log g of 4.5, while the remaining models find a lower solution of 3.5.  \nAs with the other HR 8799 planets the retrievals favour highly enriched solutions, finding [M / H] = 1 . 27 + 0 . 05 -0 . 05 . The disequilibrium chemistry retrievals find slightly lower metallicities than the free chemistry retrievals, with the most-favoured disequilibrium retrieval finding a metallicity of 1 . 05 ± 0 . 04. The data for HR 8799 c are highly discrepant in the H band (Figure 3), with the CHARIS data and photometry being around 50% brighter than the SPHERE and GPI data. As the metallicity is highly sensitive to the amplitude of the J, H, and K band peaks, such discrepancies need to be resolved in the data to ensure reliable measurement of this parameter. All of the grid-fit solutions require high metallicity and are limited by the grid boundaries. The C / O ratios between the Bayesian fits and the χ 2 minimisation are consistent, typically favouring stellar to slightly substel-  \nTable 8. Grid-retrieval results', '5.3. HR 8799 d': 'For HR 8799 d we included all available spectroscopic data as described in Section 3, but include a scaling factor for the ALES data set as otherwise it is incompatible with NACO photometric observations. This shifts the mean L-band flux of HR 8799 d to a similar magnitude as e and c.  \nUnlike HR 8799 b or c, d is well fit by a broad selection of models, and no single retrieval dominates the Bayesian average of figure 5. All of the preferred models ( ∆ log 10 Z < 2) used disequilibrium chemistry and the M20 or Z23 temperature profiles. Using the Bayesian average, the retrieved e ff ective temperature is 1179 + 31 -28 K, compatible with the ranges found by the grid fits, which find T e ff from 1150 K - 1300 K. The planet mass (9 . 2 ± 0 . 1 M Jup) is tightly constrained by the dynamical mass prior, which in turn allowed for precise measurement of the planet radius (1 . 26 ± 0 . 07 R Jup) and log g (4 . 18 ± 0 . 05). This surface gravity is consistent with the estimates from the Bayesian grid fits, but is significantly higher than the 3.0-4.0 range found by the χ 2 fits. The self-consistent fits find marginally smaller radii than the retrievals, generally between 1.1 and 1.2 R Jup .  \nThe metallicity of HR 8799 d is consistent with HR 8799 b and c, with [M / H] = 1 . 2 ± 0 . 2. While most grid-fit solutions also favoured high-metallicity atmospheres, the Bayesian fit with the ATMO model find a solution consistent with solar metallicity, though the radius was inconsistent with evolutionary models (0 . 926 ± 0 . 004 R Jup ). The best fit Exo-Rem model find a metallicity of 0.78 ± 0 . 03, a radius of 1 . 168 ± 0 . 006 R Jup and an e ff ective temperature of 1155 ± 5 K. The C / O ratio is always found to be consistent with the stellar value, with retrievals finding C / O = 0 . 60 + 0 . 04 -0 . 06 . The grid-fits are also typically consistent with stellar, though the best fit Exo-Rem model found a substellar C / O ratio of 0.2.  \nPatchy clouds are marginally disfavoured by the Bayes factor, and the patchiness is poorly constrained, finding f cloud =  \n0 . 45 ± 0 . 3. The most favoured solutions require either amorphous or crystalline MgSiO3 or crystalline MgSiO3, with a marginal preference for the amorphous structure. Each of these compositions displays slightly di ff erent near infrared slopes, shown in Figure 16. However, such a slope can be induced by various sources of continuum opacity that may not be fully accounted for in the retrieval. Thus mid-infrared observations of the silicate absorption features are necessary to robustly constrain the composition and particle geometry. There is no preference for a free cloud base pressure compared to the equilibrium position, suggesting that the AM01 model is su ffi cient to describe the clouds in this atmosphere.  \nWhile the disequilibrium retrievals are favoured over the free chemistry retrievals, there is excellent agreement between the freely retrieved abundances and the disequilibrium chemical profiles. Water and CO are both highly abundant, and the freely retrieved abundances agree with the disequilibrium profiles to within 1 σ . No other species are both highly abundant and well constrained, so we do not perform leave-one-out retrievals to test for their presence. However, even at low abundances the freely retrieved CH4 abundance is compatible with the disequilibrium profile.', '5.4. HR 8799 e': 'The measurements of HR 8799 e largely reinforce existing literature values. The single most favoured retrieval, which also dominates the group A ∩ B , used free chemistry and the Z23 temperature profile, together with a cloud base calculated using equilibrium condensation. This lead to similar results as for HR 8799 c and d, and is consistent with the results of the grid fits. An e ff ective temperature of 1138 + 30 -22 K is retrieved, compatible to within the uncertainties of the grid-fits, which found a temperature range of 1100 K to 1200 K. The mass posterior  \nFig. 11. Mass fraction abundance profiles as a function of pressure in the atmospheres of HR 8799 b, c, d, and e. The solid lines indicate the median values of the most favoured disequilibrium retrieval for each planet, with the shaded region indicating the 90% confidence interval. The circular markers show the median values for each species from the most favoured free chemistry retrieval, with the error bars indicating the 90% confidence interval. The position along the pressure axis is arbitrary. The minimum mass fraction allowed in the free retrievals was 10 -7 .  \n<!-- image -->  \nwas determined by the dynamical mass prior, as was the planet radius, finding Mpl = 7 . 5 + 0 . 6 -0 . 6 and Rpl = 1 . 13 + 0 . 05 -0 . 05 respectively. This leads to a log g of 4 . 18 + 0 . 06 -0 . 05 , slightly higher than the selfconsistent estimates of 3.5-4.0. The overall best self-consistent model by the χ 2 was Exo-Rem, which finds an e ff ective temperature of 1100K, a radius of 1.15 R Jup , and a somewhat low log g of 3.5. ATMO is the most favoured self-consistent model when using the Bayesian framework, though it found an unphysically small radius and higher temperature than other models. HR 8799 e is the only companion for which the MIRI photometry is not convincingly fit, as seen below in figure 20. Every model underestimated the flux beyond 10 µ m relative to the measurements, though this may be due to contamination from the host star or inner disc (Boccaletti et al. 2023).  \nCompared to the other three planets, HR 8799 e is found to have an even more metal rich atmosphere, with [M / H] = 1 . 9 + 0 . 1 -0 . 2 . Metallicities > 1 were a universal feature of the retrievals for e. This was reinforced by the grid-fits, which uniformly find strong enrichment, running into the upper grid boundaries. Free chemistry retrievals are always preferred over the disequilibrium retrievals; from these we found the group A ∩ B C / O ratio is 0 . 87 + 0 . 02 -0 . 02 . As the free chemistry retrieval C / O ratio only accounts for the gas-phase abundances, there is additional oxygen  \nsequestered in the silicate clouds that could reduce the C / Oratio. However, the most favoured disequilibrium retrieval finds a similar value of 0.83 ± 0 . 02, suggesting that HR 8799 e is somewhat of an outlier compared to the other three planets. Using a similar setup to Mollière et al. (2020) (e.AB.11), we find C / O = 0 . 78 + 0 . 03 -0 . 03 . The grid fits tend to find C / O ratios compatible with the stellar value, though the overall best fit Exo-Rem model also finds a higher value of 0.8.  \nIn addition to the water and CO rich atmosphere, HCN is found to be highly abundant, with log X HCN = -2 . 26 ± 0 . 11. This detection is strongly favoured by the Bayes factor, with log 10 Z = 7 . 5. As with HR 8799 c and shown in Fig. C.3, this detection was driven by using the HCN opacity to fit the ALES spectrum, though changes to the shape in the H and K-bands also provide a slightly better fit as well. This is slightly enriched compared to equilibrium chemistry predictions, but is expected for a metal rich planet with a relatively high C / Oratio (Giacobbe et al. 2021), and can be produced through photochemical reactions (Moses et al. 2013). In contrast, the presence of CH4 is poorly constrained.  \nAs with HR 8799 c and d, the most favoured retrieval for e favours silicate clouds with a deeper iron deck, both condensing at their equilibrium locations, with no preference for patchy clouds. A broad range of disequilibrium retrievals (e.AB.7-11)  \nFig. 12. Temperature profiles for HR 8799 e. In blue are temperature profiles from disequilibrium retrievals, while in red are free chemistry retrievals. The shaded regions indicate 90% confidence intervals. Also included are the temperature profiles from the best fit self-consistent models.  \n<!-- image -->  \nfound consistent cloud properties, with condensation at the equilibrium location preferred, with a silicate abundance lower than predicted by equilibrium chemistry (between 10 × to 100 × less than equilibrium), and iron abundances consistent with equilibrium. Individual f sed parameters for the silicate and iron clouds were not required, with variations in the Bayes factor driven more by the choice of patchiness and temperature profile.', '5.5. Impacts of modelling choices': 'While we examined the primary atmospheric characteristics of each of the four HR 8799 planets, it is crucial to understand how the choice of model impacts these measurements. As described in Section 4, we performed retrievals using a broad selection of thermal structures, chemical parameterisations, and cloud models, each of which we is examined in detail below.', '5.5.1. Thermal structure': 'For each planet we performed retrievals using four di ff erent temperature profile parameterisations. The Z23 profiles is used in the most favoured retrieval for three of the four planets, while for HR 8799 d there is equal evidence for the Z23 and M20 profiles. The Guillot profile is found to be the second most preferred profile for HR 8799 c and e, while the spline profile using six nodes is strongly disfavoured by the Bayes factor. In Figure 12 we compare the retrieved and self-consistent temperature profiles for HR 8799 e. We find that there is excellent agreement between nearly all of the retrievals in the photosphere region, as well as with the best-fit Exo-Rem temperature profile. There is little variation in the photosphere region between the disequilibrium and free chemistry models. Only one model - the  \nspline profile with free chemistry - found a profile more similar to that of the clear ATMO profile, though it is strongly disfavoured by the Bayes factor compared to the other free chemistry retrievals. The spline profile is the only profile that does not explicitly assume an adiabat deep in the atmosphere or rely on assumptions from self-consistent models, and so we cannot fully rule out the diabatic profiles of Tremblin et al. (2015). The bulk atmospheric properties are also reasonably consistent across the di ff erent temperature parameterisations. While there are statistically significant variations in the C / O ratio between the di ff erent parameterisations, they remain broadly consistent between 0.7 and 0.9.  \nAlthough the spline profile is disfavoured by the retrievals, it is a useful parameterisation to determine the amount of flexibility required by the model, and to explore the known degeneracies between the atmospheric thermal structure and clouds Tremblin et al. (2015, 2016). We performed a series of retrievals on HR 8799 b, varying the number of nodes in the spline profile and observing how the retrieved profile changes with the increased flexibility. We repeated this test for both a clear atmosphere model and an model with clouds condensing at the equilibrium base pressure. We find no significant di ff erences in the temperature profiles between the clear and cloudy atmospheres. For HR 8799 b, the Bayes factor favours retrievals with three or four nodes in the spline profile. Fewer nodes mean the profile cannot be accurately modelled, while more nodes add additional parameters without improving the fit to the spectra.', '5.5.2. Chemistry': 'For all four planets we performed retrievals using a grid derived from an equilibrium chemistry solver with disequilibrium H2O-CO-CH4 quenching, as well as free chemistry retrievals where we directly retrieved the mass-fraction abundance of various species. Fig. 11 shows the abundance profiles from the best fit disequilibrium and free chemistry retrievals for each planet. Both types of retrievals produce consistent metallicities and C / O ratios for each planet: overall there is excellent agreement in the water and CO abundances, which are the primary opacity sources in these atmospheres. Only HR 8799 b shows statistically significant discrepancy between the two methods for these species, with the free retrieval finding a slightly lower metallicity than the disequilibrium retrieval. Even for trace species, the free retrievals and disequilibrium retrievals are largely compatible, though only a few species have statistically significant detections in the free retrievals.  \nThe strongest trace species detections are HCN in HR 8799 c and e, at abundances far higher than predicted by the equilibrium model. HR 8799 b also has a well constrained H2S measurement, though it is not statistically significant. CH4 is significantly detected in the atmosphere of HR 8799 c, demonstrating that with su ffi cient S / N and wavelength coverage, it is possible to constrain abundances at below 10 -4 by mass. The free chemistry detection of CH4 is at a moderately higher abundance than in the best fit disequilibrium model. While it is likely also present in the cooler atmosphere of b, additional wavelength coverage or higher S / N observations are required for a significant detection. For HR 8799 c, we include in Figures C.2 and C.3 in the appendix comparisons between the HR 8799 c data and models both with and without the contribution of CH4 and HCN opacity, demonstrating the impact of these species on the spectral shape. While the HCN detection is driven primarily by the low flux of the ALES data, there is a significant change in the H-band shape, as well as a slight change in the peak amplitude of the K-band.  \nFig. 13. Metallicities and elemental number ratios for the HR 8799 planets derived from Baysian averaged group A ∩ B free chemistry retrievals.  \n<!-- image -->  \nAs the ALES data are relatively low S / N , additional H and L band data should be obtained to confirm this detection. However, the CH4 detection is driven by modest improvements in the fit throughout the K-band. Several abundant species predicted by the equilibrium network are not confidently detected by the free chemistry retrievals, such as CO2, NH3, and H2S. Additional wavelength coverage or higher spectral resolution may allow for the characterisation of such species.  \nIf we take the averaged free retrieval results at face value, we can derive elemental abundance ratios for each of the four planets, using a similar method to calculating the metallicity. Taking the volume mixing ratios of each molecular species, we can count the total number of C, N, O, and S atoms, and calculate the ratio relative to the planetary hydrogen abundance. Thus for example  \nC / H = XCO + XCO 2 + XCH 4 + XHCN X H2 + 2 X H2O + 2 X H2S + X HCN + 4 X CH4 + 3 X NH3 , (16)  \nwhere all abundances are measured in number fraction. These ratios are then normalised to the solar values from Asplund et al. (2009).  \nIn Fig. 13, we present the elemental abundance ratios for each of the four planets. We find that most elements are enhanced relative to solar for all four planets. HR 8799 b appears depleted in nitrogen relative to the other planets, likely due to the nondetection of NH3, which will require observations of the 10 µ m feature to characterise. The HCN detections in HR 8799 c and e tightly constrain the nitrogen enhancement, though these planets still appear less enriched in nitrogen than in carbon or oxygen, though this is again likely due to additional nitrogen stored in N2 and NH3, whose opacities are inaccessible at these wavelengths. The sulphur elemental ratio is poorly constrained for all of the planets apart from HR 8799 b, which has a precise - though not statistically significant - constraint on the H2S abundance. HR 8799 b appears sulphur rich, while the remaining planets appear  \nconsistent with the solar value, or slightly depleted in sulphur, though this is largely due to a lack of measured chemical species.  \nThe C / O ratio is a key consideration for planetary atmospheres. However, measuring the atmospheric C / O ratio and linking it to the bulk planet composition is far from trivial. Lodders & Fegley (2002) and Lodders (2003) explore the chemistry and condensation of substellar atmospheres, identifying the condensation sequence of refractory species throughout these atmospheres, finding that at typical L-dwarf temperatures there will be silicate clouds in the photosphere region. Fonte et al. (2023) demonstrate how oxygen is sequestered in silicate clouds and other refractory species. This was followed by the recent work from Calamari et al. (2024), who calculate the bulk planet C / O ratio from the atmospheric ratio, finding that the median sequestration of oxygen due to this condensation is 17 . 8 + 1 . 7 -2 . 3 %. They also identify a relation between the bulk and observed C / Oratio:  \n(C / O) obs ≈ (C / O) bulk 1 -0 . 371 (C / O) bulk . (17)  \nSolving for (C / O) bulk , we find that for HR 8799 e, with an observed C / O ratio of 0 . 88 + 0 . 02 -0 . 02 , should have a bulk C / O ratio of 0.66, much closer to the stellar value of 0.54. Likewise, HR 8799 d has the lowest observed C / O ratio of 0 . 61 + 0 . 03 -0 . 04 , which translates to a modestly substellar bulk C / O ratio of 0.50. In general, this relation reduces the variation between the four planets, and brings the planetary C / Oratio more in line with the known stellar value.  \nWhile both chemistry models are compatible, they also share similar biases. The free chemistry model measures the gas phase abundance in the photosphere, and is primarily impacted by the atmosphere above the silicate clouds. Conversely, the underlying equilibrium model does remove oxygen from the gas phase due to condensation, though the additional flexibility in the cloud parameterisation means that it is only exactly correct for f sed = 1. By parameterising disequilibrium via fixing the chemical abun-  \ndances above a quench point, the model may lose this sensitivity, and therefore measures the abundances of CO, H2O, and CH4 in a similar fashion to the free chemistry model. Thus the C / O ratio as inferred by both models will be strongly impacted by the oxygen-depleted region above the silicate clouds, leading to over-estimates of the C / Oratio. Throughout this work we present these measurements, but we note that the adjustment introduced by Calamari et al. (2024) is likely a more accurate estimate of the bulk planet composition.  \nIn addition to the elemental ratios, we also computed Z pl / Z ⋆ , which allows us to directly compare our metallicities to literature values, such as those of Thorngren et al. (2016). We converted the metallicity [M / H] of each atmosphere to Z pl using the methods of Thorngren & Fortney (2019), adapting for our own notation:  \n10 [M / H] = 1 + Y / X GLYPH<16> Z -1 pl -1 GLYPH<17> GLYPH<16> µ Z µ H GLYPH<17> , (18)  \nwhere X , Y , and Z are the solar hydrogen, helium and metal mass fractions and µ is the mean molecular weight of the metal content of the atmosphere. Rearranging and substituting in the measured atmospheric metallicity [M / H], we find:  \nZ pl =          1 + 1 + Y / X 10 [M / H] GLYPH<16> µ Z µ H GLYPH<17> ( Z / H ⊙ )          -1 . (19)  \nWe take the same assumptions as Thorngren & Fortney (2019), taking µ Z to be 18, assuming most of the metal content is in water, µ H to be 1 for atomic hydrogen, and Y / X to be 0.3383 as in Asplund et al. (2009). As the metallicity of HR 8799 A is near solar, Z ⋆/ H is taken to be the solar value of Z / H ⊙ = 1 . 04 × 10 -3 . To normalise to the stellar metallicity we follow Thorngren et al. (2016) and calculate the Z ⋆ as  \nZ ⋆ = 0 . 014 × 10 Fe / H . (20)  \nFor HR 8799, we used solar metallicity to calculate Fe / H, but refer to the discussion in Section 2.1.  \nDisequilibrium chemistry has long been thought to play a key role in shaping the composition of the HR 8799 atmospheres (e.g. Marois et al. 2008). With well-constrained chemical abundances, we can start to place limits on the strength of vertical mixing that drives this disequilibrium. The quench pressure we retrieve is defined as the level below which (in pressure) the abundances of H2O, CO, and CH4 become vertically constant. This parameterises dynamical mixing that homogenises the upper layers of the atmosphere. More rigorously, the quench point is defined as the point at which the chemical timescale t chem and the mixing t mix are equal. Following Zahnle & Marley (2014), the mixing timescale is defined through the ratio of the local atmospheric scale height H to the vertical eddy di ff usion coe ffi -cient, K zz:  \nt mix = H 2 K zz . (21)  \nThe chemical timescale depends on the reaction rates involved. Considering the CO-CH4 reaction chain, Zahnle & Marley (2014) derive a timescale at the quench point ( tq ) for CO. For strong mixing, pulling from material at depths below the point where the atmosphere is 1000 K, they find the timescale well described by an Arrhenius relation for quench pressure p in bar, metallicity m , where m = 10 [M / H] , and temperature T in kelvin:  \nt q1 = 1 . 5 × 10 -6 p -1 m -0 . 7 exp (42000 / T ) s . (22)  \nTable 10. Quench pressures, vertical mixing parameters, and sedimentation fractions for the HR 8799 planets.', '5.5.3. Clouds': 'For the inner three planets, we find the most favoured solution is an optically thin silicate cloud lying above a deeper, optically thick iron cloud deck, as seen in figure 14. In our standard setup, the clouds were parameterised as in AM01, with the clouds condensing at the intersection of their condensation curve and the temperature profile, with their extent determined by f sed. The cloud mass fraction was allowed to scale from equilibrium. In totally free retrievals, the cloud abundances, locations, and vertical extents were all free parameters of the model. This decouples the clouds from both the chemistry and the atmospheric thermal structure, allowing them to fit the spectral shape, but in potentially nonphysical configurations. In this framework the cloud extent was then parameterised as in AM01, determined by f sed and K zz. Using this setup we find an optically thin silicate cloud lying above a compact iron cloud, while the AM01 setup finds an iron cloud that extends high above the silicate cloud. Depending on the choice of other parameters, either the free cloud base or the equilibrium base can be preferred by the Bayes factor. HR  \n8799 e free chemistry retrievals strongly favour the equilibrium condensed clouds, while the disequilibrium retrievals favour the free cloud base setup. However, in general the clouds condensing at equilibrium are the most favoured setup for each planet. Without broad wavelength coverage and high spectral resolution to probe a high dynamic range in pressure, it is di ffi cult to robustly distinguish between the di ff erent potential cloud structures.  \nThere is a marginal preference for clouds parameterised using a Hansen (1971) particle size distribution over a log-normal distribution (e.AB.20 over e.AB.25, ∆ log 10 Z = 0 . 7), though this was only compared for HR 8799 e using the Z23 profile, and assuming the clouds condense at their equilibrium saturation location. In this case ah is calculated from f sed and K zz, and a lower f sed for the MgSiO3 cloud was retrieved than with the log normal distribution ( f sed , Hansen = 1 . 18 ± 0 . 20). The e ff ective distribution width parameter, bh was found to be 0.016, which is narrower than the distributions found by Burningham et al. (2021).  \nIn general our clouds are comparable to those of Burningham et al. (2021) and other similar studies (e.g. Mollière et al. 2020; Vos et al. 2023; Balmer et al. 2023). Burningham et al. (2021) find a combination of MgSiO3, SiO2, and Fe clouds provides the best-fit model to an ultracool field dwarf, 2MASSW J2224438015852. Figure 14 shows that the MgSiO3 clouds in their model are located at 10 -3 bar with a maximum optical depth of τ = 0 . 3 at 1 µ m and an e ff ective particle radius of ∼ 0 . 04 µ m. The SiO2 clouds are slightly deeper, at 10 -2 bar, and are optically thick at 1 µ m. This is the same location where we find MgSiO3 clouds condensing in the atmospheres of HR 8799 c, d, and e. The iron cloud is found to be deep in the atmosphere, though with an extended structure, with some contribution at the same altitude as the silicate clouds. While we did not fit for a three cloud model, the similar locations and optical depths of these clouds suggests similar structures between the objects, even though they di ff er in e ff ective temperature by hundreds of kelvin.  \nFigure 15 highlights how the di ff erence in particle radius contribute to the di ff erence in the wavelength dependence of the cloud optical depth. The cloud particle radius in the AM01 model is a function of many atmospheric factors, including the temperature, mixing strength, f sed, and particle number density. We see that changes in the particle radius are correlated with  \nFig. 15. E ff ective particle radii as a function of altitude for silicate (blue) and iron (red) clouds. The solid lines indicate the radii for Case e.AB.13, which used a free cloud base pressure and abundance, while the dashed lines are for Case e.AB.2, which used equilibrium condensation and scaled equilibrium abundances. The horizontal lines indicate the cloud base pressure. The green line indicates the temperature profile. The shaded regions indicate unphysical particle sizes where the opacity contribution is set to 0.  \n<!-- image -->  \nchanges in both the temperature and the particle density. Luna & Morley (2021) explore the impact of particle size and composition on the spectral signatures of clouds in young brown dwarfs in the mid-infrared, finding that small particle sizes will lead to visible features in the planetary spectrum. However, we see in Figure 15 that the particle sizes in regions of the atmospheres with significant cloud mass fraction tend to be larger ( > 1 µ m), and that there are no indications of deep cloud absorption features in the mid-infrared, even though silicate clouds are present in the atmosphere. Small silicate particles should produce deeper absorption features, which are not observed in the mid-infrared, suggesting that the impact of the small mean particle size in the upper atmosphere does not contribute strongly to the cloud opacity.  \nDue to the di ff erent slopes in the near-infrared opacity as a function of wavelength, as shown in Figure 16, there is some sensitivity to di ff erent compositions and particle geometries. This can explain the mild preference for amorphous MgSiO3 clouds in certain like-for-like retrieval comparisons. In reverse, the lack of features leads to a preference for Na2S clouds in HR 8799 b. Performing additional tests on HR 8799 e, we find that Mg2SiO4 (e.AB.14) are mildly disfavoured by the retrievals, and Al2O3 (e.AB.24) clouds are strongly disfavoured, though this is again more likely from their condensation location rather than from the impact of aerosol spectral features. Observations of the silicate absorption features at 10 µ m would allow more precise measurement of this wavelength dependence, and in turn place better constraints on the cloud structure and composition.', '6.1. Highly enriched atmospheres': 'While some enrichment is expected in giant planets formed through core accretion, metallicities of nearly 100 × the stellar value are far beyond the expectations for planets with masses larger than that of Jupiter. Having transformed the atmospheric metallicities from the retrievals to Z pl / Z ⋆ using equations 19 and 20, we can compare the HR 8799 planets to the broader population. Figure 17 shows how the inferred metallicities of the HR 8799 system compare to those of other directly imaged planets,  \nFig. 16. Cloud absorption (left) and scattering (right) opacities for different condensate compositions and structure, for 1µ m particles. Dark blue indicates MgSiO3, light blue is Mg2SiO4. Solid (dashed) lines are for crystalline (amorphous) substances. The solid red (green) line is for crystalline iron (Na2S).  \n<!-- image -->  \nand to a fit from Thorngren et al. (2016) derived from a sample of transiting exoplanets. The HR 8799 planets are clear outliers amongst the directly imaged planets, whose metallicities were taken from the literature. We again used equations 19 and 20 to convert from a measurement of [M / H] to Z pl / Z ⋆ . Where stellar metallicities are not available, we assumed a metallicity [Fe / H] = 0. Only the 2 M Jup planet Af Lep b has a comparable degree of enrichment to the HR 8799 planets (Zhang et al. 2023). However, they are comparable to hot jupiters observed in transmission, such as WASP 39 b (Rustamkulov et al. 2023). Looper et al. (2008) and Stephens et al. (2009) demonstrate how high metallicity atmospheres facilitate condensation, in turn leading to the strong reddening seen in a subset of L dwarfs. The retrieval of highly metal rich and very cloudy atmospheres is consistent with this picture from brown dwarfs. Further supporting the high [M / H] retrievals are their consistency with the self-consistent grids, which always favour their upper limits.  \nTo validate these findings, we performed a series of test retrievals for HR 8799 e, fixing the metallicity to solar composition, and compared cases where [M / H] is fixed to values between 0.0 and 2.0 in steps of 0.5 dex. Figure 18 shows the best fit spectra from each of these retrievals, showing the clear di ff erences in the J, H, and K bands between the di ff erent metallicity cases that are unable to be compensated for by varying other atmospheric parameters. We find that the cases of [M / H] = 1.5 and 2.0 are strongly favoured over the cases between 0.0 and 1.0, with solar composition disfavoured at ∆ log Z > 10. The remaining atmospheric parameters also significantly varied between the di ff erent retrievals: the C / O ratio increases with increasing [M / H], maintaining a relatively constant abundance of H2O in the atmosphere while allowing the CO abundance to increase. This combines with the decreasing f sed, increasing the cloudiness of the planet to dampen the stronger molecular features at high metallicity. A more thorough treatment of the condensation and chemistry in the retrieval framework are likely necessary to accurately infer both of these parameters.  \nThe spectra shown in Figure 18 show that the retrieved metallicity is strongly dependant on both the height of the J, H, and K band peaks, as well as the shape of these features. However, the J bands is only covered by the low-resolution SPHERE data, with relatively poor S / N , and di ff erent instruments measure significantly di ff erent flux in the H band. Without compatible measurements in this spectral regime, robust conclusions  \nFig. 17. Metallicities of the exoplanet planet population. In blue are transiting planets, adapted from Thorngren et al. (2016) (dark blue). In red are directly imaged planets, with references listed below. Indicated with stars are the HR 8799 planets as measured in this work. The grey line indicates the fit by Thorngren et al. (2016).  \n<!-- image -->  \nNotes: β Pic b (Gravity Collaboration et al. 2020); PDS 70 b (Wang et al. 2021); 51 Eri b (Whiteford et al. 2023); VHS 1256 b (Hoch et al. 2022); HIP 65426 b (Petrus et al. 2021); κ And b (Bonnefoy et al. 2014; Wilcomb et al. 2020); YSES 1 b (Zhang et al. 2021) AB Pic b (Palma-Bifani et al. 2023); AF Lep b (Zhang et al. 2023); GJ 504 b (Bonnefoy et al. 2018); HD 95086 b (Desgrange et al. 2022); Ross 458 c (Burgasser et al. 2010); DH Tau b (Patience et al. 2012); HN Peg b (Leggett et al. 2008; Suárez et al. 2021); CT Cha b (Schmidt et al. 2008); GQ Lup b (Demars et al. 2023).  \nabout the metallicity are hard to draw. Additional constraints can be obtained from other wavelength regions; Lodders & Fegley (2002) find that the CO2 abundance scales proportionally to [M / H] 2 . With strong features between 3 and 4 µ m, as well as in the mid-infrared, future observations should be able to place robust constraints on this parameter. Section 6.4 discusses the potential of JWST to make such observations.  \nFinally, we performed independent comparisons to the moderate resolution OSIRIS spectra presented in Ru ffi o et al. (2021). Following their methods, we used the parameters of the single best-fit disequilibrium models for HR 8799 b, c, and d to compute a high spectral resolution using line-by-line opacity lists. We convolved this model to the OSIRIS instrumental spectral resolution, and binned the model to the OSIRIS wavelength grid. The resulting spectra was multiplied by the atmospheric transmission function provided by Ru ffi o et al. (2021). The continuum was measured by high-pass filtering the spectrum, and was subsequently subtracted from the model. We fit a scaling factor between the model and the OSIRIS data and computed the resulting χ 2 . This exercise was repeated, setting the metallicity of the model to solar. For HR 8799 c and b we find that the high metallicity model provides a better fit to the OSIRIS data than the solar model, with the caveat that fixing the metallicity during the retrieval may result in a better fit than setting it a posteriori, without changing other atmospheric parameters. For HR 8799 d the high metallicity model is only a marginally better fit than the solar metallicity model. In general, the fits from the GRAVITY and remaining archival data provide reasonable fits to the R ≈ 4000 OSIRIS data, of similar quality to the fits displayed in Ru ffi o et al. (2021), and are included in the appendix in Figure C.1.', '6.2. Impacts of data selection': 'Given the inhomogeneity in the data in terms of spectral resolution, S / N , observing strategy and more, we performed a series of retrievals to examine the impact of di ff erent datasets on the retrieval results. We first performed retrievals using only the GRAVITY data to determine the constraining power of this new dataset. We found that using only the GRAVITY data we could rule out clear atmosphere solutions at ∆ log Z > 7 (e.g. e.A.31 over e.A.32), using the Z23 profile and disequilibrium chemistry. Using only the GRAVITY observations for HR 8799 e, we could obtain estimates of the metallicity (2 ± 0 . 3), e ff ective temperature (1143 + 38 -32 K), and C / O ratio (0 . 71 + 0 . 08 -0 . 2 ), which are broadly consistent with the results from the combined dataset. We obtain similarly reliable estimates for b (b.A.32), but using only the GRAVITY data for d (d.A.12) and c (c.A.11), resulted in significantly lower metallicities. All of these retrievals use the same setup of the Z23 profile, disequilibrium chemistry, and silicate and iron clouds condensing at their equilibrium saturation point. Given the higher spectral resolution and S / N of the GRAVITY data, this leads to the conclusion that solutions to the full retrievals are largely driven by the fits to the GRAVITY spectra. The C / O ratios for d and b are also incompatible when using only the GRAVITY data, finding substellar values for both planets.  \nConversely, we also performed retrievals that exclude the GRAVITY spectra. For HR 8799 e (e.A.33) and d (d.A.13), we find that the retrieved parameters again broadly agree with the retrieval including the GRAVITY data. Thus even if the retrievals are dominated by the GRAVITY data, the conclusions we draw are robust even when excluding the GRAVITY spectra.  \nFig. 18. Comparison retrievals of HR 8799 e where the metallicity is fixed to a constant value. Black points denote the equivalent retrieval where the metallicity is a free parameter. The top panel shows how the best fit spectra di ff er across the set of retrievals, while the bottom row shows how di ff erent atmospheric parameters vary as a function of metallicity.  \n<!-- image -->', '6.3. Formation': 'The formation mechanism of the HR 8799 planets has seen much debate since their discovery: simply put, how can one form four such massive planets in the same system? The C / O ratio and metallicity are the best formation tracers observed to date in these planets. Nevertheless, neither gravitational instability nor core accretion scenarios have been ruled out. The high degree of enrichment in these objects would seem to suggest a core accretion formation scenario. However, Wang (2023) finds that even for [M / H] = 0.5, approximately 100 earth masses of solids are required to enrich the atmosphere of HR 8799 e. With the even higher metallicities measured in this work, this number would increase, and it is unclear if it is possible to have nearly 1000 earth masses of metals available in a protoplanetary disc and accreted with high e ffi ciency. Wang (2023) additionally find that late accretion of planetesimals would require less material to result in a similar degree of enrichment, potentially only requiring a few hundred earth masses of solids to achieve the high metallicities of all four planets, though this material should quickly settle out of the atmosphere.  \nAll four companions have stellar C / O ratios or higher, and the C / O ratio in this system varies with separation, decreasing from b to d, before a sharp increase in the C / O ratio for the innermost planet. Super-stellar C / O ratios have been tied to core-accretion formation together with pebble drift and evaporation (Schneider & Bitsch 2021b,a; Mollière et al. 2022). A pathway to significant metal enrichment was found by Bitsch & Mah (2023), though it predicts that the CO rich pebbles should evaporate near the CO iceline, which is outside the radius of even HR 8799 b. Planetesimal accretion cannot be ruled out either: if large amounts of solids, with near-stellar composition,  \nare accreted, more metal-rich planets will have C / O ratios that approach the stellar value. This is consistent with the trends in metallicity and C / O between the b, c, and d planets, though e remains an exception. Such a transition could be explained by the outer three planets, particularly d, trapping water ice and preventing these solids from reaching the innermost planet. Some combination of these mechanisms could explain both the atmospheric enrichment and the trends in the C / O ratio: for example, early enrichment from evaporating pebbles could lead to the high planetary metallicities, while late accretion of planetesimals could then drive the C / O ratio down towards the solar value. Alternatively, Chen et al. (2024) demonstrate that the opening of gaps in a protoplanetary disc can significantly alter the composition of the gas and ices available to accrete onto forming planets, and it seems likely that substructure induced by the four HR 8799 planets would strongly impact their eventual composition.  \nWith e ff ective temperatures well over 1000 K and radii significantly larger than that of jupiter, the coldest initial condition scenarios of Marley et al. (2007) can be excluded. Beyond this constraint, the masses, luminosities, and radii seem largely consistent with a broad range of potential evolutionary tracks (e.g. Bara ff e et al. 2003; Saumon & Marley 2008; Mordasini et al. 2017). Further work to measure more formation tracers is clearly necessary to unravel this system. Mid-infrared spectroscopic observations could characterise the NH3 abundance, at least in HR 8799 b, as well as the debris disc observed in Boccaletti et al. (2023). Higher spectral resolution could enable the measurement of carbon isotopes (e.g. as well as place better constraints on the metallicity of each of these atmospheres). Such measurements, combined with dedicated formation models of giant planets outside the water iceline, are necessary to  \nFig. 19. Predictions for NIRSpec / G395H based on most favoured disequilibrium (solid) and free chemistry (dashed) retrievals, together with the best-fit self-consistent models from each grid.  \n<!-- image -->  \ndetermine whether these four planets share a formation pathway, or whether there were di ff erent mechanisms impacting di ff erent regions of the protoplanetary disc.', '6.4. Predictions for JWST': 'The high spectral resolution modes of JWST may allow us to verify the measurements made in this work, particularly through observation of CO2, CO, and CH4 features in the near infrared and the silicate absorption features near 10 µ m. To this end, we present the range of model predictions in these wavelength regions at the spectral resolution of the JWST instruments, for NIRSpec in figure 19 and for the MIRI / MRS in figure 20. The comparisons to the MIRI photometry shown in figure 20 demonstrate the typical degree of compatibility between the models and the data in these mid-infrared wavelengths. In the NIRSpec / G395H wavelength range we see significant discrepancy between models for the same planet in the amplitude of the CO2 feature at 3.8 µ m, as well as the CH4 feature at 3.3 µ m and the CO lines between 4.5 and 5 µ m. These observations will also be able to confirm the presence of HCN in the atmospheres of HR 8799 c and e. Precise measurement of these features should provide robust constraints on the metallicity of these objects, verifying the degree of enrichment found via the ground-based observations. While silicate clouds are preferred in the retrieval comparison, none of the models show signs of deep silicate absorption features near 10 µ m, but spectroscopic observations are required to validate these models. Mid-infrared observations will be particularly valuable for HR 8799 b, and will allow the clear detection of ammonia. If combined with a chemical model to determine the ratios of NH3:HCN:N2, this will allow for the measurement of the N / O ratio, which can also be used as a formation diagnostic (Turrini et al. 2021; Pacetti et al. 2022). Recently  \nFig. 20. Predictions for MIRI / MRSbased on the most favoured disequilibrium (solid) and free chemistry (dashed) retrievals, together with the best-fit self-consistent models from each grid.  \n<!-- image -->  \nRu ffi o et al. (2023) demonstrated the potential for high-contrast imaging with NIRSpec; even without the use of a coronagraph it should be possible to obtain flux calibrated, moderate resolution spectroscopy of HR 8799 b, c, and d using the NIRSpec IFU through a combination of forward modelling and reference differential imaging. In the case where the planet signal is unable to be separated from that of the host star, Patapis et al. (2022) demonstrated that it will be possible to at least identify trace species through molecular mapping in the mid-infrared, though this will be unable to characterise the broad wavelength features of the silicate clouds.', '7. Conclusions': "After more than 15 years of study, the HR 8799 planets remain mysterious, though increasing data quality is allowing us to peer deeper into these atmospheres than ever before. We present new K-band spectra from the VLTI / GRAVITY, which together with a large set of archival data form the basis of the our atmospheric analysis. Using petitRADTRANS retrievals and fits to self-consistent grids, we inferred the atmospheric properties of all four companions, with reasonable agreement between the two methods. Our results are broadly consistent with the literature in terms of e ff ective temperature, mass, surface gravity, and radius for all four planets. The use the dynamical mass as a prior in the retrievals when determining log g allows us to reliably retrieve physically reasonable planet radii.  \nWe find that all four planets are strongly enriched in metals, though there is still discrepancy between di ff erent models in constraining the precise value. This was validated by running retrievals using di ff erent temperature profiles and chemical models, and comparing to self-consistent grids. Further selfconsistent modelling is necessary, particularly to extend model  \ngrids out to high metallicities. The C / O ratio is stellar to superstellar for all four planets. It decreases from the outermost planet to HR 8799 d, while HR 8799 e has a higher C / O ratio than the other companions. We confidently detect HCN in HR 8799 c and e, at abundances far higher than predicted by equilibrium chemistry; though this detection is largely driven by lowS / N data from LBT / ALES. CH4 is also confidently detected in HR 8799 c for the first time. From the disequilibrium chemistry retrievals, H2S appears to be a highly abundant species in all of the planets, but higher S / N and spectral resolution are required for a confident detection in a free retrieval framework. Using our retrieved quench pressure and chemical abundances, we are able to derive a vertical mixing strength, finding K zz values compatible with high-metallicity predictions from Soni & Acharyya (2023). The mixing strength is stronger for the warmer planets, at log K zz ≈ 6, and is lower for HR 8799 b with log K zz ≈ 2.  \nAll of the planets are highly cloudy. For the inner three planets, these clouds are composed of silicate clouds lying above the photosphere, and deep, dense iron clouds forming the base of the photosphere. Cooler than the other three planets, the most favoured model for HR 8799 b requires Na2S clouds. All of the planets have e ff ective temperatures consistent with literature values, with HR 8799 b still unique in its lower temperature and mass compared to its siblings.  \nWe emphasise the use of robust model comparison in this work: while it may be di ffi cult to present precise measurements of certain properties, the use of multiple methods and models allows us to draw a robust portrait of each of these atmospheres. We also note that our conclusions rely on data with significant incompatibilities, particularly in the H-band flux. While we performed extensive analysis to mitigate the influence of any individual dataset, further observations are required to obtain reliable spectroscopic measurements in the near-infrared. While the HR8799 planets share many similarities, much like our own Solar System there are di ff erences in their atmospheric properties, which require further study.  \nAcknowledgements. We would like to thank Quinn Konopacky, Alice Zurlo, Alex Greenbaum, Beth Biller, Olivier Flasseur, David Doelman, and Pengyu Lui for providing the archival datasets used in this work. We are also grateful to Gilles Loupe and Malavika Vasist for providing helpful suggestions on how to implement the spline temperature profile. Thanks as well to Daniel Thorngren for providing the model for calculating the planet metallicity Z . Finally, we are grateful to our anonymous reviewer for their thorough and insightful report. Software used: petitRADTRANS, pyKLIP, species, VIP-HCI, pyMultiNest, phot\\_utils, Python, numpy, matplotlib, astropy, sympy, Aspro . SL acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-21-CE31-0017 (project ExoVLTI). J.J.W., A.C., and S.B. acknowledge the support of NASA XRP award 80NSSC23K0280. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 'Exploring the Diversity of Extrasolar Planets' (MA 9185 / 1) amd from the Swiss National Science Foundation under grant 200021\\_204847 'PlanetsInTime'. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. This work is based on observations collected at the European Southern Observatory under ESO programme 1104.C-0651. This publication makes use of VOSA, developed under the Spanish Virtual Observatory project supported by the Spanish MINECO through grant AyA2017-84089. VOSA has been partially updated by using funding from the European Union's Horizon 2020 Research and Innovation Programme, under Grant Agreement nº 776403 (EXOPLANETS-A)", 'References': "Ackerman, A. S. & Marley, M. S. 2001, ApJ, 556, 872 Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A.  \n- 2001, ApJ, 556, 357\n- Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765\n- Allard, N. F., Spiegelman, F., & Kielkopf, J. F. 2016, A&A, 589, A21\n- Allard, N. F., Spiegelman, F., Leininger, T., & Mollière, P. 2019, A&A, 628, A120\n- Amara, A. & Quanz, S. P. 2012, MNRAS, 427, 948\n- Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40\n- Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481\n- Azzam, A. A. A., Yurchenko, S. N., Tennyson, J., & Naumenko, O. V. 2016, MNRAS, 460, 4063\n- Baines, E. K., White, R. J., Huber, D., et al. 2012, ApJ, 761, 57\n- Balmer, W. O., Pueyo, L., Stolker, T., et al. 2023, ApJ, 956, 99\n- Bara ff e, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701\n- Barber, R. J., Strange, J. K., Hill, C., et al. 2014, MNRAS, 437, 1828\n- Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61\n- Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65\n- Barton, E. J., Yurchenko, S. N., & Tennyson, J. 2013, MNRAS, 434, 1469\n- Baudino, J. L., Bézard, B., Boccaletti, A., et al. 2015, A&A, 582, A83\n- Baudino, J.-L., Mollière, P., Venot, O., et al. 2017, ApJ, 850, 150\n- Bayo, A., Rodrigo, C., Barrado Y Navascués, D., et al. 2008, A&A, 492, 277\n- Bell, T. J., Welbanks, L., Schlawin, E., et al. 2023, Nature, 623, 709\n- Benneke, B. & Seager, S. 2012, ApJ, 753, 100\n- Benneke, B. & Seager, S. 2013, ApJ, 778, 153\n- Bergfors, C., Brandner, W., Janson, M., Köhler, R., & Henning, T. 2011, A&A, 528, A134\n- Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Proc. SPIE, 701418\n- Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155\n- Biller, B. A., Apai, D., Bonnefoy, M., et al. 2021, MNRAS, 503, 743\n- Bitsch, B. & Mah, J. 2023, A&A, 679, A11\n- Blain, D., Charnay, B., & Bézard, B. 2021, A&A, 646, A15\n- Boccaletti, A., M., M., P., B., et al. 2023, arXiv e-prints, arXiv:2310.13414 Bodenheimer, P. 1974, Icarus, 23, 319\n- Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2\n- Bonnefoy, M., Currie, T., Marleau, G. D., et al. 2014, A&A, 562, A111\n- Bonnefoy, M., Perraut, K., Lagrange, A. M., et al. 2018, A&A, 618, A63  \nBonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58  \nBoss, A. P. 1997, Science, 276, 1836  \n- Bowler, B. P., Liu, M. C., Dupuy, T. J., & Cushing, M. C. 2010, ApJ, 723, 850 Brandt, G. M., Brandt, T. D., Dupuy, T. J., Michalik, D., & Marleau, G.-D. 2021, ApJ, 915, L16  \nBrown-Sevilla, S. B., Maire, A. L., Mollière, P., et al. 2023, A&A, 673, A98 Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125 Burgasser, A. J., Simcoe, R. A., Bochanski, J. J., et al. 2010, ApJ, 725, 1405 Burningham, B., Faherty, J. K., Gonzales, E. C., et al. 2021, MNRAS, 506, 1944 Burningham, B., Marley, M. S., Line, M. R., et al. 2017, MNRAS, 470, 1177 Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47 Burrows, A. & Sharp, C. M. 1999, ApJ, 512, 843 Burrows, A., Sudarsky, D., & Hubeny, I. 2006, ApJ, 640, 1063 Calamari, E., Faherty, J. K., Visscher, C., et al. 2024, ApJ, 963, 67 Cameron, A. G. W. 1978, Moon and Planets, 18, 5 Cannon, A. J. & Pickering, E. C. 1993, VizieR Online Data Catalog, III / 135A Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 Chabrier, G., Bara ff e, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464 Charnay, B., Bézard, B., Baudino, J.-L., et al. 2018, ApJ, 854, 172 Chen, K., Kama, M., Pinilla, P., & Keyte, L. 2024, MNRAS, 527, 2049 Chubb, K. L., Rocchetto, M., Yurchenko, S. N., et al. 2021, A&A, 646, A21 Coles, P. A., , Yurchenko, S. N., & Tennyson, J. 2019, MNRAS, 490, 4638 Contro, B., Wittenmyer, R. A., Horner, J., & Marshall, J. P. 2015, Origins of Life and Evolution of the Biosphere, 45, 41  \n- Covarrubias, S., Blunt, S., & Wang, J. J. 2022, Research Notes of the American Astronomical Society, 6, 66  \nCridland, A. J., Eistrup, C., & van Dishoeck, E. F. 2019, A&A, 627, A127  \nCridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2020, A&A,  \n642, A229  \n- Currie, T., Burrows, A., Girard, J. H., et al. 2014, apj, 795, 133\n- Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128\n- Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., & Plavchan, P. 2012, ApJ, 755, L34  \nCushing, M. C., Marley, M. S., Saumon, D., et al. 2008, ApJ, 678, 1372  \nCushing, M. C., Roellig, T. L., Marley, M. S., et al. 2006, ApJ, 648, 614  \n- Demars, D., Bonnefoy, M., Dougados, C., et al. 2023, A&A, 676, A123\n- Desgrange, C., Chauvin, G., Christiaens, V., et al. 2022, A&A, 664, A139\n- Dodson-Robinson, S. E., Veras, D., Ford, E. B., & Beichman, C. A. 2009, ApJ, 707, 79\n- Doelman, D. S., Stone, J. M., Briesemeister, Z. W., et al. 2022, AJ, 163, 217 Eddington, A. S. 1930, MNRAS, 90, 279\n- Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14  \nEsposito, S., Mesa, D., Skemer, A., et al. 2013, A&A, 549, A52  \nFaherty, J. K. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Springer, Cham), 188  \nFaramaz, V., Marino, S., Booth, M., et al. 2021, AJ, 161, 271  \nFegley, Bruce, J. & Lodders, K. 1996, ApJ, 472, L37  \n- Feroz, F. & Hobson, M. P. 2008, MNRAS, 384, 449\n- Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10\n- Flasseur, O., Denis, L., Thiébaut, E., & Langlois, M. 2018, A&A, 618, A138  \nFlasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020, A&A, 637, A9  \n- Fonte, S., Turrini, D., Pacetti, E., et al. 2023, MNRAS, 520, 4683\n- Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., & Marley, M. S. 2011, ApJ, 729, 32\n- Fragoso, T. M., Bertoli, W., & Louzada, F. 2018, International Statistical Review, 86, 1\n- Fukagawa, M., Itoh, Y., Tamura, M., et al. 2009, ApJ, 696, L1  \nGaia Collaboration. 2020, VizieR Online Data Catalog, I  \n/  \n350  \n- Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJ, 739, L41  \nGehren, T. 1977, A&A, 59, 303  \nGeiler, F., Krivov, A. V., Booth, M., & Löhne, T. 2018, MNRAS, 483, 332  \nGiacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205  \nGo'zdziewski, K. & Migaszewski, C. 2018, ApJS, 238, 6  \nGravity Collaboration, Abuter, R., Accardo, M., et al. 2017, A&A, 602, A94 GRAVITY Collaboration, Bauböck, M., Dexter, J., et al. 2020, A&A, 635, A143 Gravity Collaboration, Lacour, S., Nowak, M., et al. 2019, A&A, 623, L11 Gravity Collaboration, Nowak, M., Lacour, S., et al. 2020, A&A, 633, A110 Gray, R. O., Corbally, C. J., Garrison, R. F., McFadden, M. T., & Robinson, P. E.  \n2003, AJ, 126, 2048  \nGray, R. O. & Kaye, A. B. 1999, AJ, 118, 2993  \nGreco, J. P. & Brandt, T. D. 2016, ApJ, 833, 134  \nGreenbaum, A. Z., Pueyo, L., Ru  \nffi  \no, J.-B., et al. 2018, AJ, 155, 226  \nGro ff , T., Chilcote, J., Brandt, T., et al. 2017, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10400, Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, ed. S. Shaklan, 1040016  \n- Gro ff , T. D., Kasdin, N. J., Limbach, M. A., et al. 2015, in Techniques and Instrumentation for Detection of Exoplanets VII, ed. S. Shaklan, Vol. 9605 (SPIE), 96051C\n- Guillot, T. 2005, Annual Review of Earth and Planetary Sciences, 33, 493 Guillot, T. 2010, A&A, 520, A27\n- Hansen, J. E. 1971, Journal of Atmospheric Sciences, 28, 1400  \nHelled, R. & Bodenheimer, P. 2010, Icarus, 207, 503  \nHoch, K. K. W., Konopacky, Q. M., Barman, T. S., et al. 2022, AJ, 164, 155  \nHoch, K. K. W., Konopacky, Q. M., Theissen, C. A., et al. 2023, AJ, 166, 85  \nIngraham, P., Marley, M. S., Saumon, D., et al. 2014, ApJ, 794, L15  \n- Janson, M., Bergfors, C., Goto, M., Brandner, W., & Lafrenière, D. 2010, ApJ, 710, L35\n- Jullion, A. & Lambert, P. 2007, Computational Statistics & Data Analysis, 51, 2542  \nKaralidi, T., Marley, M., Fortney, J. J., et al. 2021, ApJ, 923, 269  \n- Kass, R. E. & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 773\n- Kaye, A. B., Handler, G., Krisciunas, K., Poretti, E., & Zerbi, F. M. 1999, PASP, 111, 840\n- Kirkpatrick, J. D., Reid, I. N., Liebert, J., et al. 1999, ApJ, 519, 802  \nKitzmann, D., Heng, K., Oreshenko, M., et al. 2020, ApJ, 890, 174  \n- Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398\n- Lacour, S., Wang, J. J., Nowak, M., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11446, Optical and Infrared Interferometry and Imaging VII, ed. P. G. Tuthill, A. Mérand, & S. Sallum, 114460O  \nLafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148 Lang, S. & Brezger, A. 2004, Journal of Computational and Graphical Statistics, 13, 183  \n- Lapeyrère, V., Kervella, P., Lacour, S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9146, Optical and Infrared Interferometry IV, ed. J. K. Rajagopal, M. J. Creech-Eakman, & F. Malbet, 91462D\n- Larkin, J., Barczys, M., Krabbe, A., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6269, Ground-based and Airborne Instrumentation for Astronomy, ed. I. S. McLean & M. Iye, 62691A  \nLaughlin, G. & Rozyczka, M. 1996, ApJ, 456, 279 Lavie, B., Mendonça, J. M., Mordasini, C., et al. 2017, AJ, 154, 91 Lee, J.-M., Heng, K., & Irwin, P. G. J. 2013, ApJ, 778, 97 Leggett, S. K., Saumon, D., Albert, L., et al. 2008, ApJ, 682, 1256 Linder, E. F., Mordasini, C., Mollière, P., et al. 2019, A&A, 623, A85  \n- Line, M. R., Teske, J., Burningham, B., Fortney, J. J., & Marley, M. S. 2015, ApJ, 807, 183\n- Liu, P., Bohn, A. J., Doelman, D. S., et al. 2023, A&A, 674, A115 Lodders, K. 2003, ApJ, 591, 1220\n- Lodders, K. & Fegley, B. 2002, Icarus, 155, 393\n- Lodders, K., Palme, H., & Gail, H. P. 2009, Landolt B&ouml;rnstein, 4B, 712  \nLooper, D. L., Kirkpatrick, J. D., Cutri, R. M., et al. 2008, ApJ, 686, 528  \nLuna, J. L. & Morley, C. V. 2021, ApJ, 920, 146  \n- Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proceedings of the National Academy of Science, 111, 12661\n- Madhusudhan, N. 2019, Annu. Rev. Astron. Astrophys., 57, 617\n- Madhusudhan, N., Burrows, A., & Currie, T. 2011, ApJ, 737, 34  \nMadhusudhan, N., Sarkar, S., Constantinou, S., et al. 2023, ApJ, 956, L13  \n- Madhusudhan, N. & Seager, S. 2009, ApJ, 707, 24\n- Maire, A. L., Skemer, A. J., Hinz, P. M., et al. 2015, A&A, 576, A133  \nMarleau, G.-D., Klahr, H., Kuiper, R., & Mordasini, C. 2017, ApJ, 836, 221  \nMarleau, G.-D., Mordasini, C., & Kuiper, R. 2019, ApJ, 881, 144  \n- Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541\n- Marley, M. S., Saumon, D., Cushing, M., et al. 2012, ApJ, 754, 135\n- Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJ, 723, L117\n- Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85\n- Marocco, F., Day-Jones, A. C., Lucas, P. W., et al. 2014, MNRAS, 439, 372\n- Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348\n- Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080\n- Matthews, B., Kennedy, G., Sibthorpe, B., et al. 2014, ApJ, 780, 97\n- McCarthy, A. M., Muirhead, P. S., Tamburo, P., et al. 2024, arXiv e-prints, arXiv:2402.15001\n- McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019, MNRAS, 488, 2836\n- McKemmish, L. K., Yurchenko, S. N., & Tennyson, J. 2016, MNRAS, 463, 771 Metchev, S., Marois, C., & Zuckerman, B. 2009, ApJ, 705, L204  \nMiles, B. E., Biller, B. A., Patapis, P., et al. 2023, ApJ, 946, L6  \n- Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74\n- Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131\n- Mollière, P., van Boekel, R., Bouwman, J., et al. 2017, A&A, 600, A10\n- Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47  \nMollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67 Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111 Mordasini, C., Marleau, G.-D., & Mollière, P. 2017, A&A, 608, A72 Morley, C. V., Mukherjee, S., Marley, M. S., et al. 2024, arXiv e-prints, arXiv:2402.00758  \n- Moses, J. I., Madhusudhan, N., Visscher, C., & Freedman, R. S. 2013, ApJ, 763, 25\n- Moses, J. I., Marley, M. S., Zahnle, K., et al. 2016, ApJ, 829, 66  \nMoya, A., Amado, P. J., Barrado, D., et al. 2010, MNRAS, 405, L81 Moya, A., Amado, P. J., Barrado, D., et al. 2010, MNRAS, 406, 566 Nasedkin, E., Mollière, P., Wang, J., et al. 2023, A&A, 678, A41 Nasedkin, E., Mollière, P., & Blain, D. 2024, Journal of Open Source Software,  \n- 9, 5875\n- Nero, D. & Bjorkman, J. E. 2009, ApJ, 702, L163\n- Nixon, M. C., Welbanks, L., McGill, P., & Kempton, E. M. R. 2023, arXiv eprints, arXiv:2310.03713  \nNowak, M., Lacour, S., Lagrange, A. M., et al. 2020, A&A, 642, L2 Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 Oppenheimer, B. R., Baranec, C., Beichman, C., et al. 2013, ApJ, 768, 24 Pacetti, E., Turrini, D., Schisano, E., et al. 2022, ApJ, 937, 36 Palma-Bifani, P., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, A90 Patapis, P., Nasedkin, E., Cugno, G., et al. 2022, A&A, 658, A72 Patience, J., King, R. R., De Rosa, R. J., et al. 2012, A&A, 540, A85 Perri, F. & Cameron, A. G. W. 1974, Icarus, 22, 416 Petit dit de la Roche, D. J. M., van den Ancker, M. E., Kissler-Patig, M., Ivanov,  \nV. D., & Fedele, D. 2019, MNRAS, 491, 1795  \nPetrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59  \nPetrus, S., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, L9  \nPhillips, M. W., Tremblin, P., Bara  \nff  \ne, I., et al. 2020, A&A, 637, A38  \nPodolak, M. 2003, Icarus, 165, 428  \n- Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, ıcarus, 124, 62  \nPolyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597  \nPueyo, L. 2016, ApJ, 824, 117  \n- Radigan, J. 2014, ApJ, 797, 120\n- Radigan, J., Lafrenière, D., Jayawardhana, R., & Artigau, E. 2014, ApJ, 793, 75 Rahman, M. M. 2005, PhD thesis, University of British Columbia\n- Rajan, A., Barman, T., Soummer, R., et al. 2015, ApJ, 809, L33\n- Rhee, J. H., Song, I., Zuckerman, B., & McElwain, M. 2007, ApJ, 660, 1556 Rossow, W. B. 1978, Icarus, 36, 1\n- Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spectr. Rad. Transf., 111, 2139  \nRowland, M. J., Morley, C. V., & Line, M. R. 2023, ApJ, 947, 6 Ru ffi o, J.-B., Konopacky, Q. M., Barman, T., et al. 2021, AJ, 162, 290 Ru ffi o, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200 Ru ffi o, J.-B., Perrin, M. D., Hoch, K. K. W., et al. 2023, arXiv e-prints, arXiv:2310.09902  \n- Rustamkulov, Z., Sing, D. K., Mukherjee, S., et al. 2023, Nature, 614, 659 Sadakane, K. 2006, PASJ, 58, 1023\n- Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57\n- Saumon, D. & Marley, M. S. 2008, ApJ, 689, 1327\n- Schmidt, T. O. B., Neuhäuser, R., Seifahrt, A., et al. 2008, A&A, 491, 311\n- Schneider, A. D. & Bitsch, B. 2021a, A&A, 654, A71\n- Schneider, A. D. & Bitsch, B. 2021b, A&A, 654, A72  \nSepulveda, A. G. & Bowler, B. P. 2022, AJ, 163, 52  \n- Sepulveda, A. G., Huber, D., Li, G., et al. 2023, Research Notes of the American  \nAstronomical Society, 7, 2  \nShowman, A. P. & Kaspi, Y. 2013, ApJ, 776, 85  \n- Skemer, A. J., Hinz, P., Montoya, M., et al. 2015, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9605, Techniques and Instrumentation for Detection of Exoplanets VII, ed. S. Shaklan, 96051D  \nSkemer, A. J., Hinz, P. M., Esposito, S., et al. 2012, ApJ, 753, 14  \n- Skemer, A. J., Marley, M. S., Hinz, P. M., et al. 2014, ApJ, 792, 17  \nSkrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163  \nSoni, V. & Acharyya, K. 2023, ApJ, 946, 29  \nSoummer, R., Hagan, J. B., Pueyo, L., et al. 2011, ApJ, 741, 55  \n- Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28\n- Sousa-Silva, C., Al-Refaie, A. F., Tennyson, J., & Yurchenko, S. N. 2015, MNRAS, 446, 2337  \nStephens, D. C., Leggett, S. K., Cushing, M. C., et al. 2009, ApJ, 702, 154  \nStolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59  \nStolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182  \nSu, K. Y. L., Rieke, G. H., Stapelfeldt, K. R., et al. 2009, ApJ, 705, 314  \n- Suárez, G., Vos, J. M., Metchev, S., Faherty, J. K., & Cruz, K. 2023, ApJ, 954, L6\n- Suárez, G. & Metchev, S. 2022, MNRAS, 513, 5701\n- Suárez, G. & Metchev, S. 2023, MNRAS, 523, 4739\n- Suárez, G., Metchev, S., Leggett, S. K., Saumon, D., & Marley, M. S. 2021, ApJ, 920, 99  \nSzulágyi, J. & Mordasini, C. 2017, MNRAS, 465, L64  \n- Tan, X. & Showman, A. P. 2021a, MNRAS, 502, 678\n- Tan, X. & Showman, A. P. 2021b, MNRAS, 502, 2198\n- Tennyson, J. & Yurchenko, S. N. 2012, MNRAS, 425, 21\n- Thompson, W., Marois, C., Do Ó, C. R., et al. 2023, AJ, 165, 29\n- Thorngren, D. & Fortney, J. J. 2019, ApJ, 874, L31\n- Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64\n- Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19\n- Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17\n- Tremblin, P., Chabrier, G., Bara ff e, I., et al. 2017, ApJ, 850, 46  \nTremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144  \n- Turrini, D., Schisano, E., Fonte, S., et al. 2021, ApJ, 909, 40\n- Vos, J. M., Biller, B. A., Bonavita, M., et al. 2019, MNRAS, 483, 480\n- Vos, J. M., Burningham, B., Faherty, J. K., et al. 2023, ApJ, 944, 138  \nWahhaj, Z., Milli, J., Romero, C., et al. 2021, A&A, 648, A26  \nWaldmann, I. P., Rocchetto, M., Tinetti, G., et al. 2015, apj, 813, 13  \n- Wang, J. 2023, arXiv e-prints, arXiv:2310.00088\n- Wang, J., Mawet, D., Fortney, J. J., et al. 2018a, AJ, 156, 272\n- Wang, J., Wang, J. J., Ma, B., et al. 2020, AJ, 160, 150\n- Wang, J., Wang, J. J., Ru ffi o, J.-B., et al. 2023, AJ, 165, 4\n- Wang, J. J., Gao, P., Chilcote, J., et al. 2022, AJ, 164, 143\n- Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263\n- Wang, J. J., Graham, J. R., Dawson, R., et al. 2018b, AJ, 156, 192\n- Wang, J. J., Ru ffi o, J.-B., Morris, E., et al. 2021, AJ, 162, 148\n- Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148\n- Wende, S., Reiners, A., Seifahrt, A., & Bernath, P. F. 2010, A&A, 523, A58 Whiteford, N., Glasse, A., Chubb, K. L., et al. 2023, MNRAS\n- Wilcomb, K. K., Konopacky, Q. M., Barman, T. S., et al. 2020, AJ, 160, 207  \nWilner, D. J., MacGregor, M. A., Andrews, S. M., et al. 2018, ApJ, 855, 56  \nWright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868  \n- Yurchenko, S. N., Amundsen, D. S., Tennyson, J., & Waldmann, I. P. 2017, A&A, 605, A95  \nYurchenko, S. N., Mellor, T. M., Freedman, R. S., & Tennyson, J. 2020, MNRAS, 496, 5282  \nZahnle, K. J. & Marley, M. S. 2014, ApJ, 797, 41 Zerbi, F. M., Rodríguez, E., Garrido, R., et al. 1999, MNRAS, 303, 275 Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021, Nature, 595, 370 Zhang, Z., Mollière, P., Hawkins, K., et al. 2023, AJ, 166, 198 Zuckerman, B., Rhee, J. H., Song, I., & Bessell, M. S. 2011, ApJ, 732, 61 Zuckerman, B. & Song, I. 2004, ARA&A, 42, 685 Zurlo, A., Go'zdziewski, K., Lazzoni, C., et al. 2022, A&A, 666, A133 Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57  \n1  \nMax-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany  \n- 2 LESIA, Observatoire de Paris, PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Janssen, 92195 Meudon, France\n- 3 European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany\n- 4 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom\n- 5 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands\n- 6 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA\n- 7 Department of Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA\n- 8 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 9 Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstraße 1, 85748 Garching, Germany\n- 10 Universidade de Lisboa - Faculdade de Ciências, Campo Grande, 1749-016 Lisboa, Portugal\n- 11 CENTRA - Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049-001 Lisboa, Portugal\n- 12 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France\n- 13 Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France\n- 14 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France\n- 15 STAR Institute, Université de Liège, Allée du Six Août 19c, 4000 Liège, Belgium\n- 16 Department of Astrophysical & Planetary Sciences, JILA, Duane Physics Bldg., 2000 Colorado Ave, University of Colorado, Boulder, CO 80309, USA\n- 17 1. Institute of Physics, University of Cologne, Zülpicher Straße 77, 50937 Cologne, Germany\n- 18 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany\n- 19 Universidade do Porto, Faculdade de Engenharia, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal\n- 20 School of Physics, University College Dublin, Belfield, Dublin 4, Ireland\n- 21 Departments of Physics and Astronomy, Le Conte Hall, University of California, Berkeley, CA 94720, USA\n- 22 European Southern Observatory, Casilla 19001, Santiago 19, Chile\n- 23 Advanced Concepts Team, European Space Agency, TEC-SF, ESTEC, Keplerlaan 1, NL-2201, AZ Noordwijk, The Netherlands\n- 24 Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ, USA\n- 25 University of Exeter, Physics Building, Stocker Road, Exeter EX4 4QL, United Kingdom\n- 26 Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany\n- 27 Institüt für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany\n- 28 Physikalisches Institut, Universität Bern, Gesellschaftsstr. 6, 3012 Bern, Switzerland\n- 29 Astronomy Department, University of Michigan, Ann Arbor, MI 48109 USA\n- 30 Academia Sinica, Institute of Astronomy and Astrophysics, 11F AstronomyMathematics Building, NTU / AScampus, No. 1, Section 4, Roosevelt Rd., Taipei 10617, Taiwan\n- 31 European Space Agency, ESA O ffi ce, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 32 Department of Astronomy & Astrophysics, University of California, San Diego, La Jolla, CA 92093, USA\n- 33 Department of Earth & Planetary Sciences, Johns Hopkins University, Baltimore, MD 21218, USA\n- 34 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany\n- 35 Excellence Cluster ORIGINS, Boltzmannstraße 2, 85748 Garching bei München, Germany", 'Appendix A: Data logs': 'Table A.1. Near infrared stellar Photometry of HR8799, using apparent flux normalised to 10 pc, retrieved from the Spanish Virtual Observatory (Bayo et al. 2008).', 'Appendix B: Reprocessing SPHERE and GPI datasets': 'To resolve the known discrepancies in the H-band flux between the archival SPHERE and GPI datasets, we reprocessed each using KLIP (Soummer et al. 2012; Pueyo 2016), ANDROMEDA (Cantalloube et al. 2015) and PynPoint (Amara & Quanz 2012; Stolker et al. 2019). We optimised the choice of algorithm parameters through a series of injection / extraction tests into each dataset, as in Nasedkin et al. (2023). Using two di ff erent goodness-of-fit metrics on injections representative of the true companion contrast and separation, we choose the number of principal components used in the PSF subtraction in order to extract the companion spectra with minimal bias. Figures B.1 and B.2 shows the extracted spectra for each of HR 8799 c, d, and e, for SPHERE and GPI respectively, compared to published literature spectra from Zurlo et al. (2016), Greenbaum et al. (2018), and Flasseur et al. (2018). As the goodness-of-fit metrics favoured the KLIP extractions, we used these as the basis of our retrieval analysis.  \nFig. B.1. KLIP and ANDROMEDA extractions from SPHERE for HR 8799 c, d, and e compared to the spectra published in Zurlo et al. (2016) and Flasseur et al. (2018).  \n<!-- image -->', 'Appendix C: Retrieval validation': 'Extensive validation of the pRT retrieval module was performed as part of this work. Following updates described in Nasedkin  \nTable A.2. Photometric data for HR 8799 bcde.', 'Notes:': 'References: C12 Currie et al. (2012); M10 Marois et al. (2010); G11 Galicher et al. (2011); C14 Currie et al. (2014); Z16 Zurlo et al. (2016); C11 Currie et al. (2011); B23 Boccaletti et al. (2023).  \nTable A.3. Spectroscopic Observation Log.', 'Appendix D: Using the Hansen distribution with EDDYSED': 'The E ddy S ed cloud model from Ackerman & Marley (2001) is implemented in pRT , and is the most physically motivated model incorporated to date. Typically, it assumes a log-normal particle size distribution, where the geometric particle radius will vary throughout the atmosphere as a function of the vertical di ff usion coe ffi cient K zz and the sedimentation fraction f sed. Here, we substitute the log-normal particle size distribution with the Hansen distribution, originally introduced in Hansen (1971), and rederive the calculation for the particle radius as a function of K zz and f sed.  \nWe begin with a review of the E ddy S ed model: the distribution of the number of particles as a function of particle radius, n ( r ) is approximated as a log-normal distribution with width σ g and characteristic geometric radius rg .  \nn ( r ) = N r √ 2 π log σ g exp         -log 2 GLYPH<16> r / rg GLYPH<17> 2 log 2 σ g         , (D.1)  \nN is the total number of cloud particles.  \nThe goal of the E ddy S ed model is to calculate rg for each layer in the atmosphere, given K zz and f sed. It balances the upwards vertical mixing, parameterised by K zz and the particle settling velocity, vf  \nvf = w ∗ r rw ! α . (D.2)  \nHere w ∗ is the convective velocity scale. We note that rw , rg . rw is the radius at which the particle settling velocity equals the  \nconvective velocity scale:  \nw ∗ = K zz L , (D.3)  \nwhere L is the convective mixing length. Since w ∗ is known, and vf can be found analytically as in Ackerman & Marley (2001); Podolak (2003), a linear fit can be used to find both α and rw .  \nWith both of these quantities known, we follow AM01 and define f sed as:  \nf sed = R ∞ 0 r 3 + α n ( r ) dr r α w R ∞ 0 r 3 n ( r ) dr (D.4)  \nFor the log-normal distribution, one finds:  \nZ ∞ 0 r β n ( r ) dr = Nr β g exp 1 2 β 2 log 2 σ g ! (D.5)  \nWhich we can then use to solve for rg :  \nrg = rw f 1 /α sed exp -α + 6 2 log 2 σ g ! (D.6)  \nIn order to use the Hansen distribution, we must recalculate the total number of particles N , and integrate the distribution for f sed. We note here that the Hansen distribution is parameterised by the e ff ective radius, ¯ r , rather than the geometric mean radius. In this derivation we do not correct for this di ff erence in definition, as both act as nuisance parameters in the context of an atmospheric retrieval.  \nWe start by giving the Hansen distribution in full:  \nn Hansen( r ) = N (¯ rve ) (2 ve -1) / ve Γ [(1 -2 ve ) / ve ] r (1 -3 ve ) / ve exp -r ¯ rve ! (D.7)  \nIn Hansen (1971), the authors use the parameters a and b to denote the mean e ff ective radius and e ff ective variance, which we write as ¯ r and ve respectively. These di ff er from the simple mean radius and variance by weighting them by the particle area, as the cloud particle scatters an amount of light proportional to its area. Thus:  \n¯ r = R ∞ 0 r π r 2 n ( r ) dr R ∞ 0 π r 2 n ( r ) dr (D.8)  \nand  \nve = R ∞ 0 ( r -¯ r ) 2 r 2 n ( r ) dr ¯ r 2 R ∞ 0 π r 2 n ( r ) dr (D.9)  \nAs in E ddy S ed , we fit for the settling velocity, which will provide us with α and rw , which we can use to find f sed, as in D.4. However, we must now integrate the Hansen distribution. We find that:  \nZ ∞ 0 r β n Hansen( r ) dr = v β e ( ve β + 2 ve + 1) GLYPH<16> 1 ¯ r GLYPH<17> -β Γ GLYPH<16> β + 1 + 1 ve GLYPH<17> GLYPH<16> -ve + v β + 3 e + 1 GLYPH<17> Γ GLYPH<16> 1 + 1 ve GLYPH<17> (D.10)  \nWe can then use Eqns. D.4 and D.10 to solve for ¯ r :  \n¯ r =         f sed r α w v -α e GLYPH<16> v 3 + α e -ve + 1 GLYPH<17> Γ GLYPH<16> 1 + 1 ve GLYPH<17> ( ve α + 2 ve + 1) Γ GLYPH<16> α + 1 + 1 ve GLYPH<17>         1 α . (D.11)  \nFig. C.1. Comparison of best-fit disequilibrium models to OSIRIS data from Ru ffi o et al. (2021). From top to bottom is HR 8799 b, c, and d. In blue are the best-fit disequilibrium models, with the spectra generated using high-resolution line-by-line opacities, before being convolved, binned, and normalised for comparison. In orange is the same, but with the metallicity set to 0.  \n<!-- image -->  \nFig. C.2. Comparison of best-fit disequilibrium models (black) of HR 8799 c to the data, with residuals shown in the bottom panel. In blue are the same spectra, but without opacity contributions from CH4 .  \n<!-- image -->  \nFig. C.3. Comparison of best-fit disequilibrium models (black) of HR 8799 c to the data, with residuals shown in the bottom panel. In blue are the same spectra, but without opacity contributions from HCN.  \n<!-- image -->  \nThus for a given K zz, f sed, and ve , we can find the e ff ective particle radius for every layer in the atmosphere.  \nHowever, in order to compute the cloud opacity, we still require the total particle count. For a volume mixing ratio of a given species, χ i , we can integrate n ( r ) to find N :  \nN = χ i GLYPH<0> ¯ r 3 ve -1 GLYPH<1> (2 ve -1) (D.12)', 'Appendix E: Complete retrieval results': "We include in the text abridged tables that present key parameters of interest. The complete set of inferred parameters for every retrieval is available on Zenodo.  \nLegend: Chemistry / Profile / Clouds / Data / Info.  \n- -Chemistry: (D)isequilibrium or (F)ree.\n- -Profile: (M)olliere, (Z)hang, (G)uillot or (S)pline(NNodes).\n- -Clouds: Clear (CLR); (f)ree or (eq)uilibrium condensation location, (species)\\_(cd / am)\\_(P)atchy\\_(h)ansen. '*' indicates f sed was retrieved independently for each cloud species.\n- -Data: -(not included) or (only included). O indicates OSIRIS data was used in place of GRAVITY data.\n- -Info: '-' indicates not included. 'mr' indicates mass and radius were used as parameters instead of log g and radius.  \nIn the text, models will be referred to as planet.group.index .  \nb  \n8799  \nHR  \nresults  \nal  \nv  \nretrie  \nAbridged  \nE.1.  \nable  \nT  \n<!-- image -->  \nNotes", 'T able E.2. Abridged retrie v al results HR 8799 c': '<!-- image -->  \nNotes  \npercentiles.  \n1%  \n.  \n34  \n±  \nthe  \nas  \nen  \nv  \ngi  \nuncertainties  \nwith  \nfits,  \nthe  \nfrom  \nalues  \nv  \nmedian  \nthe  \nare  \npresented  \nalues  \nv  \nAll  \nT able E.3. Abridged retrie v al results HR 8799 d  \n<!-- image -->  \nNotes All v alues presented are the median v alues from the fits, with uncertainties gi v en as the ± 34 . 1% percentiles.', 'T able E.4. Abridged retrie v al results HR 8799 e': '<!-- image -->  \nNotes'}