# Four-of-a-kind? Comprehensive atmospheric characterisation of the HR 8799 planets with VLTI/GRAVITY

E. Nasedkin, P. Mollière, S. Lacour, M. Nowak, L. Kreidberg, T. Stolker, J. J. Wang, W. O. Balmer, J. Kammerer, J. Shangguan, et al.

► To cite this version:

E. Nasedkin, P. Mollière, S. Lacour, M. Nowak, L. Kreidberg, et al., Four-of-a-kind? Comprehensive atmospheric characterisation of the HR 8799 planets with VLTI/GRAVITY. Astronomy & Astrophysics - A&A, 2024, 687, 10.1051/0004-6361/20249328. insu-04726457

HAL Id: insu-04726457

https://insu.hal.science/insu-04726457v1

Submitted on 11 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Four-of-a-kind? Comprehensive atmospheric characterisation of the HR 8799 planets with VLTI/GRAVITY

E. Nasedkin¹, P. Mollière¹, S. Lacour²,³, M. Nowak⁴, L. Kreidberg¹, T. Stolker⁵, J. J. Wang⁶, W. O. Balmer⁷,⁸, J. Kammerer³, J. Shangguan⁹, R. Abuter³, A. Amorim¹⁰,¹¹, R. Asensio-Torres¹, M. Benisty¹², J.-P. Berger¹², H. Beust¹², S. Blunt⁶, A. Boccaletti², M. Bonnefoy¹², H. Bonnet³, M. S. Bordini⁹, G. Bourdarot⁹, W. Brandner¹, F. Cantalloube¹³, P. Caselli⁹, B. Charnay², G. Chauvin¹⁴, A. Chavez⁶, E. Choquet¹³, V. Christiaens¹⁵, Y. Clénet², V. Coudé du Foresto², A. Cridland⁵, R. Davies⁹, R. Dembet², J. Dexter¹⁶, A. Drescher⁹, G. Duvert¹², A. Eckart¹⁷,¹⁸, F. Eisenhauer⁹, N. M. Förster Schreiber⁹, P. Garcia¹¹,¹⁹, R. Garcia Lopez²⁰,¹, E. Gendron², R. Genzel⁹,²¹, S. Gillessen⁹, J. H. Girard⁸, S. Grant⁹, X. Haubois²², G. Heißel²³,², Th. Henning¹, S. Hinkley²⁵, S. Hippler¹, M. Houllé¹⁴, Z. Hubert¹², L. Jocou¹², M. Keppler²⁴, P. Kervella², N. T. Kurtovic⁹, A.-M. Lagrange¹²,², V. Lapeyrère², J.-B. Le Bouquin¹², D. Lutz⁹, A.-L. Maire¹², F. Mang⁹, G.-D. Marleau²⁶,²⁷,²⁸,¹, A. Mérand³, J. D. Monnier²⁹, C. Mordasini²⁸, T. Ott⁹, G. P. P. L. Otten³⁰, C. Paladini²², T. Paumard², K. Perraut¹², G. Perrin², O. Pfuhl³, N. Pourré¹², L. Pueyo⁸, D. C. Ribeiro⁹, E. Rickman³¹, J. B. Ruffio³², Z. Rustamkulov³³, T. Shimizu⁹, D. Sing⁷,³³, J. Stadler³⁴,³⁵, O. Straub³⁵, C. Straubmeier¹⁷, E. Sturm⁹, L. J. Tacconi⁹, E. F. van Dishoeck⁵,⁹, A. Vigan¹³, F. Vincent², S. D. von Fellenberg¹⁸, F. Widmann⁹, T. O. Winterhalder³, J. Woillez³, Ş. Yazıcı⁹, and the GRAVITY Collaboration

(Affiliations can be found after the references)

Received 24 January 2024 / Accepted 25 April 2024

## 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 different levels of model flexibility to understand the temperature structure, chemistry, and clouds of each planet using both petitRADTRANS atmospheric retrievals and fits to self-consistent 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<sup>+0.12</sup><sub>-0.10</sub> for d to 0.78<sup>+0.03</sup><sub>-0.04</sub> 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 CH₄ 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 Na₂S. 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.

**Key words.** radiative transfer – instrumentation: interferometers – methods: observational – planets and satellites: atmospheres – planets and satellites: composition

## 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-surface-gravity exoplanets display unique spectral shapes and colours (Faherty 2018), indicative of differences in their chemistry, clouds, and formation history (Marley et al. 2010, 2012; Charnay et al. 2018).

The HR 8799 system is a benchmark target for directly imaged exoplanets, containing four planets (Marois et al. 2008, 2010), an inner debris disc (Boccaletti et al. 2024), and an outer Kuiper-belt-like disc (Su et al. 2009). This is among the best-studied systems of exoplanets, with a wide range of photometric and spectroscopic data. The spectroscopic data cover the near-infrared 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. 2024). Most of these studies, together with extensive modelling work, have tried to answer the following main questions:

1. How did the HR 8799 system form?
2. What are the dynamics of the system? Is it stable, and how do the planets and disc interact?

3. What are the atmospheres of each planet made of, and how have they evolved through time?

In the present work, we attempt to directly answer question (3), which has implications for question (1). Using Bayesian atmospheric 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 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.

With effective 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 CH₄ in exoplanets remain elusive: there have been tentative detections in the atmosphere of HR 8799 b (Barman et al. 2015; Ruffio 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 CH₄ would allow better constraints on this vertical mixing, which is typically parameterised by the vertical diffusion coefficient $K_{zz}$.

The 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 effective 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 blue-ward shift characteristic of T dwarfs as the effective 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 µm as a function of temperature. The strength of this absorption feature was found to correlate positively with the near-infrared 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.

The mechanism through which four super-Jupiter planets can form in a single system is unclear. The presence of both an inner and 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/O ratio 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.

In 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.

While 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 direct-imaging 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/O ratios than when disequilibrium chemistry models are considered (Lavie et al. 2017; Wang et al. 2020b), possibly due to additional oxygen sequestered in refractory clouds (Fonte et al. 2023). The inferred effective temperatures ($T_{eff}$) 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 H-band; 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
<figure><img src="image_2.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

Many 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 CH₄ 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⁻⁵ contrast threshold at 300 mas obtained in Ruffio et al. (2023).

The present study provides a comprehensive examination of the atmospheres of the HR 8799 companions, making use of new, high-S/N observations obtained with VLTI/GRAVITY, together with a combination of retrieval methods and self-consistent modelling. We present further context and background information on the HR 8799 system in Sect. 2. The data used in this work are described in Sect. 3, while the details of the petitRADTRANS (pRT) forward model are described in Sect. 4, with the self-consistent models introduced in Sect. 4.7. The results of the atmospheric retrievals and self-consistent grid fits are presented in Sect. 5. We discuss the limitations of this study, additional work to validate our results, and the implications of our findings in Sect. 6. The appendices contain details of the data and data analysis (Appendices A and B), model validation (Appendix C), implementation details (Appendix D), and tables of the complete set of retrieval results (Appendix E).

## 2. The planetary system of HR 8799

While HR 8799 is one of the most well studied exoplanetary systems (as seen in Fig. 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 classified 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 H₂O 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–4 for planets b, c, d, and e respectively.

### 2.1. HR 8799 A: The host star

HR 8799 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. 2010b) due to its depletion of heavy elements in the atmosphere. Due to this depletion, the [Fe/H] of HR 8799 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. 2010b). 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. (2010b). Using high resolution spectroscopy from the LBT/PEPSI and HARPS instruments, Wang et al. (2020b) 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. For this work, we use a BT-Nextgen stellar model as fitted in Nasedkin et al. (2023) with best-fit parameters of T<sub>eff</sub> = 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.

Most 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. (2010a) 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
<table><caption>**Table 1.** Summary of literature and derived planet properties for HR 8799 b.</caption><thead><tr><th>Planet</th><th>Ref.</th><th>Clouds</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th><th>log <i>g</i><br>(cgs)</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>R<br>(<i>R</i><sub>Jup</sub>)</th><th>[M/H]</th><th>C/O</th><th>log <i>L</i><sub>bol</sub>/<i>L</i><sub>☉</sub></th></tr></thead><tbody><tr><td rowspan="15">b</td><td>B11a</td><td>Slab</td><td>0.1–3.3</td><td>3.5 ± 0.5</td><td>1100 ± 100</td><td>0.63–0.92</td><td>...</td><td>...</td><td>-5.1 ± 0.1</td></tr><tr><td>C11</td><td>Thick</td><td>5–15</td><td>4–4.5</td><td>800–1000</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>G11</td><td>Slab</td><td>1.8</td><td>4</td><td>1100</td><td>0.69</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M11</td><td>Power law</td><td>2–12</td><td>3.5–4.3</td><td>750–850</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M12</td><td>AM01</td><td>26</td><td>4.75</td><td>1000</td><td>1.11</td><td>0</td><td>...</td><td>-4.95 ± 0.06</td></tr><tr><td>L13</td><td>Slab</td><td>16<sup>+5</sup><sub>-4</sub></td><td>5.0<sup>+0.1</sup><sub>-0.2</sub></td><td>900<sup>+30</sup><sub>-60</sub></td><td>0.66<sup>+0.07</sup><sub>-0.04</sub></td><td>...</td><td>0.96</td><td>-5.1 ± 0.1</td></tr><tr><td>B15</td><td>Slab</td><td>...</td><td>3.5</td><td>1000</td><td>...</td><td>...</td><td>0.55–0.7</td><td>...</td></tr><tr><td>B16</td><td>ER4</td><td>...</td><td>3.4–3.8</td><td>1100–1200</td><td>0.6–0.7</td><td>0.5</td><td>...</td><td>...</td></tr><tr><td>L17</td><td>Mie</td><td>...</td><td>4 ± 0.1</td><td>320 ± 20<sup>(a)</sup></td><td>1.08 ± 0.02</td><td>...</td><td>0.92 ± 0.01</td><td>...</td></tr><tr><td>W21<sup>(b)</sup></td><td>BT-Settl</td><td>...</td><td>4.8<sup>+0.4</sup><sub>-0.8</sub></td><td>1423.3<sup>+212.6</sup><sub>-278.4</sub></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>R21<sup>(c)</sup></td><td>Slab</td><td>...</td><td>3.1<sup>+0.03</sup><sub>-0.03</sub></td><td>1180<sup>+14</sup><sub>-14</sub></td><td>...</td><td>0</td><td>0.578<sup>+0.004</sup><sub>-0.005</sub></td><td>...</td></tr><tr><td>Best A ∩ B<sup>(d)</sup></td><td>AM01</td><td>6.0<sup>+0.3</sup><sub>-0.3</sub></td><td>4.10<sup>+0.03</sup><sub>-0.04</sub></td><td>942<sup>+12</sup><sub>-13</sub></td><td>1.11<sup>+0.03</sup><sub>-0.03</sub></td><td>0.96<sup>+0.08</sup><sub>-0.08</sub></td><td>0.78<sup>+0.03</sup><sub>-0.04</sub></td><td>-5.08<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>BMA A ∩ B<sup>(e)</sup></td><td>AM01</td><td>6.0<sup>+0.4</sup><sub>-0.3</sub></td><td>4.10<sup>+0.03</sup><sub>-0.04</sub></td><td>942<sup>+12</sup><sub>-15</sub></td><td>1.10<sup>+0.03</sup><sub>-0.03</sub></td><td>0.96<sup>+0.08</sup><sub>-0.08</sub></td><td>0.78<sup>+0.03</sup><sub>-0.04</sub></td><td>-5.08<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>A ∩ B<sup>(f)</sup></td><td>AM01</td><td>6.0<sup>+0.3</sup><sub>-0.3</sub></td><td>4.10<sup>+0.06</sup><sub>-0.06</sub></td><td>936<sup>+22</sup><sub>-34</sub></td><td>1.11<sup>+0.08</sup><sub>-0.08</sub></td><td>1.1<sup>+0.1</sup><sub>-0.2</sub></td><td>0.73<sup>+0.04</sup><sub>-0.04</sub></td><td>-5.08<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td>Grids</td><td>Various</td><td>4.6–6.0</td><td>3.5–4.5</td><td>850–1100</td><td>0.73–1.2</td><td>> 0.7</td><td>0.3–0.55</td><td>-5.14 to -5.28</td></tr></tbody></table>

