© The Authors 2024

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 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 Sec. 2. The data used in this work are described in Sec. 3, while the details of the petitRADTRANS (pRT) forward model are described in Sec. 4, with the self-consistent models introduced in Sec. 4.7. The results of the atmospheric retrievals and self-consistent grid fits are presented in Sec. 5. We discuss the limitations of this study, additional work to validate our results, and the implications of our findings in Sec. 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).

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.

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 υ 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.09}^{+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_eff = 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 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; Soummer 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 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_⊙. This model produced results typical of the range of models explored; the mass estimates for each companion are: e = 7.6 ± 0.9 M_Jup, d = 9.2 ± 0.1 M_Jup, c = 7.7 ± 0.7 M_Jup, and b = 5.8 ± 0.4 M_Jup. These dynamical mass estimates, as well as those of Brandt et al. (2021), are broadly consistent with mass estimates from evolutionary models, assuming hot start conditions (Marley et al. 2012). Further astrometric analysis of the GRAVITY data will be examined in a forthcoming paper from Chavez et al. (in prep.).

2.3 Spectroscopic characterisation

In addition to the multitude of photometric observing campaigns, the spectroscopic characterisation of the HR 8799 planets has traced the development of dedicated exoplanet instrumentation, from long-slit spectrographs (Janson et al. 2010) to high-contrast integral field spectrographs (IFS) (Ingraham et al. 2014; Zurlo et al. 2016) to fibre-fed high resolution spectrometers (Wang et al. 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, low surface gravity, and low CH₄ 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₂O, CO, and tentatively CH₄. 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.05}^{+0.10}$. The Gemini/GPI instrument provided the first spectra obtained using a corona-graphic 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⁻¹ for c, and measurements of ${10.1}_{-2.7}^{+2.8}$ km s⁻¹ for d and ${15.0}_{-2.6}^{+2.3}$ km s⁻¹ for e. Wang et al. (2023) combine several of these datasets and perform a retrieval analysis to constrain the composition of HR 8799 c, finding [C/H] = ${0.55}_{-0.39}^{+0.36}$ [O/H] = ${0.47}_{-0.32}^{+0.31}$, and C/O = ${0.67}_{-0.15}^{+0.12}$. 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.13}^{+0.12}$, [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⁻¹, −11.1 ± 0.4 km s⁻¹, and −11.6 ± 0.8 km s⁻¹ respectively, placing important constraints on the allowed orbits for the planets. They also confirm the presence of water and CO, but are unable to significantly detect CH₄, 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 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 longstanding 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 enrichement of [M/H]=${0.48}_{-0.29}^{+0.25}$ and a C/O ratio of ${0.60}_{-0.08}^{+0.07}$. 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₄, and NH₃ abundances as a function of vertical mixing and effective temperature, identifying the key transition between CH₄ 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₄, assuming a solar composition. Soni & Acharyya (2023) extend this to superstellar metallicities and vertical mixing strengths, using the constraints on the CO and CH₄ abundances from Barman et al. (2015) to infer a vertical mixing strength of log K_zz ∈ [7, 10] for the 10× solar metallicity case. To better understand the planet structure, Thorngren et al. (2016) derive a mass-metallicity relationship. As the mass of the object increases, the metallicity tends to decrease, consistent with predictions of core accretion formation, as heavier objects accrete and retain more H₂ and He relative to a lower mass object. From their relationship, they predict that a 6 M_Jup planet should have a Z_pl/Z_* ratio of between 3 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 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_Jup on the allowed mass of companions interior to HR 8799 e (Goździewski & Migaszewski 2018).

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 × 10−15 W m−2 (μ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.

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

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². 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. 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_{{\lambda }_1},c_{{\lambda }_2},\dots ,c_{{\lambda }_n})}^T$ and its associated covariance matrix W arre given by: $$W={\left[{\displaystyle \sum _t}{W}_t^{-1}\right]}^{-1}$$(1) $$C=W\cdot [{\displaystyle \sum _t}{W}_t^{-1}\cdot C_t]\text{, }$$(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 Sec. 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.

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

3.2 Archival data

In addition to the new GRAVITY spectra, we also include archival data covering a broad wavelength range, presented in Fig. 3. Mollière et al. (2020) noted that the SPHERE (Zurlo et al. 2016) and GPI (Greenbaum et al. 2018) data are inconsistent with each other in the H-band. In order to reduce systematic variation and to account for correlations, we rereduce the data with up-to-date pipelines, and reprocesses the datasets optimally as described in Nasedkin et al. (2023) using KLIP(Soummer et al. 2012; Pueyo 2016). However, in order to best extract the planet signal we use 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 K1, K2 and H bands respectively. These were reduced using the standard GPI reduction pipeline (version 1.4.0), and reprocessed with KLIP using the same methods as the SPHERE data. As the new GRAVITY observations supersede the GPI data in the K-band, we only consider the GPI H-band data for this work.

3.2.3 CHARIS

Subaru/CHARIS (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 ≈ 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≈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.2925pm0.15 pc (Gaia Collaboration 2020).

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³. 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_pl with uniform priors or R_pl and M_pl, with Gaussian priors set by the dynamical mass estimates of Zurlo et al. (2022) and broad Gaussian priors centred at 1.1 R_Jup, in line with estimates from evolutionary models (Marley et al. 2012).

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 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. PchipInterpolator function. Following the prescription of Line et al. (2015), we penalised curvature in the temperature profile by adding an additional term to the likelihood function, $$\log p(\mathbf{T})=\frac{1}{2\gamma }{\displaystyle \sum _{i=1}^{N-1}}{(T_{i+1}-2T+T_{i-1})}^2-\text{log\hspace{0.17em}}(2\pi \gamma ).$$(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): $${\mathrm{\Gamma }}^{-1}(\gamma )=\frac{{\beta }^{\alpha }}{\mathrm{\Gamma }(\alpha )}{\left(\frac{1}{\gamma }\right)}^{\alpha +1}\exp (-\frac{\beta }{\gamma }),$$(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: $$\begin{array}{rl}{T}_{\text{Guillot }}^4=& \frac{3T_{\text{int }}^4}{4}(\frac{2}{3}+\tau )\\ & +\frac{3T_{\text{irr }}^4}{4}(\frac{2}{3}+\frac{1}{\gamma \sqrt{3}}+(\frac{\gamma }{\sqrt{3}}-\frac{1}{\gamma \sqrt{3}})e^{-\gamma \tau \sqrt{3}}),\end{array}$$(5)

where $T_{\text{irr }}=\sqrt{2}{T}_{\text{equ }}\text{ and }\tau =P\times {\kappa }_{\mathrm{I}\mathrm{R}}/g.T_{\text{equ }}$ is the equilibrium temperature of an irradiated body, T_int is the intrinsic temperature of the planet, and g is the surface gravity. κ_ir is the mean infrared opacity, and γ is the ratio between the optical and infrared opacities. While these parameters can be physically interpreted, we treat them as nuisance parameters that control the shape of the profile, rather than self-consistently linking them to the chemical opacities. As the HR 8799 planets are widely separated, T_irr is small, and thus the profile reduces to an Eddington (1930) profile, which corresponds to keeping only the first term on the 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 τ as a function of pressure (P): $$\tau =\delta P^{\alpha }$$(6)

and retrieve parameters of log δ and α, together with T_int, as in the G10 profile.

The upper atmosphere is defined as the region above τ = 0.1. Above this level, four pressure points are defined, equidistant in log P. The deepest pressure point, at τ = 0.1 is fixed to the temperature of the Eddington profile of the middle atmospheric region, while the remaining temperature points are freely retrieved parameters, subject to the constraint that the temperature decreases with altitude (Kitzmann et al. 2020), as inversions are not expected in self-luminous objects. The temperature profile is then interpolated from a cubic spline between the three points. Combined with the Eddington profile parameters, this results in a total of 6 parameters to describe the temperature profile.

The base of the atmosphere is defined as a moist adia-bat, up to the radiative–convective boundary. This boundary occurs when the temperature gradient of the Eddington profile is Schwarzchild unstable: $$\frac{\mathrm{d}T}{\mathrm{d}r}<\frac{T}{P}\frac{\mathrm{d}P}{\mathrm{d}r}(1-\frac{1}{{\gamma }_{\mathrm{a}\mathrm{d}}}).$$(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³ bar and 10⁻³ 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 }}$$(8) $$T_{i+1}=\exp (\log T_i+(\log P_{i+1}-\log P_i){\left(\frac{\mathrm{d}\text{\hspace{0.17em}log\hspace{0.17em}}T}{\mathrm{d}\text{\hspace{0.17em}log\hspace{0.17em}}P}\right)}_i)$$(9)

The atmosphere was isothermal above 10⁻³ 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₂O (Polyansky et al. 2018), CO (Rothman et al. 2010), CH₄ (Yurchenko et al. 2017), CO₂ (Yurchenko et al. 2020), NH₃ (Coles et al. 2019), HCN (Barber et al. 2014), H₂S (Azzam et al. 2016), PH₃ (Sousa-Silva et al. 2015), FeH (Wende et al. 2010), Na (Allard et al. 2019), K (Allard et al. 2016), SiO (Barton et al. 2013), TiO (McKemmish et al. 2019), and VO (McKemmish et al. 2016).

4.2.1 Free chemistry

In the free chemistry retrievals, we assumed a vertically constant mass fraction for each species, and retrieved the log mass fraction abundance (log Xi for each of H₂O, CO, CH₄, CO₂, HCN, H₂S, NH₃, FeH, Na, and K), subject to the constraint that the sum of the mass fractions is less than one. Due to the lack of spectroscopic data in in the Y and J bands, we did not retrieve FeH, Na or K for HR 8799 b in the free retrievals to reduce the number of free parameters. The hydrogen and helium mass fractions were calculated by using the solar abundances (0.766 for H₂, 0.234 for He), and multiplying by one minus the sum of the retrieved mass fractions (1 − ∑_iX_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₂O, CO, and CH₄, 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₂O, CO, CH₄, CO₂, NH₃, and H₂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₄, H₂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 Sec. 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 σ_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 σ_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₃ (both crystalline and amorphous particle shapes), Mg₂SiO₄, Fe, Al₂O₃, KCl, and Na₂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₀, and decreases with altitude to the power of f_sed: $$X_i(P)=X_{i,0}{\left(\frac{P}{P_0}\right)}^{f_{\text{sed }}}$$(10)

At pressures higher than P₀, 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: $$\mathit{S}=f_{\text{c}}{\mathit{S}}_{\text{cd}}+(1-f_{\text{c}}){\mathit{S}}_{\text{cl}}$$(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 $\mathcal{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 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₁₀ $\mathcal{Z}$ > 1 is considered substantial evidence, and Δ log₁₀ $\mathcal{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₁₀ $\mathcal{Z}$ = 2.1 for ‘strong’ evidence in favour of one hypothesis over another, equivalent to 3.6σ significance.

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.

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_i for data D, with parameters θ_i: $$P({\mathit{\theta }}_i\mid \mathit{D},M_i)=\frac{P(\mathit{D}\mid {\mathit{\theta }}_i,M_i)P({\mathit{\theta }}_i\mid M_i)}{P(\mathit{D}\mid M_i)}$$(12)

We are interested in obtain a joint posterior probability distribution P(θ|D, M) for the subset of parameters θ that are shared between the set of models. From each model we require posterior probability distribution P(θ_i|D, M_i), the likelihood P(D|θ_i, M_i), and the evidence P(D|M_i) ≡ $\mathcal{Z}$ = ∫ P(D|θ_i, M_i)P(θ_i|M_i)dθ_i. We then assume a prior probability distribution over the full set of models under consideration, and therefore each model has an associated prior probability P(M_i). 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\left(M_i\right)=\frac{1}{N}.$$(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\mid \mathit{D})=\frac{P(\mathit{D}\mid M_i)P\left({\mathit{M}}_{\mathit{i}}\right)}{\sum _{j=1}^{N}P(\mathit{D}\mid M_j)P\left(M_j\right)}.$$(14)

The marginal posterior distribution for a single parameter θ present in all of the models is thus $$P(\theta \mid \mathit{D})={\displaystyle \sum _{j=1}^N}P(\theta \mid \mathit{D},M_j)P(M_j\mid \mathit{D}).$$(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.

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 λ_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 Na2S and KCl clouds, and the hot grid adds Mg₂SiO₄ 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.

In addition to the Bayesian fits, we performed a simple χ² 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.

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

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 χ² 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 χ² 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 χ²|ν < 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₂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₄ 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₂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 produced masses consistent with dynamical mass measurements and radii consistent with evolutionary models. In the group A ∩ B retrievals that used dynamical mass priors, the retrieved mass estimate is within 1σ of the dynamical mass for b, d, and e, with the posterior width consistent with the prior width which. For HR 8799 c the retrieved mass is moderately higher than the dynamical mass estimate at 8.5 ± 0.5 M_Jup compared to the dynamical mass of 7.7 ± 0.7 M_Jup. However, if the dynamical mass is not used as a constraint and log g is freely retrieved, the resulting mass estimate is found to be much larger than the dynamical mass estimate, highlighting the importance of including 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 brown dwarfs, consistent with their inferred effective temperature. 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 Sec. 5.5 we present a detailed discussion of the results and challenges of comparing different thermal structures, chemistry and cloud parameterisations, and data inclusion.

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.

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.

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 temperature of ${942}_{-16}^{+12}$ K and a mass of ${6.0}_{-0.3}^{+0.4}$ M_J, 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_pl = ${1.10}_{-0.03}^{+0.03}$); combining the mass and radius estimates leads to log g = ${4.10}_{-0.04}^{+0.03}$. This is consistent with the Bayesian grid fits, though the single-best χ² 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₁₀ $\mathcal{Z}$ of 4 relative to the next best retrieval. This single best retrieval uses the Z23 temperature profile and free chemistry, finding a metallicity of ${0.96}_{-0.08}^{+0.08}$ and a C/O of ${0.78}_{-0.04}^{+0.03}$.

The grid fits find temperatures between 850 and 1100 K, with the single best fit, found using Exo-REM, finding T_eff = 850 K. Using the Bayesian fit, Exo-REM is again the most favoured model by the Bayes factor, and finds T_eff = ${862}_{-6}^{+21}$ K and a radius of ${1.19}_{-0.03}^{+0.03}$ R_Jup, consistent with expectations from evolutionary models (e.g. Marley et al. 2012). The ATMO and Sonora Diamondback models finds a somewhat higher temperature and a smaller radius, but are disfavoured by the Bayes factor. petitCODE is divided into cool-cloudy and hot-cloudy grids, with the cool grid extending up to 850 K, and the hot grid beginning at 1000 K. As the 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_sed are compatible with the other self-consistent models.

Figure 7 shows the variation of the metallicity across the different models considered. We see that strong enrichment solution is almost always found, particularly by models preferred by the Bayes factor. The degree of enrichment does not systematically vary between disequilibrium and free chemistry retrievals. All of the self-consistent models favour high metallicity solutions, reaching the upper bounds of the grid in all cases.

Similarly to the metallicity, the C/O ratio is shown for the different retrievals in Fig. 8 and T_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 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_{{\mathrm{C}\mathrm{H}}_4}$ = −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 Na2S clouds are preferred over silicate, iron, and KCl clouds (Δ log₁₀ $\mathcal{Z}$ ≥. 2). Patchy clouds were also explored, but no evidence was found for patchiness, regardless of cloud composition. At the temperature of HR 8799 b silicate clouds are expected to occur below the photosphere, with 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.

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.

Fig. 9 Teff 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.

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 χ2 model from the grid.

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 Δ log₁₀ $\mathcal{Z}$ = 2 compared to the next best retrieval. This retrieval used the Z23 temperature profile and free chemistry, and requires high-altitude silicate clouds (log P_base = −3.4 ± 1.8 bar). The inferred 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_Jup) is slightly higher than the dynamical mass estimate, though radius (1.10 ± 0.01 R_Jup) is compatible with evolutionary models, and from these we derive a log g of ${4.26}_{-0.03}^{+0.02}$. The grid-fits also found plausible radii, favouring values slightly larger than 1 R_Jup. All of the self-consistent models found temperatures of 1100–1200 K, and log g between 3.5 and 4.5. From the single-best χ² 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 [M/H] = ${1.27}_{-0.05}^{+0.05}$. The disequilibrium chemistry retrievals find slightly lower metallicities than the free chemistry retrievals, with the most-favoured disequilibrium retrieval finding a metallicity of 1.05 ± 0.04. The data for HR 8799 c are highly discrepant in the H band (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 χ² 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 ∩ B C/O of ${0.60}_{-0.01}^{+0.02}$. 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_HCN = −2.54 ± 0.05, orders of magnitude higher than the predictions from equilibrium chemistry. This finding is strongly favoured by the Bayes factor, with log₁₀ $\mathcal{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 Sec. 4.2. If HCN is excluded, the overall metallicity is also increased, mostly due to a 3 dex increase in the H₂S abundance to log $X_{{\mathrm{H}}_2\mathrm{S}}$ = −2.38 ± 0.06. While this solution is disfavoured, the higher H₂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₄ abundance is relatively low, with log $X_{{\mathrm{C}\mathrm{H}}_4}$ = −4.3 ± 0.06, it is precisely constrained and detected with high confidence, Δ log₁₀ $\mathcal{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_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}$ = −3.4 ± 1.8 bar) MgSiO₃ cloud with a mass fraction of log $X_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}$ = −4.7 ± 1.2, together with a deep iron deck. The vertical mixing strength for the clouds is log K_zz = 8.0 ± 0.9, while the f_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 χ², 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 ~ 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.

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.

5.3 HR 8799 d

For HR 8799 d, we included all available spectroscopic data as described in Sec. 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 (Δ log₁₀ $\mathcal{Z}$ < 2) used disequilibrium chemistry and the M20 or Z23 temperature profiles. Using the Bayesian average, the retrieved effective temperature is ${1179}_{-28}^{+31}$ K, compatible with the ranges found by the grid fits, which find T_eff from 1150–1300 K. The planet mass (9.2 ± 0.1 M_Jup) is tightly constrained by the dynamical mass prior, which in turn allowed for precise measurement of the planet radius (1.26 ± 0.07 R_Jup) and log g (4.18 ± 0.05). This surface gravity is consistent with the estimates from the Bayesian grid fits, but is significantly higher than the 3.0–4.0 range found by the χ² fits. The self-consistent fits find marginally smaller radii than the retrievals, generally between 1.1 and 1.2 R_Jup.

The metallicity of HR 8799 d is consistent with HR 8799 b and c, with [M/H] = 1.2 ± 0.2. While most grid-fit solutions also favoured high-metallicity atmospheres, the Bayesian fit with the ATMO model find a solution consistent with solar metallicity, though the radius was inconsistent with evolutionary models (0.926 ± 0.004 R_Jup). The best fit Exo-Rem model find a metallicity of 0.78±0.03, a radius of 1.168 ± 0.006 R_Jup and an effective temperature of 1155±5 K. The C/O ratio is always found to be consistent with the stellar value, with retrievals finding C/O = ${0.60}_{-0.06}^{+0.04}$. 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_cloud = 0 45 ± 0 3. The most favoured solutions require eitheramorphous 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σ. 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 ∩ 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}_{-20}^{+30}$ 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_pl = ${7.5}_{-0.6}^{+0.6}$ and R_pl = ${1.13}_{-0.05}^{+0.05}$ respectively. This leads to a log g of ${4.18}_{-0.05}^{+0.06}$, slightly higher than the self-consistent estimates of 3.5–4.0. The overall best self-consistent model by the χ² was Exo-Rem, which finds an effective temperature of 1100K, a radius of 1.15 R_Jup, and a somewhat low log g of 3.5. ATMO is the most favoured self-consistent model when using the Bayesian framework, though it found an unphysically small radius and higher temperature than other models. HR 8799 e is the only companion for which the MIRI photometry is not convincingly fit, as seen below in 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 [M/H]=${1.9}_{-0.2}^{+0.1}$. Metallicities >1 were a universal feature of the retrievals for e. This was reinforced by the grid-fits, which uniformly find strong enrichment, running into the upper grid boundaries. Free chemistry retrievals are always preferred over the disequilibrium retrievals; from these we found the group A ∩ B C/O ratio is ${0.87}_{-0.02}^{+0.02}$. As the free chemistry retrieval C/O ratio only accounts for the gas-phase abundances, there is additional oxygen 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±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 $\mathrm{C}/\mathrm{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 ± 0.11. This detection is strongly favoured by the Bayes factor, with log₁₀ $\mathcal{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₄ 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_sed parameters for the silicate and iron clouds were not required, with variations in the Bayes factor driven more by the choice of patchiness and temperature profile.

5.5 Impacts of modelling choices

While we examined the primary atmospheric characteristics of each of the four HR 8799 planets, it is crucial to understand how the choice of model impacts these measurements. As described in Sec. 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 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.

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.

5.5.2 Chemistry

For all four planets, we performed retrievals using a grid derived from an equilibrium chemistry solver with disequilibrium H₂O-CO-CH₄ 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 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 $$\mathrm{C}/\mathrm{H}=\frac{X_{\mathrm{C}\mathrm{O}}+X_{{\mathrm{C}\mathrm{O}}_2}+X_{{\mathrm{C}\mathrm{H}}_4}+X_{\mathrm{H}\mathrm{C}\mathrm{N}}}{X_{{\mathrm{H}}_2}+2X_{{\mathrm{H}}_2\mathrm{O}}+2X_{{\mathrm{H}}_2\mathrm{S}}+X_{\mathrm{H}\mathrm{C}\mathrm{N}}+4X_{{\mathrm{C}\mathrm{H}}_4}+3X_{{\mathrm{N}\mathrm{H}}_3}},$$(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 nondetection 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}_{-2.3}^{+1.7}\mathrm{\%}$. They also identify a relation between the bulk and observed C/O ratio: $$(\mathrm{C}/\mathrm{O}{)}_{\text{obs }}\approx \frac{(\mathrm{C}/\mathrm{O}{)}_{\text{bulk }}}{1-0.371(\mathrm{C}/\mathrm{O}{)}_{\text{bulk }}}\text{. }$$(17)

Solving for (C/O)_bulk, we find that for HR 8799 e, with an observed C/O ratio of ${0.88}_{-0.02}^{+0.02}$, should have a bulk C/O ratio of 0.66, much closer to the stellar value of 0.54. Likewise, HR 8799 d has the lowest observed C/O ratio of ${0.61}_{-0.04}^{+0.03}$, 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_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₂O, and CH₄ 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_pl/Z_*, which allows us to directly compare our metallicities to literature values, such as those of Thorngren et al. (2016). We converted the metallicity [M/H] of each atmosphere to Z_pl using the methods of Thorngren & Fortney (2019), adapting for our own notation: $${10}^{[\mathrm{M}/\mathrm{H}]}=\frac{1+Y/X}{(Z_{\mathrm{p}\mathrm{l}}^{-1}-1)\left(\frac{{\mu }_Z}{{\mu }_H}\right)},$$(18)

where X, Y, and Z are the solar hydrogen, helium and metal mass fractions and μ is the mean molecular weight of the metal content of the atmosphere. Rearranging and substituting in the measured atmospheric metallicity [M/H], we find: $$Z_{\mathrm{p}\mathrm{l}}={(1+\frac{1+Y/X}{{10}^{[\mathrm{M}/\mathrm{H}]}\left(\frac{{\mu }_{\mathrm{Z}}}{{\mu }_H}\right)\left(Z/H_{\odot }\right)})}^{-1}.$$(19)

We take the same assumptions as Thorngren & Fortney (2019), taking μ_Z to be 18, assuming most of the metal content is in water, μ_H to be 1 for atomic hydrogen, and Y/X to be 0.3383 as in Asplund et al. (2009). As the metallicity of HR 8799 A is near solar, Z_*/H is taken to be the solar value of Z/H_⊙ = 1.04 × 10⁻³. To normalise to the stellar metallicity we follow Thorngren et al. (2016) and calculate the Z_* as $$Z_{\star }=0.014\times {10}^{\mathrm{F}\mathrm{e}/\mathrm{H}}\text{. }$$(20)

For HR 8799, we used solar metallicity to calculate Fe/H, but refer to the discussion in Sec. 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₂O, CO, and CH₄ become vertically constant. This parameterises dynamical mixing that homogenises the upper layers of the atmosphere. More rigorously, the quench point is defined as the point at which the chemical timescale t_chem and the mixing t_mix are equal. Following Zahnle & Marley (2014), the mixing timescale is defined through the ratio of the local atmospheric scale height H to the vertical eddy diffusion coefficient, K_zz: $$t_{\text{mix }}=\frac{H^2}{K_{zz}}.$$(21)

The chemical timescale depends on the reaction rates involved. Considering the CO-CH₄ 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^[M/H], and temperature T in kelvin: $$t_{\mathrm{q}1}=1.5\times {10}^{-6}{p}^{-1}{\mathrm{m}}^{-0.7}\exp (42000/T)\mathrm{s}.$$(22)

For weak mixing, and therefore drawing from low temperatures with little CO, the timescale is found to be $$t_{\mathrm{q}2}=40p^{-2}\exp (25000/T)\mathrm{s},$$(23)

Combining the two, the total chemical timescale is defined as: $$t_{\mathrm{C}\mathrm{O}}={(\frac{1}{t_{\mathrm{q}1}}+\frac{1}{t_{\mathrm{q}2}})}^{-1}$$(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_{\mathrm{z}\mathrm{z}}=\frac{H^2}{t_{\mathrm{C}\mathrm{O}}}$$(25)

where the scale height is defined as $$H=\frac{k_{\mathrm{B}}T}{\mu g},$$(26)

for temperature T, surface gravity g, mean molecular mass μ, and the Boltzmann constant k_B. In order to calculate these quantities for the HR 8799 planets, we take T_eff as a representative temperature to calculate the timescales and scale height, g as the measured surface gravity, and μ as the average mean molecular weight of the atmosphere. As K_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_zz ≈ 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 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 ∩ 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 ∩ B, we derived K_zz from the quench pressure. We then turn to Soni & Acharyya (2023) for a comparison, who provide predictions for CH₄ and CO abundances for varying K_zz, T_eff, and log g across a range of metallicities. For our measured T_eff and CH₄ 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₄ 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.7}^{+0.6}$ 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 ~ 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.

Fig. 13 Metallicities and elemental number ratios for the HR 8799 planets derived from Baysian averaged group A ∩ B free chemistry retrievals.

Fig. 14 Cloud properties for HR 8799 e, showing the optical depth due to clouds as a function of pressure (colour map), with the τ = 0.3 and τ = 1.0 contours highlighted by the dashed lines. The solid lines indicate the mass fraction abundance of the MgSiO3 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 SFe = 0.0 ± 1.1, log $S_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}$ = −0.8 ± 1.0).

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, Δ log₁₀ $\mathcal{Z}$ = 0.7), though this was only compared for HR 8799 e using the Z23 profile, and assuming the clouds condense at their equilibrium saturation location. In this case a_h is calculated from f_sed and K_zz, and a lower f_sed for the MgSiO₃ cloud was retrieved than with the log normal distribution (f_{sed, Hansen} = 1.18 ± 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₃, SiO₂, and Fe clouds provides the best-fit model to an ultracool field dwarf, 2MASSW J2224438-015852. Figure 14 sh₃ows that the MgSiO₃ clouds in their model are located at 10⁻³ bar with a maximum optical depth of τ = 0.3 at 1 μm and an effective particle radius of ~0.04 μm. The SiO₂ clouds are slightly deeper, at 10⁻² bar, and are optically thick at 1 μm. This is the same location where we find MgSiO₃ 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 AM01 model is a function of many atmospheric factors, including the temperature, mixing strength, f_sed, and particle number density. We see that changes in the particle radius are correlated with 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, 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₃ clouds in certain like-for-like retrieval comparisons. In reverse, the lack of features leads to a preference for Na₂S clouds in HR 8799 b. Performing additional tests on HR 8799 e, we find that Mg₂SiO₄ (e.AB.14) are mildly disfavoured by the retrievals, and Al₂O₃ (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.

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.

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

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_* 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_*. Where stellar metallicities are not available, we assumed a metallicity [Fe/H]=0. Only the 2 M_Jupplanet Af Lep b has a comparable degree of enrichment to the HR 8799 planets (Zhang et al. 2023). However, they are comparable to hot jupiters observed in transmission, such as WASP 39 b (Rustamkulov et al. 2023). Looper et al. (2008) and Stephens et al. (2009) demonstrate how high metallicity atmospheres facilitate condensation, in turn leading to the strong reddening seen in a subset of L dwarfs. The retrieval of highly metal rich and very cloudy atmospheres is consistent with this picture from brown dwarfs. Further supporting the high [M/H] retrievals are their consistency with the self-consistent grids, which always favour their upper limits.

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 Δ log $\mathcal{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₂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 dependant on both the height of the J, H, and K band peaks, as well as the shape of these features. However, the J bands is only covered by the low-resolution SPHERE data, with relatively poor S/N, and 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₂ abundance scales proportionally to [M/H]². 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 χ². This exercise was repeated, setting the metallicity of the model to solar. For HR 8799 c and b, we find that the high metallicity model provides a better fit to the OSIRIS data than the solar model, with the caveat that fixing the metallicity during the retrieval may result in a better fit than setting it a posteriori, without changing other atmospheric parameters. For HR 8799 d, the high metallicity model is only a marginally better fit than the solar metallicity model. In general, the fits from the GRAVITY and remaining archival data provide reasonable fits to the R ≈ 4000 OSIRIS data, of similar quality to the fits displayed in Ruffio et al. (2021), and are included in the appendix in Fig. C.1.

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: β 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); κ 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).

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 GRAVITY data to determine the constraining power of this new dataset. We found that using only the GRAVITY data we could rule out clear atmosphere solutions at Δ log $\mathcal{Z}$ > 7 (e.g. e.A.31 over e.A.32), using the Z23 profile and disequilibrium chemistry. Using only the GRAVITY observations for HR 8799 e, we could obtain estimates of the metallicity (2 ± 0.3), 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.

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.

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

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.

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.

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 superstellar for all four planets. It decreases from the outermost planet to HR 8799 d, while HR 8799 e has a higher C/O ratio than the other companions. We confidently detect HCN in HR 8799 c and e, at abundances far higher than predicted by equilibrium chemistry; though this detection is largely driven by 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_zz values compatible with high-metallicity predictions from Soni & Acharyya (2023). The mixing strength is stronger for the warmer planets, at log K_zz ≈ 6, and is lower for HR 8799 b with log K_zz ≈ 2.

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 Na2S 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. While 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) amd from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime”. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. This work is based on observations collected at the European Southern Observatory under ESO programme 1104.C-0651. This publication makes use of VOSA, developed under the Spanish Virtual Observatory project supported by the Spanish MINECO through grant AyA2017-84089. VOSA has been partially updated by using funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement n° 776403 (EXOPLANETS-A).

Appendix A Data logs

Appendix B Reprocessing SPHERE and GPI datasets

To resolve the known discrepancies in the H-band flux between the archival SPHERE and GPI datasets, we reprocessed each using KLIP (Soummer et al. 2012; Pueyo 2016), ANDROMEDA (Cantalloube et al. 2015) and PynPoint (Amara & Quanz 2012; Stolker et al. 2019). We optimised the choice of algorithm parameters through a series of injection/extraction tests into each dataset, as in Nasedkin et al. (2023). Using two different goodness-of-fit metrics on injections representative of the true companion contrast and separation, we choose the number of principal components used in the PSF subtraction in order to extract the companion spectra with minimal bias. Figures B.1 and B.2 shows the extracted spectra for each of HR 8799 c, d, and e, for SPHERE and GPI respectively, compared to published literature spectra from Zurlo et al. (2016), Greenbaum et al. (2018), and Flasseur et al. (2018). As the goodness-of-fit metrics favoured the KLIP extractions, we used these as the basis of our retrieval analysis.

Fig. B.1 KLIP and ANDROMEDA extractions from SPHERE for HR 8799 c, d, and e compared to the spectra published in Zurlo et al. (2016) and Flasseur et al. (2018).

Fig. B.2 KLIP and ANDROMEDA extractions from GPI for HR 8799 c, d, and e compared to the spectra published in Greenbaum et al. (2018)

Appendix C Retrieval validation

Extensive validation of the pRT retrieval module was performed as part of this work. Following updates described in Nasedkin et al. (2024), we verified that the results of Mollière et al. (2020) could be reproduced. We independently tested updates to the c-k mixing implementation, the adaptive mesh refinement implementation, updated opacity sources for H₂O and CO, bug fixes for convergence on multiple scattering in the clouds, the inclusion of photometric data, the inclusion of scaling factors on the SPHERE and GPI datasets, including or excluding the GPI K-band spectra, updates to each of the SPHERE, GPI, and GRAVITY datasets, different prior widths and the number of live points used in the retrieval. All of the posterior distributions were fully consistent to within 2σ with most falling well within 1σ of the published results, apart from the inclusion of new epochs of GRAVITY data, which led to a significantly higher retrieved metallicity ([M/H]=1.1 ± 0.32) and an f_sed of 5 ± 2.6.

We verified several model assumptions through retrievals that only include the GRAVITY datasets, or the GRAVITY data and photometry. We find that the GRAVITY data alone could not distinguish between clear and cloudy models (Δ log₁₀ $\mathcal{Z}$ < 1), while cloudy models were strongly favoured once the broad wavelength coverage of the photometry was included (Δ log₁₀ $\mathcal{Z}$ > 10).

Models that use the dynamical mass estimates as priors for calculating the surface gravity were marginally favoured over those that freely retrieve log g and R_pl, but that the posteriors parameter distributions were generally consistent, with T_int and R_pl showing the greatest discrepancy. Using the dynamical mass as a prior and setting a Gaussian prior on the radius led to more reasonable estimates of the radius of HR 8799 e (0.97 ± 0.04) compared to the free retrieval (0.79 ± 0.05). However, the composition of the planet and the degree of cloudiness did not vary significantly between the two models.

Using the full dataset for HR 8799 b, we verified that retrievals including scattering clouds are strongly favoured over those without scattering (Δ log₁₀ $\mathcal{Z}$ > 10). Without scattering, both the temperature and composition ([M/H] and C/O) are significantly discrepant from retrievals that include scattering clouds.

Fig. C.1 Comparison of best-fit disequilibrium models to OSIRIS data from Ruffio et al. (2021). From top to bottom is HR 8799 b, c, and d. In blue are the best-fit disequilibrium models, with the spectra generated using high-resolution line-by-line opacities, before being convolved, binned, and normalised for comparison. In orange is the same, but with the metallicity set to 0.

Fig. C.2 Comparison of best-fit disequilibrium models (black) of HR 8799 c to the data, with residuals shown in the bottom panel. In blue are the same spectra, but without opacity contributions from CH4.

Fig. C.3 Comparison of best-fit disequilibrium models (black) of HR 8799 c to the data, with residuals shown in the bottom panel. In blue are the same spectra, but without opacity contributions from HCN.

Appendix D Using the Hansen distribution with EDDYSED

The EDDYSEDcloud model from Ackerman & Marley (2001) is implemented in pRT, and is the most physically motivated model incorporated to date. Typically, it assumes a log-normal particle size distribution, where the geometric particle radius will vary throughout the atmosphere as a function of the vertical diffusion coefficient K_zz and the sedimentation fraction f_sed. Here, we substitute the log-normal particle size distribution with the Hansen distribution, originally introduced in Hansen (1971), and rederive the calculation for the particle radius as a function of K_zz and f_sed.

We begin with a review of the EDDYSEDmodel: the distribution of the number of particles as a function of particle radius, n(r) is approximated as a log-normal distribution with width σ_g and characteristic geometric radius r_g. $$n(r)=\frac{N}{r\sqrt{2\pi }\text{\hspace{0.17em}log\hspace{0.17em}}{\sigma }_g}\exp (-\frac{{\log }^2\left(r/r_g\right)}{2{\text{\hspace{0.17em}log\hspace{0.17em}}}^2{\sigma }_g}),$$(D.1)

N is the total number of cloud particles.

The goal of the EDDYSEDmodel is to calculate r_g for each layer in the atmosphere, given K_zz and f_sed. It balances the upwards vertical mixing, parameterised by K_zz and the particle settling velocity, υ_f $$v_f=w_{\ast }{\left(\frac{r}{r_w}\right)}^{\alpha }.$$(D.2)

Here w_* is the convective velocity scale. We note that r_w ≠ r_g. r_w is the radius at which the particle settling velocity equals the convective velocity scale: $$w_{\ast }=\frac{K_{\mathrm{z}\mathrm{z}}}{L}\text{, }$$(D.3)

where L is the convective mixing length. Since w_* is known, and υ_f can be found analytically as in Ackerman & Marley (2001); Podolak (2003), a linear fit can be used to find both α and r_w.

With both of these quantities known, we follow AM01 and define f_sed as: $$f_{\text{sed }}=\frac{{\int }_0^{\mathrm{\infty }}{r}^{3+\alpha }n(r)dr}{r_w^{\alpha }{\int }_0^{\mathrm{\infty }}{r}^{3}n(r)dr}$$(D.4)

For the log-normal distribution, one finds: $${\int }_0^{\mathrm{\infty }}{r}^{\beta }n(r)dr=N{r}_g^{\beta }\exp (\frac{1}{2}{\beta }^2{\log }^2{\sigma }_g)$$(D.5)

Which we can then use to solve for r_g: $$r_g=r_w{f}_{\mathrm{s}\mathrm{e}\mathrm{d}}^{1/\alpha }\exp (-\frac{\alpha +6}{2}{\log }^2{\sigma }_g)$$(D.6)

In order to use the Hansen distribution, we must recalculate the total number of particles N, and integrate the distribution for f_sed. We note here that the Hansen distribution is parameterised by the effective radius, $\overline{r}$, rather than the geometric mean radius. In this derivation we do not correct for this difference in definition, as both act as nuisance parameters in the context of an atmospheric retrieval.

We start by giving the Hansen distribution in full: $$n_{\text{Hansen }}(r)=\frac{N{\left(\overline{r}{v}_e\right)}^{(2v_e-1)/v_e}}{\mathrm{\Gamma }\left[(1-2v_e)/v_e\right]}{r}^{(1-3v_e)/v_e}\exp (-\frac{r}{\overline{r}{v}_e})$$(D.7)

In Hansen (1971), the authors use the parameters a and b to denote the mean effective radius and effective variance, which we write as $\overline{r}$ and υ_e respectively. These differ from the simple mean radius and variance by weighting them by the particle area, as the cloud particle scatters an amount of light proportional to its area. Thus: $$\overline{r}=\frac{{\int }_0^{\mathrm{\infty }}r\pi r^{2}n(r)dr}{{\int }_0^{\mathrm{\infty }}\pi r^{2}n(r)dr}$$(D.8)

and $$v_e=\frac{{\int }_0^{\mathrm{\infty }}(r-\overline{r}{)}^2{r}^{2}n(r)dr}{{\overline{r}}^2{\int }_0^{\mathrm{\infty }}\pi r^{2}n(r)dr}$$(D.9)

As in EDDYSED, we fit for the settling velocity, which will provide us with α and r_w, which we can use to find f_sed, as in D.4. However, we must now integrate the Hansen distribution. We find that: $${\int }_0^{\mathrm{\infty }}{r}^{\beta }{n}_{\text{Hansen }}(r)dr=\frac{v_e^{\beta }(v_e\beta +2v_e+1){\left(\frac{1}{\overline{r}}\right)}^{-\beta }\mathrm{\Gamma }(\beta +1+\frac{1}{v_e})}{(-v_e+v_e^{\beta +3}+1)\mathrm{\Gamma }(1+\frac{1}{v_e})}$$(D.10)

We can then use Eqns. D.4 and D.10 to solve for $\overline{r}$: $$\overline{r}={\left(\frac{f_{\mathrm{sed}}{r}_w^{\alpha }{v}_e^{-\alpha }(v_e^{3+\alpha }-v_e+1)\mathrm{\Gamma }(1+\frac{1}{v_e})}{(v_e\alpha +2v_e+1)\mathrm{\Gamma }(\alpha +1+\frac{1}{v_e})}\right)}^{\frac{1}{\alpha }}.$$(D.11)

Thus for a given K_zz, f_sed, and υ_e, we can find the effective particle radius for every layer in the atmosphere.

However, in order to compute the cloud opacity, we still require the total particle count. For a volume mixing ratio of a given species, χ_i, we can integrate n(r) to find N: $$N=\frac{{\chi }_i}{({\overline{r}}^3{v}_e-1)(2v_e-1)}$$(D.12)

Appendix E Complete retrieval results

We include in the text abridged tables that present key parameters of interest. The complete set of inferred parameters for every retrieval is available on Zenodo.

Legend: Chemistry/Profile/Clouds/Data/Info.

• Chemistry: (D)isequilibrium or (F)ree.

• Profile: (M)olliere, (Z)hang, (G)uillot or (S)pline(NNodes).

• Clouds: Clear (CLR); (f)ree or (eq)uilibrium condensation location, (species)_(cd/am)_(P)atchy_(h)ansen. ‘*’ indicates f_sed was retrieved independently for each cloud species.

• Data: -(not included) or (only included). O indicates OSIRIS data was used in place of GRAVITY data.

• Info: ‘-’ indicates not included. ‘mr’ indicates mass and radius were used as parameters instead of log g and radius.

In the text, models will be referred to as planet.group.index.

References

All Tables

All Figures

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

	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 × 10−15 W m−2 (μ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.
In the text

	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).
In the text

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

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

In the text

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

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

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

	Fig. 9 Teff 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.
In the text

	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 χ2 model from the grid.
In the text

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.

In the text

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

	Fig. 13 Metallicities and elemental number ratios for the HR 8799 planets derived from Baysian averaged group A ∩ B free chemistry retrievals.
In the text

Fig. 14 Cloud properties for HR 8799 e, showing the optical depth due to clouds as a function of pressure (colour map), with the τ = 0.3 and τ = 1.0 contours highlighted by the dashed lines. The solid lines indicate the mass fraction abundance of the MgSiO3 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 SFe = 0.0 ± 1.1, log $S_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}$ = −0.8 ± 1.0).

In the text

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.

In the text

	Fig. 16 Cloud absorption (left) and scattering (right) opacities for different condensate compositions and structure, for 1-μm particles. Dark blue indicates MgSiO3, light blue is Mg2SiO4. Solid (dashed) lines are for crystalline (amorphous) substances. The solid red (green) line is for crystalline iron (Na2S).
In the text

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: β 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); κ 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).

In the text

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

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

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

	Fig. B.1 KLIP and ANDROMEDA extractions from SPHERE for HR 8799 c, d, and e compared to the spectra published in Zurlo et al. (2016) and Flasseur et al. (2018).
In the text

	Fig. B.2 KLIP and ANDROMEDA extractions from GPI for HR 8799 c, d, and e compared to the spectra published in Greenbaum et al. (2018)
In the text

	Fig. C.1 Comparison of best-fit disequilibrium models to OSIRIS data from Ruffio et al. (2021). From top to bottom is HR 8799 b, c, and d. In blue are the best-fit disequilibrium models, with the spectra generated using high-resolution line-by-line opacities, before being convolved, binned, and normalised for comparison. In orange is the same, but with the metallicity set to 0.
In the text

	Fig. C.2 Comparison of best-fit disequilibrium models (black) of HR 8799 c to the data, with residuals shown in the bottom panel. In blue are the same spectra, but without opacity contributions from CH4.
In the text

	Fig. C.3 Comparison of best-fit disequilibrium models (black) of HR 8799 c to the data, with residuals shown in the bottom panel. In blue are the same spectra, but without opacity contributions from HCN.
In the text