Retrieving scattering clouds and disequilibrium chemistry in the atmosphere of HR 8799e

Clouds are ubiquitous in exoplanet atmospheres and represent a challenge for the model interpretation of their spectra. Complex cloud models are too numerically costly for generating a large number of spectra, while more efficient models may be too strongly simplified. We aim to constrain the atmospheric properties of the directly imaged planet HR 8799e with a free retrieval approach. We use our radiative transfer code petitRADTRANS for generating spectra, which we couple to the PyMultiNest tool. We added the effect of multiple scattering which is important for treating clouds. Two cloud model parameterizations are tested: the first incorporates the mixing and settling of condensates, the second simply parameterizes the functional form of the opacity. In mock retrievals, using an inadequate cloud model may result in atmospheres that are more isothermal and less cloudy than the input. Applying our framework on observations of HR 8799e made with the GPI, SPHERE and GRAVITY, we find a cloudy atmosphere governed by disequilibrium chemistry, confirming previous analyses. We retrieve that ${\rm C/O}=0.60_{-0.08}^{+0.07}$. Other models have not yet produced a well constrained C/O value for this planet. The retrieved C/O values of both cloud models are consistent, while leading to different atmospheric structures: cloudy, or more isothermal and less cloudy. Fitting the observations with the self-consistent Exo-REM model leads to comparable results, while not constraining C/O. With data from the most sensitive instruments, retrieval analyses of directly imaged planets are possible. The inferred C/O ratio of HR 8799e is independent of the cloud model and thus appears to be a robust. This C/O is consistent with stellar, which could indicate that the HR 8799e formed outside the CO$_2$ or CO iceline. As it is the innermost planet of the system, this constraint could apply to all HR 8799 planets.