**Notes.** 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 diffusion and mixing of aerosol particles (Allard et al. 2012). (a) Only *T*<sub>int</sub>, a model parameter, is reported. (b) W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7.2 ± 0.7 *M*<sub>Jup</sub> for the inner three planets and 5.8 ± 0.5 *M*<sub>Jup</sub> for HR 8799 b. A radius of 1.2 ± 0.1 *R*<sub>Jup</sub> was used for all planets. (c) R21 uncertainties were limited by the coarseness and boundaries of their model grid. (d) Single best retrieval parameters. (e) Bayesian model averaged parameters from group A ∩ B. (f) Unweighted average parameters from group A ∩ B.

**References.** 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. (2021b); R21: Ruffio et al. (2021).

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 et al. 2011, 2012, 2014; Bergfors et al. 2011; Galicher et al. 2011; Soummier et al. 2011; Skemer et al. 2012; Esposito et al. 2013; Skemer 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. 2024) 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).

In 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.

Bonnefoy et al. (2016) explores the implications of the near-infrared photometry for all four of the companions,

comparing them to spectrally similar field objects from the SpeX PRISM library. 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 coefficients for corundum (Al₂O₃), iron (Fe), enstatite (MgSiO₃), and forsterite (Mg₂SiO₄) 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).

Marley 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 Baraffe et al. (2003) and Saumon & Marley (2008). With self-consistent, radiative equilibrium models, the planet radius is often difficult 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 this amplitude is rare outside of the transition regime (Radigan 2014). However, as we view the HR 8799 system nearly
<table><caption>**Table 2.** Summary of literature and derived planet properties for HR 8799 c.</caption><thead><tr><th>Planet</th><th>Ref.</th><th>Clouds</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th><th>log <i>g</i><br>(cgs)</th><th><i>T</i><sub>eff</sub><br>(K)</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th><th>[M/H]</th><th>C/O</th><th>log <i>L</i><sub>bol</sub>/<i>L</i><sub>☉</sub></th></tr></thead><tbody><tr><td rowspan="15">c</td><td>C11</td><td>Thick</td><td>7–17.5</td><td>4–4.5</td><td>1000–1200</td><td>...</td><td>0</td><td>...</td><td>−4.7 ± 0.1</td></tr><tr><td>G11</td><td>Slab</td><td>1.1</td><td>3.5</td><td>1200</td><td>0.97</td><td>...</td><td>...</td><td>−4.7 ± 0.1</td></tr><tr><td>M11</td><td>Power law</td><td>3–11</td><td>4.0–4.3</td><td>950–1025</td><td>...</td><td>...</td><td>...</td><td>−4.7 ± 0.1</td></tr><tr><td>M12</td><td>AM01</td><td>8–11</td><td>4.1 ± 0.1</td><td>950 ± 60</td><td>1.32–1.39</td><td>0</td><td>...</td><td>−4.90 ± 0.1</td></tr><tr><td>K13</td><td>...</td><td>3–7</td><td>3.5–4.0</td><td>1100 ± 100</td><td>1–1.5</td><td>...</td><td>0.65 ± 0.1</td><td>...</td></tr><tr><td>B16</td><td>ER4</td><td>...</td><td>3.8–3.9</td><td>1200</td><td>1.0</td><td>0.5</td><td>...</td><td>...</td></tr><tr><td>L17</td><td>Mie</td><td>...</td><td>4.5 ± 0.1</td><td>960 ± 20<sup>(a)</sup></td><td>1.25 ± 0.02</td><td>...</td><td>0.55 ± 0.01</td><td>...</td></tr><tr><td>G18</td><td>Various</td><td>...</td><td>3.5–4.0</td><td>1100–1350</td><td>0.7–1.2</td><td>...</td><td>...</td><td>−4.58 to −4.82</td></tr><tr><td>W20</td><td>Deck</td><td>...</td><td>3.97<sup>+0.03</sup><sub>−0.03</sub></td><td>1054<sup>+7</sup><sub>−5</sub></td><td>1.47 ± 0.02</td><td>...</td><td>0.58<sup>+0.06</sup><sub>−0.06</sub></td><td>−4.59 ± 0.004</td></tr><tr><td>W20<sup>(b)</sup></td><td>Deck</td><td>...</td><td>3.95<sup>+0.04</sup><sub>−0.12</sub></td><td>1102 ± 2</td><td>1.20</td><td>...</td><td>0.39<sup>+0.06</sup><sub>−0.06</sub></td><td>−4.69 ± 0.0002</td></tr><tr><td>W21<sup>(c)</sup></td><td>BT-Settl</td><td>...</td><td>5.4<sup>+0.1</sup><sub>−0.2</sub></td><td>1474.4<sup>+24.4</sup><sub>−36.3</sub></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>R21<sup>(d)</sup></td><td>Slab</td><td>...</td><td>3.63<sup>+0.03</sup><sub>−0.02</sub></td><td>1200<sup>+*</sup><sub>−14</sub></td><td>...</td><td>0</td><td>0.562<sup>+0.004</sup><sub>−0.005</sub></td><td>...</td></tr><tr><td>W23</td><td>AM01</td><td>...</td><td>4.17<sup>+0.41</sup><sub>−0.47</sub></td><td>1421<sup>+92<sup>(e)</sup></sup><sub>−72</sub></td><td>1.01<sup>+0.09</sup><sub>−0.08</sub></td><td>0.51<sup>+0.40<sup>(f)</sup></sup><sub>−0.43</sub></td><td>0.67<sup>+0.12</sup><sub>−0.15</sub></td><td>...</td></tr><tr><td>Best A ∩ B</td><td>AM01</td><td>8.5<sup>+0.4</sup><sub>−0.4</sub></td><td>4.26<sup>+0.02</sup><sub>−0.03</sub></td><td>1158<sup>+11</sup><sub>−12</sub></td><td>1.10<sup>+0.01</sup><sub>−0.01</sub></td><td>1.27<sup>+0.06</sup><sub>−0.06</sub></td><td>0.66<sup>+0.01</sup><sub>−0.01</sub></td><td>−4.71<sup>+0.02</sup><sub>−0.02</sub></td></tr><tr><td>BMA A ∩ B</td><td>AM01</td><td>8.5<sup>+0.4</sup><sub>−0.5</sub></td><td>4.26<sup>+0.02</sup><sub>−0.03</sub></td><td>1159<sup>+11</sup><sub>−12</sub></td><td>1.10<sup>+0.01</sup><sub>−0.01</sub></td><td>1.27<sup>+0.05</sup><sub>−0.06</sub></td><td>0.66<sup>+0.01</sup><sub>−0.01</sub></td><td>−4.71<sup>+0.02</sup><sub>−0.02</sub></td></tr><tr><td>A ∩ B</td><td>AM01</td><td>8.6<sup>+0.4</sup><sub>−0.4</sub></td><td>4.25<sup>+0.04</sup><sub>−0.14</sub></td><td>1159<sup>+24</sup><sub>−76</sub></td><td>1.10<sup>+0.23</sup><sub>−0.03</sub></td><td>1.2<sup>+0.1</sup><sub>−0.1</sub></td><td>0.63<sup>+0.05</sup><sub>−0.02</sub></td><td>−4.70<sup>+0.03</sup><sub>−0.03</sub></td></tr><tr><td>Grids</td><td>Various</td><td>1.24–10.3</td><td>3.5–4.5</td><td>1100–1200</td><td>0.8–1.31</td><td>> 1.0</td><td>0.3–0.8</td><td>−4.65 to −4.72</td></tr></tbody></table>

**Notes.** (a) Only $T_{\text{int}}$, a model parameter, is reported. (b) W20 compared using strong and weak radius priors to enforce physicality. (c) W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7.2 ± 0.7 $M_{\text{Jup}}$ for the inner three planets and 5.8 ± 0.5 $M_{\text{Jup}}$ for HR 8799 b. A radius of 1.2 ± 0.1 $R_{\text{Jup}}$ was used for all planets. (d) R21 uncertainties were limited by the coarseness and boundaries of their model grid. (e) W23 report the temperature at 3.3 bar rather than the effective temperature. (f) The metallicity of W23 was found by averaging the retrieved C/H and O/H ratios.

**References.** 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. (2020b); W21: Wang et al. (2021b); R21: Ruffio et al. (2021); W23: Wang et al. (2023).

