HELIOSRetrieval: An Opensource, Nested Sampling Atmospheric Retrieval Code, Application to the HR 8799 Exoplanets and Inferred Constraints for Planet Formation
Abstract
We present an opensource retrieval code named HELIOSRETRIEVAL, designed to obtain chemical abundances and temperaturepressure profiles from inverting the measured spectra of exoplanetary atmospheres. We use an exact solution of the radiative transfer equation, in the pure absorption limit, in our forward model, which allows us to analytically integrate over all of the outgoing rays. Two chemistry models are considered: unconstrained chemistry and equilibrium chemistry (enforced via analytical formulae). The nested sampling algorithm allows us to formally implement Occam’s Razor based on a comparison of the Bayesian evidence between models. We perform a retrieval analysis on the measured spectra of the four HR 8799 directly imaged exoplanets. Chemical equilibrium is disfavored for HR 8799b and c. We find supersolar C/H and O/H values for the outer HR 8799b and c exoplanets, while the inner HR 8799d and e exoplanets have a range of C/H and O/H values. The C/O values range from being superstellar for HR 8799b to being consistent with stellar for HR 8799c and being substellar for HR8799d and e. If these retrieved properties are representative of the bulk compositions of the exoplanets, then they are inconsistent with formation via gravitational instability (without latetime accretion) and consistent with a core accretion scenario in which latetime accretion of ices occurred differently for the inner and outer exoplanets. For HR 8799e, we find that spectroscopy in the K band is crucial for constraining C/O and C/H. HELIOSRETRIEVAL is publicly available as part of the Exoclimes Simulation Platform (ESP; www.exoclime.org).
1 Introduction
1.1 Motivation
Traditionally, the masses and radii of brown dwarfs and substellar objects have been inferred from applying evolutionary tracks to measurements of their luminosities and ages (e.g., Burrows et al. 1997; Chabrier et al. 2000; Baraffe et al. 2002). On rare occasions, brown dwarfs and lowmass stars may transit their binary companions and allow for their other properties to be studied (see Burrows, Heng & Nampaisarn 2011 and references therein). A particularly important study was conducted by Konopacky et al. (2010), who were able to obtain dynamical masses for 15 brown dwarfs residing in binaries. By comparing the dynamical and photometric masses, Konopacky et al. (2010) showed that both the Burrows et al. (1997) and Chabrier et al. (2000) models underpredicted the masses of M and L dwarfs and overpredicted the mass of the lone T dwarf in their sample by (tens of percent). From studying a sample of 46 L dwarfs, Hiranaka et al. (2016) suggested that a dust haze of submicronsized particles exist in their upper atmospheres, which are neglected by the standard evolutionary tracks.
Taken together, these results suggest that the traditional approach of using selfconsistent evolutionary tracks may be incomplete and motivates alternative and complementary ways of interpreting the spectra of brown dwarfs and substellar objects. We expect this train of thought to apply to the recently discovered directly imaged exoplanets as well, since the interpretation of their photometry and spectroscopy is typically performed using the evolutionary tracks computed for brown dwarfs (e.g., Bonnefoy et al. 2016).
1.2 Theoretical Improvements
Selfconsistent forward modeling starts with a set of assumptions and computes forward to predict the temperaturepressure profile and synthetic spectrum of an object. Atmospheric retrieval is a complementary approach borrowed from the Earth remote sensing community, where one applies an inversion method to obtain the temperaturepressure profile and chemical abundances from finding the bestfit solution to the measured spectrum (e.g., Madhusudhan & Seager 2009; Benneke & Seager 2013; Lee, Heng & Irwin 2013; Line et al. 2013, 2016; Barstow et al. 2015; Waldmann et al. 2015). It sacrifices selfconsistency and sophistication for simplicity, which allows for a more thorough exploration of parameter space. Atmospheric retrieval is particularly wellsuited for addressing questions regarding planet formation, since it allows for the posterior distributions of the carbontooxygen ratio (C/O) and the elemental abundances of carbon (C/H) and oxygen (O/H) to be computed.
The HR 8799 system hosts four exoplanets (Marois et al., 2008, 2010), whose formation mechanisms remain an enigma (Kratter, MurrayClay & Youdin, 2010). Spectra with resolutions of about 30 to 4000 have been obtained by, e.g., Barman et al. (2011, 2015), Konopacky et al. (2013), Oppenheimer et al. (2013), Ingraham et al. (2014) and Zurlo et al. (2016). Since these spectra have resolutions that are considerably higher than those obtained for hot Jupiters using WFC3 on the Hubble Space Telescope (e.g., Deming et al. 2013; Mandell et al. 2013; Kreidberg et al. 2014; Stevenson et al. 2014), they present an opportunity for performing remote sensing of exoplanetary atmospheres that is similar to what planetary scientists had to work with a few decades ago, before the advent of probes. A key difference is that the radii and masses of these directly imaged exoplanets are unknown^{1}^{1}1Meaning they are typically not directly measured, but rather inferred using evolutionary models, which means the radii and masses are modeldependent., unlike for transiting exoplanets. A recent review of directly imaged exoplanets, which includes the HR 8799 system, may be found in Bowler (2016).
The first atmospheric retrieval analysis of directly imaged exoplanets was performed by Lee, Heng & Irwin (2013), who studied only the HR 8799b exoplanet. In the current study, we collect all of the published spectra of the HR 8799b, c, d and e exoplanets and subject them to the same retrieval method with the intention of using the retrieved chemistry to constrain planet formation scenarios.
Besides the novelty of our analysis, the current study is also a method paper for our new atmospheric retrieval code named HELIOSR, which we constructed from scratch to study exoplanetary atmospheres. HELIOSR is part of the HELIOS radiation package of the Exoclimes Simulation Platform^{2}^{2}2http://www.exoclime.org and has the following features (Figure 1):

We have implemented a nested sampling algorithm to explore the multidimensional parameter space (Skilling, 2006; Feroz, Hobson & Bridges, 2009; Benneke & Seager, 2013; Waldmann et al., 2015; Line et al., 2016). Unlike other approaches (e.g., Markov Chain Monte Carlo, nonlinear optimal estimation), nested sampling allows for the Bayesian evidence to be directly calculated, which in turn allows for models with different parametrizations (and number of parameters) to be compared on an equal footing. Models with extra complexity are penalized, which allows for Occam’s Razor^{3}^{3}3Whether Occam’s Razor always yields the correct answer is another matter. In the current study, we are guided by Occam’s Razor in the limit of sparse data. to be formally enforced. For example, our retrieval analysis allows us to formally determine if chemical equilibrium is favored or disfavored in an atmosphere in a completely datadriven manner. As another example, it allows us to determine the number and types of molecules to be included in the retrieval.

Our temperaturepressure profile is taken from Heng, Mendonça & Lee (2014), who generalized the work of Guillot (2010) and Heng et al. (2012) to include nonisotropic scattering and nonconstant opacities. When stellar irradiation and scattering are omitted, the temperaturepressure profile reduces to the classical solution of Milne for selfluminous objects (Mihalas, 1970). By construction, it conserves energy in an analytical and exact sense.

Our atmospheric cross sections are computed using our customized opacity calculator named HELIOSK, which was previously published by Grimm & Heng (2015).

To combine the cross sections of different molecules, one needs to have a chemistry model that calculates their relative abundances. We use the analytical solutions of Heng, Lyons & Tsai (2016), Heng & Lyons (2016) and Heng & Tsai (2016), which have been shown to be accurate at the level (or better) when benchmarked against numerical solutions using Gibbs free energy minimization. These analytical solutions allow for fast computation if one wishes to enforce chemical equilibrium.

Our radiative transfer scheme, which translates cross sections and temperatures into fluxes (and hence allows us to compute the synthetic spectrum), uses the exact analytical solution in the limit of isothermal model layers and pure absorption (Heng, Mendonça & Lee, 2014). It allows us to analytically integrate over all of the incoming and outgoing angles associated with every ray.