Introduction
The description of clouds in exoplanets and brown dwarfs is one of the major uncertainties when modeling the structures and spectra of self-luminous atmospheres (e.g., Marley et al. 2013).Fully modeling the microphysics of clouds is difficult, due to modeling uncertainties and long computational timescales, as the chemistry, nucleation process, condensation, particle coalescence, settling and mixing need to be accurately described.Moreover, even if all these processes are taken into account, it is not straightforward which values to prescribe for the remaining free parameters.Some of the "free parameters" of such elaborate cloud models are likely not free at all, but are determined by the full solution of the (multi-dimensional) atmospheric structure, which is a function of the cloud properties itself, due to the radiative feedback of the clouds.At the same time such complicated cloud models are useful and necessary, because they allow us to understand the interplay of physical processes during cloud formation, and will ultimately have to explain the cloud properties of all exoplanets and brown dwarfs, irradiated or self-luminous.Two examples for such complete cloud models are described in Woitke & Helling (2004); Helling et al. (2008b); Woitke et al. (2020) and Gao et al. (2018); Powell et al. (2018).When comparing synthetic cloudy spectra to observations, often one-dimensional self-consistent models are used, where the cloud opacity is radiatively coupled to the atmospheric temperature structure.Here iterating the structure is necessary.This requires a faster cloud modeling approach, which parametrizes parts of the cloud formation process.Examples are models based on timescale comparisons as in Allard et al. (2001Allard et al. ( , 2003)), implementing the approach of Rossow (1978), or Ackerman & Marley (2001), which uses the ratio of the cloud particle settling and mixing velocities ( f sed ) as a free parameter, and solves for the particle size assuming a log-normal particle size distribution and a vertical diffusion coefficient K zz .The model of Charnay et al. (2018) is again different, and mixes the two previous approaches: as in Ackerman & Marley (2001), the vertical distribution of the cloud mass is determined assuming a steady state between mixing, settling and cloud condensation in every layer, while the average particle size is found using the timescale approach of Rossow (1978).Thus no f sed needs to be specified in Charnay et al. (2018) and this model can reproduce the L-T transition, including the effect of low gravity, which moves the transition to lower effective temperatures.Another interesting approach is the recent model by Ormel & Min (2019), which determines the cloud mass fraction and average particle size by solving the steady-state differential equations including cloud settling, mixing, nucleation, condensation and coagulation as a function of K zz and the nucleation rate.A summary of other models can be found in Helling et al. (2008a); Helling & Casewell (2014).
The two modeling philosophies described above (full microphysical or simplified for increased speed) are invaluable for understanding both cloud (micro)physics and the self-consistent radiative feedback of clouds.However, they are challenging when one aims at fitting cloudy spectra.The former, more complete models may take prohibitively long when calculating a large number of atmospheric structures and spectra.The latter, more parametrized models, allow for the calculation of larger self-consistent model grids.However, an important question is whether the simplification steps during model construction were all justified, and if all explicit (and implicit) free parameters of the model have been varied sufficiently.Moreover, an update of the model requires the calculation of a new grid, the models of which are demanding to produce and converge especially with clouds (see, e.g., Morley et al. 2014;Mollière et al. 2017), whereas considering additional parameters increases the grid size by orders of magnitude.
In the work presented here, we use a different approach.Namely we attempt to retrieve the characteristics of cloudy, selfluminous atmospheres by means of free retrievals.This is done by parameterizing the temperature profile, as well as the cloud properties, while using chemical equilibrium abundances with a simple quench pressure treatment to account for atmospheric mixing.While clouds in various parameterizations have been included in retrievals of transmission spectra of exoplanets (see MacDonald & Madhusudhan 2017;Fisher & Heng 2018;Tsiaras et al. 2018;Pinhas et al. 2019;Barstow 2020, for some recent examples), cloudy retrievals are still a comparatively novel approach for self-luminous targets.The use of a free retrieval approach for fitting the spectra of brown dwarfs and directly imaged planets is motivated by Line et al. (2015Line et al. ( , 2017)); Zalesky et al. (2019), who retrieved the atmospheric properties of clear T-and Y-dwarfs, Burningham et al. (2017), who studied cloudy L-dwarfs, and Lee et al. (2013); Lavie et al. (2017) who attempted to retrieve the properties of the cloudy HR 8799 planets for the first time.These pioneering works show the power of free retrievals for constraining condensation physics and clouds in cloud-free and cloudy brown dwarfs, and how retrieved planetary abundances may be connected to planet formation.
The radiative transfer tool used in the retrievals here is pe-titRADTRANS (Mollière et al. 2019), which we update to include the effect of scattering, which can no longer be ignored for cloudy atmospheres.By parametrizing the clouds, rather than making assumptions on how to simplify the cloud modelling process, we let the data constrain basic cloud characteristics such as cloud mass, location, and particle size, provided that the signal-to-noise of the data is high enough.Our model has the advantage that any changes in the (cloud) parametrization approach can be quickly implemented and tested, without the need of recalculating cloudy model grids.
The capabilities of our retrieval model are demonstrated by analyzing new and archival spectra of the cloudy planet HR 8799e, taken with the GRAVITY (K band, see Gravity Collaboration et al. 2019), SPHERE (YJH bands, see Zurlo et al. 2016), and GPI (H band, see Greenbaum et al. 2018) instruments.While not included in the fit, we also compare to archival mid-infrared (MIR) photometry. 1 The HR 8799 system is especially interesting because it hosts four directly imaged planets that orbit their star within a massive debris disk (Marois et al. 2008(Marois et al. , 2010;;Currie et al. 2011;Su et al. 2009) at distances from 15 to 70 au (e.g., Wang et al. 2018).This allows for the comparative characterization of the planets' spectral properties.In particular the planets' atmospheric abundances may shed light on how they formed from the circumstellar disk.Consequently, the HR 8799 planets have been extensively studied in the literature, and have been classified to bear the hallmarks of thick clouds and disequilibrium chemistry, placing them in the low-gravity, cool end of L spectral sequence (see Section 5.2 for a more detailed discussion of the literature).In addition to confirming these findings for HR 8799e, we also derive the planet's metallicity and, for the first time, carbon-to-oxygen number ratio (C/O), to study possible formation pathways.
Our retrieval model is detailed in Section 2, including our description of the scattering implementation, the temperature and chemistry parameterization, as well as our cloud model parameterizations.Section 3 describes tests to verify our model setup with mock retrievals.Section 4 contains our retrieval study of the planet HR8799e.Section 5 discusses the implications of our results for the formation of HR 8699e, and compares to the extensive body of literature on the HR 8799 planets.We summarize our findings in Section 6 and provide an outlook for the further development and application of our new retrieval model.

Adding scattering to petitRADTRANS
The objects we seek to model are expected to be inherently cloudy.Hence, scattering is an important process that needs to be considered during the radiative transfer calculations.To calculate quantities such as the photon destruction probability, it is necessary to compare the scattering and absorption opacities.For this, the total atmospheric opacity, combined from the individual absorber opacitites, needs to be calculated.Thus the correlated-k treatment of petitRADTRANS needs to be adapted to combine the k-tables (opacity tables) of individual atmospheric absorbers.2This is in contrast to the case treating purely emission, which can be handled by using the products of the transmissions of individual species (see, e.g., Irwin et al. 2008;Mollière et al. 2019).Because our goal is to run retrievals, the k-table combination has to be as fast and as accurate as possible.
Here we present a newly developed method to quickly combine k-tables of different absorber species, which works by sampling the opacity distribution functions of individual absorbers.Computational time is saved by sampling the indices of the ktable entries, instead of interpolating the k-tables to sampled values of the cumulative probability.The k-table mixing process is described and tested in Appendix A.
For solving the radiative transfer equation we then use the same treatment as described in our self-consistent petitCODE (Mollière et al. 2015(Mollière et al. , 2017)), namely by using the Feautrier method (Feautrier 1964), converging the scattering source function with local Accelerated Lambda Iteration (ALI) (Olson et al. 1986) and Ng acceleration (Ng 1974).The scattering process is assumed to be isotropic, with a (1 − g a ) correction factor applied to the scattering cross-sections, where g a is the scattering anisotropy.The scattering implementation is further described in Appendix A.6 of Mollière et al. (2017).We show a verification of the cloudy spectra of petitRADTRANS, including scattering, in Appendix B. For this verification, we used petitCODE to calculate a self-consistent HR 8799e model in radiative-convective and chemical equilibrium, which included clouds of MgSiO 3 and Fe.The spectra of petitRADTRANS and petitCODE agree excellently.It takes petitRADTRANS a few seconds to calculate a cloudy emission spectrum in the YJHK bands, which is fast enough for retrievals on computational clusters.

Temperature model
Our goal is to parameterize the vertical temperature profile of the atmosphere in a way that imposes as few prior constraints on the solution as possible.It would hence appear to be ideal to retrieve the temperatures in every layer of the discretized atmosphere independently, which is an approach commonly followed in the planetary science community (e.g., Rodgers 2000;Irwin et al. 2008), and which has also been applied to cloud-free brown dwarfs (Line et al. 2014a) and exoplanets (e.g., Lee et al. 2012).
However, if the data are sparse or of low signal-to-noise, a level-by-level retrieval of the temperature can lead to overfitting and thus unphysical oscillations in the inferred temperature profile.One way of reducing such oscillations is by smoothing the resulting P-T profile (as done in Irwin et al. 2008;Lee et al. 2012).Another way to circumvent this problem is to retrieve temperatures at a limited number of altitudes in the atmosphere, which are then connected via (spline) interpolation to yield a temperature at all layers, thereby reducing the number of free parameters (Line et al. 2015;Kitzmann et al. 2019).In addition, Line et al. (2015) included a penalty term on the spatial sum of second derivatives of the temperature profile which further discouraged oscillatory solutions.Because the weight of said penalty can bias the results, its optimal value is also fitted in the inference process.This method has been used to retrieve the temperature profiles of cloud-free T-and Y-dwarfs in Line et al. (2015Line et al. ( , 2017)); Zalesky et al. (2019).
When investigating different models suitable for retrieving the atmospheres of cloudy self-luminous exoplanets, we settled on a model that uses both freely variable and physically motivated parameterizations, based on the atmospheric altitude.This temperature model, which allows us to retrieve the synthetic structures of cloudy atmospheres, is split into three parts, going from high, to middle, to low altitude.The spatial coordinate of the temperature model is an optical depth τ4 , which we relate to the pressure P by τ = δP α , where δ and α are free parameters, and P is the atmospheric pressure in units of dyn cm −2 .This mapping is required because P is the vertical coordinate of petitRADTRANS.We then setup the atmospheric temperature profile, starting with the middle altitudes, that is, the 'photosphere'.
'Photosphere' (middle altitudes) This region stretches from τ = 0.1 to the radiative-convective boundary.Here we set the temperature according to the Eddington approximation where T 0 is a free parameter.The optical depth τ is obtained from Equation 1 above.In the original Eddington solution, from which we take the functional form of the temperature profile, this corresponds to the internal temperature.Likewise, we note that the 'photospheric' region does not necessarily have to correspond to the true photosphere of the planet.Because the Eddington solution will always lead to an isothermal upper atmosphere, which is not expected to occur in reality, the highaltitude region of the atmosphere is treated separately, described immediately below.

High altitude
This region extends from the top of the atmosphere (P = 10 −6 bar) to τ = 0.1.Here we split the atmosphere into four equidistant locations in log(P) space and treat the temperature at the three upper locations as free parameters.The temperature at the lowest altitude, which is at τ = 0.1, is taken from the Eddington approximation of the 'photosphere'.The temperatures in this atmospheric region are then found from a cubic spline interpolation.

Troposphere (low altitudes)
This region starts from the radiative-convective boundary and extends to the bottom of the atmosphere.The radiative convective boundary is found by comparing the atmospheric temperature gradient of the Eddington approximation with the moist adiabatic temperature gradient of the atmosphere.This is done by interpolating the moist adiabatic gradient in the T -P-[Fe/H]-C/O space of the chemistry table (see Section 2.3).As soon as the atmosphere is found to be Schwarzschild-unstable, the atmosphere is forced onto the moist adiabat.

Priors
We restrict α (see Equation 1) to vary between 1 and 2, following Robinson & Catling (2012).Moreover, to prevent the formation of temperature inversions, which are not expected in selfluminous objects, we required that the three free temperature points in the high altitude region of the atmosphere are colder than the highest point of the 'photosphere' (middle altitude region), and that they decrease in temperature monotonically with increasing altitude.This prior is enforced by setting the upper boundary of the allowed temperature range of such a free temperature point equal to the temperature of the underlying temperature point, identical to the treatment in Kitzmann et al. (2019).
In Gravity Collaboration et al. (2020) we had found that we had to impose further priors on the temperature parameterization, based on the structure of the atmospheric opacity of a given forward modeling realization.For example, the power law index of the optical depth, α, was not allowed to deviate too far from the power law index measured from the opacity structure in the forward model, within the spectral range of the retrieved data.We found that this was necessary because of the high dimensionality of our retrieval model, and our inability to make the MCMC sampler find the global maximum of the log-probability in a finite amount of time otherwise.These opacity priors restricted the parameter space for the MCMC sufficiently.In this publication we use nested sampling (Skilling 2004), and using a sufficiently large number of live points made such opacity priors unnecessary.

Chemistry model
In the retrievals presented in this work, the chemical abundances within the atmosphere are determined by means of interpolation in a chemical equilibrium table, with a simple quench layer approximation used to account for atmospheric mixing.The abundance tables are calculated with easyCHEM, our CEA (Gordon & McBride 1994;McBride & Gordon 1996) clone described in Mollière et al. (2017).The equilibrium condensation of the following species is included in the abundance calculations: Al 2 O 3 , Fe, FeO, Fe 2 O 3 , Fe 2 SiO 4 , H 2 O, H 3 PO 4 , KCl, MgSiO 3 , Mg 2 SiO 4 , Na 2 S, SiC, TiO, TiO 2 , VO.Because rainout is expected to remove Si from the upper layers of the atmosphere (see, e.g.Lodders 2010), we do not include feldspars, thereby inhibiting the sequestration of Na and K at high temperatures.This is consistent with abundance constraints of the alkalis from systematic retrieval analyses of T-and Y-dwarfs (Line et al. 2017;Zalesky et al. 2019).
The chemical abundances (mass fractions) are tabulated as a function of pressure P, temperature T , carbon-to-oxygen number ratio C/O, and metallicity [Fe/H].The pressure ranges from 10 −8 to 1000 bar, in 100 points spaced equidistantly in log(P) space.The temperature ranges from 60 to 4000 K, in 100 equidistant points.The C/O values go from 0.1 to 1.6, in 20 equidistant points and the metallicity is tabulated for [Fe/H] values going from -2 to 1.84, in 31 equidistant points.Fourdimensional linear interpolation is used to interpolate the logabundances of all absorbers.If a T -P-[Fe/H]-C/O coordinate falls outside of the grid, the abundances interpolated to the closest boundary point are used.The C/O ratio is varied by varying the oxygen abundance.For minimizing the Gibbs free energy we use the thermodynamic data of the CEA code, or the references detailed in Mollière et al. (2017).The data for FeH were obtained from M. Line (priv. comm.).
We approximate the effect of disequilibrium chemistry by setting the quench pressure P quench as a free parameter.For atmospheric pressures P < P quench we take the abundances of CO, H 2 O, and CH 4 to be constant, and equal to the abundances at P = P quench .This follows the result from, for example, Zahnle & Marley (2014), namely that the abundances of CO, H 2 O, and CH 4 can be taken to be constant above the quenching point.This treatment has been further verified by comparing to the results of the reaction network by Venot et al. (2012Venot et al. ( , 2015) ) in Baudino et al. (2017).
The chemical abundance table also contains the value of the adiabatic temperature gradient, ∇ ad , which was calculated with easyCHEM, using Equations 2.50, 2.59, and 2.75 of Gordon & McBride (1994).The derivatives used for the calculation of the specific heat of the mixture are so-called equilibrium derivatives (Gordon & McBride 1994), meaning that the ∇ ad value used in this work accounts for any change in the abundances (and thus heat release) during the adiabatic temperature change.Hence our adiabats are moist adiabats.The interpolation in the ∇ ad table is used when constructing the temperature profile in the troposphere of our temperature model.

Cloud model 1
As discussed in Section 1, our goal is to impose as few prior assumptions on the cloud properties as possible.Hence the ideal setup would be to freely retrieve the altitude-dependent distribution of cloud particle radii, as well as the vertical cloud density for every layer, independently.Such an approach would require many free parameters, which in turn requires enough data points of sufficiently high signal-to-noise to prevent over-fitting.Instead, we start with a more modest approach.As the Ackerman & Marley (2001) model is already implemented in petitRAD-TRANS, we use its three free parameters to control the mean particle size, cloud mass fraction and particle size distribution independently.
Hence, for the cloud to be retrieved, we set the settling parameter f sed as a free parameter, which controls the altitudedependent cloud mass fraction X c above the cloud base via where the cloud mass fraction at the cloud base X c 0 is an additional free parameter.The pressure at the cloud base P base is found by intersecting the saturation vapor pressure curve of the considered cloud species with the temperature profile of the atmosphere.In our retrievals below, we only consider MgSiO 3 and Fe clouds, because silicates and iron likely dominate the atmosphere at the temperatures and surface gravities reported for HR 8799 planets (see, e.g., Figure 7 in Morley et al. 2012, andMarley et al. 2012;Charnay et al. 2018).The recent findings of Gao et al. (2020), based on micro-physical cloud modeling, corroborate the importance of especially silicate clouds at these temperatures.
We also let the vertical eddy diffusion coefficient K zz vary as a free parameter, which effectively sets the particle size, given an f sed value.In self-consistent calculations, the K zz parameter is usually set by mixing length theory with some lower limit assumption (see, e.g., Ackerman & Marley 2001;Morley et al. 2014;Mollière et al. 2017;Charnay et al. 2018), or fixed to a constant value (e.g., Marley et al. 2012;Samland et al. 2017).
Here we let it float as a free parameter (and take it to be vertically constant), so as to determine the average particle size independently from f sed .We note that a retrieved K zz value could be inconsistent with the derived chemical quench pressure (see Section 2.3).This could imply either a true shortcoming of our retrieval model, a shortcoming in a disequilibrium kinetic network, or a deviation from how K zz actually sets the average particle size to how it is implemented in the Ackerman & Marley (2001) model.
Lastly, the particle size distribution is fitted by letting the width of the log-normal size distribution, σ g , be a free parameter.This parameter is usually not varied if the Ackerman & Marley (2001) cloud model is used.We note that it has been shown that a log-normal particle size distribution can be a poor choice when compared to the often bi-modal particle size distributions found from microphysics (Gao et al. 2018;Powell et al. 2018).However, in Gao et al. (2018), the K zz and σ g values were both fixed when fitting f sed to the cloud structure of the microphysics result.Here we let f sed , K zz , and σ g vary independently, so the retrieval should be flexible enough to determine the values of these three parameters that describe the cloud mass fraction, effective particle size, and dispersion of sizes around that value, independently.In principle, this model can also describe monodisperse particle distributions, if a retrieval were to favor cases with σ g close to 1.In the limit σ g → 1 a log-normal particle size distribution approaches a delta function.In general, the retrieved cloud parameter values are expected to describe those visible atmospheric layers which are most affected by the clouds.

Cloud model 2
The choice of f sed , K zz and σ g in Cloud Model 1 may just be a glorified way of setting the cloud spectral slope and single scattering albedo.It is also questionable, for example, whether a retrieved K zz does actually correspond to the true vertical diffusion coefficient of the atmosphere.Rather, it could effectively be a nuisance parameter of the retrieval, varied to mimic the true properties of the cloud opacity, as alluded to above.
To test for a less physically motivated treatment, we constructed Cloud model 2. This approach is motivated by the cloud parameterization of Burningham et al. (2017).Our model retrieves the spectral slope ξ (which we take to be vertically constant) of the cloud opacity directly.Specifically, we set where κ tot is the total (scattering + absorption) cloud opacity, κ(P) is its value at 1 µm, at pressure P, and λ the wavelength.
The κ(P) we describe as and set it to zero for pressures larger than the cloud base pressure P base .The f sed again describes the power law decrease of the cloud with altitude.Additionally, we set the single-scattering albedo, ω, as a free parameter, which we also take to be vertically constant.For the absorption opacity it then holds that In summary, the five free parameters of Cloud model 2 are κ 0 , ξ, f sed , P base , and ω.Because there is a degeneracy between κ 0 and P base (see Equation 5) if the atmosphere below the cloud deck cannot be probed by the observations, we put a prior on Table 1.Parameters for generating the synthetic observations of the cloudy exoplanet spectrum retrieved for verification purposes in Section 3. X eq is the mass fraction predicted for the cloud species when assuming equilibrium condensation at the cloud base location.
where P Fe and P MgSiO 3 are the cloud base positions of Fe and MgSiO 3 , respectively, obtained from intersecting their saturation vapor pressure curves with the atmospheric temperature profile.
Clearly, this will also tend to favor clouds that lie close to the expected cloud base positions for these condensate species.We note that this prior choice was made to not introduce parameters which are degenerate by construction.As mentioned, silicate and iron clouds may be the dominant cloud opacity carriers for the temperatures and surface gravities generally inferred for the HR 8799 planets.Removing this prior would allow to better probe situations where different cloud species dominate in the planetary atmosphere.This will be tested in our future work.
To test how well such a cloud description may be suited to describe more "physically consistent" clouds we carried out the following test.The cloud properties of the atmospheric model we used to create the synthetic observation for the verification retrieval (see Section 3) were generated with Cloud model 1.Fitting these cloud opacities with Equation 5, and the spectral slope with Equation 4, we found that the cloud parameters could indeed be well represented with Cloud model 2. In addition, the single-scattering albedo, taking the spectral average over the 0.9 to 2.5 µm range, was consistent with a large (ω ∼ 0.85), vertically constant value.

Retrieval with Cloud model 1
Here we present our retrieval tests when using Cloud model 1 (see Section 2.4), that is, MgSiO 3 and Fe clouds parameterized using the Ackerman & Marley (2001) model, varying all of its three free parameters.We generated a synthetic observation as follows: for the input parameters of the temperature profile, we fitted our temperature model to the self-consistent atmospheric structure used for verifying our scattering implementation (see Section 2.1 and Appendix B).The values of all input parameters of the model are shown in Table 1.
To derive the posterior abundances of our fit we used the PyMultiNest5 package (Buchner et al. 2014), which is a Python wrapper of the MultiNest method (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2013) ) for nested sampling (Skilling 2004).Nested sampling has the benefit of being able to approximate model evidences (i.e., the probability of the model, given the data), which allows for the pair-wise vetting of different models.Moreover, it is sampling the parameter space more thoroughly.This minimizes the problem of sampling the posterior distribution around Wavelength (micron) < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h j L Z b t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 a local, but not the global, maximum of the log-probability.It does not fully alleviate this problem, however (see discussion below).To ensure a high sample acceptance fraction of our highdimensional model, we ran MultiNest in the constant efficiency mode, with a prescribed sampling efficiency of 5 %.When using MultiNest in Importance Nested Sampling mode, evidences can still be calculated, even when prescribing a sampling efficiency6 .We used 4000 live points in our retrievals, which we found necessary to cover the parameter space sufficiently.
To initially test our retrieval framework under idealized conditions, a synthetic observation was created assuming a continuous wavelength spacing of λ/∆λ = 400 between 0.95 and 2.45 µm.We focused on this spectral region because it overlaps with the YJHK-bands of the SPHERE, GPI and GRAVITY instruments.The flux error was chosen to be constant across this wavelength range, with a mean S/N value of 10 per wavelength step.For comparison, the S/N per wavelength step of HR 8799e is about 4 (λ/∆λ ≈ 70), 7 (λ/∆λ ≈ 200), and 11 (λ/∆λ ≈ 1000) for the SPHERE YJH (Zurlo et al. 2016), GPI H (Greenbaum et al. 2018) and GRAVITY K (Gravity Collaboration et al. 2019) band data, respectively.In order to not be affected by a given noise instantiation during the verification retrieval, we took the observational errors into account for calculating the loglikelihood, but did not perturb the mock observations using these error bars.
The results of this verification retrieval are shown in In Panel (c) the retrieved pressure temperature structure of the atmosphere is shown, with the percentiles setting the boundaries of the uncertainty envelopes corresponding to the 1-, 2-and 3-σ ranges of a Gaussian distribution.In order to illustrate which altitudes of the atmosphere can actually be probed by the observation, the opaqueness of the temperature uncertainty envelopes has been scaled by the atmospheric contribution function, with a minimum value of 10 %.For this the contribution function was flux-averaged.As can be seen, the uncertainty envelopes follow the input P-T profile.The input profile lies within the 1-σ envelope.
Finally, Panel (d) shows the corner plot of the remaining parameters.The planetary radius, surface gravity, metallicity and C/O ratio can all be well retrieved.Only an upper limit is found for the quench pressure, which is as expected, as no quenching was considered.The cloud parameters are well constrained.We note that the retrieval appears to have found a bi-modal parameter distribution, but the one-dimensional posteriors constrain the input parameters well.A clear positive correlation can be seen between the metallicity and C/O ratio.This is attributed to the fact that we vary the oxygen abundance when changing C/O, such that a higher C/O corresponds to less oxygen, and hence water.
In summary, the retrieval verification test presented here can be regarded as successful.We were able to retrieve the temperature, composition, and cloud properties of the atmosphere.Nonetheless, the following challenges can be identified: even though an excellent fit to the spectrum has been achieved (see Panel (a) of Figure 1) the median values of the retrieved parameters were not exactly at the input values (although within the 1 σ envelope).This is unexpected because the synthetic observation was not perturbed by the assumed error bars, in order to not be sensitive to stochastic noise of a given noise instantiation.Another point is the bi-modality of the posterior.
In our initial tests we found that using 400 live points in PyMultiNest resulted in biased retrieval results, such that median parameter values could be more than 1 σ away from the input values, and residuals in the spectral fit that were larger.Increasing the number of live points to 4000 ameliorated this issue: the input parameters were retrieved at higher accuracy (within 1σ), and the residuals to the best-fit spectrum became smaller.We deduce from this that for input models of high dimensionality a sufficient number of live points has to be used, so as to increase chances that the positions of the live points sampled during the early stages of the nested sampling run will fall into the vicinity of the global maximum of the likelihood.Because the nested sampling method will zero-in on the highest likelihood regions during the retrieval, the danger exists that the true maximum location in parameter space will be missed.This problem is especially pressing for observations of high S/N, such as used in our example presented here, because the high-likelihood volume of the parameter space will shrink.This means that higher quality observations require a larger number of live points, and thus more computational time.This is especially important for the large spectral coverage, high S/N data to be taken with JWST.
Further exploring the bi-modality of the posterior shown in Figure 1, we ran a second retrieval with smaller prior ranges.They were restricted by the high-likelihood regions of the posterior from the initial fit, enclosing its bi-modal posterior distribution.The resulting posterior is shown in the corner plot in Figure C.1.It is unimodal and consistent with the input parameters.In general, the bimodality and the offset between the median and input parameters of the posteriors indicate that multiple parameter combinations can lead to excellent spectral fits, while within 1 σ of the true parameters.This may indicate that a retrieval model with fewer free parameters could be favored, leading to unique solutions, which we will explore in future studies.