pole on, it is difficult to see the effects 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 post-processing 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_{\odot}$. This model produced results typical of the range of models explored; the mass estimates for each companion are: $e = 7.6 \pm 0.9 M_{\text{Jup}}$, $d = 9.2 \pm 0.1 M_{\text{Jup}}$, $c = 7.7 \pm 0.7 M_{\text{Jup}}$, and $b = 5.8 \pm 0.4 M_{\text{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. 2021b). 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 different 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.

Many 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 additional H-band OSIRIS observations, and inferred the low temperature,
<table><caption>**Table 3.** Summary of literature and derived planet properties for HR 8799 d.</caption><thead><tr><th>Planet</th><th>Ref.</th><th>Clouds</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th><th>log <i>g</i><br>(cgs)</th><th><i>T</i><sub>eff</sub><br>(K)</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th><th>[M/H]</th><th>C/O</th><th>log <i>L</i><sub>bol</sub>/<i>L</i><sub>☉</sub></th></tr></thead><tbody><tr><td rowspan="12">d</td><td>C11</td><td>Thick</td><td>5–17.5</td><td>3.75–4.5</td><td>1000–1200</td><td>...</td><td>0</td><td>...</td><td>−4.7 ± 0.1</td></tr><tr><td>G11</td><td>Slab</td><td>6</td><td>4.0</td><td>1100</td><td>1.25</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M11</td><td>Power law</td><td>3–11</td><td>3.5–4.2</td><td>850–1000</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M12</td><td>AM01</td><td>8–11</td><td>4.1 ± 0.1</td><td>1000 ± 75</td><td>1.33–1.41</td><td>0</td><td>...</td><td>−4.80 ± 0.09</td></tr><tr><td>B16</td><td>ER4</td><td>...</td><td>4.4–4.5</td><td>1200–1300</td><td>0.9–1.1</td><td>0.5</td><td>...</td><td>...</td></tr><tr><td>L17</td><td>Mie</td><td>...</td><td>4.2 ± 0.2</td><td>1420 ± 10<sup>(a)</sup></td><td>0.96 ± 0.05</td><td>...</td><td>0</td><td>...</td></tr><tr><td>G18</td><td>Various</td><td>...</td><td>3.5–4.0</td><td>1100–1600</td><td>0.65–1.4</td><td>...</td><td>...</td><td>−4.58 to −4.82</td></tr><tr><td>W21<sup>(b)</sup></td><td>BT-Settl</td><td>...</td><td>5.1<sup>+0.3</sup><sub>−0.4</sub></td><td>1558.8<sup>+50.9</sup><sub>−91.4</sub></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>R21<sup>(c)</sup></td><td>Slab</td><td>3.7<sup>+0.03</sup><sub>−0.03</sub></td><td>1200<sup>+*</sup><sub>−14</sub></td><td>...</td><td>...</td><td>0</td><td>0.551<sup>+0.004</sup><sub>−0.005</sub></td><td>...</td></tr><tr><td>Best <i>A</i> ∩ <i>B</i></td><td>AM01</td><td>9.19<sup>+0.08</sup><sub>−0.07</sub></td><td>4.18<sup>+0.04</sup><sub>−0.03</sub></td><td>1177<sup>+21</sup><sub>−21</sub></td><td>1.26<sup>+0.05</sup><sub>−0.06</sub></td><td>1.2<sup>+0.2</sup><sub>−0.1</sub></td><td>0.61<sup>+0.03</sup><sub>−0.04</sub></td><td>−4.63<sup>+0.04</sup><sub>−0.04</sub></td></tr><tr><td>BMA <i>A</i> ∩ <i>B</i></td><td>AM01</td><td>9.19<sup>+0.08</sup><sub>−0.07</sub></td><td>4.18<sup>+0.06</sup><sub>−0.04</sub></td><td>1179<sup>+31</sup><sub>−28</sub></td><td>1.26<sup>+0.06</sup><sub>−0.08</sub></td><td>1.2<sup>+0.2</sup><sub>−0.2</sub></td><td>0.60<sup>+0.04</sup><sub>−0.06</sub></td><td>−4.62<sup>+0.05</sup><sub>−0.04</sub></td></tr><tr><td><i>A</i> ∩ <i>B</i></td><td>AM01</td><td>9.20<sup>+0.09</sup><sub>−0.08</sub></td><td>4.19<sup>+0.07</sup><sub>−0.04</sub></td><td>1179<sup>+38</sup><sub>−36</sub></td><td>1.25<sup>+0.06</sup><sub>−0.09</sub></td><td>1.2<sup>+0.4</sup><sub>−0.3</sub></td><td>0.55<sup>+0.12</sup><sub>−0.10</sub></td><td>−4.61<sup>+0.05</sup><sub>−0.05</sub></td></tr><tr><td>Grids</td><td>Various</td><td>8.2–9.9</td><td>3.5–4.5</td><td>1200–1300</td><td>0.97–1.21</td><td>&gt; 0.0</td><td>0.2–0.55</td><td>−4.59 to −4.65</td></tr></tbody></table>

**Notes.** (a) Only $T_{\text{int}}$, a model parameter, is reported. (b) W21 used high resolution spectroscopy, and did not infer masses or radii, using masses of 7.2 ± 0.7 $M_{\text{Jup}}$ for the inner three planets and 5.8 ± 0.5 $M_{\text{Jup}}$ for HR 8799 b. A radius of 1.2 ± 0.1 $R_{\text{Jup}}$ was used for all planets. (c) R21 uncertainties were limited by the coarseness and boundaries of their model grid.

**References.** 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. (2021b); R21: Ruffio et al. (2021).

low surface gravity, and low CH$_4$ 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 H$_2$O, CO, and tentatively CH$_4$. 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.

Janson 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 1300 K, with a log g around 4.0. As with

the photometry, the self-consistent models they used to infer the planet properties struggled to obtain radius estimates consistent with predictions of evolutionary models. Wang et al. (2021b) use high resolution spectroscopy to measure the rotation of c, d, and e, finding an upper limit of 14 km s$^{-1}$ for c, and measurements of 10.1$^{+2.8}_{-2.7}$ km s$^{-1}$ for d and 15.0$^{+2.3}_{-2.6}$ km s$^{-1}$ 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. (2020b), 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. Ruffio et al. (2019, 2021) and Wang et al. (2021b) explore HR 8799 c using moderate and high resolution spectroscopy respectively. Both works characterise the dynamics of the planets, with Ruffio et al. (2021) measuring the radial velocities for planets b, c, and d, finding them to be −9.1 ± 0.4 km s$^{-1}$, −11.1 ± 0.4 km s$^{-1}$, and −11.6 ± 0.8 km s$^{-1}$ 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 CH$_4$, which was consistent with the results of Wang et al. (2018a).

The 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 L6-L8 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 effective temperatures
<table><caption>**Table 4.** Summary of literature and derived planet properties for HR 8799 e.</caption><thead><tr><th>Planet</th><th>Ref.</th><th>Clouds</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th><th>log <i>g</i><br>(cgs)</th><th><i>T</i><sub>eff</sub><br>(K)</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th><th>[M/H]</th><th>C/O</th><th>log <i>L</i><sub>bol</sub>/<i>L</i><sub>☉</sub></th></tr></thead><tbody><tr><td rowspan="10">e</td><td>B16</td><td>ER4</td><td>...</td><td>3.7–4.1</td><td>1200–1300</td><td>0.9–1.0</td><td>0.5</td><td>...</td><td>...</td></tr><tr><td>L17</td><td>Mie</td><td>...</td><td>3.8 ± 0.3</td><td>1230 ± 30<sup>(a)</sup></td><td>1.2 ± 0.1</td><td>...</td><td>0</td><td>...</td></tr><tr><td>G18</td><td>Various</td><td>...</td><td>3.5–4.0</td><td>1100–1650</td><td>0.6–1.4</td><td>...</td><td>...</td><td>-4.58 to -4.75</td></tr><tr><td>M20</td><td>AM01</td><td>4.81<sup>+8.78</sup><sub>-3.33</sub></td><td>4.00<sup>+0.46</sup><sub>-0.52</sub></td><td>1154<sup>+49</sup><sub>-48</sub></td><td>1.12<sup>+0.09</sup><sub>-0.09</sub></td><td>0.48<sup>+0.25</sup><sub>-0.29</sub></td><td>0.60<sup>+0.07</sup><sub>-0.08</sub></td><td>...</td></tr><tr><td>W21<sup>(b)</sup></td><td>BT-Settl</td><td>...</td><td>3.7<sup>+0.3</sup><sub>-0.1</sub></td><td>1345.6<sup>+57.0</sup><sub>-53.3</sub></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Best A ∩ B</td><td>AM01</td><td>7.5<sup>+0.6</sup><sub>-0.6</sub></td><td>4.20<sup>+0.06</sup><sub>-0.06</sub></td><td>1172<sup>+29</sup><sub>-27</sub></td><td>1.14<sup>+0.05</sup><sub>-0.05</sub></td><td>1.9<sup>+0.1</sup><sub>-0.1</sub></td><td>0.88<sup>+0.02</sup><sub>-0.02</sub></td><td>-4.71<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>BMA A ∩ B</td><td>AM01</td><td>7.5<sup>+0.7</sup><sub>-0.7</sub></td><td>4.20<sup>+0.06</sup><sub>-0.06</sub></td><td>1161<sup>+33</sup><sub>-34</sub></td><td>1.12<sup>+0.05</sup><sub>-0.05</sub></td><td>1.9<sup>+0.1</sup><sub>-0.2</sub></td><td>0.88<sup>+0.02</sup><sub>-0.02</sub></td><td>-4.72<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td>A ∩ B</td><td>AM01</td><td>7.5<sup>+0.7</sup><sub>-0.7</sub></td><td>4.3<sup>+0.1</sup><sub>-0.1</sub></td><td>1198<sup>+41</sup><sub>-77</sub></td><td>1.05<sup>+0.15</sup><sub>-0.08</sub></td><td>1.8<sup>+0.3</sup><sub>-0.4</sub></td><td>0.84<sup>+0.06</sup><sub>-0.07</sub></td><td>-4.71<sup>+0.07</sup><sub>-0.08</sub></td></tr><tr><td>Grids</td><td>Various</td><td>1.07–8.8</td><td>3.5–4.5</td><td>1100–1400</td><td>0.75–1.24</td><td>> 1.0</td><td>> 0.55</td><td>-4.70 to -4.78</td></tr></tbody></table>

**Notes.** (a) Only $T_{int}$, a model parameter, is reported. (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.

**References.** B16: Bonnefoy et al. (2016); L17: Lavie et al. (2017); G18: Greenbaum et al. (2018); M20: Mollière et al. (2020); W21: Wang et al. (2021b).

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 provided 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 (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 different shape than the other two planets, but that all three shared supersolar metallicities and effective 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 difficult to distinguish between different 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 (2019), together 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 enrichment of [M/H]=0.48<sup>+0.25</sup><sub>-0.29</sub> and a C/O ratio of 0.60<sup>+0.07</sup><sub>-0.08</sub>. Finally, Wang et al. (2020b, 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 efforts (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 effective 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. 2018), 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).

Zahnle & Marley (2014) provide an in-depth exploration of the impacts of disequilibrium chemistry on cool, self-luminous atmospheres, providing predictions for the CO, CH$_4$, and NH$_3$ abundances as a function of vertical mixing and effective temperature, identifying the key transition between CH$_4$ 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 effective temperatures, and provides column abundance predictions for HR 8799 b, finding that the CO abundance should dominate over CH$_4$, assuming a solar composition. Soni & Acharyya (2023) extend this to superstellar metallicities and vertical mixing strengths, using the constraints on the CO and CH$_4$ abundances from Barman et al. (2015) to infer a vertical mixing strength of $\log K_{zz} \in [7, 10]$ for the $10\times$ 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 H$_2$ 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_\star$ ratio of between 3 and 5 (in a 68% confidence interval).

In addition to the 1D modelling efforts, 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 low-gravity 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 formed through gravitational instability, it is unlikely that disc fragmentation could have formed the inner three companions.

In 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 effect 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.

The different formation scenarios can lead to dramatically different 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 different radiative efficiencies can lead to different 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 (Baraffe 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 Baraffe 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.

Finally, 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. 2024), 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
<figure><img src="image_3.png" /><figcaption>**Fig. 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 \times 10^{-15} \text{ W m}^{-2} (\mu\text{m})^{-1}$ offset. 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.</figcaption></figure>

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. 2024). Modelling efforts have placed tentative limits of ~1 $M_{\oplus}$ on the allowed mass of companions interior to HR 8799 e (Goździewski & 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<fnref n="1" />.

### 3.1. GRAVITY data

In Fig. 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 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 offers a resolution of *R*~500 over a nominal wavelength range of 2.0–2.4 µm.

Two different 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 using 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/off-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 (2020b).

The data-reduction was performed using the tools developed for the ExoGRAVITY large program<fnref n="2" />. The main steps of the reductions are as follows:

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.

2. For the on-axis observations, the phase reference is extracted from the on-star observations and subtracted from the on-planet observations. For the off-axis observations, this phase-reference is extracted from the observations of the binary-calibrator HD 25535. In both cases, the amplitude reference is taken using the on-star observations.

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.

<footnote n="1">¹ https://zenodo.org/records/10914429</footnote>

<footnote n="2">² https://gitlab.obspm.fr/mnowak/exogravity</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 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. (2024). References: OSIRIS (Barman et al. 2011); SPHERE (Zurlo et al. 2016); GPI (Greenbaum et al. 2018); CHARIS (Wang et al. 2020b, 2022); ALES (Doelman et al. 2022).</figcaption></figure>

The astrometry of the planet is then extracted from the observations.

4. The contrast spectrum is then extracted using a model that also takes into account the residual starlight, and the planet astrometry previously extracted.

This procedure yields a contrast spectrum for each planet, at each observation date. The spectrum extraction, which consists entirely 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.

For 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$ are given by:

$$W = \left[ \sum_{t} W_{t}^{-1} \right]^{-1} \quad (1)$$

$$C = W \cdot \left[ \sum_{t} W_{t}^{-1} \cdot C_{t} \right], \quad (2)$$

where $C_t$ and $W_t$ represent the contrast spectrum and its associated covariance matrix on a given observing date $t$.

The 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 Sect. 2.1 and are based on Nasedkin et al. (2023).

The 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 yr 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 reread 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 KLIPin 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: $N_b = 297$, $N_c = 391$, $N_d = 387$, $N_e = 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 KLIPin 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.
Additional 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 insufficient S/N. While we attempted to calibrate the companion spectra using the 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 *K*1, *K*2 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 (Groff et al. 2015, 2017) observations of HR 8799 c, d, and e were presented in Wang et al. (2020b, 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 $\sqrt{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 GRAVITY in measuring the near infrared spectrum of the planet. With an unbinned spectral resolution of $R \approx 4000$, and an integration time of 2700 s in the *H*-band and 1800 s in the *K*-band, when binned to a spectral resolution of $R \approx 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).

<table><caption>**Table 5.** Retrieval prior; temperature profile priors are included in Table 6. $N(\mu, \sigma)$, $U(\text{low, high})$.</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td>log g</td><td>U(2.5, 5.5)</td></tr><tr><td>Radius (R<sub>Jup</sub>)</td><td>U(0.7, 2.0)<br>N(1.1, 0.1)</td></tr><tr><td>Mass (M<sub>Jup</sub>)</td><td>N(μ<sub>M, dyn</sub>, σ<sub>M, dyn</sub>)</td></tr><tr><td>log P<sub>Quench</sub> (log bar)</td><td>U(-6.0, 3.0)</td></tr><tr><td>[M/H]</td><td>U(0.5, 2.5)</td></tr><tr><td>C/O</td><td>U(0.1, 1.6)</td></tr><tr><td>σ<sub>LN</sub></td><td>U(1.05, 3.0)</td></tr><tr><td>f<sub>sed</sub></td><td>U(0.0, 10.0)</td></tr><tr><td>log K<sub>zz</sub></td><td>U(5.0, 13.0)</td></tr><tr><td>log Eq. Cloud Scaling</td><td>U(-2.5, 2.5)</td></tr><tr><td>log Cloud Mass Fracs.</td><td>U(-6.5, 0.0)</td></tr><tr><td>log Cloud P<sub>Base</sub></td><td>U(-6.0, 3.0)</td></tr><tr><td>log Mass Fracs.</td><td>U(-7.0, 0.3)</td></tr></tbody></table>

We also include the *K*-band spectra of Konopacky et al. (2013; Fig. 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 different measurements on the retrieved atmospheric parameters, and to determine if our methods can reproduce the results of earlier work.

More recent observations of the HR 8799 planets using OSIRIS have been explored in Ruffio et al. (2021), but these spectra are continuum subtracted, requiring a somewhat different 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<fnref n="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.

To allow for both a data-driven and physically motivated approach, we retrieved either log g and R<sub>pl</sub> with uniform priors or R<sub>pl</sub> and M<sub>pl</sub>, with Gaussian priors set by the dynamical mass estimates of Zurlo et al. (2022) and broad Gaussian priors centred at 1.1 R<sub>Jup</sub>, in line with estimates from evolutionary models (Marley et al. 2012).

As 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

<footnote n="3">³ https://petitradtrans.readthedocs.io/</footnote>
<table><caption>**Table 6.** Priors for temperature profiles.</caption><thead><tr><th colspan="2">Spline</th><th colspan="2">Guillot (G10)</th><th colspan="2">Mollière (M20)</th><th colspan="2">Zhang (Z23)</th></tr><tr><th>Name</th><th>Prior</th><th>Name</th><th>Prior</th><th>Name</th><th>Prior</th><th>Name</th><th>Prior</th></tr></thead><tbody><tr><td>T<sub>0</sub></td><td>U(0 K, 300 K)</td><td>T<sub>int</sub></td><td>U(300 K, 2500 K)</td><td>T<sub>int</sub></td><td>U(300 K, 2000 K)</td><td>T<sub>bot</sub></td><td>U(2000 K, 12 000 K)</td></tr><tr><td>T<sub>i</sub></td><td>U(300 K, 11 900 K)</td><td>T<sub>equ</sub></td><td>U(10 K, 100 K)</td><td>T<sub>0</sub></td><td>U(0, T<sub>Edd</sub>)</td><td>d log T/d log P<sub>0</sub></td><td>N(0.25, 0.025)</td></tr><tr><td>γ</td><td>Γ<sup>-1</sup>(1, 5 × 10<sup>-5</sup>)</td><td>γ</td><td>N(0, 2)</td><td>T<sub>1</sub></td><td>U(0, T<sub>0</sub>)</td><td>d log T/d log P<sub>1</sub></td><td>N(0.25, 0.045)</td></tr><tr><td>...</td><td>...</td><td>log κ<sub>IR</sub></td><td>U(-4, 1)</td><td>T<sub>2</sub></td><td>U(0, T<sub>1</sub>)</td><td>d log T/d log P<sub>2</sub></td><td>N(0.26, 0.05)</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>α</td><td>U(1, 2)</td><td>d log T/d log P<sub>3</sub></td><td>N(0.2, 0.05)</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>log δ</td><td>U(0, 1)</td><td>d log T/d log P<sub>4</sub></td><td>N(0.12, 0.045)</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>d log T/d log P<sub>5</sub></td><td>N(0.07, 0.07)</td></tr></tbody></table>

**Notes.** N(μ, σ): Gaussian prior. U(low, high): uniform prior. Γ<sup>-1</sup>(α, β): inverse Gamma function prior.

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 different cloud parameterisations.

## 4.1. Thermal structure

We 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 different temperature profiles we computed an effective 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 effective 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.Pchip Interpolator` 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,

$$ \log p(\mathbf{T}) = \frac{1}{2\gamma} \sum_{i=1}^{N-1} (T_{i+1} - 2T + T_{i-1})^2 - \log(2\pi\gamma). \quad (3) $$

This 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):

$$ \Gamma^{-1}(\gamma) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} \left(\frac{1}{\gamma}\right)^{\alpha+1} \exp\left(-\frac{\beta}{\gamma}\right), \quad (4) $$

for 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 planets:

$$ T_{\text{Guillot}}^4 = \frac{3T_{\text{int}}^4}{4} \left(\frac{2}{3} + \tau\right) + \frac{3T_{\text{irr}}^4}{4} \left(\frac{2}{3} + \frac{1}{\gamma\sqrt{3}} + \left(\frac{\gamma}{\sqrt{3}} - \frac{1}{\gamma\sqrt{3}}\right) e^{-\gamma\tau\sqrt{3}}\right), \quad (5) $$

where $T_{\text{irr}} = \sqrt{2}T_{\text{equ}}$ and $\tau = P \times \kappa_{\text{IR}}/g$. $T_{\text{equ}}$ is the equilibrium temperature of an irradiated body, $T_{\text{int}}$ is the intrinsic temperature of the planet, and g is the surface gravity. $\kappa_{\text{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_{\text{irr}}$ is small, and thus the profile reduces to an Eddington (1930) profile, which corresponds to keeping only the first term on the right-hand side of Eq. (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.

The 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 Eq. (5). However,
for this profile we parameterise the opacity $\tau$ as a function of pressure (P):

$$ \tau = \delta P^{\alpha} \qquad (6) $$

and retrieve parameters of log $\delta$ and $\alpha$, together with $T_{int}$, as in the G10 profile.

The upper atmosphere is defined as the region above $\tau = 0.1$. Above this level, four pressure points are defined, equidistant in log P. The deepest pressure point, at $\tau = 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.

The 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:

$$ \frac{dT}{dr} < \frac{T}{P} \frac{dP}{dr} \left( 1 - \frac{1}{\gamma_{\text{ad}}} \right). \qquad (7) $$

The 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 Sect. 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 accomplished by fitting for the gradient of the temperature with respect to pressure, as opposed to directly retrieving the temperature as in the spline profile. The prior locations and widths of the gradients were determined by empirically measuring the temperature gradients in self-consistent radiative-convective models, thus providing a means to enforce the physics of these models in a retrieval framework. The atmosphere between $10^3$ bar and $10^{-3}$ bar was divided up into six layers, equidistant in log pressure. The temperature at the bottom of the atmosphere ($T_{bot}$) was freely retrieved. For the remaining layers, $d \log T / d \log P|_{i}$ were retrieved as free parameters. The temperature profile was then found by interpolating the gradient to the full pressure grid, and integrating to find the temperature at each pressure.

$$ T_0 = T_{\text{Bot}} \qquad (8) $$

$$ T_{i+1} = \exp \left( \log T_i + (\log P_{i+1} - \log P_i) \left( \frac{d \log T}{d \log P} \right)_i \right) \qquad (9) $$

The atmosphere was isothermal above $10^{-3}$ bar.

## 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. 2020), and included H$_2$O (Polyansky et al. 2018), CO (Rothman et al. 2010), CH$_4$ (Yurchenko et al. 2017), CO$_2$ (Yurchenko et al. 2020), NH$_3$ (Coles et al. 2019), HCN (Barber et al. 2014), H$_2$S (Azzam et al. 2016), PH$_3$ (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 X_i$ for each of H$_2$O, CO, CH$_4$, CO$_2$, HCN, H$_2$S, NH$_3$, 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 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 H$_2$, 0.234 for He), and multiplying by one minus the sum of the retrieved mass fractions ($1 - \Sigma_i X_i$). 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 H$_2$O, CO, and CH$_4$, 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 different setups, this setup is representative of typical retrievals, and ensures consistent comparisons.

To 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 H$_2$O, CO, CH$_4$, CO$_2$, NH$_3$, and H$_2$S 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, CH$_4$, H$_2$O 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 Sect. 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.3. Clouds

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 diffusion ($K_{zz}$), and settle back down with a speed proportional to the parameter $f_{sed}$. We retrieved each of these parameters independently, together with $\sigma_{LN}$, which is the width of the log-normal particle size distribution.

While these parameters determine the structure of the clouds, we also determined the cloud opacity through the use of cloud optical constants for different cloud compositions, allowing us to differentiate 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 $\sigma_{ln}$, we retrieved the mean effective 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.

We tested a range of cloud compositions, including MgSiO$_3$ (both crystalline and amorphous particle shapes), Mg$_2$SiO$_4$, Fe, Al$_2$O$_3$, KCl, and Na$_2$S, 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 different vertical extents for different cloud compositions.

Nominally, 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 $X_i$ was defined at this base pressure $P_0$, and decreases with altitude to the power of $f_{sed}$:

$$X_i(P) = X_{i,0} \left(\frac{P}{P_0}\right)^{f_{sed}}. \qquad (10)$$

At pressures higher than $P_0$, the cloud is not condensed, and thus $X_i = 0$.

Finally, 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.

For 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 off the cloud opacity sources, and calculated another clear atmosphere spectrum, $S_{cl}$. We then combined the two spectra, weighted by $f_c$:

$$S = f_c S_{cd} + (1 - f_c) S_{cl}. \qquad (11)$$

This 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 different patchiness fractions for

different cloud layers, allowing for different 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.

### 4.4. Retrieval setup

We 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 efficiency to 0.05, and used constant efficiency mode in order to reduce computation time. Comparisons to retrievals using 4000 live points and a sampling efficiency of 0.8 without constant efficiency showed that this choice does not bias the posterior estimates, and that the importance nested sampling evidence estimate is of sufficient 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 different 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.

**Group 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 than 1 $R_{Jup}$. Even high-mass, cold brown dwarfs over 1 Gyr in age are found to have minimum radii of ~0.88 $R_{Jup}$, 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 $R_{Jup}$. 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 $M_{Jup}$ and less than 22 $M_{Jup}$, 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.

**Group 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
<figure><img src="image_5.png" /><figcaption>Fig. 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.</figcaption></figure>

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 effectively reduces to that of Group A.

### Group 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.

The 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) for an interpretation of the Bayes factor in terms of frequentist statistical significance. Thus a Δlog₁₀ Z > 1 is considered substantial evidence, and Δlog₁₀ 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₁₀ 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.

Consider Bayes theorem for the ith model Mᵢ for data D, with parameters θᵢ:

$$P(\theta_i|D, M_i) = \frac{P(D|\theta_i, M_i)P(\theta_i|M_i)}{P(D|M_i)} \quad (12)$$

We are interested in obtaining 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(θᵢ|D, Mᵢ), the likelihood P(D|θᵢ, Mᵢ), and the evidence P(D|Mᵢ) = Z = ∫ P(D|θᵢ, Mᵢ)P(θᵢ|Mᵢ)dθᵢ. 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(Mᵢ). 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:

$$P(M_i) = \frac{1}{N} \quad (13)$$

This allows the data to determine which models should be favoured based on the evidence.

Considering all models in the range 1 to N, the posterior model probabilities given the data are

$$P(M_i|D) = \frac{P(D|M_i)P(M_i)}{\sum_{j=1}^{N} P(D|M_j)P(M_j)} \quad (14)$$

The marginal posterior distribution for a single parameter θ present in all of the models is thus

$$P(\theta|D) = \sum_{j=1}^{N} P(\theta|D, M_j)P(M_j|D) \quad (15)$$

This 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.
<table><caption>**Table 7.** Self-consistent grid boundaries and step sizes.</caption><thead><tr><th>Parameter</th><th>ATMO</th><th>Diamondback</th><th>Exo-REM</th><th>petitCODE, Cool</th><th>petitCODE, Hot</th></tr></thead><tbody><tr><td>T<sub>eff</sub> (K)</td><td>[800, 3000]<br>50</td><td>[900, 2400]<br>100</td><td>[400, 2000]<br>50</td><td>[500, 850]<br>50</td><td>[1000, 1800]<br>100</td></tr><tr><td>log g (dex)</td><td>[2.5, 5.5]<br>0.1</td><td>[3.5, 5.5]<br>0.5</td><td>[3.0, 5.0]<br>0.5</td><td>[3.0, 5.0]<br>0.5</td><td>[3.5, 5.5]<br>0.5</td></tr><tr><td>[M/H]</td><td>[-0.6, 0.6]<br>0.3</td><td>[-0.5, 0.5]<br>0.5</td><td>[-0.5, 1.0]<br>0.5</td><td>[0.0, 1.4]<br>0.2</td><td>[-0.3, 0.3]<br>0.3</td></tr><tr><td>C/O</td><td>[0.3, 0.7]<br>0.24</td><td>0.458<br>...</td><td>[0.1, 0.8]<br>0.05</td><td>0.55<br>...</td><td>[0.55, 0.75, 0.90]<br>...</td></tr><tr><td>f<sub>sed</sub></td><td>...</td><td>[1, 2, 3, 4, 8]</td><td>...</td><td>[0.5, 2.0]</td><td>[1.5, 4.5]</td></tr><tr><td>γ<sub>ad</sub></td><td>[1.01, 1.05]<br>0.02</td><td>...</td><td>...</td><td>0.5</td><td>1.5</td></tr></tbody></table>

### 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 effective adiabatic index $\gamma_{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 (Exo-REM, Baudino et al. 2015, 2017) is a self-consistent 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 effective 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 effective 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.

Sonora 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/O ratio 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 Na<sub>2</sub>S and KCl clouds, and the hot grid adds Mg<sub>2</sub>SiO<sub>4</sub> 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. (2020a), 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.
<figure><img src="image_6.png" /><figcaption>**Fig. 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 different 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).</figcaption></figure>

In addition to the Bayesian fits, we performed a simple $\chi^2$ minimisation over each grid to avoid potential issues with interpolating the spectra along the different 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, with 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 $\chi^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 $\chi^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.

Nearly 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 $\chi^2/\nu < 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.

Both 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 the carbon and oxygen abundances as measured through H$_2$O and CO, and particularly through the CO absorption feature at 2.3 µm. The C/O ratios 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 CH$_4$ is detected in HR 8799 c.

We 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 Na$_2$S clouds.

The 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
<table><caption>**Table 8.** Grid-retrieval results.</caption><thead><tr><th>Planet</th><th>Model</th><th>Δ log<sub>10</sub> Z</th><th>T<sub>eff</sub><br>(K)</th><th>log g<br>(cgs)</th><th>[M/H]</th><th>f<sub>sed</sub></th><th>C/O</th><th>Radius<br>(R<sub>Jup</sub>)</th><th>log L/L<sub>☉</sub></th><th>Mass<br>(M<sub>Jup</sub>)</th></tr></thead><tbody><tr><td rowspan="4">b</td><td>ATMO</td><td>-63</td><td>1020<sup>+9</sup><sub>-7</sub></td><td>4.42<sup>+0.04</sup><sub>-0.04</sub></td><td>>0.6</td><td>...</td><td>&lt;0.3</td><td>0.73<sup>+0.02</sup><sub>-0.02</sub></td><td>-5.278<sup>+0.008</sup><sub>-0.009</sub></td><td>5.6<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>Diamondback</td><td>-98</td><td>959<sup>+10</sup><sub>-8</sub></td><td>4.18<sup>+0.03</sup><sub>-0.03</sub></td><td>>0.5</td><td>&lt;1.0</td><td>0.458</td><td>0.95<sup>+0.02</sup><sub>-0.02</sub></td><td>-5.155<sup>+0.007</sup><sub>-0.007</sub></td><td>5.6<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>Exo-Rem</td><td>0</td><td>862<sup>+21</sup><sub>-6</sub></td><td>4.0<sup>+0.02</sup><sub>-0.03</sub></td><td>0.69<sup>+0.11</sup><sub>-0.03</sub></td><td>...</td><td>0.50<sup>+0.10</sup><sub>-0.02</sub></td><td>1.19<sup>+0.03</sup><sub>-0.03</sub></td><td>-5.14<sup>+0.02</sup><sub>-0.02</sub></td><td>5.5<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>petitCODE<sup>(a)</sup></td><td>-105</td><td>>850</td><td>3.48<sup>+0.01</sup><sub>-0.02</sub></td><td>>1.4</td><td>2.53<sup>+0.01</sup><sub>-0.02</sub></td><td>0.55</td><td>2.00<sup>+0.02</sup><sub>-0.04</sub></td><td>-4.73<sup>+0.01</sup><sub>-0.02</sub></td><td>4.8<sup>+0.8</sup><sub>-0.2</sub></td></tr><tr><td rowspan="4">c</td><td>ATMO</td><td>0</td><td>1195<sup>+5</sup><sub>-6</sub></td><td>4.34<sup>+0.03</sup><sub>-0.03</sub></td><td>>0.6</td><td>...</td><td>0.425<sup>+0.012</sup><sub>-0.014</sub></td><td>1.01<sup>+0.01</sup><sub>-0.01</sub></td><td>-4.723<sup>+0.003</sup><sub>-0.003</sub></td><td>9.1<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>Diamondback</td><td>-245</td><td>1237<sup>+2</sup><sub>-2</sub></td><td>&lt;3.5</td><td>>0.5</td><td>2.32<sup>+0.02</sup><sub>-0.02</sub></td><td>0.458</td><td>0.99<sup>+0.01</sup><sub>-0.01</sub></td><td>-4.683<sup>+0.002</sup><sub>-0.002</sub></td><td>1.25<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>Exo-Rem</td><td>-30</td><td>1065<sup>+1</sup><sub>-1</sub></td><td>4.15<sup>+0.02</sup><sub>-0.02</sub></td><td>>1.0</td><td>...</td><td>0.501<sup>+0.003</sup><sub>-0.002</sub></td><td>1.314<sup>+0.005</sup><sub>-0.005</sub></td><td>-4.694<sup>+0.002</sup><sub>-0.002</sub></td><td>9.8<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>petitCODE</td><td>-622</td><td>1400<sup>+0.2</sup><sub>-0.5</sub></td><td>4.00<sup>+0.01</sup><sub>-0.01</sub></td><td>>0.3</td><td>>4.5</td><td>0.55<sup>+0.00</sup><sub>-0.00</sub></td><td>0.80<sup>+0.00</sup><sub>-0.00</sub></td><td>-4.652<sup>+0.002</sup><sub>-0.002</sub></td><td>2.57<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td rowspan="4">d</td><td>ATMO</td><td>-38</td><td>1300<sup>+2</sup><sub>-2</sub></td><td>4.39<sup>+0.01</sup><sub>-0.01</sub></td><td>0.07<sup>+0.02</sup><sub>-0.02</sub></td><td>...</td><td>&lt;0.3</td><td>0.926<sup>+0.004</sup><sub>-0.004</sub></td><td>-4.652<sup>+0.003</sup><sub>-0.003</sub></td><td>8.4<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Diamondback</td><td>-76</td><td>1234<sup>+8</sup><sub>-7</sub></td><td>4.27<sup>+0.02</sup><sub>-0.02</sub></td><td>>0.5</td><td>1.95<sup>+0.04</sup><sub>-0.06</sub></td><td>0.458</td><td>1.05<sup>+0.01</sup><sub>-0.02</sub></td><td>-4.632<sup>+0.002</sup><sub>-0.002</sub></td><td>8.4<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Exo-Rem</td><td>0</td><td>1155<sup>+5</sup><sub>-4</sub></td><td>4.23<sup>+0.01</sup><sub>-0.01</sub></td><td>0.78<sup>+0.03</sup><sub>-0.03</sub></td><td>...</td><td>0.200<sup>+0.002</sup><sub>-0.002</sub></td><td>1.168<sup>+0.006</sup><sub>-0.005</sub></td><td>-4.656<sup>+0.007</sup><sub>-0.005</sub></td><td>9.3<sup>+0.6</sup><sub>-0.2</sub></td></tr><tr><td>petitCODE</td><td>-185</td><td>1400<sup>+1</sup><sub>-1</sub></td><td>4.47<sup>+0.01</sup><sub>-0.01</sub></td><td>>0.3</td><td>2.99<sup>+0.06</sup><sub>-0.07</sub></td><td>0.55<sup>+0.000</sup><sub>-0.000</sub></td><td>0.859<sup>+0.003</sup><sub>-0.003</sub></td><td>-4.587<sup>+0.003</sup><sub>-0.002</sub></td><td>8.9<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td rowspan="4">e</td><td>ATMO</td><td>0</td><td>1251<sup>+22</sup><sub>-20</sub></td><td>4.38<sup>+0.06</sup><sub>-0.07</sub></td><td>>0.6</td><td>...</td><td>0.52<sup>+0.03</sup><sub>-0.03</sub></td><td>0.87<sup>+0.03</sup><sub>-0.03</sub></td><td>-4.778<sup>+0.006</sup><sub>-0.006</sub></td><td>7.2<sup>+0.8</sup><sub>-0.8</sub></td></tr><tr><td>Diamondback</td><td>-16</td><td>1244<sup>+6</sup><sub>-6</sub></td><td>&lt;3.5</td><td>>0.5</td><td>2.31<sup>+0.05</sup><sub>-0.05</sub></td><td>0.458</td><td>0.92<sup>+0.01</sup><sub>-0.01</sub></td><td>-4.738<sup>+0.005</sup><sub>-0.005</sub></td><td>1.12<sup>+0.07</sup><sub>-0.05</sub></td></tr><tr><td>Exo-Rem</td><td>-14</td><td>1064<sup>+3</sup><sub>-3</sub></td><td>4.13<sup>+0.04</sup><sub>-0.04</sub></td><td>>1.0</td><td>...</td><td>0.50<sup>+0.01</sup><sub>-0.01</sub></td><td>1.24<sup>+0.01</sup><sub>-0.01</sub></td><td>-4.747<sup>+0.006</sup><sub>-0.005</sub></td><td>8.0<sup>+0.8</sup><sub>-0.8</sub></td></tr><tr><td>petitCODE</td><td>-56</td><td>1399<sup>+2</sup><sub>-4</sub></td><td>4.00<sup>+0.02</sup><sub>-0.03</sub></td><td>>0.3</td><td>4.48<sup>+0.02</sup><sub>-0.03</sub></td><td>0.55<sup>+0.00</sup><sub>-0.00</sub></td><td>0.755<sup>+0.005</sup><sub>-0.004</sub></td><td>-4.703<sup>+0.004</sup><sub>-0.004</sub></td><td>2.3<sup>+0.1</sup><sub>-0.1</sub></td></tr></tbody></table>

**Notes.** All values presented are the median values from the fits, with uncertainties given as the ±34.1% percentiles. Values without uncertainties were fixed during the fit. (a) the petitcode-cool-cloudy grid was used to fit HR 8799 b.

<table><caption>**Table 9.** Grid-fit $\chi^2$ results.</caption><thead><tr><th>Planet</th><th>Model</th><th>χ<sup>2</sup></th><th>T<sub>eff</sub><br>(K)</th><th>log g<br>(cgs)</th><th>[M/H]</th><th>f<sub>sed</sub></th><th>C/O</th><th>γ<sub>ad</sub></th><th>Radius<br>(R<sub>Jup</sub>)</th></tr></thead><tbody><tr><td rowspan="4">b</td><td>ATMO</td><td>744</td><td>1000</td><td>4.5</td><td>0.6</td><td>...</td><td>0.3</td><td>1.05</td><td>1.00</td></tr><tr><td>Diamondback</td><td>879</td><td>1100</td><td>3.5</td><td>0.5</td><td>2.0</td><td>0.458</td><td>...</td><td>1.07</td></tr><tr><td>Exo-REM</td><td>510</td><td>850</td><td>3.5</td><td>0.5</td><td>...</td><td>0.55</td><td>...</td><td>1.05</td></tr><tr><td>petitCODE<sup>(a)</sup></td><td>968</td><td>850</td><td>3.5</td><td>1.4</td><td>2.0</td><td>...</td><td>...</td><td>0.79</td></tr><tr><td rowspan="4">c</td><td>ATMO</td><td>2162</td><td>1200</td><td>4.5</td><td>0.6</td><td>...</td><td>0.3</td><td>1.01</td><td>0.97</td></tr><tr><td>Diamondback</td><td>3900</td><td>1200</td><td>3.5</td><td>0.5</td><td>2.0</td><td>0.458</td><td>...</td><td>1.08</td></tr><tr><td>Exo-REM</td><td>1868</td><td>1100</td><td>3.5</td><td>1.0</td><td>...</td><td>0.8</td><td>...</td><td>1.22</td></tr><tr><td>petitCODE</td><td>7374</td><td>1200</td><td>3.5</td><td>0.3</td><td>3.0</td><td>0.55</td><td>...</td><td>1.09</td></tr><tr><td rowspan="4">d<sup>(b)</sup></td><td>ATMO</td><td>943</td><td>1300</td><td>3.0</td><td>0.6</td><td>...</td><td>0.55</td><td>1.01</td><td>0.94</td></tr><tr><td>Diamondback</td><td>1336</td><td>1200</td><td>3.5</td><td>0.5</td><td>1.0</td><td>0.458</td><td>...</td><td>1.12</td></tr><tr><td>Exo-REM</td><td>936</td><td>1200</td><td>3.0</td><td>1.0</td><td>...</td><td>0.55</td><td>...</td><td>1.12</td></tr><tr><td>petitCODE</td><td>2254</td><td>1200</td><td>4.0</td><td>0.0</td><td>1.5</td><td>0.55</td><td>...</td><td>1.21</td></tr><tr><td rowspan="4">e</td><td>ATMO</td><td>669</td><td>1300</td><td>4.0</td><td>0.6</td><td>...</td><td>0.55</td><td>1.03</td><td>0.81</td></tr><tr><td>Diamondback</td><td>819</td><td>1200</td><td>3.5</td><td>0.5</td><td>2.0</td><td>0.458</td><td>...</td><td>1.02</td></tr><tr><td>Exo-REM</td><td>623</td><td>1100</td><td>3.5</td><td>1.0</td><td>...</td><td>0.8</td><td>...</td><td>1.15</td></tr><tr><td>petitCODE</td><td>1136</td><td>1400</td><td>4.0</td><td>0.3</td><td>4.5</td><td>0.55</td><td>...</td><td>0.74</td></tr></tbody></table>

**Notes.** All values presented are the single best-fit value according to the $\chi^2$. All of the models share a similar number of parameters, which are listed in Table 7. Additional parameters for data scaling and covariance fitting are shared between all models. (a) for HR 8799 b, the petitcode-cool-cloudy grid was used rather than the petitcode-hot-cloudy grid. (b) the fits for HR 8799 d included a scaling parameter for the ALES data.
<figure><img src="image_7.png" /><figcaption>**Fig. 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.</figcaption></figure>

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<sub>Jup</sub> compared to the
dynamical mass of 7.7 ± 0.7 M<sub>Jup</sub>. 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 mass estimate, highlighting the importance of includ-
ing additional constraints in a retrieval framework. The estimate
of the effective 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
<figure><img src="image_8.png" /><figcaption>**Fig. 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.</figcaption></figure>

brown dwarfs, consistent with their inferred effective tempera-
ture. HR 8799 c is more likely fit by an early-T spectral type as
it is a few kelvin cooler than e and d and is beginning to show
spectral features from CH₄. At 950 K HR 8799 b is solidly in the
T dwarf regime.

In this section, we present the measured properties for each
of the four planets individually. Following that, in Sect. 5.5
we present a detailed discussion of the results and challenges
of comparing different 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

data 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.

HR 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 effective temper-
ature of 942<sup>+12</sup><sub>-16</sub> K and a mass of 6.0<sup>+0.4</sup><sub>-0.3</sub> M<sub>J</sub>, 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
(R<sub>pl</sub> = 1.10<sup>+0.03</sup><sub>-0.03</sub>); combining the mass and radius estimates leads
to log *g* = 4.10<sup>+0.03</sup><sub>-0.04</sub>. This is consistent with the Bayesian grid fits,
though the single-best χ<sup>2</sup> 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<sub>10</sub> Z
of 4 relative to the next best retrieval. This single best retrieval
<figure><img src="image_9.png" /><figcaption>**Fig. 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.</figcaption></figure>

uses the Z23 temperature profile and free chemistry, finding a metallicity of 0.96<sup>+0.08</sup><sub>-0.08</sub> and a C/O of 0.78<sup>+0.03</sup><sub>-0.04</sub>.

The grid fits find temperatures between 850 and 1100 K, with the single best fit, found using Exo-REM, finding $T_{\text{eff}} = 850$ K. Using the Bayesian fit, Exo-REM is again the most favoured model by the Bayes factor, and finds $T_{\text{eff}} = 862_{-6}^{+21}$ K and a radius of $1.19_{-0.03}^{+0.03} R_{\text{Jup}}$, consistent with expectations from evolutionary models (e.g. Marley et al. 2012). The ATMO and Sonora Dia-mondback models finds a somewhat higher temperature and a smaller radius, but are disfavoured by the Bayes factor. petit-CODE is divided into **cool-cloudy** and **hot-cloudy** grids, with the cool grid extending up to 850 K, and the hot grid begin-ning at 1000 K. As the effective temperature of HR 8799 b likely falls between these grids, it poorly fit the data, although the remaining parameters of log g and $f_{\text{sed}}$ are compatible with the other self-consistent models.

Figure 7 shows the variation of the metallicity across the dif-ferent 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.

Similarly to the metallicity, the C/O ratio is shown for the different retrievals in Fig. 8 and $T_{\text{eff}}$ in Fig. 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 a lower C/O ratio. 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
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** *T*<sub>eff</sub> 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.</figcaption></figure>

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. H₂S and CH₄ 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 posterior 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 CH₄ and H₂S respectively. While the inferred CH₄ abundance of log X<sub>CH4</sub> = -5.0 ± 0.4 is compatible with Barman et al. (2015), there is no statistical evidence to support the detection.

In 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 Na₂S clouds are preferred over silicate, iron, and KCl clouds (Δlog₁₀ Z ≥ 2). Patchy clouds were also explored, but no evidence was found for patchiness, regardless of cloud composition.
<figure><img src="image_11.png" /><figcaption>**Fig. 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 $\chi^2$ model from the grid.</figcaption></figure>

At the temperature of HR 8799 b silicate clouds are expected to occur below the photosphere, with Na₂S or KCl clouds becoming the primary aerosol opacity source. This result is likely driven by the condensation temperature of Na₂S rather than the optical properties; unlike silicate clouds which have strong absorption features in the mid-infrared, crystalline Na₂S 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 binning 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.

Like HR 8799 b, the group A ∩ B retrievals of HR 8799 c are dominated by a single retrieval, with $\Delta \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_{\text{base}} = -3.4 \pm 1.8$ bar). The inferred effective temperature of $1159_{-12}^{+11}$ 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.5}^{+0.4} M_{\text{Jup}}$) is slightly higher than the dynamical mass estimate, though radius ($1.10 \pm 0.01 R_{\text{Jup}}$) is compatible with evolutionary models, and from these we derive a log g of $4.26_{-0.03}^{+0.02}$. The grid-fits also found plausible radii, favouring values slightly larger than $1 R_{\text{Jup}}$. All of the self-consistent models found temperatures of 1100–1200 K,
<figure><img src="image_12.png" /><figcaption>**Fig. 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}$.</figcaption></figure>

and $\log g$ between 3.5 and 4.5. From the single-best $\chi^2$ fits, the ATMO model found a $\log g$ of 4.5, while the remaining models find a lower solution of 3.5.

As with the other HR 8799 planets the retrievals favour highly enriched solutions, finding $[\text{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 \pm 0.04$. The data for HR 8799 c are highly discrepant in the H band (Fig. 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 $\chi^2$ minimisation are consistent, typically favouring stellar to slightly substellar C/O. The retrievals present a more consistent picture, with most retrievals favouring a mildly super-stellar C/O ratio, with the average group $A \cap B$ C/O of $0.60^{+0.02}_{-0.01}$. The most favoured free chemistry retrieval for HR 8799 c find water and CO abundances consistent with most favoured disequilibrium retrieval. The free retrieval also finds an extremely high HCN mass fraction of

$\log X_{\text{HCN}} = -2.54 \pm 0.05$, orders of magnitude higher than the predictions from equilibrium chemistry. This finding is strongly favoured by the Bayes factor, with $\log_{10} Z = 30$ in favour of including HCN (c.AB.3 over c.A.8), with both retrievals using the Z23 temperature profile and clouds condensing at their equilibrium locations. The detection of HCN was largely driven by the ALES data; the wavelength dependence of the detection is discussed further in Sect. 4.2. If HCN is excluded, the overall metallicity is also increased, mostly due to a 3 dex increase in the H$_2$S abundance to $\log X_{\text{H}_2\text{S}} = -2.38 \pm 0.06$. While this solution is disfavoured, the higher H$_2$S abundance is more compatible with the equilibrium chemistry predictions. High resolution spectroscopy is likely required to precisely characterise the sulphur and nitrogen elemental abundances, and determine reliable abundances for these trace species. Although the CH$_4$ abundance is relatively low, with $\log X_{\text{CH}_4} = -4.3 \pm 0.06$, it is precisely constrained and detected with high confidence, $\Delta \log_{10} Z = 11.5$ (c.AB.3 over c.A.6).

HR 8799 c is host to a highly cloudy atmosphere. The most favoured retrieval finds high altitude ($\log P_{\text{MgSiO}_3} = -3.4 \pm 1.8$ bar) MgSiO$_3$ cloud with a mass fraction of $\log X_{\text{MgSiO}_3} = -4.7 \pm 1.2$, together with a deep iron deck. The vertical mixing strength for the clouds is $\log K_{zz} = 8.0 \pm 0.9$, while the $f_{\text{sed}}$
for both the silicate and iron clouds are compatible, with values between 5 and 6. Cloud composition could not be robustly determined for HR 8799 c, due to difficulties in retrieval convergence. However, we find that crystalline MgSiO₃ clouds (c.AB.5) provide a better fit by the $\chi^2$, and are favoured by the Bayes factor over patchy amorphous MgSiO₃ (c.AB.8), both using the M20 profile and disequilibrium chemistry. The use of patchy cloud layers may improve this fit, allowing individual layers to impact the spectrum independently, but this would come at the cost of substantially increasing the number of parameters to fit the continuum shape. The crystalline morphology provide a marginally better fit to the MIRI photometric data, but spectroscopic characterisation of the silicate feature is likely necessary to robustly distinguish these cases.

We ran a disequilibrium and free chemistry retrieval using the OSIRIS data in place of the GRAVITY data to check for consistency and to determine if we could reproduce the findings of Konopacky et al. (2013). In order to use the correlated-k method of pRT, we binned the OSIRIS data by a factor of 4 to a spectral resolving power of $R \sim 1000$. We find that a high metallicity solution is still found using these data, with effective temperatures, and surface gravities consistent with the GRAVITY retrievals. For the disequilibrium retrieval, the C/O ratio is found to be significantly higher than any of the GRAVITY retrievals, as well as the free chemistry OSIRIS retrieval. While the free retrievals using the GRAVITY data find a slightly higher metallicity overall, the OSIRIS data finds a slightly higher CH₄ abundance of -3.81 ± 0.09, as well as a much higher H₂S abundance. Overall, we find that the main findings of metal-rich atmospheres are reproducible regardless of whether we use the GRAVITY or OSIRIS data, and confirm the detection of water and CO in the atmosphere of HR 8799 c.

### 5.3. HR 8799 d

For HR 8799 d, we included all available spectroscopic data as described in Sect. 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.

Unlike HR 8799 b or c, d is well fit by a broad selection of models, and no single retrieval dominates the Bayesian average of Fig. 5. All of the preferred models ($\Delta \log_{10} Z < 2$) used disequilibrium chemistry and the M20 or Z23 temperature profiles. Using the Bayesian average, the retrieved effective temperature is $1179^{+31}_{-28}$ K, compatible with the ranges found by the grid fits, which find $T_{\text{eff}}$ from 1150–1300 K. The planet mass ($9.2 \pm 0.1 \ M_{\text{Jup}}$) is tightly constrained by the dynamical mass prior, which in turn allowed for precise measurement of the planet radius ($1.26 \pm 0.07 \ R_{\text{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 $\chi^2$ fits. The self-consistent fits find marginally smaller radii than the retrievals, generally between 1.1 and 1.2 $R_{\text{Jup}}$.

The metallicity of HR 8799 d is consistent with HR 8799 b and c, with $[\text{M/H}] = 1.2 \pm 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 \pm 0.004 \ R_{\text{Jup}}$). The best fit Exo-Rem model find a metallicity of $0.78 \pm 0.03$, a radius of $1.168 \pm 0.006 \ R_{\text{Jup}}$ and an effective temperature of $1155 \pm 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.

Patchy clouds are marginally disfavoured by the Bayes factor, and the patchiness is poorly constrained, finding $f_{\text{cloud}} = 0.45 \pm 0.3$. The most favoured solutions require either amorphous or crystalline MgSiO₃ or crystalline MgSiO₃, with a marginal preference for the amorphous structure. Each of these compositions displays slightly different near infrared slopes, shown in Fig. 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 sufficient to describe the clouds in this atmosphere.

While 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\sigma$. 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 CH₄ 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 \cap 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 effective 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 was determined by the dynamical mass prior, as was the planet radius, finding $M_{\text{pl}} = 7.5^{+0.6}_{-0.6}$ and $R_{\text{pl}} = 1.13^{+0.05}_{-0.05}$ respectively. This leads to a $\log g$ of $4.18^{+0.06}_{-0.05}$, slightly higher than the self-consistent estimates of $3.5-4.0$. The overall best self-consistent model by the $\chi^2$ was Exo-Rem, which finds an effective temperature of 1100K, a radius of 1.15 $R_{\text{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 Fig. 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. 2024).

Compared to the other three planets, HR 8799 e is found to have an even more metal rich atmosphere, with $[\text{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 \cap 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 sequestered in the silicate clouds that could reduce the C/O ratio. However, the most favoured disequilibrium retrieval finds a similar value of $0.83 \pm 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.

In addition to the water and CO rich atmosphere, HCN is found to be highly abundant, with $log X_{HCN} = -2.26 \pm 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/O ratio (Giacobbe et al. 2021), and can be produced through photochemical reactions (Moses et al. 2013). In contrast, the presence of CH$_4$ is poorly constrained.

As 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) found 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*<sub>sed</sub> 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 Sect. 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 different 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 Fig. 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 spline 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 different temperature parameterisations. While there are statistically significant variations in the C/O ratio between the different

<figure><img src="image_13.png" /><figcaption>**Fig. 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.</figcaption></figure>

parameterisations, they remain broadly consistent between 0.7 and 0.9.

Although 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 a model with clouds condensing at the equilibrium base pressure. We find no significant differences 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 H$_2$O–CO–CH$_4$ quenching, as well as free chemistry retrievals where we directly retrieved the mass-fraction abundance of various species. Figure 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
<figure><img src="image_14.png" /><figcaption>Fig. 13. Metallicities and elemental number ratios for the HR 8799 planets derived from Baysian averaged group A ∩ B free chemistry retrievals.</figcaption></figure>

a few species have statistically significant detections in the free retrievals.

The 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 H₂S measurement, though it is not statistically significant. CH₄ is significantly detected in the atmosphere of HR 8799 c, demonstrating that with sufficient S/N and wavelength coverage, it is possible to constrain abundances at below 10⁻⁴ by mass. The free chemistry detection of CH₄ 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 Figs. C.2 and C.3 comparisons between the HR 8799 c data and models both with and without the contribution of CH₄ 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. As the ALES data are relatively low S/N, additional H and L band data should be obtained to confirm this detection. However, the CH₄ 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 CO₂, NH₃, and H₂S. Additional wavelength coverage or higher spectral resolution may allow for the characterisation of such species.

If 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

$$C/H = \frac{X_{CO} + X_{CO_2} + X_{CH_4} + X_{HCN}}{X_{H_2} + 2X_{H_2O} + 2X_{H_2S} + X_{HCN} + 4X_{CH_4} + 3X_{NH_3}}, \quad (16)$$

where all abundances are measured in number fraction. These ratios are then normalised to the solar values from Asplund et al. (2009).

In 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 non-detection of NH₃, 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 N₂ and NH₃, 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 H₂S abundance. HR 8799 b appears sulphur rich, while the remaining planets appear consistent with the solar value, or slightly depleted in sulphur, though this is largely due to a lack of measured chemical species.

The 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}_{-0.3}\%$. They also identify a relation between the bulk and observed C/O ratio:

$$ (C/O)_{\text{obs}} \approx \frac{(C/O)_{\text{bulk}}}{1 - 0.371 (C/O)_{\text{bulk}}}. \qquad (17) $$

Solving for $(C/O)_{\text{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/O ratio more in line with the known stellar value.

While 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_{\text{sed}} = 1$. By parameterising disequilibrium via fixing the chemical abundances above a quench point, the model may lose this sensitivity, and therefore measures the abundances of CO, H$_2$O, and CH$_4$ 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/O ratio. 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.

In addition to the elemental ratios, we also computed $Z_{\text{pl}}/Z_{\star}$, 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_{\text{pl}}$ using the methods of Thorngren & Fortney (2019), adapting for our own notation:

$$ 10^{[\text{M/H}]} = \frac{1 + Y/X}{(Z_{\text{pl}}^{-1} - 1)\left(\frac{\mu_Z}{\mu_H}\right)}, \qquad (18) $$

where X, Y, and Z are the solar hydrogen, helium and metal mass fractions and $\mu$ is the mean molecular weight of the metal content of the atmosphere. Rearranging and substituting in the measured atmospheric metallicity [M/H], we find:

$$ Z_{\text{pl}} = \left( 1 + \frac{1 + Y/X}{10^{[\text{M/H}]} \left( \frac{\mu_Z}{\mu_H} \right) (Z/H_{\odot})} \right)^{-1}. \qquad (19) $$

We take the same assumptions as Thorngren & Fortney (2019), taking $\mu_Z$ to be 18, assuming most of the metal content is in water, $\mu_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_{\star}/H$ is taken to be the solar value of $Z/H_{\odot} = 1.04 \times 10^{-3}$. To normalise to the stellar metallicity we follow Thorngren et al. (2016) and calculate the $Z_{\star}$ as

$$ Z_{\star} = 0.014 \times 10^{\text{Fe/H}}. \qquad (20) $$

For HR 8799, we used solar metallicity to calculate Fe/H, but refer to the discussion in Sect. 2.1.

Disequilibrium 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 H$_2$O, CO, and CH$_4$ 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_{\text{chem}}$ and the mixing $t_{\text{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 diffusion coefficient, $K_{\text{zz}}$:

$$ t_{\text{mix}} = \frac{H^2}{K_{\text{zz}}}. \qquad (21) $$

The chemical timescale depends on the reaction rates involved. Considering the CO-CH$_4$ reaction chain, Zahnle & Marley (2014) derive a timescale at the quench point ($t_q$) 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^{[\text{M/H}]}$, and temperature $T$ in kelvin:

$$ t_{q1} = 1.5 \times 10^{-6} p^{-1} m^{-0.7} \exp(42000/T) \text{ s.} \qquad (22) $$

For weak mixing, and therefore drawing from low temperatures with little CO, the timescale is found to be

$$ t_{q2} = 40p^{-2} \exp(25000/T) \text{ s,} \qquad (23) $$

Combining the two, the total chemical timescale is defined as:

$$ t_{\text{CO}} = \left( \frac{1}{t_{q1}} + \frac{1}{t_{q2}} \right)^{-1}, \qquad (24) $$

which will favour the lower of the two values $t_{q1}$ and $t_{q2}$. Equating the mixing and CO reaction timescales, we can infer the strength of vertical mixing in the atmospheres of the HR 8799 planets:

$$ K_{\text{zz}} = \frac{H^2}{t_{\text{CO}}}, \qquad (25) $$

where the scale height is defined as

$$ H = \frac{k_{\text{B}}T}{\mu g}, \qquad (26) $$

for temperature $T$, surface gravity $g$, mean molecular mass $\mu$, and the Boltzmann constant $k_B$. In order to calculate these quantities for the HR 8799 planets, we take $T_{\text{eff}}$ as a representative temperature to calculate the timescales and scale height, $g$ as the measured surface gravity, and $\mu$ as the average mean molecular weight of the atmosphere. As $K_{\text{zz}}$ is exponential in temperature, the choice of what temperature to use strongly influences the measured value. Using the temperature at the quench pressure, typically deep in the atmosphere, results in unphysically strong vertical mixing, with $\log K_{\text{zz}} \approx 20$. A more thorough analysis could try to measure the vertical mixing as a function of temperature throughout the atmosphere, but the current data quality is of insufficient resolution or S/N for such measurements. Thus
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Cloud properties for HR 8799 e, showing the optical depth due to clouds as a function of pressure (colour map), with the $\tau = 0.3$ and $\tau = 1.0$ contours highlighted by the dashed lines. The solid lines indicate the mass fraction abundance of the MgSiO$_3$ and Fe clouds in blue and red respectively. *Left*: cloud properties e.AB.13, using a freely retrieved cloud base pressure and abundances. *Right*: the same, but for Case e.AB.2, which uses the equilibrium condensation to determine the location of the cloud base. The abundances are determined using equilibrium chemistry, and retrieving a scaling factor, ($\log S_{Fe} = 0.0 \pm 1.1$, $\log S_{MgSiO_3} = -0.8 \pm 1.0$).</figcaption></figure>

<table><caption>**Table 10.** Quench pressures, vertical mixing parameters, and sedimentation fractions for the HR 8799 planets.</caption><thead><tr><th>Planet</th><th>log P<sub>q</sub><br>(bar)</th><th>K<sub>zz,q</sub><br>(cm<sup>2</sup> s<sup>-1</sup>)</th><th>K<sub>zz,AM01</sub><br>(cm<sup>2</sup> s<sup>-1</sup>)</th><th colspan="2">f<sub>sed</sub></th></tr><tr><th></th><th></th><th></th><th></th><th>MgSiO<sub>3</sub></th><th>Fe</th></tr></thead><tbody><tr><td>b</td><td>1.7<sup>+0.2</sup><sub>-0.2</sub></td><td>2.9<sup>+0.6</sup><sub>-0.7</sub></td><td>8.6<sup>+0.8</sup><sub>-1.0</sub></td><td>1.1<sup>+0.3</sup><sub>-0.3</sub></td><td>4<sup>+2</sup><sub>-1</sub></td></tr><tr><td>c</td><td>2.3<sup>+0.1</sup><sub>-0.1</sub></td><td>6.3<sup>+0.3</sup><sub>-0.5</sub></td><td>9.2<sup>+0.4</sup><sub>-0.4</sub></td><td>3.3<sup>+0.3</sup><sub>-0.3</sub></td><td>6<sup>+2</sup><sub>-2</sub></td></tr><tr><td>d</td><td>1.0<sup>+0.6</sup><sub>-0.7</sub></td><td>4.8<sup>+0.9</sup><sub>-1.3</sub></td><td>9.2<sup>+0.7</sup><sub>-0.9</sub></td><td>2<sup>+4</sup><sub>-1</sub></td><td>6<sup>+2</sup><sub>-2</sub></td></tr><tr><td>e</td><td>2.1<sup>+0.4</sup><sub>-1.2</sub></td><td>5.9<sup>+0.9</sup><sub>-0.6</sub></td><td>8.6<sup>+0.7</sup><sub>-0.8</sub></td><td>1.4<sup>+0.4</sup><sub>-0.3</sub></td><td>6<sup>+2</sup><sub>-3</sub></td></tr></tbody></table>

**Notes.** Measured from the group $A \cap B$ retrievals, providing median and $\pm 34\%$ confidence regions.

we treat $T_{eff}$ as a representative temperature with which to determine the vertical mixing strength. We include the results of these calculations, together with the retrieved $K_{zz}$ used to parameterise the AM01 clouds in Table 10.

The quench pressure is well constrained in the Bayesian average of group $A \cap B$ retrievals for all four planets; the values of which are listed in Table 10. All of the planets quench below the photosphere, with d quenching at the highest altitude, around 10 bar. Including only the disequilibrium retrievals in the Bayesian average of group $A \cap B$, we derived $K_{zz}$ from the quench pressure. We then turn to Soni & Acharyya (2023) for a comparison, who provide predictions for CH$_4$ and CO abundances for varying $K_{zz}$, $T_{eff}$, and $\log g$ across a range of metallicities. For our measured $T_{eff}$ and CH$_4$ abundances, we should expect $\log K_{zz}$ of around 6 for the warmer three companions, regardless of whether we use the measured CO or CH$_4$ abundance. For HR 8799 b a much lower value (less than ~2) is expected, assuming a 10× solar metallicity. Our inferred $K_{zz}$ values for c, d, and e are compatible with this prediction, finding $\log K_{zz}$ between 5 and 6. For HR 8799 b we also measure weak mixing, with $\log K_{zz} = 2.9^{+0.6}_{-0.7}$, which is again compatible with the predictions of Soni & Acharyya (2023). These measurements from the quench pressure are also inconsistent with the parameter used in the AM01 clouds, which require stronger vertical mixing of $\log K_{zz} \sim 9$. This discrepancy is perhaps not surprising: 3D modelling predicts that $K_{zz}$ should vary with altitude

throughout the atmosphere, and the larger cloud particles likely respond to the atmospheric motion differently than the gas phase constituents. Ultimately, more precise constraints on the thermal structure and chemical abundances, as well as trace species detections are necessary to derive a more precise vertical mixing strength. Further modelling work is also necessary to provide a more physically motivated transport model than a vertically constant eddy diffusion coefficient.

### 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 Fig. 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 8799 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 difficult to robustly distinguish between the different potential cloud structures.

There 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, $\Delta \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
<figure><img src="image_16.png" /></figure>

**Fig. 15.** Effective 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.

saturation location. In this case $a_h$ is calculated from $f_{sed}$ and $K_{zz}$, and a lower $f_{sed}$ for the $MgSiO_3$ cloud was retrieved than with the log normal distribution ($f_{sed, Hansen} = 1.18 \pm 0.20$). The effective distribution width parameter, $b_h$ was found to be 0.016, which is narrower than the distributions found by *Burningham et al. (2021)*.

In 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 $MgSiO_3$, $SiO_2$, and Fe clouds provides the best-fit model to an ultracool field dwarf, 2MASS J2224438-015852. Figure 14 shows that the $MgSiO_3$ clouds in their model are located at $10^{-3}$ bar with a maximum optical depth of $\tau = 0.3$ at 1 µm and an effective particle radius of ~0.04 µm. The $SiO_2$ clouds are slightly deeper, at $10^{-2}$ bar, and are optically thick at 1 µm. This is the same location where we find $MgSiO_3$ 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 differ in effective temperature by hundreds of kelvin.

Figure 15 highlights how the difference in particle radius contribute to the difference in the wavelength dependence of the cloud optical depth. The cloud particle radius in the AMOI 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 changes 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 Fig. 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,

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

**Fig. 16.** Cloud absorption (left) and scattering (right) opacities for different condensate compositions and structure, for 1-µm particles. Dark blue indicates $MgSiO_3$, light blue is $Mg_2SiO_4$. Solid (dashed) lines are for crystalline (amorphous) substances. The solid red (green) line is for crystalline iron ($Na_2S$).

suggesting that the impact of the small mean particle size in the upper atmosphere does not contribute strongly to the cloud opacity.

Due to the different slopes in the near-infrared opacity as a function of wavelength, as shown in Fig. 16, there is some sensitivity to different compositions and particle geometries. This can explain the mild preference for amorphous $MgSiO_3$ clouds in certain like-for-like retrieval comparisons. In reverse, the lack of features leads to a preference for $Na_2S$ clouds in HR 8799 b. Performing additional tests on HR 8799 e, we find that $Mg_2SiO_4$ (e.AB.14) are mildly disfavoured by the retrievals, and $Al_2O_3$ (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. Discussion

## 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_\star$ using Eqs. (18) and (19), 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, and 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 Eqs. (18) and (19) to convert from a measurement of [M/H] to $Z_{pl}/Z_\star$. Where stellar metallicities are not available, we assumed a metallicity [Fe/H]=0. Only the 2 $M_{sup}$ 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
<figure><img src="image_18.png" /><figcaption>**Fig. 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). Notes: $\beta$ Pic b (GRAVITY Collaboration 2020a); PDS 70 b (Wang et al. 2021a); 51 Eri b (Whiteford et al. 2023); VHS 1256 b (Hoch et al. 2022); HIP 65426 b (Petrus et al. 2021); $\kappa$ 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).</figcaption></figure>

[M/H] retrievals are their consistency with the self-consistent grids, which always favour their upper limits.

To 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 differences in the *J*, *H*, and *K* bands between the different 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 $\Delta \log Z > 10$. The remaining atmospheric parameters also significantly varied between the different retrievals: the C/O ratio increases with increasing [M/H], maintaining a relatively constant abundance of H$_2$O 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.

The spectra shown in Fig. 18 show that the retrieved metallicity is strongly dependent 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 different instruments measure significantly different flux in the *H* band. Without compatible measurements in this spectral regime, robust conclusions about the metallicity are hard to draw. Additional constraints can be obtained from other wavelength regions; Lodders & Fegley (2002) find that the CO$_2$ 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.

Finally, we performed independent comparisons to the moderate resolution OSIRIS spectra presented in Ruffio 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 Ruffio 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 $\chi^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 \approx 4000$ OSIRIS data, of similar quality to the fits displayed in Ruffio et al. (2021), and are included in the appendix in Fig. 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 different datasets on the retrieval results. We first performed retrievals using only the
<figure><img src="image_19.png" /><figcaption>**Fig. 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 differ across the set of retrievals, while the bottom row shows how different atmospheric parameters vary as a function of metallicity.</figcaption></figure>

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 $\Delta \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 \pm 0.3$), effective temperature ($1143_{-32}^{+38}$ K), and C/O ratio ($0.71_{-0.2}^{+0.08}$), 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.

Conversely, 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.

## 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 efficiency. 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.

All 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, are 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
<figure><img src="image_20.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

With effective 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. Baraffe 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 NH₃ abundance, at least in HR 8799 b, as well as the debris disc observed in Boccaletti et al. (2024). 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 determine whether these four planets share a formation pathway, or whether there were different mechanisms impacting different regions of the protoplanetary disc.

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Predictions for MIRI/MRS based on the most favoured disequilibrium (solid) and free chemistry (dashed) retrievals, together with the best-fit self-consistent models from each grid.</figcaption></figure>

## 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 CO₂, CO, and CH₄ 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 Fig. 19 and for the MIRI/MRS in Fig. 20. The comparisons to the MIRI photometry shown in Fig. 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 CO₂ feature at 3.8 µm, as well as the CH₄ 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 NH₃:HCN:N₂, 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 Ruffio 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 effective 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.

We find that all four planets are strongly enriched in metals, though there is still discrepancy between different models in constraining the precise value. This was validated by running retrievals using different temperature profiles and chemical models, and comparing to self-consistent grids. Further self-consistent modelling is necessary, particularly to extend model grids out to high metallicities. The C/O ratio is stellar to stellar 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 low-S/N data from LBT/ALES. CH₄ is also confidently detected in HR 8799 c for the first time. From the disequilibrium chemistry retrievals, H₂S 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<sub>zz</sub> values compatible with high-metallicity predictions from Soni & Acharyya (2023). The mixing strength is stronger for the warmer planets, at log K<sub>zz</sub> ≈ 6, and is lower for HR 8799 b with log K<sub>zz</sub> ≈ 2.

All 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 Na₂S clouds. All of the planets have effective temperatures consistent with literature values, with HR 8799 b still unique in its lower temperature and mass compared to its siblings.

We emphasise the use of robust model comparison in this work: while it may be difficult 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.

the HR 8799 planets share many similarities, much like our own Solar System there are differences in their atmospheric properties, which require further study.

Acknowledgements. 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. S.L. 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) and 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. 2001, ApJ, 556, 357

Allard, F., Homeier, D., & Freytag, B. 2012, Philos. Trans. Roy. Soc. Lond. Ser. A, 370, 2765

Allard, N. F., Spiegelman, F., & Kielkopf, J. F. 2016, A&A, 589, A21

Allard, N. F., Spiegelman, F., Leininger, T., & Mollière, P. 2019, A&A, 628, A120

Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948

Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40

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

Azzam, A. A., Yurchenko, S. N., Tennyson, J., & Naumenko, O. V. 2016, MNRAS, 460, 4063

Baines, E. K., White, R. J., Huber, D., et al. 2012, ApJ, 761, 57

Balmer, W. O., Pueyo, L., Stolker, T., et al. 2023, ApJ, 956, 99

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Barber, R. J., Strange, J. K., Hill, C., et al. 2014, MNRAS, 437, 1828

Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65

Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61

Barton, E. J., Yurchenko, S. N., & Tennyson, J. 2013, MNRAS, 434, 1469

Baudino, J. L., Bézard, B., Boccaletti, A., et al. 2015, A&A, 582, A83

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

Bayo, A., Rodrigo, C., Barrado y Navascués, D., et al. 2008, A&A, 492, 277

Bell, T. J., Welbanks, L., Schlawin, E., et al. 2023, Nature, 623, 709

Benneke, B., & Seager, S. 2012, ApJ, 753, 100

Benneke, B., & Seager, S. 2013, ApJ, 778, 153

Bergfors, C., Brandner, W., Janson, M., Köhler, R., & Henning, T. 2011, A&A, 528, A134

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, SPIE Conf. Ser., 7014, 701418

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Biller, B. A., Apai, D., Bonnefoy, M., et al. 2021, MNRAS, 503, 743

Bitsch, B., & Mah, J. 2023, A&A, 679, A11

Blain, D., Charnay, B., & Bézard, B. 2021, A&A, 646, A15

Boccaletti, A., Mâlin, M., Baudoz, P., et al. 2024, A&A, 686, A33

Bodenheimer, P. 1974, Icarus, 23, 319

Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2

Bonnefoy, M., Currie, T., Marleau, G. D., et al. 2014, A&A, 562, A111

Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58

Bonnefoy, M., Perraut, K., Lagrange, A. M., et al. 2018, A&A, 618, A63

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

Bowler, B. P., Liu, M. C., Dupuy, T. J., & Cushing, M. C. 2010, ApJ, 723, 850