Our cloud model is based on the basic principles of Mie theory (e.g., Pierrehumbert 2010). It assumes a monodisperse set of particles, which may be interpreted as the dominant size in a size distribution of particles (e.g., Burrows, Heng & Nampaisarn 2011). It includes a dimensionless parameter that is a proxy for the cloud composition. When the particles are small compared to the wavelength, it reproduces Rayleigh scattering. By contrast, models that implement a constant cloudtop pressure implicitly assume the cloud particles to be large (compared to the wavelength observed) and preclude Rayleigh scattering by construction.
While each component of HELIOSR may not be novel by itself, the assembly of all of these components into a single code and retrieval tool is a novel endeavor. Furthermore, we have designed HELIOSR to run on Graphics Processing Units (GPUs), which affords speedups of at least a factor of several compared to the CPU version. With a UCrg model (see Table 1) retrieval performed on the HR 8799b dataset, the GPU version is 5 times faster than the CPU version on a Macbook Pro laptop equipped with a NVIDIA GeForce GT 750M GPU card and an Intel Core i7 2.5 GHz CPU. For this analysis, we used our GPU cluster of NVIDIA K20 cards; it takes seconds to evaluate one likelihood of this UCrg model.
In §2, we provide a detailed description of each component or ingredient of HELIOSR. In §3, we subject HELIOSR to several tests before applying it to the measured spectra of the HR 8799b, c, d and e directly imaged exoplanets. In §4, we present our retrieval results of the HR 8799 system. In §5, we compare our study to previous work and describe opportunities for future work. Table 1 shows the suite of models tested in the current study. Table 2 states the priors used for our fitting parameters. Table 3 summarizes our retrieval results. Appendix A states our fast analytical formulae for evaluating the exponential integral of the first order. Appendix B includes, for completeness, the full posterior distributions of the best models for the atmospheres of HR 8799b, c, d and e.
Notation  Meaning 

U  Unconstrained chemistry 
E  Equilibrium chemistry 
B  Cloudfree (“blue sky”) 
C  Cloudy 
1  HO is included in retrieval 
2  CO is included in retrieval 
5  CO is included in retrieval 
6  CH is included in retrieval 
r  planet radius is included in retrieval 
g  planet surface gravity is included in retrieval 
d  distance of the system is included in retrieval 
Note: “1”, “2”, “5” and “6” refer to the HITRAN/HITEMP labels for these molecules. When no number is specified, it means that all four molecules are included in the retrieval. Example: UBrg16 is a cloudfree model with unconstrained chemistry, where the mixing ratios of water and methane, as well as the planetary radius and surface gravity, are included as fitting parameters. By contrast, the UB model includes all four molecules in the retrieval, but fixes the planetary radius and surface gravity to userspecified values.
Symbol  Prior Used  Value 

Gaussian  
Gaussian  (cgs)  
Loguniform  to  
Loguniform  to 10 (mks)  
Uniform  10 to 1500 K  
Uniform  1 to 100  
Loguniform  to m  
Loguniform  to  
Gaussian  pc 
cgs: centimeters, grams and seconds.
mks: meters, kilograms and seconds.
2 Methodology
The executive summary is that each model of the retrieval contains up to 11 parameters: the radius, the surface gravity, 2 for the temperaturepressure profile, 2 or 4 for the chemistry (depending on whether one adopts equilibrium or unconstrained chemistry) and 3 for the cloud model. The mean molecular weight is not a parameter and is constructed from the mixing ratios. Each HR 8799 exoplanet typically has between 40 and 120 data points for its measured spectrum: 68 for b, 105 for c, 115 for d and 48 for e.
To construct an atmospheric retrieval model, we need a forward model. By “forward model”, we refer to the temperaturepressure profile, atmospheric opacities, chemistry model, radiative transfer scheme and cloud model. We also need a method to scan the vast multidimensional parameter space of our forward model to locate the highest likelihood region, i.e., the best solution that fits the data (e.g., for a review, see Press et al. 2007).
2.1 Nested Sampling
We use a nested sampling algorithm (Skilling, 2006) to scan the diverse, multidimensional parameter space describing our onedimensional model atmospheres. Benneke & Seager (2013) previously gave a detailed overview of the nested sampling method. Waldmann et al. (2015) and Line et al. (2016) also used nested sampling. Here, we provide a concise description of our implementation.
Consider a model with a set of parameters , where is the number of parameters. Consider a set of models labeled by the index : . The probability density function (PDF) on the parameters for a given model is , which is also known as the “prior”.
Discussions of any Bayesian method necessarily start with Bayes’s rule, which states that the PDF of a model given the data (denoted by ) is (e.g., Skilling 2006)
(1) 
The quantity is the “likelihood”. We assume to be the same Gaussian function as equation (5) of Benneke & Seager (2013).
We will term the “posterior”. Since it normalizes to unity, the Bayesian evidence is given by the multidimensional integral,
(2) 
Fitting a model to a measured spectrum is an exercise in which a better fit is obtained when more free parameters (e.g., more molecules) are introduced. Model selection is essentially the enforcing of Occam’s Razor, meaning that we select the model that has a level of sophistication or complexity that is commensurate with the quality of data available. It prevents the overfitting of data by a model that is too complex. For example, Hansen, Schwartz & Cowan (2014) find that for some of the exoplanets the photometric data of Spitzer alone may be fitted with a Planck function and a more complex model is unnecessary. As the data quality improves, so does the complexity of the best model.
The essence of nested sampling is to reduce the computation of the Bayesian evidence to a onedimensional integral (Skilling, 2006),
(3) 
where the likelihood now only depends on a single variable and is denoted by . This variable is termed the “prior mass” and is bounded between 0 and 1. A visualization of the prior mass and its relationship to the Bayesian evidence is given in Figure 3 of Skilling (2006) and Figure 1 of Benneke & Seager (2013). Numerically, we use the trapezoid rule to compute the Bayesian evidence as a finite sum,
(4) 
We begin by randomly drawing points from the parameter space() subjected to the constraint of the chosen prior. We use either Gaussian (radius, logarithm of gravity, distance), loguniform (mixing ratios, mean opacity, cloud particle radius, cloud mixing ratio) or uniform (temperature, cloud composition parameter) priors. For a set of points drawn, we compute their likelihood values. At each step of the algorithm, we discard the worst point and replace it with a newly drawn point until the convergence criteria is met (see Skilling 2006). This newly drawn point needs to have a higher likelihood than the worst point just discarded. Specifically, we use the opensource software named PyMultiNest^{4}^{4}4https://github.com/JohannesBuchner/PyMultiNest/ (Buchner et al., 2014), which is a Python wrapper for the opensource MultiNest^{5}^{5}5https://ccpforge.cse.rl.ac.uk/gf/project/multinest/ program written in Fortran 90 (Feroz & Hobson, 2008; Feroz, Hobson & Bridges, 2009; Feroz et al., 2013). For each model, we run the nested sampling algorithm using 40 000 living points parallelised into 20 runs of 2000 “living points” each. For comparison, Waldmann et al. (2015) uses living points. Benneke & Seager (2013) use between and 10,000 living points. Line et al. (2016) do not specify the number of living points used. Equation (4) is used to compute the Bayesian evidence. As a byproduct of this procedure, one also obtains posteriordistribution samples of the model parameters.
For the purpose of comparing two models, which we denote by and , it is useful to define a quantity known as the Bayes factor, which is the ratio of the Bayesian evidences (Trotta, 2008),
(5) 
The Bayes factor is equal to the posterior odds when both models are considered equally likely. As shown in Table 2 of Trotta (2008), which is reproduced in Table 2 of Benneke & Seager (2013), there is a relationship between the Bayes factor, the value of the frequentists and the significance in terms of the number of standard deviations. We use the Jeffreys scale (Kass & Raftery, 1995) to evaluate model significances. Weak, moderate and strong evidence for favoring the th model over the th model correspond to , 2.5 and 5, respectively.
2.2 TemperaturePressure Profile
For the temperaturepressure profile, we assume a onedimensional, planeparallel model atmosphere. Its layers are evenly spaced in the logarithm of pressure between 1 bar and 1 kbar. We implement equation (126) of Heng, Mendonça & Lee (2014), who previously generalized the work of Guillot (2010) (pure absorption limit and constant opacities) and Heng et al. (2012) (isotropic scattering, constant shortwave/optical opacity) to include nonisotropic scattering and a nonconstant shortwave/optical opacity. Since the HR 8799 exoplanets are nonirradiated, we essentially use a reduced version of equation (126) of Heng, Mendonça & Lee (2014),
(6) 
where is the internal/interior temperature, is the longwave/infrared scattering parameter, is the constant component of the longwave/infrared opacity and is the opacity associated with collisioninduced absorption (CIA). The column mass is denoted by , while is the column mass referenced to the bottom of the model atmosphere. We set kbar, where is the surface gravity.
Equation (6) is essentially a generalization of the classical Milne’s solution (Mihalas, 1970) to include scattering and CIA. In the limit of pure absorption () and in the absence of CIA, we obtain when , somewhat different from the classical Milne value of 2/3. It is worth emphasizing that equation (6) is, by construction, a temperaturepressure profile in radiative equilibrium, which implies that both local and global energy conservation are guaranteed in an exact, analytical sense (Heng, Mendonça & Lee, 2014; Heng & Lyons, 2016). By contrast, the versatile fitting function used by Madhusudhan & Seager (2009) does not, by construction, obey energy conservation and this has to be enforced as a separate numerical condition. However, in using a mean opacity equation (6) sacrifices accuracy for simplicity, which makes the temperaturepressure profile more isothermal, at high altitudes, than if a more realistic radiative transfer calculation was performed.
In principle, and are mean opacities that may be calculated directly from the spectroscopic line lists. However, in deriving these analytical temperaturepressure profiles Guillot (2010), Heng et al. (2012) and Heng, Mendonça & Lee (2014) have assumed that the absorption, flux, Planck and Rosseland mean opacities are equal, which makes it unclear how to exactly compute and . Therefore, we opt to use and as fitting parameters instead. In other words, our temperaturepressure profile is not selfconsistent with the atmospheric opacities used.
We find that using and as fitting parameters have a negligible effect on our results (not shown). In practice, the use of equation (6) with only and as fitting parameters (i.e., setting and ) is sufficient for our retrieval calculations.
We use a constant value of the surface gravity, as we are sensing at most 6 orders of magnitude in pressure, which corresponds to 13.8 scale heights. This means that the region of the atmosphere being sensed is only several percent of the radius of the exoplanet. A constant surface gravity is thus not unreasonable.
2.3 Atmospheric Cross Sections
Property  Exoplanet  Value 
HR8799b  
HR8799c  
HR8799d  
HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  
C/O 
HR8799b  
C/O  HR8799c  
C/O  HR8799d  
C/O  HR8799e  
C/H 
HR8799b  
C/H  HR8799c  
C/H  HR8799d  
C/H  HR8799e  
O/H 
HR8799b  
O/H  HR8799c  
O/H  HR8799d  
O/H  HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  