Retrieving Cloud Model 1 with Cloud model 2
In this section we describe what happens when retrieving a mock observation made with Cloud Model 1 using Cloud Model 2. The retrieval was thus set up identically to the one described in the section immediately above, but used Cloud Model 2, while the mock observation was identical to the test above, that is made with Cloud Model 1.We find that the spectral fit is again very good, although there are a few regions of systematic residuals of 1 σ.Hence, Cloud model 2 appears to satisfactorily describe the properties of the synthetic observations generated with Cloud model 1.However, we find significant differences in the retrieved atmospheric properties.C/O is constrained to 0.58 +0.01  −0.01 (input was 0.55), [Fe/H] is retrieved to be 0.14 +0.08  −0.08 (input was 0), and log(g) = 4.01 (input was 3.75).Thus, we find values close to, but offset from the true input values.Moreover, instead of probing down to 2 bar at most, the atmosphere can now be probed down to 10 bar, and is more isothermal than the input temperature profile: here the retrieval mimics the effect of a thick cloud.Instead of the cloud hiding the deep hot regions from view, these regions are erroneously constrained to be less hot by the retrieval.At the same time the cloud is thus too optically thin and deep, with the retrieved cloud position at 8 bar (which is well constrained).The emission contribution consequently shows that the emission stems from a more extended region than in the retrieval described in the section immediately above.Hence Cloud Model 2 was not able to describe the clouds made with Cloud Model 1 accurately enough, such that the retrieval modified the P-T profile instead.
We conclude that a good fit to the spectrum alone is a dangerous measure when assessing whether or not a fit result is reasonable.All retrieved parameter values need to be carefully vetted, the retrieved P-T profile should also be compared to that of a self-consistent atmospheric code, when running the latter using the best-fit or median parameters of the free retrieval.We note that with free retrievals alone it may also be challenging to determine whether the atmospheric temperature profile is truly shallower than expected, while being less cloudy at the same time.This has been suggested by Tremblin et al. (2015Tremblin et al. ( , 2016Tremblin et al. ( , 2017)), challenged in Leconte (2018), and defended in Tremblin et al. (2019).As they suggested, we find here that a shallow P-T profile can indeed result in an excellent fit to the observations, even though we know that in this case the input model was actually more cloudy than retrieved.At the same time it is somewhat reassuring that in our example shown here the absolute deviation between the input and retrieved parameters such as C/O, [Fe/H], log(g), etc. is not large.The values, however, are biased, and we stress that a more detailed study needs to be performed on how strongly cloud model assumptions can affect the retrieved bestfit parameters.

Retrieving HR 8799e
In this section we describe how we used GRAVITY and archival SPHERE and GPI data to retrieve the atmospheric properties of HR 8799e, which is located at 15 au from its host star (e.g., Wang et al. 2018).We also compare our results to archival photometry for this planet, while not including these data in the retrieval itself.

GRAVITY K band spectroscopy
We use three separate GRAVITY (Gravity Collaboration et al. 2017) observations of HR8799e.First, the observation presented in Gravity Collaboration et al. (2019).In addition, we here report on two new observations, taken November 9th, 2019 and November 11th.2019, as part of the ExoGRAVITY Large Program.The log of the observations is presented in Table 2.The observations on the 9th were made from a short observing block, with a total integration time of 180 s.The phase referencing was done on the star using the fringe tracker (Lacour et al. 2019).The zero point of the metrology was obtained by directly observing the star on the spectrometer.The observations on the 11th were done using the roof mirror as a field splitter: 100% of the planetary flux could be used, but such an observation needs a binary to calibrate the zero point.This zero point was obtained on the binary system HD 25335.
The data were reduced analogously to the data reduction presented for β Pic b in Gravity Collaboration et al. (2020).The flux of stellar origin is removed, and the spectra were obtained from the ratio between the coherent flux on the planet from the coherent flux on the star.This ratio is then multiplied by a theoretical BT-NextGen spectra of the star (Allard et al. 2012a).The spectra is therefore calibrated from the telluric absorption.For the retrieval, the full spectral covariance is considered when deriving log-likelihoods.

SPHERE and GPI archival data
We use the YJH-band spectroscopy of SPHERE reported in Zurlo et al. (2016).In addition, we use the GPI H-band spectroscopy reported in Greenbaum et al. (2018).We do not take the spectral covariance into account for GPI or SPHERE.For both SPHERE and GPI we fit for a scaling factor with respect to the GRAVITY observation during the retrieval, as was also done in Gravity Collaboration et al. (2020).This also appears necessary given the noticeable shift between the SPHERE and GPI observations in their overlapping region at ∼ 1.6 micron.Similar shifts between GPI and SPHERE observations have been reported in Samland et al. (2017).However, this difference could also be caused by variability, because spectral template brown dwarfs that reproduce the spectral properties of HR 8799e well have been found to exhibit considerable variability, possibly up to a 20-30% peak-to-peak amplitude (see the discussion of Mace et al. 2013;Biller et al. 2015in Bonnefoy et al. 2016).

Archival photometry
Although not included during the fit, we compare our results with archival photometry in the mid-infrared.We consider the 3.3 µm LBT photometry reported in Skemer et al. (2012), the L' band and [4.05]-Brα photometry reported in Currie et al. (2014), as well as the M' band upper limit of Galicher et al. (2011).The photometry was converted from magnitudes to flux using the species7 toolkit, which has been described in Stolker et al. (2020).We decided against including the photometric fluxes in the fit because their relatively low signal-to-noise would add little constraining power to the retrieval, when compared to the spectra, but would double the run-time of our retrievals due to the increased spectral range.

Retrieval model setup
We set up our nominal retrieval model with 18 free parameters, using Cloud Model 1, which we describe in the following.Additionally, we will compare to retrievals using Cloud Model 2, see Section 4.4.The free parameters and prior ranges are listed in Table 3. Priors of the HR 8799e retrieval.U stands for a uniform distribution, with the two parameters being the range boundaries.
The units for the parameters are the same as the ones used for Table 1.f SPHERE and f GPI are the scaling factors retrieved for the SPHERE and GPI data, respectively.(a) and (b): please see Section 4.2 for a definition of P phot and T connect .It holds that Xi = X i 0 /X i eq , where the latter quantity has been defined in Table 1.
Table 3.As in Gravity Collaboration et al. ( 2020), we fitted multiplicative scaling factors f SPHERE and f GPI to account for systematic biases in the flux normalization of these datasets.The T connect quantity referenced in the prior range of T 3 is the uppermost temperature of the 'photospheric' layer, and was calculated by setting τ = 0.1 in Equation 2. Like the priors for T 2 and T 1 , this ensures a temperature profile that is monotonically decreasing with altitude, also see Section 2.2.δ was sampled from the prior by assuming a log-uniform prior on P phot , where we defined P phot as the pressure where τ = 1 in Equation 1.This allowed to solve for δ for a given P phot value.The following absorber species were included: CO, CO  2019) for the full reference list and description of the opacity sources.The FeH opacity has been multiplied by a factor 1/2 due to the partition function correction described in Charnay et al. (2018).We convolved the synthetic spectra using a Gaussian kernel, in order to approximate the line spread function of the SPHERE, GPI, and GRAVITY instruments.The instrumental resolving power was assumed to be 30, 45 and 500 for SPHERE, GPI and GRAVITY, respectively.The resolution element ∆λ of the spectrograph was assumed to be the FWHM of the line spread function.This means that the standard deviation σ LSF of the Gaussian kernel used for convolution is defined by ∆λ = 2 √ 2 ln2 σ LSF .The retrieval was run with PyMultiNest, using 4000 live points.

Pressure grid
We use an adaptive spacing for the atmosphere's pressure grid.For cloud-free calculations the atmosphere would be separated into 60 points, spaced equidistantly in log-pressure between 10 −6 and 1000 bar.For cloudy calculations the retrieval code considers all cloud base pressures P base and increases the spatial resolution for P ∈ [0.5 P base , 1.12 P base ] (corresponding to a pressure range of -0.3 and 0.05 dex) by a factor of 12.This better resolves the placement of the cloud base within the atmosphere, in case a P base does not fall into the immediate proximity of a grid point of the coarse pressure grid.In addition, the abrupt increase in atmospheric opacity at the cloud deck position, and its decline ∝ P f sed for lower pressures, is resolved better.For two spatially separated cloud decks this leads to 104 grid points.We found that this treatment is as accurate as running the whole calculation at a 12 times-increased resolution, which would result in 720 grid points.For this we compared to a baseline calculation made at a 24 times higher resolution, using 1440 grid points.

Testing the retrieval model
Because all test retrievals mentioned in Section 3 were carried out on data sets of homogeneous wavelength coverage in YJHK bands, we verified our retrieval setup by running a mock retrieval that had the same wavelength spacing and error properties as the actual HR8799e data sets of SPHERE, GPI, and GRAVITY.As input for the synthetic observation we used a posterior sample of the actual HR8799e retrieval, the result is shown in Figure E.1.We find that we can retrieve all parameters well, except for [Fe/H] and K zz , which are biased, as are the scaling parameters.Running a second fit that neglected the scaling parameters lead to a well retrieved [Fe/H], but slightly too small radius.This behavior could be due to the random noise instantiation used in the retrieval, and the fact that especially the scaling may introduce biases in the retrieved atmospheric parameters in cases where differences are introduced between model and observation.This has also been described in Kitzmann et al. (2019).In their case the differences arose from using two different models for generating the mock observation and the retrieval, here the differences arise from the random noise properties.We note that the scaling value we retrieve for the actual HR8799e data below is consistent with unity.
To test the impact of the random noise instantiation further we also ran test retrievals for the same synthetic observation, but neglecting the random perturbation of the data due to the noise.Similar to our noise-free test retrieval presented in Section 3.1, we found that the noise-free test led to a bi-modal posterior, with the modes bracketing the input values.An analogous approach (zooming in on the prior volume populated by the bi-modal posterior) lead to a uni-modal posterior, retrieving the input parameters.Thus we reconfirm our observation that noiseless test retrievals can lead to multi-modal posteriors if the ratio of the prior volume and the number of live points is large.This indicates that multiple parameter combinations can lead to excellent spectral fits, while within 1 σ of the true parameters.As stated in the manuscript before this may indicate that a retrieval model with fewer free parameters could be favored, leading to more unique solutions.We note that in these test retrievals described here, the atmospheric C/O ratio was a robustly retrieved parameter in all retrieval setups, and that our two retrievals for the real HR 8799e data presented below, using either Cloud Model 1 or 2, lead to consistent C/O, [Fe/H] and log(g) constraints.

Free retrieval results of HR8799e
In this section we describe the retrieval results.The results will also be discussed in view of the possible planet formation history, and compared to existing literature studies of HR8799e in sections 5.1 and 5.2, respectively.The spectral fit for HR8799e is presented in Figure 2. In general, the retrieval model is able to explain the observations well.The residuals scatter around zero, with some systematic differences visible at 1.325, 1.525, 1.725, 2.07, 2.11, 2.18, 2.275 and 2.425 micron.These differences could be due to the model missing absorbers, or not being flexible enough to fit intricacies in the atmospheric temperature, abundance, or cloud structure.Another likely possibility are remaining systematics in the observations: the difference between the SPHERE and GPI observations in their overlap region (1.525 micron), as well as the overall wiggly appearance of the GRAVITY observation may indicate this.We note that also the photometric flux measurements in the MIR are fit well, except for the narrow [4.05] band measurement by Currie et al. (2014). 8This is especially interesting as these points were not included in the retrieval.Especially the 3.3 µm LBT -L' and L'-[4.05]colors have been noted to be difficult to explain with self-consistent models, see discussion in Section 5.2.
From the spectral appearance of the H and K band observations it is already clear that CH 4 is not an important absorber in the atmosphere: the flux decrease expected from CH 4 absorption at 1.6 and 2.2 micron is absent.The fact that only an upper limit is found by observations in the M-band, together with the comparatively high flux in the 3.3 µm LBT and L' bands, also speaks for a CH 4 -poor atmosphere.
The retrieved pressure-temperature uncertainty envelopes of HR8799e are shown in the left panel of Figure 3.As in the analogous plot shown in Figure 1, the opacity of the uncertainty envelopes is scaled by the flux-averaged emission contribution function of the best-fit model, to show where in the atmosphere the observations are probing.This contribution function is also shown in the right panel of Figure 3.For comparison, we also show the P-T curves derived with our self-consistent petitCODE using the best-fit parameters of the retrieval.This is further discussed in Section 4.5.
The corner plot of the one-and two-dimensional projections of the 18-dimensional posterior distribution of the retrieval is shown in Figure 4. We summarize a few of the most striking results here, while the implications of the retrieved parameter values will be discussed in sections 5.1 and 5.2.In general, we note that none of the retrieved parameters ran into its prior boundaries and that all parameters (except for the Fe mass fraction at the cloud base) are well constrained.In addition, we derive that the atmospheric C/O ratio is C/O = 0.60 +0.07 −0.08 and the retrieved metallicity is [Fe/H] = 0.48 +0.25  −0.29 .Together with the planet's mass, which we derive to be 4.81 +8.78  −3.33 M J , this has important implications for how the planet could have formed, see Section 5.1.We note that the surface gravity and radius retrieved for HR8799e, log(g) = 4.00 +0.46  −0.52 and R P = 1.12 +0.09 −0.09R J , are constrained with a symmetric, uni-modal peak.Hence also the logarithm of the planet mass is constrained with a symmetric, uni-modal peak, but the distribution of the mass itself is skewed towards lower masses, with large mass uncertainties due to the large uncertainty on log(g).The atmosphere is clearly affected by disequilibrium chemistry, with a large quench pressure of log(P quench /1 bar) = 1.35 +0.50  −0.56 .We find that disabling quenching at the best-fit parameters leads to strong CH 4 absorption features in the spectrum, which are inconsistent with the data.Moreover, the scaling value retrieved for SPHERE is consistent with one, 1.01 +0.08  −0.07 .The scaling retrieved for GPI is significantly smaller than 1, namely 0.90 +0.06 −0.05 .Only when applying this scaling on GPI do the SPHERE and GPI datasets agree in their overlapping region, see Figure 2. From sampling the posterior distribution 300 times, and calculating the spectra between 0.5 and 28 micron we derive an effective temperature of T eff = 1154 +49 −48 K. Using this derived temperature, our retrieved surface gravity, and This plot shows the one and two-dimensional projections of the 18-dimensional posterior distribution for all but the six nuisance parameters of the temperature structure (see the left panel of Figure 3 for the P-T confidence envelopes).M P,spec is the mass of HR 8799e in units of Jupiter masses, derived from the retrieval posterior distributions of log(g) and R P .
Equation (4) of Zahnle & Marley (2014), we find that the upper limit for the atmospheric mixing is log(K zz,max ) = 10.2, which assumes that all flux is transported by convection.Our derived value, which is used solely for determining the particle size for a given f sed , is log(K zz ) = 9.80 +1.15 −1.39 , so below the theoretical upper limit (while still large).

Retrieval with the non-nominal Cloud Model 2
To test the robustness of the constrained atmospheric properties, we also retrieved HR8799e with Cloud Model 2. Like before, we find that Cloud Model 2 leads to retrieved temperature gradients which are too shallow when compared to physically consistent solutions.The retrieved atmospheric solution is bi-modal, with one solution corresponding to a P-T structure with a shallow temperature gradient and intermediate cloudiness, and a second solution corresponding to an even more isothermal, cloud-free atmosphere.The spectral fit, full posterior distribution, and P-T uncertainty envelopes are shown in Appendix F, where we also describe the prior setup of this retrieval.
Focusing on the bulk atmospheric properties (log(g), C/O, [Fe/H]), we find that the cloudy mode of the retrieval with Cloud Model 2 is fully consistent with the results from our nominal retrieval with Cloud Model 1.The one-dimensional posterior dis- tributions look almost identical.Thus, even though the atmospheric temperature and cloud structure retrieved with different cloud models can differ, quantities such as C/O, [Fe/H] and the atmospheric gravity may be very robust.The clear, isothermal solution of Cloud Model 2 leads to different values for these retrieved parameters, but their 1-σ uncertainty regions overlap.Figure 5 shows the marginalized one-dimensional posterior distributions of HR8799e's gravity, metallicity, and C/O, derived with Cloud Model 1 and Cloud Model 2.