HR8799b  
HR8799c  
HR8799d  
HR8799e  
d 
HR8799b  
d  HR8799c  
d  HR8799d  
d  HR8799e 
Note: we have listed the uncertainties, which were computed by locating the 15.87th and 84.13th percentile points on the horizontal axis. In the limit of a symmetric Gaussian function, these would yield the fullwidth at halfmaximum of the Gaussian. For planet d and e, the molecules abundances and the mean molecular weight are given at 1 bar. Values are in (except for C/O) and dimensionless (except for which is in meters).
We first distinguish between our use of the terms “cross section” and “opacity”. The former has units of area. The latter is the cross section per unit mass. We previously designed and wrote an opensource opacity calculator (Grimm & Heng, 2015), based on implementing Algorithm 916 (Zaghloul & Ali, 2012) to perform fast computations of the Voigt profile by recasting it as a Faddeeva function. Typically, HELIOSK is able to compute an opacity or cross section function with spectral lines in s on a NVIDIA K20 GPU. In principle, it is agnostic about the spectroscopic line list being used and is able to take any line list as an input. The details of how to take the inputs of a line list and use them to compute the integrated line strengths and line shapes have previously been summarized in Grimm & Heng (2015) and we will not repeat them here.
We restrict ourselves to only four molecules: carbon monoxide (CO), carbon dioxide (CO), water (HO) and methane (CH). For CO and CO, we use the HITEMP database (Rothman et al., 2010). For HO and CH, we use the ExoMol line list (Barber et al., 2006; Yurchenko & Tennyson, 2014). Acetylene, ammonia, ethylene and hydrogen cyanide have been omitted, because they are subdominant at the photospheric temperatures of the HR 8799 exoplanets (Madhusudhan, 2012; Heng & Tsai, 2016; Moses et al., 2016). In particular, see Figure 10 of Moses et al. (2016).
In the current study, we choose to deal with cross sections instead of opacities. For our HELIOS selfconsistent radiative transfer code, we chose to use opacities instead (Malik et al., 2017). There are various strategies to construct the cross section function of the atmosphere. By “cross section function”, we refer to the function that depends on temperature, pressure, wavenumber and type of molecule. The cross section function is a theoretical construction: it may be defined continuously or be sampled at an arbitrary number of discrete points. We consider the way in which the cross section function is sampled as an issue of implementation, which we will now discuss. Regardless of the approach used to construct and sample the cross section function, the end goal is the same: to use them to construct transmission functions and ultimately integrate fluxes over a waveband.
The first approach is to use the “kdistribution method”, which resamples the highly erratic cross section function into a monotonically increasing cumulative distribution function (Lacis & Oinas, 1991; Fu & Liou, 1992; Grimm & Heng, 2015). Since the kdistribution method is only exact for a homogeneous atmosphere with one molecule (Grimm & Heng, 2015), one has to apply the “correlatedk approximation” as well, which assumes that the spectral lines are perfectly correlated (see Chapter 4.4.5 of Pierrehumbert 2010).
The second approach is to use “opacity sampling”, which is to discretely sample the opacity function, typically at a smaller number of points than there are lines. In our context, it is perhaps more accurate to term it “cross section sampling”.
The “linebyline” limit occurs when the integrated fluxes over a waveband is exact (to machine precision). It is essentially the second approach, but where the cross section function is sampled at more wavenumber points than there are lines. Since there are (or more) lines for the water molecule alone, this is a formidable computational challenge and is currently infeasible for any retrieval code dealing with hot exoplanetary atmospheres. We note that a cross section function that includes all of the lines of a given line list does not qualify it as being “line by line”, if the sampling is not fine enough to resolve each line profile.
In the current study, we adopt the second approach, which is also used by Madhusudhan & Seager (2009), Benneke & Seager (2013), Line et al. (2013) and Waldmann et al. (2015). Our spectral resolution used is 1 cm, evenly sampled across wavenumber. We note that Line et al. (2013, 2015) and Waldmann et al. (2015) also used a spectral resolution of 1 cm. Some authors do not specify the spectral resolution of their atmospheric cross section function (e.g., Madhusudhan & Seager 2009; Benneke & Seager 2013; Barstow et al. 2015; Line et al. 2016). We precompute our cross sections on a grid across wavenumber, pressure and temperature: 100 to 2900 K (in increments of 200 K) and 1 bar to 1 kbar (with two points per dex in pressure) for CO, CO, CH and HO. The grid is then interpolated to obtain values of the cross sections for any temperature and pressure within the stated ranges.
A lingering issue, which stems from an unsolved physics problem, is that the far line wings of Voigt profiles do not accurately represent the wings of real lines. Various groups have adopted different ad hoc approaches to truncating the Voigt profiles (see Grimm & Heng 2015 and references therein). Hedges & Madhusudhan (2016) discuss this issue, but do not provide any solution for it. In the current study, we adopt a 100 cm cutoff.
2.4 Chemistry
Once the cross sections have been computed, they may be used to compute the optical depth, of each model layer, for all of the molecules,
(7) 
where and are the mixing ratio and cross section of the th molecule, respectively. is the thickness of the layer in terms of the difference in pressure. The mean molecular mass is given by , where is the mean molecular weight and is the atomic mass unit. The preceding expression assumes hydrostatic equilibrium, isothermal layers and that the surface gravity is constant throughout our model atmosphere.
Generally, the mixing ratios of molecular hydrogen and helium only show up via CIA as a contribution to the continuum of a spectrum, which implies that they cannot be as definitively determined as that of the molecules. Our CIA opacities are obtained from Richard et al. (2012). In the range of –0.9, the effect on the continuum of the synthetic spectrum is very similar (not shown). The effects of HHe CIA are even more subtle. As such, we adjust to render the sum of the mixing ratios unity,
(8) 
where we have assumed that to reflect cosmic abundance. By denoting the mass of the th molecule by , the corresponding mean molecular weight is calculated using
(9) 
For example, if we have , and , then we have . In models with equilibrium chemistry, the mean molecular weight changes slightly for each layer, because the molecular abundances vary from layer to layer even for the same metallicity.
In the current study, we consider two chemistry models. “Unconstrained chemistry” refers to using each as a fitting parameter in the retrieval. “Equilibrium chemistry” means that the may be determined using only the elemental abundances of carbon () and oxygen (), if CHO gaseous chemistry is considered. In this case, the 4parameter system of unconstrained chemistry reduces to 2 parameters. To compute the four values given and , we use the validated analytical formulae of Heng, Lyons & Tsai (2016), Heng & Lyons (2016) and Heng & Tsai (2016). Specifically, we implement equations (12), (20) and (21) of Heng & Lyons (2016) for gaseous CHO chemistry. The benchmarking of these formulae against calculations of Gibbs free energy minimization was previously performed by Heng & Tsai (2016), who showed that they are accurate at the level or better. Further validation of these formulae comes from matching the trends found by Madhusudhan (2012) and Moses et al. (2013).
For unconstrained chemistry, the carbontooxygen ratio is computed using
(10) 
The elemental abundances are inferred using
(11) 
Each mixing ratio is assumed to be constant over the entire model atmosphere. The alternative, which is to have a different value of the mixing ratio for each of the 100 model layers we assume, would result in 400 free parameters. This is unwarranted given the sparseness of the data, i.e., we have less than 400 data points.
For equilibrium chemistry, the carbontooxygen ratio is simply
(12) 
and and are directly the fitting parameters of the retrieval. Since the mixing ratios of all of the molecules can be exactly specified for each layer, which has its own temperature and pressure, the assumption of constant mixing ratios across height/pressure is unnecessary for the models with equilibrium chemistry. The 400 values of the mixing ratios are specified by just two free parameters: and .
Unlike in previous studies, we do not manually decide whether to pick unconstrained or equilibrium chemistry. Rather, we compute both of these models and select between them based on the Bayesian evidence computed.
We note that, as part of the ESP, we have previously developed a chemical kinetics solver named VULCAN (Tsai et al., 2017).
2.5 Radiative Transfer Scheme
With the cross sections and temperaturepressure profiles in hand, one may compute the optical depth and hence the transmission function for each layer of the model atmosphere. To propagate fluxes through the atmosphere and thus obtain the synthetic spectrum, we need a radiative transfer scheme. Beer’s law^{6}^{6}6Also known as the BeerLambertBouguer law. is the simplest example of such a scheme, where incident radiation through a passive medium is exponentially attenuated. A more sophisticated radiative transfer scheme needs to account for both the fluxes incident upon a layer and the thermal emission associated with the layer itself, since each layer has a finite temperature. To this end, we use equation (B4) of Heng, Mendonça & Lee (2014),
(13) 
where the fluxes are computed at the  and th interfaces. The Planck function () is evaluated within each layer. The transmission function is given by equation (B5) of Heng, Mendonça & Lee (2014)
(14) 
with being the exponential integral of the first order. The optical thickness of each layer is given by . Appendix A describes an analytical fitting formula for that is highly accurate and allows for the computation to be significantly sped up.
We use equation (13) to propagate the boundary condition at the bottom of the atmosphere (i.e., the internal/interior heat flux), which is the Planck function with a temperature given by the temperaturepressure profile at the bottom boundary. The outgoing flux at the top of the atmosphere is then the synthetic spectrum.
We emphasize equation (13) is an exact solution of the radiative transfer solution in the limit of isothermal layers and pure absorption. It is an improvement over using approximate solutions (e.g., twostream solutions) and allows us to implement a radiative transfer scheme without taking any approximations besides assuming pure absorption. Equation (13) is equivalent to the approach of Line et al. (2016), who used fourpoint Gaussian quadrature to account for angledependent flux propagation. In our exact solution, the integration over angle has been performed analytically and is encapsulated in the exponential integral of the first order. We gain computational efficiency both by bypassing the need for performing Gaussian quadrature and also by evaluating using an analytical approach (Appendix A). The overall accuracy is relegated to the number of discrete layers used.
The radius of the exoplanet () only appears as a scaling factor between the observed flux () and the flux escaping from the top of the atmosphere (),
(15) 
where is the distance between the observer and the object. The HR 8799 system is located at pc (van Leeuwen, 2007), but the measured fluxes are usually reported as if it were located at pc (i.e., absolute fluxes).
2.6 Cloud Model
The need for a cloud model is motivated by previous suggestions that the atmospheres of the HR 8799 exoplanets are cloudy (Barman et al., 2011; Madhusudhan, Burrows & Currie, 2011; Marley et al., 2012), and also by the finding that each cloud configuration essentially corresponds to a different massradius relationship (Burrows, Heng & Nampaisarn, 2011; Lee, Heng & Irwin, 2013). Our cloud model is based on the notion that, while cloud formation is challenging to model from first principles (e.g., Helling & Woitke 2006), once clouds do form it is somewhat easier to describe their effects on the synthetic spectrum, since this derives from our knowledge of classical optics and Mie theory (Pierrehumbert, 2010).
Following Lee, Heng & Irwin (2013), we consider the presence of clouds to add an extra contribution to the optical depth,
(16) 
where is the extinction efficiency, is the radius of the (spherical) particles, is the number density of clouds and is the spatial thickness of the layer. The cloud mixing ratio is and it is this quantity that we set a prior on (see Table 2 for its range of values). We assume the cloud to be uniformly distributed throughout the atmosphere.
In a departure from the approach of Lee, Heng & Irwin (2013), we do not use a specific composition of cloud (e.g., enstatite). Specifically, we adopt their approximate fitting formula (listed in the appendix of Lee, Heng & Irwin 2013 but not used in their analysis),
(17) 
where and is the wavelength. When the particles are small (), we recover Rayleigh scattering: . Large particles () produce a roughly constant . By contrast, Benneke & Seager (2013) assume their clouds to be described by only one number, which is the cloudtop pressure. Their model carries the implicit assumption that the cloud particles are large compared to the range of wavelengths examined.
The dimensionless quantity serves as a proxy for the cloud composition. Refractory species (e.g., silicates) have , while volatile species (e.g., ammonia, methane, water) have –80. By using as a fitting parameter in the retrieval, we can constrain the composition of the clouds. The other fitting parameters in our cloud model are and .
Since we do not selfconsistently treat the cloud physics and gaseous chemistry, the caveat is that our retrieved C/O values are representative of only the gaseous component of the atmosphere. It is conceivable that the true C/O values, which must account for the material sequestered in the cloud particles, are different.
2.7 Data Selection: Spectra of HR 8799b, c, d and e
The spectra and photometric data points of the HR 8799b, c, d and e exoplanets have been taken from Bonnefoy et al. (2016) and Zurlo et al. (2016). The new SPHERE data were presented in Zurlo et al. (2016), while Bonnefoy et al. (2016) unified all of the previous data of the four exoplanets. Specifically, we use the data from Figure 4 of Bonnefoy et al. (2016).
To compute the flux in a photometric waveband, we simply integrate the synthetic spectrum over the range of wavelengths of the filter and assume a Heaviside function with a value of unity throughout. Unlike Lee, Heng & Irwin (2013), we do not apply filter functions with nonunity values to our synthetic spectrum, because this correction has already been done en route to reporting the observed fluxes in Bonnefoy et al. (2016) and Zurlo et al. (2016). It is unclear what has been done in previous studies. Madhusudhan, Burrows & Currie (2011) display filter functions in their Figure 1, but do not describe whether these filter functions were applied to their synthetic spectra. Line et al. (2013) state that, “For the broadband points we simply integrate the flux from the highresolution model spectrum with the appropriate filter function for that point,” but do not provide quantitative descriptions of their filter functions. It cannot be ruled out that these filter functions have values of unity throughout. For the spectroscopic data points, we do not convolve the synthetic spectrum with the instrument’s response function, because the impact is minor for lowresolution spectra.
3 Tests
Before analyzing the measured spectra of the HR 8799b, c, d and e directly imaged exoplanets, we subject HELIOSR to a battery of tests.
3.1 Number of Atmospheric Layers
The number of layers used in a onedimensional model atmosphere is a critical but often overlooked or unexplored detail. We wish to quantify the mean and maximum errors associated with assuming a specific number of model layers. We use the measured spectrum of HR 8799b as an illustration. We consider an ensemble of cloudfree models with unconstrained chemistry. For each model, we randomly select our parameter values: 2 parameters for the temperaturepressure profile, 4 parameters for the mixing ratios and 1 for the surface gravity. The range of parameter values used is listed in Table 2. No model selection is performed for this test. We consider forward models with both isothermal and nonisothermal layers. For the latter, we use equation (B6) of Heng, Mendonça & Lee (2014).
For each of the models, the spectrum computed with 10,000 nonisothermal layers is used as a reference. We then compute coarser models with between 10 and 8000 isothermal or nonisothermal layers and calculate the fractional error on the synthetic spectrum compared to the reference model. In Figure 2, we show both the mean and maximum errors associated with the synthetic spectrum. With 100 layers, we see that models with isothermal and nonisothermal layers have the same mean and maximum errors of about 2.5% and 8%, respectively. For the rest of the paper, we will use 100 isothermal layers. For comparison, Madhusudhan & Seager (2009), Lee, Heng & Irwin (2013) and Line et al. (2013) used 100, 43 and 90 layers, respectively.
3.2 Validating the Forward Model
We previously developed a selfconsistent radiative transfer code named HELIOS, which solves the radiative transfer equation in tandem with the first law of thermodynamics to obtain onedimensional model atmospheres in radiative equilibrium (Malik et al., 2017). HELIOS was validated against the radiative transfer model of MillerRicci & Fortney (2010). In the limit of pure absorption, we also demonstrated that the twostream and exact solutions produce excellent agreement if the diffusivity factor is set to 2 (Heng, Mendonça & Lee, 2014; Malik et al., 2017).
The forward model of HELIOS uses the same equation as HELIOSR, but was implemented independently by the first author of each study. Here, we compare the forward models of HELIOSR and HELIOS to verify that our implementation is bugfree. In Figure 3, we constructed a cross section function consisting purely of water and used the kdistribution method to compute the fluxes. Malik et al. (2017) used instead an opacity function, but also consisting purely of water and using the kdistribution method as well. The kdistribution tables were constructed using a resolution of cm evenly distributed across wavenumber (not shown). Other assumptions include a hydrogendominated atmosphere (), a water mixing ratio of and a surface gravity of in cgs units ( 19.5 m s). We then assumed an input temperaturepressure profile, as shown in the insert of Figure 3, in tandem with the kdistribution tables to compute the synthetic spectrum using both HELIOSR and HELIOS. The excellent agreement validates our implementation of the forward model.
3.3 Retrieval on a Mock Dataset
A useful test is to create a mock dataset, where we know what the “ground truth” is concerning the synthetic spectrum, temperaturepressure profile, molecular abundances, surface gravity, etc. We assume a cloudfree model with unconstrained chemistry, which has the following input parameters:
(18) 
where is the radius of Jupiter. Using this setup, we create 3 mock datasets: a full mock spectrum from 0.7–5m with 0.01 m resolution, HR8799blike and HR8799elike data coverage. We assume this mock object to be located pc away.
Such a test serves three purposes. First, if and are fixed to their input values (and excluded from being fitting parameters in the retrieval), then it is a test of the ability of our nested sampling algorithm to correctly recover the molecular abundances and temperaturepressure profile. Second, if we now include and as fitting parameters, it allows us to study the degeneracies associated with our ignorance of the surface gravity and/or radius. Third, by adapting and degrading the mock spectrum to the data resolution and spectral coverage of HR 8799b and HR 8799e, we may study the effects of incomplete or sparse data on the retrieved molecular abundances. The key difference between the currently available data for HR 8799b and HR 8799e is that the latter does not have K band spectroscopy.
Figure 4 shows the outcomes of these tests. When and are fixed to their input values, HELIOSR correctly recovers the input values of the mixing ratios and  profile parameters from the full mock spectrum (first row, first column). Surprisingly, our ability to recover these input values appear to be insensitive to whether the mock spectrum is degraded or not (first row, second and third columns), if and are known.
When the radius is implemented as a uniform prior, its value is correctly recovered, although the posterior distributions of the other fitting parameters become a little broader (third row of Figure 4). With HR 8799elike data coverage, we see clear signs of degeneracies being introduced into the posterior distributions. It suggests that the K band spectrum contains important information on the molecular abundances, an issue we will explore further in §4.4.
Allowing the surface gravity to be a fitting parameter has more serious consequences, as it introduces degeneracies into all of the other fitting parameters (second row of Figure 4). Even full data coverage does not lift these degeneracies (second row, first column). It suggests that an informative prior needs to be set on the surface gravity.
Surprisingly, the retrieved posterior distribution of C/O appears to be robust to the different model assumptions (Figure 5). It suggests that the C/O is a robust outcome of the retrieval.
Overall, these exercises teach us that wavelength coverage and spectral resolution are generally not as important as knowledge of the surface gravity, although the K band spectrum appears to encode crucial information on the molecular abundances. In §4.1, we will argue for setting Gaussian priors on as well as when analyzing real data from the HR 8799 exoplanets.
4 Results
4.1 Setting Priors on Radius and Surface Gravity
The strongest demonstration of why our assumptions for the prior distributions of input parameters are important comes from examining a model where the radius and surface gravity are implemented as uniform priors in the retrieval. Specifically, we perform a retrieval on the measured spectrum of HR 8799b using model UBrg in Figure 6, where and are specified as uniform priors. We see that the retrieved solution is , which is physically unreasonable. The surface gravity takes on unphysical values of –6. As we have learned from the mockretrieval exercises in §3.3, these difficulties stem from specifying the radius and surface gravity as unconstrained fitting parameters.
We now discuss why the values for are physically unreasonable. There are indirect arguments for why retrieved solutions with radii well below a Jupiter radius should be rejected. First, brown dwarfs and lowmass stars with masses between 20 and 100 Jupiter masses have transit radii that are at least (see Burrows, Heng & Nampaisarn 2011 and references therein), including CoRoT3b, which is a lowmass brown dwarf with a dynamical mass of and a transit radius of (Deleuil et al., 2008). Second, a review of the data for all of the transiting Jupiterlike exoplanets reveals also that objects with radii below do not exist (Figure 7). When a cut is made to only include objects with zeroalbedo equilibrium temperatures below 1000 K (to exclude objects that are “inflated” by some unknown mechanism related to stellar heating, e.g., Demory & Seager 2011), we find that the radii are bound between 0.8 and . The single outlier with is Kep447b, which has an extremely grazing transit (LilloBox et al., 2015) that may render its radius measurement unreliable. Third, objects with a mass of Jupiter (or higher) are partially degenerate and it is theoretically challenging to get their radius to be less than that of Jupiter’s (Burrows & Liebert, 1993).
While it may be tempting to fix our model radius at between 0.8 and , we should be reminded of the fact that these radii are measured for Gyrold objects, whereas the HR 8799 exoplanets are estimated to be –100 Myr old. Guided by evolutionary models (Mordasini et al., 2012; Spiegel & Burrows, 2012), we set as a Gaussian prior of our retrievals. The uncertainty of is the fullwidth at halfmaximum of the Gaussian. We note that Moses et al. (2016) assume a fixed value of for their selfconsistent model of HR 8799b.
The bottom panel of Figure 7 is also revealing, as it shows the measured surface gravities of transiting Jupitersized exoplanets to be hovering around for objects with masses of , where is the mass of Jupiter. Since we expect the HR 8799 exoplanets to have radii that are slightly larger than Jupiter’s, we expect their surface gravities to also be . Surface gravities of –5.0 are only appropriate when one crosses over into the brown dwarf regime (), e.g., CoRoT3b has . The photometric masses of HR 8799b, c, d and e are less than half that of CoRoT3b (Marois et al., 2008, 2010). Based on the evolutionary calculations of Marleau & Cumming (2014), who estimated –13 for the HR 8799 exoplanets, we set a Gaussian prior of on the surface gravity (taking into account ). This range of surface gravities is somewhat higher than the values considered by Barman et al. (2015).
In summary, we find that what we assume for the prior distributions of the input quantities is critical to the outcome of the retrieval. Uniform or loguniform priors may not always be the best choice as they may lead to unphysical or even nonsensical outcomes. Gaussian priors are better choices in these instances, but only when they are guided by physics. We find our retrievals to be physically meaningful only when Gaussian priors are set on the radius and surface gravity, which is a departure from the HR 8799b analysis of, e.g., Lee, Heng & Irwin (2013).
4.2 Model Selection Using Bayesian Evidence
Traditionally, model selection is performed manually by the modeler or theorist. One starts with a set of assumptions, computes forward and arrives at a prediction for the thermal structure and synthetic spectrum. These assumptions include chemical equilibrium or disequilibrium, a value for the strength of atmospheric mixing, the number of atoms and molecules included in the model, the metallicity and C/O, etc. Other assumptions are more closely related to technique, e.g., the approximate or limiting form of the radiative transfer equation being solved.
Like all of the other previous studies involving both forward modeling and retrieval, we inevitably make a set of both physical and technical assumptions. However, we use our nested sampling approach to go a step further: we compute the Bayesian evidence for models with and without equilibrium chemistry. We then compare them in order to formally quantify whether equilibrium chemistry is a warranted assumption. Instead of assuming a fixed set of cloud parameters for each retrieval, as was done by Lee, Heng & Irwin (2013), we allow our cloud model to be part of the retrieval and also compare its Bayesian evidence to a retrieval that assumes a cloudfree atmosphere. In these ways, we allow model selection based on the Bayesian evidence to inform us on whether the atmosphere is cloudy or cloudfree and in chemical equilibrium or disequilibrium.
Figure 8 shows a montage of all of the models tested for all four HR 8799 exoplanets. Table 1 explains what the labels of the models correspond to. For HR 8799b, c and d, we see that the Bayesian evidence favors model atmospheres that are not in chemical equilibrium and are cloudy. For HR 8799e, the relative lack of data, compared to the other HR 8799 exoplanets, means that we are unable to strongly select between the different models.
Figure 9 shows the bestfit spectra. Our retrieval procedure generally manages to find good fits to the data, except for the bandhead near 1 m for HR 8799c. We speculate that this mismatch could be due to the influence of an additional molecule we have not included in our analysis, but we deem it beyond the scope of the present paper to identify it. In the Appendix, Figure 21 elucidates the effects of using ExoMol methane and water versus HITRAN methane and HITEMP water.
For the rest of this paper, we will discuss the retrieved properties of the HR 8799b, c and d exoplanets based on the bestfit models only. For HR 8799e, we will discuss results from the model with equilibrium chemistry and that includes all four molecules in the retrieval. In Figures 17, 18, 19 and 20 of Appendix B, we provide the full posterior distributions of the best models for all exoplanets for completeness.
4.3 Retrieving the Cloud Properties and Inferring
Figure 10 shows the retrieved posterior distributions of the cloud particle radius () and composition parameter (). Unsurprisingly, the retrieved values of span a broad enough range (3 to 4 orders of magnitude) that they are uninformative with regards to distinguishing between different compositions, consistent with the expectation that the absorption and scattering properties of the cloud are mainly determined by the particle size and less by the composition (Heng & Demory, 2013).
The inferred values of span a broad range and lie between about 1 and 100 m. The presence of these cloud particles implies that they are being held aloft by atmospheric motion. Since these exoplanets are not being heavily irradiated (unlike for hot Jupiters), we can safely assume that the underlying mechanism driving this motion is convection (Burrows et al., 1997; Chabrier et al., 2000; Baraffe et al., 2002) and estimate approximate values for the associated “eddy diffusion coefficient”, which we denote by . We use equations (15) and (17) of Spiegel, Silverio & Burrows (2009), as well as equations (6) and (8) of Heng & Demory (2013), to calculate the terminal speed associated with a particle of radius , which we denote by . The eddy diffusion coefficient is roughly
(19) 
where the pressure scale height is and is the Boltzmann constant. We follow the prescription of Smith (1998) and use as the characteristic length scale, which is more conservative than what was assumed in Lee, Heng & Irwin (2013). We note that the preceding expression for has no dependence on , as it appears in the numerator of and the denominator of . We assume the intrinsic density of the particles to be 3 g cm.
In Figure 10, we see that spans a broad range of values from cm s to cm s as increases from 1 m to 1 mm. The deviation in the curves between and 1 bar arises from the CunninghamMillikanDavies “slip factor correction” kicking in when the mean free path for collisions between the hydrogen molecules becomes comparable to the cloud particle radius. If we place the retrieved values of corresponding to the peak of each posterior distribution on the plot, we infer – cm s, in agreement with Barman et al. (2015). Madhusudhan, Burrows & Currie (2011) assume – cm s, while Barman et al. (2011) and Marley et al. (2012) assume cm s.
4.4 Retrieving C/O, C/H and O/H for the HR 8799b, c, d and e Exoplanets and Implications for Planet Formation
4.4.1 The Star of HR 8799
We refer to the “metallicity” as the set of elemental abundances with atomic mass numbers that are larger than that of hydrogen and helium. In our current study, these would be and . For comparison, their values in the solar photosphere are and , such that (Lodders, 2003). For the star of the HR 8799 system, Sadakane (2006) has found that
(20) 
4.4.2 Retrieved C/O, C/H and O/H Values
Given the interest in the possibility of carbonrich exoplanets (Gaidos, 2000; Kuchner & Seager, 2005), our retrieval analysis yields the posterior distributions of C/O, C/H and O/H for the atmospheres of HR 8799b, c, d and e in Figures 11, 12 and 13, respectively, which we then compare to the values for the star listed in equation (20). A caveat is that the retrieved values are only for the gaseous phase and the true C/O ratio may be hidden in a condensed phase such as graphite (Moses et al., 2013). The retrieved posterior distributions of C/O and C/H for HR 8799e are not as definitive as for the other three exoplanets, because its K band spectrum has not been measured.
4.4.3 Locations of Snowlines/Icelines
Konopacky et al. (2013) have previously estimated that the HO, CO and CO snowlines or icelines are located at about 10, 90 and 600 AU, respectively. We wish to point out that the iceline locations depend on the formation history of the HR 8799 exoplanets.
In Figure 14, we show calculations of the locations of the CO, CO and HO icelines as functions of the age of the HR 8799 system. We consider two scenarios: an optically thin disk and a vertically isothermal, passively irradiated disk. For the optically thin disk, the temperatures are simply the zeroalbedo equilibrium temperatures at a given distance from the star informed by the Pisa stellar evolution models (Tognelli, Prada Moroni & Degl’Innocenti, 2011). By “passively irradiated”, we mean that viscous heating associated with turbulence is neglected (Chiang & Goldreich, 1997). Both models consider the evolution of stellar heating as the star ages. We expect more sophisticated calculations that involve temperature gradients, photoevaporation and viscous heating to produce iceline curves that are intermediate between these two scenarios. The calculations are shown for to years, because this encompasses the gasclearing phase of the protoplanetary disk. Curiously, the CO iceline sits between different pairs of HR 8799 exoplanets as its location evolves during the gasclearing phase ( years), implying that a variation in the C/O, C/H and O/H values of these exoplanets may be a natural outcome of the planet formation process.
4.4.4 Implications for Planet Formation
Our findings have implications for planet formation, if we assume the retrieved C/O, C/H and O/H values to be representative of the bulk composition of each exoplanet. Öberg, MurrayClay & Bergin (2011) have previously elucidated the chemical signatures associated with the planet formation mechanism and history of an exoplanet. If an exoplanet forms by gravitational instability, the zerothorder expectation is that its C/O, C/H and O/H values mirror that of the star, unless latetime accretion occurred. This is clearly at odds with our inferred values of C/O, C/H and O/H for the HR 8799b, c, d and e exoplanets.
In the context of the core accretion formation mechanism, all four exoplanets should have C/O values that are enhanced above stellar, but below unity, if they formed insitu and in between the water and carbon dioxide snowlines/icelines (Öberg, MurrayClay & Bergin, 2011). Our retrieved values of C/O for HR 8799b and c are consistent with this scenario, whereas HR 8799d and e have subsolar C/O values. Öberg, MurrayClay & Bergin (2011) have suggested that substellar C/O values are still consistent with core accretion if the latetime accretion of planetesimals has occurred to pollute the atmospheres. The link between latetime planetesimal accretion and atmospheric composition has been emphasized by Mordasini et al. (2016). The HR 8799b and c exoplanets have superstellar C/H and O/H values, which suggests that they accreted both carbon and oxygenrich ices. The HR 8799d and e exoplanets, which reside closer to the star, have substellar C/H values but stellar to superstellar O/H values, which suggest the accretion only of oxygenrich ices.
Overall, our retrieved values of C/O, C/H and O/H appear to be consistent with the core accretion formation mechanism and inconsistent with gravitational instability without latetime accretion, as has been suggested by, e.g., Kratter, MurrayClay & Youdin (2010).
4.4.5 Why Spectroscopy in the K Band is Crucial
A lesson we have learned from our analysis is that spectroscopy in the K band is crucial for obtaining meaningful constraints on C/H and C/O, as it affects the ability of the retrieval approach to constrain the abundances of CO and/or CH. The lack of K band spectroscopy for HR 8799e hampers our ability to make stronger statements on its C/H and C/O values. These findings have implications for the design of future instruments on the European Extremely Large Telescope (ELT). Furthermore, multiple wavebands should be monitored simultaneously in order to detect variability (Apai et al., 2016).
5 Discussion
5.1 Summary and Comparison to Previous Work
We have presented the complete methodology for a nested sampling atmospheric retrieval code named HELIOSR, which allows us to insert arbitrary prior distributions of parameters and also compute the full posterior distributions of the retrieved quantities. In its current implementation, we used analytical formulae for the forward model, temperaturepressure profile and equilibrium chemistry, as well as a customized opacity calculator (HELIOSK). By computing the Bayesian evidence, we can compare models that assume equilibrium versus unconstrained chemistry and determine which scenario is favored by the data.
We apply HELIOSR to the measured spectra of the HR 8799b, c, d and e directly imaged exoplanets. We find that the outer HR 8799b and c exoplanets are enriched in carbon and have superstellar and stellar C/O values, respectively. The inner HR 8799d and e exoplanets are diminished in carbon and C/O. All four exoplanets are possibly enriched in oxygen relative to the star, which is a clear signature of latetime accretion of waterrich planetesimals. Figure 15 provides a summary of our findings. We note that our retrieved water abundances are about 2 to 3 orders of magnitude higher than what was found by Madhusudhan et al. (2014) for three hot Jupiters, although it should be noted that these authors do not include a cloud model in their retrievals. The inclusion of a cloud model should worsen the discrepancy between these outcomes. Our retrieved molecular abundances and C/O for HR 8799b are in broad agreement with Lee, Heng & Irwin (2013), despite differences in our retrieval techniques. Table 3 summarizes the properties of the four exoplanets inferred from the retrieval.
Our conclusions differ somewhat from previous studies, which reach a diversity of conclusions. Barman et al. (2011) used selfconsistent models to interpret the H and K band spectra of HR 8799b. They infer and . We deem this radius value to be unphysical for the reasons described in §4.1. Marley et al. (2012) also used selfconsistent models and found that if the theoretical interpretation is made of the photometry alone, then the inferred radius for the HR 8799b exoplanet is 1.11 but with a surface gravity of , considerably higher than the value of Barman et al. (2011). Madhusudhan, Burrows & Currie (2011) used selfconsistent models^{7}^{7}7Strictly speaking, these are parametric models, because the cloud physics is not treated selfconsistently with the gaseous chemistry and is instead parametrized. with various cloud configurations to conclude that the HR 8799b, c and d exoplanets have masses of 2–12, 6–13 and 3–11 , respectively, and surface gravities . In these studies, solar abundance is assumed. The diversity of reported results from these studies already hint at the difficulty of using photometry and spectroscopy to infer the radius and mass of a directly imaged exoplanet from the traditional use of forward modeling.
Barman et al. (2015) performed a manual fitting of the H and K band spectra of HR 8799b and HR 8799c. They first held the CO and CH abundances fixed to their solar values, then fitted for the abundance of HO. The bandheads involving CO and CH are masked or excluded from the fit. Next, the HO abundance is held at its bestfit value (and CH is again held fixed at its solar value) and the abundance of CO is inferred. The final step involves fitting for CH. Such an approach is plausible as a first step, but does not explore the model degeneracies. It is likely that the reported value of for HR 8799b has uncertainties that are underestimated. Barman et al. (2015) themselves remark that, “The various sources of uncertainty in the models (are) not accounted for in the formal mole fraction errorbars.” Building on the work of Barman et al. (2015), Moses et al. (2016) assumed fixed values for the equilibrium temperature, surface gravity, radius, C/O, metallicity and , as well as a fixed temperaturepressure profile. They explored thermo and photochemical models of HR 8799b and produced synthetic spectra that somewhat match the measured spectrum (see their Figure 14).
Lee, Heng & Irwin (2013) analyzed the HR 8799b exoplanet and reported supersolar metallicities for their bestfits, consistent with the present study. They considered two cloud models, where the monodisperse cloud particle radius is fixed manually and not formally included as part of the retrieval. The cloud composition is also assumed to be enstatite, whereas we have allowed the cloud composition to be part of the retrieval. The models of Lee, Heng & Irwin (2013) allowed for and to be uniform or loguniform priors, whereas in the current study we have chosen and to be Gaussian priors. Somewhat surprisingly, despite these differences, they retrieve a C/O value that is similar to what we find (see Figure 11). On the technical side, Lee, Heng & Irwin (2013) used the NEMESIS code, which implements nonlinear optimal estimation (versus the nested sampling algorithm we implemented). This technique, which is also used by Barstow et al. (2015), assumes that the priors and posteriors are Gaussian and is unable to formally perform model selection via Bayesian evidence comparison. Lee, Heng & Irwin (2013) also do not consider equilibrium chemistry in their comparison of models. (See Line et al. 2013 for a comparison of these optimization methods.) Overall, HELIOSR implements a number of improvements over NEMESIS that are more appropriate for the sparse data regime of exoplanetary atmospheres (compared to the remote sensing data of Solar System objects) and is able to more rigorously explore a broader range of parameter space.
5.2 Opportunities for Future Work
There are ample opportunities for future work. Instead of unconstrained chemistry, disequilibrium chemistry may be described by some form of atmospheric mixing (e.g., eddy diffusion). More molecules may be added to the analysis, including acetylene, ethylene and hydrogen cyanide, which are known to be spectroscopically active in the infrared at temperatures higher than for the photospheres of the HR 8799 exoplanets. Ultimately, it is our hope that the collective body of work on atmospheric retrieval will stimulate and connect to work on disk chemistry (e.g., Cridland, Pudritz & Alessi 2016). It will also be insightful to train HELIOSR on a large sample of brown dwarf photometry and spectra, as Line et al. (2015) have done for two T dwarfs.
Appendix A Analytical Formula for the Exponential Integral of the First Order
We may avoid the numerical integration of the exponential integral of the first order by using the approximate, but highly accurate, analytical formulae presented in Abramowitz & Stegun (1970),
(A1) 
The fitting coefficients , and are given in equations (5.1.53) and (5.1.56) of Abramowitz & Stegun (1970), but we reproduce them here for convenience: , , , , and ; , , , , , , , and . As originally stated by Abramowitz & Stegun (1970), the formula involving has a precision better than , while that involving and is precise to better than . In Figure 16, we check these claims by evaluating using a canned routine (expint in IDL) and computing the diffusivity factor using
(A2) 
We label these calculations as “exact”. The calculations labeled “approximate” were performed using the fitting formulae in equation (A1). We see that the error is better than . By contrast, a 13th order polynomial fit to the exact solution incurs large errors ().
Appendix B Full Posterior Distributions for best models of HR 8799b, c, d and e
References
 Abramowitz & Stegun (1970) Abramowitz, M., & Stegun, I.A. 1970, Handbook of Mathematical Functions, 9th printing (New York: Dover Publications)
 Apai et al. (2016) Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40
 Baraffe et al. (2002) Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P.H. 2002, A&A, 382, 563
 Barber et al. (2006) Barber, R.J., Tennyson, J., Harris, G.J., & Tolchenov, R.N. 2006, MNRAS, 368, 1087
 Barman et al. (2011) Barman, T.S., Macintosh, B., Konopacky, Q.M., & Marois, C. 2011, ApJ, 733, 65
 Barman et al. (2015) Barman, T.S., Konopacky, Q.M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61
 Barstow et al. (2015) Barstow, J.K., Aigrain, S., Irwin, P.G.J., Kendrew, S., & Fletcher, L.N. 2015, MNRAS, 448, 2546
 Benneke & Seager (2013) Benneke, B., & Seager, S. 2013, ApJ, 778, 153
 Bonnefoy et al. (2016) Bonnefoy, M., Zurlo, A., Baudino, J.L., et al. 2016, A&A, 587, A58
 Bowler (2016) Bowler, B.P. 2016, PASP, 128, 102001
 Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
 Burrows & Liebert (1993) Burrows, A., & Liebert, J. 1993, Reviews of Modern Physics, 65, 301
 Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W.B., et al. 1997, ApJ, 491, 856
 Burrows, Heng & Nampaisarn (2011) Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47
 Chabrier et al. (2000) Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464
 Chiang & Goldreich (1997) Chiang, E.I., & Goldreich, P. 1997, ApJ, 490, 368
 Cridland, Pudritz & Alessi (2016) Cridland, A.J., Pudritz, R.E., & Alessi, M. 2016, MNRAS, 461, 3274
 Deleuil et al. (2008) Deleuil, M., Deeg, H.J., Alonso, R., et al. 2008, A&A, 491, 889
 Deming et al. (2013) Deming, D., Wilkins, A., McCullough, P., et al. 2013, ApJ, 774, 95
 Demory & Seager (2011) Demory, B.O., & Seager, S. 2011, ApJS, 197, 12
 ForemanMackey et al. (2016) ForemanMackey , D., Vousden, W. , PriceWhelan , A., et al. 2016, 10.5281/zenodo.53155
 Feroz & Hobson (2008) Feroz, F., & Hobson, M.P. 2008, MNRAS, 384, 449
 Feroz, Hobson & Bridges (2009) Feroz, F., Hobson, M.P., & Bridges, M. 2009, MNRAS, 398, 1601
 Feroz et al. (2013) Feroz, F., Hobson, M.P., Cameron, E., Pettitt, A.N. 2013, arXiv:1306.2144
 Fu & Liou (1992) Fu, Q., & Liou, K.N. 1992, Journal of the Atmospheric Sciences, 49, 2139
 Gaidos (2000) Gaidos, E.J. 2000, Icarus, 145, 637
 Grimm & Heng (2015) Grimm, S.L., & Heng, K. 2015, ApJ, 808, 182
 Guillot (2010) Guillot, T. 2010, A&A, 520, A27
 Han et al. (2014) Han, E., Wang, S.X., Wright, J.T., et al. 2014, PASP, 943, 827
 Hansen, Schwartz & Cowan (2014) Hansen, C.J., Schwartz, J.C., Cowan, N.B. 2014, MNRAS, 444, 3632
 Hedges & Madhusudhan (2016) Hedges, C., & Madhusudhan, N. 2016, MNRAS, 458, 1427
 Helling & Woitke (2006) Helling, Ch., & Woitke, P. 2006, A&A, 455, 325
 Heng et al. (2012) Heng, K., Hayek, W., Pont, F., & Sing, D.K. 2012, MNRAS, 420, 20
 Heng & Demory (2013) Heng, K., & Demory, B.O. 2013, ApJ, 777, 100
 Heng, Mendonça & Lee (2014) Heng, K., Mendonça, J.M., & Lee, J.M. 2014, ApJS, 215, 4
 Heng, Lyons & Tsai (2016) Heng, K., Lyons, J.R., & Tsai, S.M. 2016, ApJ, 816, 96
 Heng & Lyons (2016) Heng, K., & Lyons, J.R. 2016, ApJ, 817, 149
 Heng & Tsai (2016) Heng, K., & Tsai, S.M. 2016, ApJ, 829, 104
 Hiranaka et al. (2016) Hiranaka, K., Cruz, K.L., Douglas, S.T., Marley, M.S., & Baldassare, V.F. 2016, arXiv:1606.09485
 Ingraham et al. (2014) Ingraham, P., Marley, M.S., Saumon, D., et al. 2014, ApJL, 794, L15
 Kass & Raftery (1995) Kass, R.E., & Raftery, A.E. 1995, JASA, 90, 773
 Konopacky et al. (2010) Konopacky, Q.M., Ghez, A.M., Barman, T.S., et al. 2010, ApJ, 711, 1087
 Konopacky et al. (2013) Konopacky, Q.M., Barman, T.S., Macintosh, B.A., & Marois, C. 2013, Science, 339, 1398
 Kuchner & Seager (2005) Kuchner, M.J., & Seager, S. 2005, arXiv:astroph/0504214
 Kratter, MurrayClay & Youdin (2010) Kratter, K.M., MurrayClay, R.A., & Youdin, A.N. 2010, ApJ, 710, 1375
 Kreidberg et al. (2014) Kreidberg, L., Bean, J.L., Désert, J.M., et al. 2014, ApJL, 793, L27
 Lacis & Oinas (1991) Lacis, A.A., & Oinas, V. 1991, Journal of Geophysical Research, 96, 9027
 Lee, Heng & Irwin (2013) Lee, J.M., Heng, K., & Irwin, P.G.J. 2013, ApJ, 778, 97
 LilloBox et al. (2015) LilloBox, J., Barrado, D., Santos, N.C., et al. 2015, A&A, 577, A105
 Line et al. (2013) Line, M.R., Wolf, A.S., Zhang, X., et al. 2013, ApJ, 775, 137
 Line et al. (2015) Line, M.R., Teske, J., Burningham, B., Fortney, J.J., & Marley, M.S. 2015, ApJ, 807, 183
 Line et al. (2016) Line, M.R., Stevenson, K.B., Bean, J., et al. 2016, arXiv:1605.08810
 Lodders (2003) Lodders, K. 2003, ApJ, 591, 1220
 Madhusudhan & Seager (2009) Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24
 Madhusudhan, Burrows & Currie (2011) Madhusudhan, N., Burrows, A., & Currie, T. 2011, ApJ, 737, 34
 Madhusudhan (2012) Madhusudhan, N. 2012, ApJ, 758, 36
 Madhusudhan et al. (2014) Madhusudhan, N., Crouzet, N., McCullough, P.R., Deming, D., & Hedges, C. 2014, ApJL, 791, L9
 Mandell et al. (2013) Mandell, A.M., Haynes, K., Sinukoff, E., et al. 2013, ApJ, 779, 128
 Marleau & Cumming (2014) Marleau, G.D., & Cumming, A. 2014, MNRAS, 437, 1378
 Marley et al. (2012) Marley, M.S., Saumon, D., Cushing, M., et al. 2012, ApJ, 754, 135
 Malik et al. (2017) Malik, M., Grosheintz, L., Mendonça, J.M., et al. 2017, AJ, 153, 56
 Marois et al. (2008) Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348
 Marois et al. (2010) Marois, C., Zuckerman, B., Konopacky, Q.M., et al. 2010, Nature, 468, 1080
 Mihalas (1970) Mihalas, D. 1970, Stellar Atmospheres (San Francisco: Freeman and Company)
 MillerRicci & Fortney (2010) MillerRicci, E., & Fortney, J.J. 2010, ApJL, 716, L74
 Mordasini et al. (2012) Mordasini, C., Alibert, Y., Georgy, C., Dittkrist, K.M., Klahr, H., & Henning, T. 2012, A&A, 547, A112
 Mordasini et al. (2016) Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, in press (arXiv:1609.03019)
 Moses et al. (2013) Moses, J.I., Line, M.R., Visscher, C., et al. 2013, ApJ, 777, 34
 Moses et al. (2016) Moses, J.I., Marley, M.S., Zahnle, K., et al. 2016, ApJ, in press (arXiv:1608.08643)
 Öberg, MurrayClay & Bergin (2011) Öberg, K.I., MurrayClay, R., & Bergin, E.A. 2011, ApJL, 743, L16
 Oppenheimer et al. (2013) Oppenheimer, B.R., Baranec, C., Beichman, C., et al. 2013, ApJ, 768, 24
 Pierrehumbert (2010) Pierrehumbert, R.T. 2010, “Principles of Planetary Climate” (New York: Cambridge University Press)
 Press et al. (2007) Press, W.H., Teukolsky, S.A., Vetterling, W.T., & Flannery, B.P. 2007, Numerical Recipes: the Art of Scientific Computing, Third Edition (New York: Cambridge University Press)
 Richard et al. (2012) Richard, C., Gordon, I.E., Rothman, L.S., et al. 2012, Journal of Quantitative Spectroscopy & Radiative Transfer, 113, 1276
 Rothman et al. (1996) Rothman, L.S., Rinsland, C.P., Goldman, A., et al. 1996, Journal of Quantitative Spectroscopy & Radiative Transfer, 60, 665
 Rothman et al. (2010) Rothman, L.S., Gordon, I.E., Barber, R.J., et al. 2010, Journal of Quantitative Spectroscopy & Radiative Transfer, 111, 2139
 Rothman et al. (2013) Rothman, L.S., Gordon, I.E., Babikov, Y., et al. 2013, Journal of Quantitative Spectroscopy & Radiative Transfer, 130, 4
 Sadakane (2006) Sadakane, K. 2006, PASJ, 58, 1023
 Skilling (2006) Skilling, J. 2006, BayAn, 1, 833
 Smith (1998) Smith, M.D. 1998, Icarus, 132, 176
 Spiegel, Silverio & Burrows (2009) Spiegel, D.S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487
 Spiegel & Burrows (2012) Spiegel, D.S., & Burrows, A. 2012, ApJ, 745, 174
 Stevenson et al. (2014) Stevenson, K.B., Désert, J.M., Line, M.R., et al. 2014, Science, 346, 838
 Tognelli, Prada Moroni & Degl’Innocenti (2011) Tognelli, E., Prada Moroni, P.G., & Degl’Innocenti, S. 2007, A&A, 533, A109
 Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71
 Tsai et al. (2017) Tsai, S.M., Lyons, J.R., Grosheintz, L., et al. 2017, ApJS, 228, 20
 van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
 Waldmann et al. (2015) Waldmann, I.P., Tinetti, G., Rocchetto, M., et al. 2015, ApJ, 802, 107
 Yurchenko & Tennyson (2014) Yurchenko, S.N., & Tennyson, J. 2014, MNRAS, 440, 1649
 Zaghloul & Ali (2012) Zaghloul, M.R., & Ali, A.N. 2012, ACM Transactions on Mathematical Software, 38, 15 (arXiv:1106.0151)
 Zurlo et al. (2016) Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587 A57