Comparison of the results with self-consistent atmospheric models
Free retrievals allow to deviate from the rigidity of (potentially imperfect) physical assumptions in self-consistent codes.This can be both boon and bane.On one hand it allows accounting for effects that influence the atmospheric structure which are not adequately captured by the physical assumptions made in the self-consistent codes.On the other hand, the free retrieval may converge on parameter results which lead to a seemingly good fit to the data, but are unphysical.Because of this, it is crucial to verify retrieval results by comparing to constraints that can be obtained from self-consistent codes, as was also done in, for example, Line et al. (2017); Gandhi & Madhusudhan (2018).
We describe such tests below: we compared our derived atmospheric structure with self-consistent results and compared our free retrieval to a grid interpolation retrieval with self-consistent atmospheric spectra.
Self-consistent P-T structures of petitCODE petitCODE (Mollière et al. 2015(Mollière et al. , 2017) is a self-consistent code for calculating atmospheric structures and spectra.It assumes radiative-convective and chemical equilibrium to calculate the atmospheric structure.The radiative transfer includes scattering, and the code can include gas line, continuum, and cloud opacities.
We carried out two tests: for the overall planet parameters (log(g), T eff , [Fe/H], C/O), we used the median of the retrieved values.We then calculated a cloud-free atmospheric structure, assuming chemical equilibrium.As a second test we included clouds, prescribed the median values of the retrieved cloud parameters ( f sed , K zz , σ g , XFe , XMgSiO3 ), and enforced that the H 2 O, CH 4 and CO abundances be held constant below the retrieved best-fit quench pressure.The resulting structures are shown as gray and black solid lines in the left panel of Figure 3 for the cloud-free chemical equilibrium and the cloudy non-equilibrium structure, respectively.
Overall, it can be seen that the self-consistent structures follow the uncertainty envelopes of the retrieved pressure temperature structure well, in terms of absolute temperature and slope.Interestingly, we notice that the cloud-free structure in chemical equilibrium falls within the 1-σ envelope at all pressures, while the self-consistent cloudy structure which included quenching moves out into the 2-σ envelope between 0.05 and 0.3 bar.This is above the region of maximum emission as measured by the contribution function of the best-fit model, however.We conclude that we do not see any clear deviation of our retrieved temperature envelopes when compared to physical expectations, both the overall shape and absolute temperatures appear to be close to what is predicted in a self-consistent model.

Spectral fit with Exo-REM
In addition to the free retrieval with petitRADTRANS described above, we carried out a grid-interpolation retrieval to obtain the atmospheric parameters of HR8799e from its spectrum.The grid of self-consistent model spectra was obtained with Exo-REM (Baudino et al. 2015(Baudino et al. , 2017;;Charnay et al. 2018), in the version by Charnay et al. (2018), which includes scattering and disequilibrium chemistry.The chemical disequilibrium and cloud scale height is determined through taking into account the vertical atmospheric mixing, which is set through the atmospheric eddy diffusion parameter K zz .In Exo-REM, K zz is determined from the atmospheric structure in the convective region consistently, using mixing length theory.Above the convective region, K zz is determined from a convective overshooting description.Exo-REM includes the cloud opacities of spherical, amorphous Fe and Mg 2 SiO 4 grains, as well as the gas opacities of Na, K, H 2 O, CH 4 , CO, CO 2 , NH 3 , PH 3 , TiO, VO, and FeH, in addition to H 2 -H 2 and H 2 -He collision-induced absorption (CIA).The Exo-REM grid used here ranged in T eff from 1000-2000 K (∆T eff = 50 K), and in C/O from 0.3-0.75(∆C/O = 0.05).The [Fe/H] grid points were −0.5, 0, 0.5, while the log(g) points were at 3.5, 4, 4.5.
Notable points of difference in the opacity treatment between Exo-REM and petitRADTRANS are the use of different alkali wing profiles (Burrows & Volobuyev 2003 for Exo-REM, Allard et al. 2003Allard et al. , 2012b for petitRADTRANS) and that pe-titRADTRANS assumed irragularly shaped, crystalline Fe and MgSiO 3 cloud particles instead of spherical, amorphous Fe and Mg 2 SiO 4 ones.Moreover, Exo-REM varies C/O by changing C, whereas petitRADTRANS varies C/O by changing O, which can make a difference (Lodders 2010).
Applying a LSF convolution and rebinning identically to the petitRADTRANS retrieval, we used species to carry out a gridbased retrieval with PyMultiNest, where the model spectra were obtained by linearly interpolating within the Exo-REM grid. 100 spectra sampled from the Exo-REM posterior are shown in Figure 6, analogous to the spectra shown for the free retrieval with petitRADTRANS in Figure 2.
In general the quality of the spectral fit is very similar to the free retrieval with petitRADTRANS, with the difference that the residuals in the near-IR are a somewhat larger.Like petitRAD-TRANS, Exo-REM can fit the photometry in the mid-IR well, Fig. 6.Same as Figure 2, but for a grid-based interpolation retrieval with the self-consistent code Exo-REM.even though it was not included in the fit.Analogous to the pe-titRADTRANS fit, also Exo-REM does not reproduce the [4.05] band photometric point well, but is somewhat more consistent than petitRADTRANS.
We show the one-dimensional marginalizations of the Exo-REM posterior distribution in Figure 7, together with the marginalized posteriors of the corresponding parameters of the free petitRADTRANS retrieval.Because the effective temperature T eff is not a free parameter of the petitRADTRANS retrieval, the T eff distribution was obtained from sampling the posterior 300 times, and calculating the effective temperature from the spectrum ranging from 0.5 to 28 µm.We find two modes in the Exo-REM retrieval.The first mode has lower temperatures, smaller GPI and SPHERE scaling factors, and larger planet radii when compared to the petitRADTRANS retrieval.For the second mode the scaling factors, T eff and radius are consistent with the petitRADTRANS results.The surface gravity agrees with the broad petitRADTRANS posterior, for both modes.The lower temperature mode of the Exo-REM fits runs into the lower grid boundary for the effective temperature.Both modes run into to lower grid boundary for C/O, and in the upper grid boundary for [Fe/H].Due to the inherent difficulty of converging selfconsistent cloudy models the Exo-REM grid is incomplete, that is, models in the grid are missing especially at lower temperatures.This could affect the fit at low effective temperatures.Therefore, as a second test, we fixed the GPI and SPHERE scaling factors to the best-fit values of the petitRADTRANS fit, which is also shown in Figure 7.This leads to effective temperatures and radii consistent with the petitRADTRANS fit, while still running into the C/O and [Fe/H] boundaries.
For both the metallicity and C/O ratio the comparison between the petitRADTRANS and Exo-REM results is difficult: the retrieved best-fit metallicity value from the Exo-REM grid is consistent with the petitRADTRANS peak metallicity.However, the Exo-REM retrieval runs into the upper boundary of the grid.For the C/O ratio the petitRADTRANS fit peaks at around 0.6, while it is driven into the lower grid boundary (0.3) for the Exo-REM fit.The behavior of the Exo-REM posterior for these two quantities therefore makes a more detailed comparison of the pe-titRADTRANS and Exo-REM retrievals difficult.It appears as if the Exo-REM grid may not be able to fully reproduce the near-IR photometry, and trends into the boundaries of the grid in search of the true probability maximum that is not contained within its grid boundaries.
We can speculate at the near-IR region being the cause for the difficulties we face in the Exo-REM retrieval: here the residuals are larger than in the petitRADTRANS retrieval.Because the YJH bands are the most strongly affected by clouds, this could hint at a difference in the description of clouds, but also the alkalis could play a role, especially in the SPHERE Y band.Indeed, when only retrieving the atmospheric properties using the GRAVITY K-band with Exo-REM, we find that the C/O is constrained at 0.43 +0.07 −0.07 , while [Fe/H] still trends into the upper grid boundary (0.5).A larger grid extend with additional [Fe/H] and C/O values, may therefore alleviate the problems we face in our analysis here.However, we expect that a grid that also varies the cloud parameters of Exo-REM may help, as we find that decreasing the C/O (or increasing the metallicity) leads to stronger water absorption features across the YJH bands.This speaks for C/O being used to counteract too strong cloud absorption.
We conclude that the petitRADTRANS and Exo-REM retrievals are consistent with each other in the parameters that can be easily compared.However, parameters trending into the grid boundary in the Exo-REM retrieval make the comparison difficult for [Fe/H] and C/O.

Implication of the retrieved C/O and [Fe/H] for the formation location of HR 8799e
Measuring a planet's C/O has been suggested as a powerful tool to trace where in the protoplanetary disk a planet may have formed (e.g., Öberg et al. 2011;Madhusudhan et al. 2014a;Mordasini et al. 2016;Cridland et al. 2016).The general idea is to compare a planet's C/O to the C/O predicted for the disk's solid and gas phases, and using the planet's bulk enrichment to determine whether the planet's metal enrichment is dominated by solid or gas accretion.From this it may be possible to infer where in the disk a planet formed.If this formation location were to be conclusively shown to be further away from the star than the planet's current location, this would be a proof for orbital migration.In what follows we attempt such an analysis using the C/O ratio inferred for HR 8799e.We neglect the, potentially very important, effect of compositional gradients in the planet, that may lead to atmospheric abundances being different from the bulk of the planet (e.g., Leconte & Chabrier 2012;Vazan et al. 2018).
The C/O value that has been reported for the host star HR 8799 is C/O star = 0.56 (Sadakane 2006).The author finds that C and O have roughly solar abundances.The star is a λ-Boötis-type star, meaning that its iron peak elements are subsolar ([Fe/H] measured for iron specifically is −0.55 ± 0.1, Sadakane 2006).Because the analysis presented in this sub-section hinges on C/O and the total metal content of the planet, we neglect this additional information for now.We note, however, that if the composition of the photosphere of HR 8799 is representative for the disk from which its planets formed, assuming solar abundance ratios for all elements except O 9 during the atmospheric modeling is problematic.
To constrain HR 8799e's formation history, we make the assumption that the disk from which the HR 8799 planets formed had a composition as specified in Table 1 of Öberg & Wordsworth (2019).The authors of the paper present this as a model for the young solar nebula.Because the C/O ratio of HR8799 (C/O = 0.56) is essentially solar (C/O = 0.55, see Asplund et al. 2009) and we are most interested in the molecular volatiles, which represent the largest mass reservoir for solid planetary building blocks, we deem this assumption acceptable, even though HR8799 is a λ-Boötis-type star.
Using the planetary mass of 4.81 +8.78 −3.33 M J , the metallicity of 0.48 +0.25  −0.29 and the C/O of 0.60 +0.07 −0.08 , inferred from our spectral retrieval with petitRADTRANS, we find that HR 8799e is likely heavily enriched in ices, accreted in forms of pebbles and planetesimals, and most likely formed outside of the CO iceline.This conclusion was obtained from fitting the O/H and C/H content of HR 8799e, derived from our spectral retrieval, with an abundance model of a planet that forms in a disk as defined in Öberg & Wordsworth (2019).For this we treated the planets mass, accreted solid mass, and the accretion locations of the solids and gas as free parameters.For the planet mass a prior based on our spectrally retrieved mass was assumed.The planetary O/H and C/H was then fitted with PyMultiNest.Specifically, due to the increased planetary metallicity of 0.48 +0.25  −0.29 , we find that the planet has accreted between 65 and 360 M ⊕ of ices (1-σ range) that are mixed into its envelope and atmosphere.The large uncertainty stems from the large mass and metallicity uncertainties from our spectral retrieval.The atmospheric metal content of HR 8799e appears to be dominated by solids due to the large inferred atmospheric enrichment (e.g., Espinoza et al. 2017).If the metal content were dominated by gas accretion the atmospheric metallicity is expected to be close to, or smaller than, stellar.As we find that the planet has a C/O ratio consistent with its host star (we find C/O = 0.60 +0.07 −0.08 for HR8799e) this could mean that the planet has formed outside of the CO iceline.In particular, we derive that a formation location outside the CO iceline is more than twice as likely compared to a formation inside the CO iceline (more details on this analysis will be published in an upcoming study).This is explained by the fact that only outside the CO iceline the solid material in the disk will have the stellar value, see Öberg et al. (2011).In contrast, C/O values of 0.3, which is the lower boundary of the atmospheric grid approached in our Exo-REM retrieval, are possible if the planet formed within the CO iceline, where H 2 O and CO 2 dominate the ice composition.If HR 8799e did form outside the CO iceline, this could mean that all HR 8799 planets formed outside the CO iceline, because HR8799e is the innermost planet of the HR8799 system.This will have to be tested by deriving C/Os and metallicities for all HR 8799 planets.The C/O analysis of the HR 8799 planets by Lavie et al. (2017) is consistent with this assessment: for the HR 8799b and c, where the authors succeeded in deriving C/O values, they found C/O ≥ C/O star .
If HR 8799e truly formed outside the CO iceline, this also allows to put constraints on the planet's possible migration.In Öberg & Wordsworth (2019), the CO iceline is situated at ∼20 au for the young solar nebula.In the disk around HR 8799, the same temperature due to the irradiation of the star would have been 9 We vary C/O by varying O.
reached at ∼45 au (Marois et al. 2010), neglecting the evolution of the stellar luminosity at young ages.Because HR 8799e resides at ∼ 15 au (Wang et al. 2018), this could imply that the planet migrated significantly.This would be consistent with the finding by Wang et al. (2018) that the HR 8799 planets needed to migrate in the gas disk after formation to get locked into a stable resonant orbit.
The model for the disk composition in Öberg et al. (2011); Öberg & Wordsworth ( 2019) is strongly simplified.The disk's properties such as temperature, surface density, and abundance profiles are assumed to be static.The iceline positions are determined from a simple thermodynamic stability analysis of the ice species.Processes such as the viscous and chemical evolution of the disk are neglected.However, chemical evolution and ionization of the disk material can be of crucial importance, as well as the initial composition of the disk at the start of the evolution, as shown by Eistrup et al. (2016Eistrup et al. ( , 2018)).Importantly, these studies describe how gas-grain chemistry may deplete CO from the gas phase within the CO iceline, condensing it in the form of CO 2 , at the expense of also H 2 O. Thus our conclusion regarding migration, based on the C/O ratio of HR 8799e, may be based on oversimplified disk chemistry assumptions.
To assess the effect of properly treating the disk chemistry we used the ANDES physical-chemical code to compute a 2D steady-state disk physical structure and time-dependent chemistry for the HR 8799 disk (Akimkin et al. 2013;Molyarova et al. 2017).The detailed setup of the model is described in Section G.For this setup of the disk chemical model we found that the CO iceline lies at around 100 au, which would indicate that HR 8799e migrated even further after formation.However, due to the above-mentioned gas-grain chemistry, CO is converted into CO 2 ice effectively starting from around 20 au.This makes the solid C/O in the disk approach stellar values already at 20 au, such that HR 8799e may have formed as close as 20 au from the star.Hence realistic disk chemistry modeling could indicate that the planet migrated much less than when compared to simplistic disk abundance models.
Finally, we note that the nitrogen content may be a better way of constraining a planets formation location in the disk ( Öberg & Wordsworth 2019;Bosman et al. 2019), where a large N-content corresponds to a formation in the outer parts of the disk.However, this would require to study planets cooler than HR 8799e, for which most of its accreted nitrogen is in the form of N 2 , and therefore invisible due to the low N 2 opacity.

Comparison of retrieval results with literature studies
Since their discovery (Marois et al. 2008(Marois et al. , 2010;;Currie et al. 2011) the HR8799 planets have been extensively studied.We provide a summary of the studies that exist on the HR8799 planets below and how they relate to our results for HR8799e.We start with the qualitative properties of the planet, before comparing our retrieved values with the ones reported by others.

Clouds
In general, studies find that all HR8799 planets have comparable surface gravities and temperatures.In addition, all studies find that the HR8799 planets are dominated by thick clouds, as indicated by their red near-infrared (NIR) colors.For example, using NIR and MIR photometry and a grid of selfconsistent atmospheric models, Madhusudhan et al. (2011a) find that HR8799bcd are dominated by clouds much thicker than expected for field brown dwarfs.This is similar to an assessment already made by Bowler et al. (2010), studying HR8799b.Marley et al. (2012) came to the conclusion that thick clouds are required when studying HR8799bcd, but emphasized that this is not peculiar, and rather a consequence of cloud formation being gravity-dependent, a result that has also been borne out by the model calculations of Charnay et al. (2018).
Considering the apparent lack of comparison brown dwarfs that resemble the HR8799 planets, it is important to remember that the earlier results in the literature are highly heterogeneous in terms of the available data (and models) that were used to arrive at a given conclusion.High quality spectra are important to truly unlock the planetary characteristics.For example, additional SPHERE spectroscopy obtained by Zurlo et al. (2016) allowed Bonnefoy et al. (2016) to show that especially HR8799de can be well fit with low-gravity cloudy brown dwarfs of the late L spectral type.For HR8799bc the picture is less clear, and Bonnefoy et al. (2016) find that good comparison objects can only be identified when reddening T-dwarf spectra with iron or silicate extinction.Given the abundance of literature reporting on the cloudiness of the HR 8799 planets, our finding that HR 8799e is cloudy therefore does not come out of the blue (sky).

Disequilibrium chemistry
Disequilibrium chemistry has also been reported in the HR8799 planets by a variety of studies.For planets such as HR 8799bcde, disequilibrium chemistry means that CO is more and CH 4 is less abundant then predicted from chemical equilibrium, due to atmospheric mixing overruling the chemical reactions in the upper atmosphere.Using OSIRIS H and K band spectra at low resolution, Barman et al. (2011) report on HR8799b exhibiting weaker CH 4 absorption than expected, and that disequilibrium chemistry may be at play in order to decrease the CH 4 abundance.Similar findings were also reported for the medium-resolution OSIRIS data for planets b and c, where CH 4 was not detectable in c (Konopacky et al. 2013;Barman et al. 2015).Madhusudhan et al. (2011a) found that they had to neglect the 3.3 µm band containing too strong CH 4 absorption when fitting their chemical equilibrium models to the HR8799bc data obtained by Currie et al. (2011) (there was only an upper limit available for d in the 3.3 µm band).Marley et al. (2012) reported disequilibrium chemistry to be important for HR8799bcd.Skemer et al. (2014) report on disequilibrium chemistry for HR8799cd, based on NIR and MIR (narrow and broad band) photometry.Based on photometry, Currie et al. (2014) report the need for disequilibrium in planets b and c, but less strongly for d and e. Lavie et al. (2017) presented the first free retrieval analysis of the HR8799 planets and found that bc require disequilibrium chemistry, but not planets de.It is important to note, however, that no K-band spectroscopy was available for planet e, and that the K-band seemed heavily affected by systematics for planet d.Because the K band is important for detecting the presence of CO, and hence for detecting disequilibrium chemistry, their finding for d and e to be in chemical equilibrium has to be taken with caution.Using SPHERE and the same GPI K band data as in Lavie et al. (2017), Bonnefoy et al. (2016) report that HR 8799de can be fit well with chemical equilibrium models when using a non-scattering-cloud Exo-REM (Baudino et al. 2015) grid.With the new Exo-REM models by Charnay et al. (2018), which include chemical disequilibrium and scattering clouds, it is found that HR8799e requires disequilibrium to explain the GRAVITY K band data reported in Gravity Collaboration et al. (2019).In summary, while the consensus in the literature regarding HR8799e is not entirely clear, we corroborate the finding of Gravity Collaboration et al. (2019), which is based on a high-S/N spectrum in the K-band, that HR8799e is affected by disequilibrium chemistry.

Reported modeling challenges
When comparing models to observations of the HR 8799 planets a few problems have been identified.First, the models can converge toward high effective temperatures and thus small radii, in order to conserve the total flux of the planet within the model.Some radii that have been reported are smaller than expected from theoretical standpoints.Due to electron degeneracy pressure, the radii of gas giant planets and brown dwarfs are always above 0.75 R J and even this lowest limit is only reached after 10 Gyr of contraction, for brown dwarfs of 70 M J (Chabrier et al. 2009).For ages up to 5 Gyr, radii are above 1 R J for masses up to 30 M J (Mordasini et al. 2012).For objects with ages below 100 Myr, the minimum radius is therefore expected to be above 1 R J (Marley et al. 2012).A radius that is too small when compared to these theoretical constraints hints at shortcomings in the atmospheric model being used.Examples for such small inferred radii are found in Barman et al. (2011) (who report R = 0.75 R J for HR8799b), Greenbaum et al. (2018) (who report that some, but not all, of the model grids they tried have radii below 1 R J for HR8799cde) and Bonnefoy et al. (2016) (who reported the same for some but not for all of the models they applied to HR8799bcde).One approach is to reject such unphysical radii right from the start, by tying the spectral model fitting to evolutionary models, which guarantees that physically consistent radii are used for the spectral analysis (Marley et al. 2012).Because this can worsen the spectral fit and thus lead to additional biases concerning inferred properties of the atmosphere, it is questionable how much is gained from such an approach (Barman et al. 2015).In this regard, limiting the radius in our HR8799e retrieval to a minimum value of 0.9 R J can be seen as a compromise between these two approaches, and our best-fit radius of 1.12 0.09 −0.09R J is above the 1 R J limit described above.In a similar vein, it is useful to check whether the retrieved values of the planet's effective temperature, surface gravity and radius (hence also its luminosity and mass) are consistent with evolutionary models.Considering the evolutionary plots presented in Marley et al. (2012) (their figures 8 and 11) and following their analysis, places our median log(g) and T eff values between their 10 and 30 Myr isochrones, although our log(g) uncertainties also allow for ages in excess of 100 Myr.These values are consistent with the ages that can be inferred for the HR8799 system (between 30 and 60 Myr, see Marley et al. 2012, for a more complete discussion).Similarly, our median log(g) and T eff values lie between the evolutionary tracks of 5 and 10 M J models, and again our large uncertainties on log(g) allow for masses well below 5 M J and in excess of 10 M J .The mass we derive from HR8799e's spectrum, M P = 4.81 +8.78  −3.33 M J is certainly consistent with this assessment.We note here that Marley et al. (2012) made the assumption of a hot start evolution in their work, and that the HR8799 system is young enough for hot and cold start differences to play a role.It was found by Spiegel & Burrows (2012); Marley et al. (2012); Marleau & Cumming (2014), however, that at least a classical cold start assumption (Marley et al. 2007) is ruled out for these planets.In any case, recent theoretical modeling of the physics of the accretion shock (Marleau et al. 2017(Marleau et al. , 2019) ) and of the structure of accreting planets (Berardo et al. 2017;Berardo & Cumming 2017;Cumming et al. 2018) Table 4. Comparison of reported properties of HR8799e, derived from spectral/photometric analyses.The second Exo-REM2 retrieval from our study fixed the SPHERE/GPI scaling factors to the best-fit values of petitRADTRANS.References: B16 (Bonnefoy et al. 2016), L17 (Lavie et al. 2017), G18 (Greenbaum et al. 2018), C14 (Currie et al. 2014).Code references: Exo-REM2 (version of Charnay et al. 2018), Exo-REM1 (version of Baudino et al. 2015), Helios-r (Lavie et al. 2017), PHOENIX (models reported in Barman et al. 2011), COOLTLUSTY (Sudarsky et al. 2003;Hubeny et al. 2003;Burrows et al. 2006).Notes: (a) these publications compared their observations to more than one model, we report the best-fit values of the models that provide the best fit to the data, either as stated by the authors or by visual inspection.with models studying the orbital evolution of the HR 8799 system, where it was found that the mass of HR 8799e has to be below 7.6 M J to ensure orbital stability (Wang et al. 2018).
In addition, models that assume a homogeneous cloud cover have been reported to have trouble at reproducing all the MIR photometry simultaneously, especially the 3.3 µm, L' and [4.05] bands.The use of patchy cloud models, or models mixing clouds of different vertical extent, has been shown to be one promising way of solving this problem (see, e.g., Currie et al. 2011;Skemer et al. 2012;Currie et al. 2014;Skemer et al. 2014).With this in mind it is interesting to see that our high likelihood retrieval models reproduce both the 3.3 µm and L' band photometry, without including these data points in the retrieval.This would speak against a heterogeneous cloud cover being necessary to explain the data.A similar result was found by Bonnefoy et al. (2016) when fitting the HR 8799 planets with the Exo-REM code (Baudino et al. 2015) in the non-scattering, chemical equilibrium version.It is important to note, however, that neither Bonnefoy et al. (2016) nor we can reproduce the [4.05] band photometry of Currie et al. (2014), which is consistently higher than our best-fit spectra, and all the best-fit models presented in Bonnefoy et al. (2016), for all of the HR8799 planets.Thus, a heterogeneous cloud cover cannot be ruled out.Given the prevalence of variability of, especially, L-T dwarf transition objects (see, e.g., Apai et al. 2013;Crossfield 2014), such a heterogeneous cloud cover could be expected for the HR 8799 planets.

Quantitative comparison to the literature
Finally, we compare our retrieved parameter values for HR8799e with those reported in the literature.All values are listed in Table 4.Our derived gravity value, log(g) = 4.00 +0.46  −0.52 , falls within the values reported in the literature, which range from 3.5 to 4. The same holds for the effective temperature, for which we retrieve 1154 +49 −48 K.This value is bracketed by the reported values, ranging between 1000 and 1200 K. Similarly, our retrieved radius value (1.12 +0.09 −0.09R J ) falls between the reported values of 1 to 1.3 R J .Including ours, four out of six data-model comparison studies also varied the metallicity.Reported values are between 0.4 and 0.5, and our value ([Fe/H] = 0.48 +0.25  −0.29 ) is consistent with their assessments.C/O deserves a more detailed discussion.
Only one other study considered the atmospheric C/O.The retrieved values range from 0 to 0.60 +0.07 −0.08 , the latter being the value we derive with petitRADTRANS in this study.As discussed in Section 4.5, the Exo-REM fit runs into the grid boundaries for both C/O and [Fe/H], making it difficult for us to assess whether this trend to low C/O ratios is actually merited by the data or whether this parameter is used to copmensate for the too stringent cloud description.The value reported by Lavie et al. (2017) (C/O = 0) suffers from the fact that K-band spectroscopy was not available at the time of their study, such that the CO abundance could not be determined.They only report an upper limit for the atmospheric C/H value, and their best-fit spectrum does not show any CO absorption in the K-band, which we detect in the GRAVITY data.Given these caveats we conclude that a comparison to the C/O ratios derived in other studies is at this point inconclusive.

Summary
We present a new version of our retrieval radiative transfer code petitRADTRANS (Mollière et al. 2019) to which we added the effect of multiple scattering.This enables us to run free retrievals on cloudy self-luminous objects such as directly imaged planets and brown dwarfs.This updated version of petitRADTRANS will be available on the petitRADTRANS website soon10 and is already available now, upon request.
Running verification retrievals on synthetic observations, we found that we can retrieve the input parameters.The high dimensionality of the input model can lead to small offsets within the observational uncertainties, however.Increasing the number of live points in our retrievals with the nested sampling method improves this, but we expect this to be a persisting problem for models with a large number of free parameters, especially for high S/N observations which lead to narrow posterior distributions and thus let a smaller fraction of the prior volume be of interest.
We tested two different cloud models.The first is the physically motivated Ackerman & Marley (2001) cloud model.Our second cloud model simply retrieves the wavelength properties of the cloud opacities.When running retrievals with Cloud Model 2 on mock observations made with Cloud model 1, we find that we can get an excellent fit but with biased atmospheric parameters.In particular the planet's photosphere is found to be more isothermal and less cloudy than the input, which mimics the shallow temperature gradients predicted by Tremblin et al. (2015Tremblin et al. ( , 2016Tremblin et al. ( , 2017Tremblin et al. ( , 2019)).Thus, retrievals alone will likely not be enough to investigate whether shallow temperature gradients indeed occur in such thought-to-be-very-cloudy atmospheres.Retrieval analyses could be aided by longer wavelength data in the mid-IR, however, which could reveal the spectral features of cloud particles at ∼10 micron (Cushing et al. 2006).
We ran our retrieval setup on archival GPI, SPHERE, and partially new GRAVITY data for the directly imaged planet HR 8799e, using Cloud Model 1. Applying such data-driven, free retrievals for directly imaged planets becomes possible with our high S/N observations, especially GRAVITY's K band spectra at a spectral resolution of R = 500.In addition, observations in the K band can probe H 2 O, CH 4 and CO features and are therefore crucial for constraining atmospheric disequilibrium chemistry and the atmospheric C/O.
We are able to fit the observations well, and confirm a cloudy atmosphere dominated by disequilibrium chemistry.We also compare our retrieved atmospheric spectra to archival photometric observations in the L', 3.3 µm, [4.05] and M band.Our spectra are consistent with all points, except for the [4.05] band point, although the photometry has not been included in the retrieval.The L'−3.3 µm and L'−[4.05]band colors have been suggested to require heterogeneous cloud coverage to explain the data, and we cannot confirm this for the L'−3.3 µm color.
The posterior parameter values we retrieve for the atmospheric log(g), T eff and [Fe/H] are consistent with previous studies, and hot start evolutionary calculations.For the first time, we successfully constrain the C/O of HR 8799e and find that it is 0.60 +0.07 −0.08 , which is consistent with stellar.Running additional retrievals on HR 8799e with Cloud Model 2, we find that the retrieved planetary C/O, [Fe/H] and log(g) are identical to the values found with the nominal Cloud Model 1, therefore independent of our cloud model choice.This is noteworthy as the Cloud Model 2 retrievals lead to less cloudy, more isothermal atmospheres.This indicates that C/O may be a quite robust outcome of the retrievals.We also fit the HR 8799e spectrum with the self-consistent code Exo-REM, which uses a state-of-the-art onedimensional cloud model, scattering, and disequilibrium chemistry.Our free retrieval results compare well to the Exo-REM fit, except for the C/O ratio.With Exo-REM the C/O ratio is driven into the grid boundary (C/O < 0.3), as is the metallicity derived from Exo-REM ([Fe/H] > 0.5).This makes the comparison between the free retrieval with petitRADTRANS and Exo-REM difficult for C/O.A larger Exo-REM grid, which also varies the free parameters of its cloud prescription, may resolve this issue.
Using our retrieved C/O and metallicity, and a highly simplified disk model, we find that HR 8799e could have formed outside the CO iceline.This would imply that the planet migrated significantly.Because HR8799e is the innermost planet of the HR8799 system, this could indicate that all HR8799 planets formed outside of the CO iceline.Similar formation distances, relative to the icelines, have been theorized for Jupiter in the Solar System ( Öberg & Wordsworth 2019;Bosman et al. 2019).Using sophisticated gas-grain chemical modeling for the protoplanetary disk we find that the planet could also have formed more closely to the star, but outside the CO 2 iceline.This would require less migration.

Outlook
Here we introduced the first version of our retrieval framework for cloudy scattering atmospheres.With this we introduce a versatile tool for interpreting the spectra of directly imaged planets and brown dwarfs.At the same time, it is clear that there are many avenues for improving and testing our method.
One could explore different pressure-temperature parameterizations, cloud model setups, or abundance models.The latter could mean, for example, retrieving absorber abundances independently, as is often done in retrieval studies.Alternatively, the assumption of chemistry could be kept, while retrieving not just C/O (which we varied by changing O) and metallicity, but by retrieving C/H, O/H, and other atomic abundance ratios instead (as was done in, e.g., Lavie et al. 2017;Spake et al. 2019).Another important addition will be to include a parameterization for heterogeneous cloud coverage.
It is also crucial to test what happens when retrieving atmospheric parameters using a model setup that is different from the one used to generate a synthetic observation.In this way one can begin to quantify the uncertainties and biases of retrieved parameters given the model choices, and how robust certain parameters are against using a wrong model.An example is the robustness of C/O that we found in our results here, when using different cloud models.In principle, the most likely among a set of models can be found using the Bayes factor, computed with the model evidences derived from nested sampling.However, a given model may be worse than another in terms of model assumptions, which could be highly unphysical, while still being favored by a Bayes factor analysis.Thus, such comparison retrieval studies offer additional insight regarding the real parameter uncertainties, including the modeling choices.
On the observational side, additional data in the mid-IR, ideally spectroscopy, is necessary to explore the properties of clouds further.This is because the L', 3.3 µm and [4.05] band may encode information about a heterogeneous cloud cover.Excitingly, the outer planets (HR 8799bcd) will be studied in the mid-IR, using NIRSpec IFU spectroscopy with the JWST11 .HR 8799bcde will also be studied with photometry in the mid-IR with JWST12 .ESO's ERIS13 instrument, to be mounted at the VLT, is a promising option, as well as KPIC, to be mounted on Keck II (Mawet et al. 2018).ESO's imminent CRIRES+ instrument may offer the possibility to study the HR 8799 planets at high spectral resolution, where the S/N of the planetary flux measurement could be boosted using the cross-correlation method (e.g., Hoeijmakers et al. 2018), which also allows for carrying out retrievals (Brogi & Line 2019).Further in the future, the METIS spectrograph (Brandl et al. 2014) of ESO's upcoming ELT telescope as well as the PSI instrument (Skemer et al. 2018) on the TMT will be excellent instruments for mid-IR observations.For studying the potential absorption feature of silicate clouds at 10 micron, JWST will again be an excellent instrument.Here one approach could be to obtain mid-IR spectra of HR8799 planet analogs such as PSO J31814 .
Lastly, also constraining basic parameters of the HR 8799 planets better could be of great help.If astrometry were to give a mass estimate for the HR 8799 planets, a prior on the planet mass would lead to a better log(g) and therefore [Fe/H] inference during the retrievals (these two parameters being correlated).This in turn would also help to understand the planet's formation better, because the planetary metallicity can be regarded as a measure for the relative importance of the solid body accretion of a planet.

Appendix A: K-table mixing
This section describes our process for obtaining the total cumulative opacity distribution function, hereafter called k-table, of the atmosphere.The resulting k-table will contain the contribution of all absorbers.This combined k-table is required when including scattering during the computation of emission spectra.The combination is achieved by sampling the k-tables of the individual absorber species.

A.1. Gaussian quadrature grid definition
For every species, petitRADTRANS stores k-tables as a function of pressure and temperature.The spectral bins for which the ktables are stored have the width ∆λ.The width varies as a function of wavelength, because it is chosen such that λ/∆λ = 1000.Within these bins the opacity of each species is stored as a function of the cumulative probability g, where g = 0 denote the lowest opacity values, and g = 1 the highest values.These are the k-tables.An introduction to correlated-k, and why such a ktable treatment is useful to speed up calculations when compared to line-by-line calculations, is given in Irwin et al. (2008), their Section 2. Especially their Figure 1 illustrates the definition of g as the spectral coordinate, when compared to the wavelength.At low pressures, where the effect of pressure broadening is small, the k-tables of a line absorber will have a low-opacity tail extending over most of the g values.In addition, the line cores will give rise to a steep increase in opacity (by orders of magnitude), when approaching g values of unity.In order to sufficiently resolve the cumulative opacity distribution for low pressures, where the rise to highest opacity values will occur over a very narrow range at high g values, we split our g grid in two parts.The first is a Gaussian quadrature grid with coordinate g low extending from g low = 0 to g low = 0.9, while the second is a Gaussian grid extending from g high = 0.9 to g high = 1.Both grids have eight points, and their weights w have been rescaled such that This weight rescaling guarantees that, for any λand gdependent function f , ) .
(A.2)  curves of the individual species (dashed lines, same color coding as in the upper panel), the total k-table obtained from premixing the opacities of all species in wavelength space (black solid line), and total k-table obtained with the method used in this work (red solid line).

A.2. Sampling
For obtaining the total k-table, we use the standard assumption for the on-the-fly combination of the opacities of different absorbers (see, e.g., Lacis & Oinas 1991;Irwin et al. 2008;Mollière et al. 2015;Amundsen et al. 2017), namely that their opacity distribution functions are independent (also called 'random overlap').Making this assumption, an opacity value of the total k-table can be sampled by drawing opacity samples from the k-tables of each individual species, scaling them according to the respective abundances of the absorbers, and adding them.
By sampling the total k-table in this fashion many times, a good approximation of the total opacity distribution function can be constructed: the sampled values are simply sorted in magnitude, with the lowest value corresponding to g = 0, and the highest value corresponding to g = 1.This approximated k-table can then be interpolated back to the g low and g high values, to be ready for use in petitRADTRANS.
The sampling of a species' k-table could be done by drawing a uniformly distributed number between 0 and 1, and obtaining the opacity value by interpolating in the k-table of the species, to the so-drawn g value.However, this step needs to be done a number of times, and a numerically less costly option is to randomly selected an index numbering the opacities of the species' 16-point k-table, and treating the tabulated opacity at this index as the sampled value.Doing this for all species, and scaling the opacities by the abundances, yields again a sample of the total k-table, when the samples are added.
One subtlety is that not every index is equally likely, because we store the opacities on two eight-point Gaussian quadrature Table C.1.Priors of the "zoomed-in" verification retrieval of Cloud Model 1. Prior ranges marked with ( * ) have been changed when compared to the initial run, so as to "zoom in" on the high likelihood region identified by the initial retrieval.The symbols have the same meaning as in Table 3.
structure of petitCODE as the input for a petitRADTRANS calculation.This is done to verify the k-table mixing (see Appendix A) and scattering implementation (see Section 2.1) of petitRAD-TRANS.
The model we use as an input is broadly motivated by the properties of HR 8799e, with the following input parameters: T eff = 1200 K, log(g) = 3.75 (with g in cm s −2 ), [Fe/H] = 0, C/O = 0.55, f sed = 3, σ g = 2, K zz = 10 8.5 cm 2 s −1 , where f sed is the settling parameter as defined in Ackerman & Marley (2001), σ g is the width of the log-normal particle size distribution, and K zz is the atmospheric eddy diffusion coeffiecient.The rest of the symbols have their usual meaning.The clouds were assumed to consist of irregularly-shaped particles.For this we used cloud opacities calculated with the code by Min et al. (2005), which also makes uses of software by Toon & Ackerman (1981).The code assumes a distribution of hollow spheres (DHS) for the particles.Moreover, the condensates were assumed to be crystalline.
The upper panel of Figure A.2 shows the comparison of the resulting spectra of petitCODE and petitRADTRANS, which agree excellently.In addition, we also show the petitRADTRANS spectrum which would result from neglecting scattering.It is evident that scattering is an important process that needs to be accounted for.For completeness, the lower panel of Figure A.2 shows the self-consistent structures of the atmospheric temperature, cloud mass fraction and particle radius, resulting from the petitCODE calculation, which is used as an input for producing the spectrum with petitRADTRANS.

Appendix G: Setup of the ANDES disk model
Here we describe the setup of the ANDES disk model (Akimkin et al. 2013;Molyarova et al. 2017), used to study the impact of sophisticated disk chemical modeling on the inferences made about the formation location of planet HR 8799e, see Section 5.1.The gas surface density distribution was described by a tapered power law with an exponent γ = 1 and a characteristic radius R c = 100 au, the total gas mass in the disk was 0.1 M .The disk midplane temperature was calculated from the stellar and accretion luminosities.The thermal structure of the disk upper layers was calculated using ray tracing in the UV-optical wavelengths.The vertical disk density structure was derived from the hydrostatic equilibrium and was iteratively made consistent with the temperature structure.The gas and dust temperatures were assumed to be equal.The chemical model is based on the gas-grain code ALCHEMIC (Semenov & Wiebe 2011) with desorption energies updated according to Cuppen et al. (2017).An MRN-like power law size distribution (that is, following Mathis, Rumpl, & Nordsieck 1977) with the maximum dust size of 25 micron was used.The average dust radius of 0.37 micron and a dust-to-gas ratio of 0.01 were used for chemical simulations.The assumed stellar parameters representative of HR 8799  Posterior distribution of the follow-up retrieval with Cloud Model 1 (see Section 3).For this retrieval the uniform prior boundaries were chosen to enclose the regions of highest likelihood inferred in the original Cloud Model 1 retrieval, so as to test the effects of using the same number of live points in a smaller prior volume.The uniform prior boundaries are described in Appendix C.
were defined using an evolutionary track model for 1.47 Msun star (Yorke & Bodenheimer 2008).Accretion on the star was described by adding an accretion region with a temperature of 15000 K, contributing to accretion luminosity and UV radiation field.The assumed accretion rate was 10 −8 M /yr.Cosmic rays, X-rays and radioactive nuclides were included as additional sources of ionization.The chemical evolution of the disk was run for 1 Myr, starting from a pre-calculated composition of a 1 Myr old molecular cloud with "low metal" initial abundances (Lee et al. 1998).

2P.
Mollière et al.:  Scattering clouds and disequilibrium chemistry in the atmosphere of HR 8799e 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p o q x J Y x G r T o C a C S 5 Z 0 3 A j W C d R D K N A s H Y w v p 3 5 7 S e m N I / l o 5 k k z I 9 w K Hn I K R o r P V T x v F + u u D V 3 D r J K v J x U I E e j X / 7 q D W K a R k w a K l D r r u c m x s 9 Q G U 4 F m 5 Z 6 q W Y J 0 j E O W d d S i R H T f j Y / d U r O r D I g Y a x s S U P m 6 u + J D C O t J 1 F g O y M 0 I 7 3 s z c T / v G 5 q w m s / 4 z J J D Z N 0 s S h M B T E x m f 1 N B l w x a s T E E q S K 2 1 s J H a F C a m w 6 J R u C t / z y K m l d 1 D y 3 5 t 1 f V u o 3 e R x F O I F T q I I H V 1 C H O 2 h A E y g M 4 Rl e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > (b)

Fig. 1 .
Fig. 1. Results of the verification retrieval using Cloud model 1.Panel (a): synthetic observation, best-fit spectrum and residuals.Panel (b): emission contribution function.Due to the clouds, pressures larger than 1-2 bar cannot be probed.Panel (c): retrieved pressure-temperature confidence envelopes.The black dashed line shows the flux average of the emission contribution function that is shown in Panel (b).The opaqueness of the temperature uncertainty envelopes has been scaled by this contribution function, with a minimum value of 10 %.Panel (d): 2-d posterior plot of the (non-nuisance) retrieved atmospheric parameters.The red dashed lines denote the input values.The values of the cloud mass fractions at the cloud base have been divided by the mass fractions predicted when assuming equilibrium condensation at the cloud base location.
Figure 1.Panel (a) shows the synthetic observation, best-fit model, and the residuals between the two, scaled by the error bars.The residuals are flat and consistent with zero.The emission contribution function of the best-fit model is shown in Panel (b) of Figure 1.Regions between 0.004 and 2 bar are accessible, with the Fe and MgSiO 3 cloud blocking the flux from the deeper regions.Methane absorption blocks most of the flux longward of 2.2 µm, probing cool regions as high as 0.004 bar.Shortward of 2.2 µm most of the flux originates from a narrow pressure region from 0.2 to 2 bar.
Figure D.1 shows the corresponding results (see Section D).

Fig. 2 .
Fig.2.Spectral fit of HR 8799e.The upper panel shows the YJH-band observations of SPHERE and GPI, the middle panel the GRAVITY K-band observations.The lowest panel shows the photometry of the planet, which was not included during the retrieval.The 16-84 and 2-98 % flux envelopes of the sampled petitRADTRANS retrieval models are shown in all panels.Because also the SPHERE and GPI scaling factors were sampled 100 times for making this plot, there are multiple points visible at every wavelength.

Fig. 3 .
Fig. 3. Left panel: temperature distribution of the atmosphere of HR8799e, retrieved with the petitRADTRANS free retrieval setup.See the caption of Figure 1 for an explanation of how to read this plot.In addition the self-consistent P-T curves derived from petitCODE, assuming chemical equilibrium and no clouds, or chemical quenching with clouds, are shown as gray and black solid lines, respectively.Right panel: emission contribution function of the best-fit model of the HR8799e retrieval.

Fig. 4 .
Fig. 4. Corner plot of the free retrieval of HR 8799e with petitRADTRANS.This plot shows the one and two-dimensional projections of the 18-dimensional posterior distribution for all but the six nuisance parameters of the temperature structure (see the left panel of Figure3for the P-T confidence envelopes).M P,spec is the mass of HR 8799e in units of Jupiter masses, derived from the retrieval posterior distributions of log(g) and R P .

Fig. 5 .
Fig. 5. Marginalized one-dimensional posterior distributions of HR8799e's gravity, metallicity, and C/O, shown for our nominal retrieval with Cloud Model 1 (black solid line), and the retrieval with Cloud Model 2, which lead to a bi-modal solution of a cloudy (orange solid line) and ∼isothermal, clear (green solid line) atmospheric state.

PFig. 7 .
Fig. 7. One-dimensional marginalization of the HR 8799e retrieval posteriors for the self-consistent Exo-REM grid when retrieving the SPHERE/GPI flux scaling (gray), or fixing it to the median values of petitRADTRANS (purple).The results of the freely parameterized petitRADTRANS retrieval are shown in orange.The vertical dashed lines denote the upper [Fe/H], and lower C/O and T eff grid boundaries of the Exo-REM models.
(b) the best-fit model was trending into the boundaries of the Exo-REM grid, and only a boundary value can be reported for C/O.(c) same as (b), but for the [Fe/H] parameter.(d) as read of by eye from their corner plot.(e) not specified.(f) no K-band spectrum was available for their analysis, such that the CO abundance could not be constrained.Their retrieved C/O ratio is pushed against the 0 boundary.(g) as derived from their stated O/H ratio, in comparison to the solar O/H ratio.(h) consistent with radii derived from evolutionary models.

20P.
Mollière et al.:  Scattering clouds and disequilibrium chemistry in the atmosphere of HR 8799e

Fig
Fig. C.1.Posterior distribution of the follow-up retrieval with Cloud Model 1 (see Section 3).For this retrieval the uniform prior boundaries were chosen to enclose the regions of highest likelihood inferred in the original Cloud Model 1 retrieval, so as to test the effects of using the same number of live points in a smaller prior volume.The uniform prior boundaries are described in Appendix C.
t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h j L Z b t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p o q x J Y x G r T o C a C S 5 Z 0 3 A j W C d R D K N A s H Y w v p 3 5 7 S e m N I / l o 5 k k z I 9 w K H n I K R o r P V T x v F + u u D V 3 D r J K v J x U I E e j X / 7 q D W K a R k w a K l D r r u c m x s 9 Q G U 4 F m 5 Z 6 q W Y J 0 j E O W d d S i R H T f j Y / d U r O r D I g Y a x s S U P m 6 u + J D C O t J 1 F g O y M 0 I 7 3 s z c T / v G 5 q w m s / 4 z JJ D Z N 0 s S h M B T E x m f 1 N B l w x a s T E E q S K 2 1 s J H a F C a m w 6 J R u C t / z y K m l d 1 D y 3 5 t 1 f V u o 3 e R x F O I F T q I I H V 1 C H O 2 h A E y g M 4 R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h j L Z b t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p o q x J Y x G r T o C a C S 5 Z 0 3 A j W C d R D K N A s H Y w v p 3 5 7 S e m N I / l o 5 k k z I 9 w K H n I K R o r P V T x v F + u u D V 3 D r J K v J x U I E e j X / 7 q D W K a R k w a K l D r r u c m x s 9 Q G U 4 F m 5 Z 6 q W Y J 0 j E O W d d S i R H T f j Y / d U r O r D I g Y a x s S U P m 6 u + J D C O t J 1 F g O y M 0 I 7 3 s z c T / v G 5 q w m s / 4 z J J D Z N 0 s S h M B T E x m f 1 N B l w x a s T E E q S K 2 1 s J H a F C a m w 6 J R u C t / z y K m l d 1 D y 3 5 t 1 f V u o 3 e R x F O I F T q I I H V 1 C H O 2 h A E y g M 4 R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h j L Z b t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p o q x J Y x G r T o C a C S 5 Z 0 3 A j W C d R D K N A s H Y w v p 3 5 7 S e m N I / l o 5 k k z I 9 w K H n I K R o r P V T x v F + u u D V 3 D r J K v J x U I E e j X / 7 q D W K a R k w a K l D r r u c m x s 9 Q G U 4 F m 5 Z 6 q W Y J 0 j E O W d d S i R H T f j Y / d U r O r D I g Y a x s S U P m 6 u + J D C O t J 1 F g O y M 0 I 7 3 s z c T / v G 5 q w m s / 4 z J J D Z N 0 s S h M B T E x m f 1 N B l w x a s T E E q S K 2 1 s J H a F C a m w 6 J R u C t / z y K m l d 1 D y 3 5 t 1 f V u o 3 e R x F O I F T q I I H V 1 C H O 2 h A E y g M 4 R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O g m m F / u T R g U Z i C 8 8 x 4 R B t P B z J J E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h j L Z b t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p o q x J Y x G r T o C a C S 5 Z 0 3 A j W C d R D K N A s H Y w v p 3 5 7 S e m N I / l o 5 k k z I 9 w K H n I K R o r P V T x v F + u u D V 3 D r J K v J x U I E e j X / 7 q D W K a R k w a K l D r r u c m x s 9 Q G U 4 F m 5 Z 6 q W Y J 0 j E O W d d S i R H T f j Y / d U r O r D I g Y a x s S U P m 6 u + J D C O t J 1 F g O y M 0 I 7 3 s z c T / v G 5 q w m s / 4 z J J D Z N 0 s S h M B T E x m f 1 N B l w x a s T E E q S K 2 1 s J H a F C a m w 6 J R u C t / z y K m l d 1 D y 3 5 t 1 f V u o 3 e R x F O I F T q I I H V 1 C H O 2 h A E y g M 4 R l e 4 c 0 R z o v z 7 n w s W g t O P n M M f + B 8 / g C I q 4 1 K < / l a t e x i t > (b)< l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q + E k + 1 I s F e c Z y E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h r L Z T t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p Y t h k s Y h V J 6 A a B Z f Y N N w I 7 C Q K a R Q I b A f j 2 5 n f f k K l e S w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q+ E k + 1 I s F e c Z y E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h r L Z T t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p Y t h k s Y h V J 6 A a B Z f Y N N w I 7 C Q K a R Q I b A f j 2 5 n f f k K l e Sw f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < /l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q + E k + 1 I s F e c Z y E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h r L Z T t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p Y t h k s Y h V J 6 A a B Z f Y N N w I 7 C Q K a R Q I b A f j 2 5 n f f k K l e S w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n v m t g W E x 7 q w P d Q+ E k + 1 I s F e c Z y E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B a h X k o i g h 6 L X j x W t B / Q h r L Z T t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p Y t h k s Y h V J 6 A a B Z f Y N N w I 7 C Q K a R Q I bA f j 2 5 n f f k K l e S w f z S R B P 6 J D y U P O q L H S Q z U 4 7 5 c r b s 2 d g 6 w S L y c V y N H o l 7 9 6 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I z q w x I G C t b 0 p C 5 + n s i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l V d K 6 q H l u z b u / r N R v 8 j i K c A K n U A U P r q A O d 9 C A J j A Y w j O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A G K M I 1 L < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " w R j F 1 Y N S + i Q E o 0 8 0 j D U C r f P 6 d d 7 q A B T W A w g G d 4 h T d H O C / O u / O x a C 0 4 + c w x / I H z + Q O N O o 1 N < / l a t e x i t > Fig. D.1.Same as Figure1, but showing the results when a synthetic spectrum made with Cloud Model 1 is retrieved using Cloud Model 2. This test is described further in Section 3.2.

Fig. E. 1 .
Fig. E.1.Same as Figure1, but showing the results when a synthetic spectrum of the same wavelength spacing and noise properties as the actual HR8799e data is retrieved.Green posteriors indicate the retrieval where the scaling parameters of the GPI and SPHERE spectra were retrieved, while gray posteriors indicate the results that neglected the scaling parameters.This test is described further in Section 4.2.In Panel (a), the black points in the K-band correspond to synthetic observations with the same data quality as the GRAVITY observations reported in Gravity Collaboration et al. (2019), while the gray points show the data at the same quality as found for the two new GRAVITY observations presented in this work.All three synthetic GRAVITY data sets were fitted simultaneously.

Table 2 .
Log of the GRAVITY observations.
suggests that hot starts are more likely.The spectroscopic mass derived in our study is also consistent