Unveiling cloudy exoplanets: the influence of cloud model choices on retrieval solutions
MMNRAS , 1–15 (2019) Preprint 27 July 2020 Compiled using MNRAS L A TEX style file v3.0
Unveiling cloudy exoplanets: the influence of cloud modelchoices on retrieval solutions
Joanna K. Barstow, , (cid:63) School of Physical Sciences, The Open University, Walton Hall, Milton Keynes, MK7 6AA, UK Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
In recent years, it has become clear that a substantial fraction of transiting exoplanetshave some form of aerosol present in their atmospheres. Transit spectroscopy – mostlyof hot Jupiters, but also of some smaller planets – has provided evidence for this,in the form of steep downward slopes from blue to red in the optical part of thespectrum, and muted gas absorption features throughout. Retrieval studies seeking toconstrain the composition of exoplanet atmospheres must therefore account for thepresence of aerosols. However, clouds and hazes are complex physical phenomena, andthe transit spectra that are currently available allow us to constrain only some of theirproperties. Therefore, representation of aerosols in retrieval models requires that theyare described by only a few parameters, and this has been done in a variety of wayswithin the literature. Here, I investigate a range of parameterisations for exoplanetaerosol and their effects on retrievals from transmission spectra of hot Jupiters HD189733b and HD 209458b. I find that results qualitatively agree for the cloud/hazeitself regardless of the parameterisation used, and indeed using multiple approachesprovides a more holistic picture; the retrieved abundance of H O is also very robustto assumptions about aerosols. I also find strong evidence that aerosol on HD 209458bcovers less than half of the terminator region, whilst the picture is less clear for HD189733b.
Key words: radiative transfer – planets and satellites: atmospheres – techniques:spectroscopic
In recent years, the characterisation of transiting exoplanetatmospheres has evolved to the extent that several compar-ative studies of the most favourable targets have been pub-lished. These targets mostly consist of hot Jupiters, whichare ideal candidates for transit spectroscopy. Their high tem-peratures and H -He dominated atmospheres result in largeatmospheric scale heights, and therefore large fluctuationsin the transit depth as a function of wavelength.Evidence for the presence of cloud or haze (or a lack ofevidence for its absence) has been found in the majority ofhot Jupiter spectra. Whilst simulations including condensa-tional processes or cloud microphysics can predict the cloudthat is expected to form under particular circumstances,we have no prior knowledge of likely cloud structures onintensely irradiated, hot worlds. Therefore, simple, param-eterised retrieval models (e.g. Lee et al. 2012; Line et al. (cid:63) E-mail: [email protected]
Hubble SpaceTelescope
Wide Field Camera 3 (WFC3) instrument only(Tsiaras et al. 2018; Fisher & Heng 2018), and two also in-corporating
Hubble
Space Telescope Imaging Spectrograph © a r X i v : . [ a s t r o - ph . E P ] J u l Joanna K. Barstow (STIS) and
Spitzer
InfraRed Array Camera (IRAC) data(Barstow et al. 2017; Pinhas et al. 2019). The general philos-ophy of these studies is the same; to use a relatively agnosticmodel set up to retrieve the basic atmospheric properties fora range of planets, in such a way that the results are directlycomparable. Retrieved properties common to all studies in-clude H O abundance; temperature structure (with varyingassumptions); radius or pressure baseline; and cloud param-eters. A detailed comparison of the retrieval results fromthese studies is included in Barstow & Heng (2020). Of par-ticular interest in this work are the differences in cloud pa-rameterisation, and the effects of this on the retrieved cloudproperties.In this work, I contrast the different approaches to pa-rameterising cloud, and test each of the cloud models in turnon a benchmark dataset. For this purpose, I have chosen thespectra of HD 189733b Bouchy et al. (2005) and HD 209458bCharbonneau et al. (2000) as presented in Sing et al. (2016).These are the most precisely measured hot Jupiter spectraso far, and both contain clear evidence for scattering parti-cles in the atmosphere of the planet, albeit of different kinds.Whilst the HD 189733b spectrum has a substantial slope inthe visible, which may be an indicator of the presence ofsmall, Rayleigh scattering particles, the HD 209458b spec-trum has a muted H O feature but very little optical slope,suggesting its spectrum is more likely to be dominated bylarger particles. Together, these spectra provide a good testof the ability of each cloud model to match a diverse rangeof observations.Previous work by Mai & Line (2019) applies a similarapproach, for more physically-motivated cloud parameteri-sations, to simulated
James Webb Space Telescope spectraof WASP-62b. They conclude that including cloud in themodel is necessary for retrieving spectra of a cloudy planet,but that all cloud models they tested provide unbiased so-lutions for other atmospheric properties, despite the differ-ences in approach.In Sections 2.1 and 2.2 I describe the basic retrievalmodel setup. In Section 2.3, I discuss the differences betweenthe cloud models used in each paper, and the way in whichthey have been implemented. Retrieval results are presentedin Section 3, and discussed further in Section 4.
I use the NEMESIS radiative transfer and retrieval code,originally developed for Solar System planets (Irwin et al.2008) and subsequently extended for use with exoplanets(Lee et al. 2012; Barstow et al. 2017; Krissansen-Totton et al.2018). NEMESIS incorporates a fast, 1D radiative transfercalculation, which makes use of the correlated-k approxima-tion allowing absorption line data to be stored in a quick lookup table (Lacis & Oinas 1991). Whilst originally NEME-SIS used an Optimal Estimation technique to converge onthe most likely atmospheric solution (Rodgers 2000), thedependence of this method on an informative prior meansthat it is less appropriate for exoplanets since in this caseour prior knowledge is severely restricted. Instead, a recentupgrade (Krissansen-Totton et al. 2018) sees NEMESIS in- terface with the PyMultiNest algorithm (Feroz & Hobson2008; Feroz et al. 2009, 2013; Buchner et al. 2014).
The basic atmospheric models for HD 189733b and HD209458b are based on those used in Barstow et al. (2017).We restrict the spectrally active gases to Na, K and H O, asthese are the only gases for which strong evidence has beenfound, with prior ranges for log (Na,K) of -13 to -3, and for log (H O) of -8 to -2. They are assumed to be well mixed,so gas abundances account for three parameters in the re-trieval. The bulk of the atmosphere is a mixture of H andHe. Line data for H O are taken from the BT2 database(Barber et al. 2006), and those for Na and K are from theVALD database (Heiter et al. 2008). Collision-induced ab-sorption for H and He is taken from Borysow & Frommhold(1989); Borysow et al. (1989); Borysow & Frommhold(1990); Borysow et al. (2001); Borysow (2002). The Na andK features are expected to be strongly pressure-broadenedin the wings of the absorption lines. To account for this,the line wing cutoff for Na and K is extended to 6000 cm − from the line centre, as opposed to the more usual 25 cm − .Welbanks et al. (2019) explore the H -broadened shape ofthese lines further, and find that at temperatures of 2000K the wings of both Na and K features do not extend be-yond 1.4 µ m; cutting at 6000 cm − for the K band means theline extends to 1.41 µ m. In any case, for cloudy planets thepressure-broadened line wings are generally obscured by thecloud, which can be seen to be the case here.The temperature profile is represented as an isothermfor P < . atm and P > . atm, and as an adiabat inbetween. The stratospheric temperature T strat (prior range:100—3000 K) is therefore the single temperature variable inthe model. I also retrieve the planetary radius at the 10 atmpressure level, because the literature value for this is takenfrom the white light transit, and the pressure that this refersto is extremely dependent on the atmospheric properties ofthe planet. For further context as to why this is necessary, seeBarstow & Heng (2020) Section 1.2 and references therein.The prior range for this value is planet specific, from 0.9—1.3R J for HD 189337b, and from 1.1—1.5 R J for HD 209458b.Therefore, before cloud is introduced, the model atmospherecontains 5 free parameters. Cloud formation is a complex process, and reducing cloudsto something that can be represented by a minimal numberof parameters is challenging. As a result, different teamshave adopted a variety of approaches, but all have somefactors in common.The key features of a cloud that are important to cap-ture when modelling transit spectra are 1) the pressure levelat which the cloud becomes opaque and 2) the wavelengthdependence of this. These effects can be represented in a va-riety of ways. Additional effects may also be important un-der certain circumstances; if a cloud layer is optically thin,the location of the cloud base may also become important.For hot Jupiters, which are tidally locked, strong winds are
MNRAS , 1–15 (2019) loud model choices and retrievals thought to result in the temperature varying considerablybetween the morning and evening terminators, which in turncould lead to variation in cloud coverage (e.g. Line & Par-mentier 2016); it may therefore also be necessary to includesome formulation for fractional cloud cover around the ter-minator. The approaches to representing these key cloudfeatures vary between the four comparative studies, and aresummarised below.Barstow et al. (2017) (hereafter B17) use a grid of cloudmodels. They assume that the cloud has uniform specificdensity for pressures above a variable cloud top pressure.The cloud is treated either as Rayleigh scattering (scatter-ing efficiency scales as 1/ λ , representing small particles) orgrey (constant scattering efficiency). Two vertical cloud dis-tributions are tested - the cloud either extends from the toppressure to the bottom of the atmosphere, or extends down-wards by a decade in pressure. The second approach allowsa detached haze layer to be simulated as well as a deep clouddeck. Finally, total optical depth is a free parameter in anOptimal Estimation retrieval for each cloud model.In this work, the B17 model has been extended fromthe grid-search approach to include four free parameters:total nadir optical depth τ ; top pressure P top ; base pressure P base ; and scattering index γ , where the wavelength depen-dence of the extinction efficiency is proportional to λ − γ .Tsiaras et al. (2018) (T18) and Fisher & Heng (2018)(F18) use effectively the same cloud model as each other,based on that presented by Kitzmann & Heng (2018), whichin turn evolved from models used by Lee et al. (2013) andLavie et al. (2017). They use an analytical model to cap-ture the functional dependence of cloud extinction on wave-length. As described in Fisher & Heng (2018), the cloudopacity is parameterized as follows: κ cloud = κ Q x − a + x . (1)where κ is an optical depth scaling factor, Q de-termines the wavelength at which the extinction efficiencypeaks, a is a scattering slope index and x is the particle sizeparameter, given by x = π r λ (2)where r is the effective particle radius and λ is the wave-length. The model has four free parameters, κ , Q , a and r .The Q parameter is of particular interest because the wave-length at which the extinction efficiency peaks is related tothe composition of the aerosol particles, so retrieval of thisparameter could provide some constraint on possible cloudspecies.The model presented by Kitzmann & Heng (2018) couldbe extended to include constraints on the vertical locationof the cloud; however, here I adopt the version of the modelused in F18, which assumes the cloud is spread verticallythroughout the atmosphere.T18 use a special case of this formalism, which fixes Q to 50 and a to 4. They fix κ to 5, but additionally retrieve acloud mixing ratio χ C which scales the total opacity in thesame way. This model is combined with a simple opaque,grey cloud for pressures P > P top .Finally, Pinhas et al. (2019) (P19) include a cloud model combining an optically thick grey cloud (as used in T18 for P > P top ) with an overlying haze layer, parameterized byan optical depth scaling and a power law dependence of ex-tinction with wavelength. This follows on from the modelintroduced in MacDonald & Madhusudhan (2017): P < P top : κ = n H a σ ( λ / λ ) − γ P > P top : κ = ∞ where n H is the number density of H molecules, σ is the molecular scattering opacity at reference wavelength λ , and a here represents the ‘Rayleigh enhancement fac-tor’ describing the magnitude of additional scattering dueto clouds. γ is the dependence of the scattering efficiency onwavelength λ .All of these four possible parameterizations have beenincorporated into the updated version of the NEMESIS spectral retrieval code that works with PyMultiNest. Thefree and fixed parameters in each of the models considered,as set in their respective papers, are summarized in Table 1and Figure 1. These were incorporated into
NEMESIS asfar as possible in the same way as in the original papers,but the results have been processed such that they can becompared across models. For example, instead of present-ing opacity scaling factors (which differ in meaning betweenthe four models), after the retrieval I combine the relevantparameters from each model to calculate the nadir opticaldepth at 0.2 µ m, and present this value. For the B17 andP19 models, the optical depth is directly equivalent to thescaling factor due to the fact that the scattering is a simplepower law and the cross section is normalised to 1 at 0.2 µ m; for T18 and F18, the scattering cross section dependson particle size and is not normalised.The prior ranges for each case, as incorporated intoNEMESIS, are also included in Table 1. In general, verywide priors have been used to limit retrieval dependence onthe prior, with a few exceptions. The opacity scaling factorrange is slightly narrower for T18 and F18 because this doesnot directly correspond to an optical depth, and keeping theupper end of the range the same as for B17 and P19 resultedin eventual optical depths that were larger than NEMESIScould cope with; however, these ranges are still extremelywide. The other exceptions are for the F18 model, whichuses prior ranges for the scattering index and shape factorsuggested by the discussion of the model in Kitzmann &Heng (2018). This model is derived from an analytical fit toMie scattering calculations, and the authors found that typ-ical scattering indices were between 3 and 7, and the shapefactor for the various possible species varied from 0.07 to64.98.I also test options for all cloud models to incorporatea fractional cloud coverage parameter. This assumes thatsome fraction f of the terminator has cloud, while the re-mainder (1- f ) is entirely cloud-free. P19 and MacDonald& Madhusudhan (2017) already include a fractional cloudparameter in their analysis; here, I investigate whether theinclusion of an extra free parameter is justified by the infor-mation content of current data. MNRAS000
NEMESIS asfar as possible in the same way as in the original papers,but the results have been processed such that they can becompared across models. For example, instead of present-ing opacity scaling factors (which differ in meaning betweenthe four models), after the retrieval I combine the relevantparameters from each model to calculate the nadir opticaldepth at 0.2 µ m, and present this value. For the B17 andP19 models, the optical depth is directly equivalent to thescaling factor due to the fact that the scattering is a simplepower law and the cross section is normalised to 1 at 0.2 µ m; for T18 and F18, the scattering cross section dependson particle size and is not normalised.The prior ranges for each case, as incorporated intoNEMESIS, are also included in Table 1. In general, verywide priors have been used to limit retrieval dependence onthe prior, with a few exceptions. The opacity scaling factorrange is slightly narrower for T18 and F18 because this doesnot directly correspond to an optical depth, and keeping theupper end of the range the same as for B17 and P19 resultedin eventual optical depths that were larger than NEMESIScould cope with; however, these ranges are still extremelywide. The other exceptions are for the F18 model, whichuses prior ranges for the scattering index and shape factorsuggested by the discussion of the model in Kitzmann &Heng (2018). This model is derived from an analytical fit toMie scattering calculations, and the authors found that typ-ical scattering indices were between 3 and 7, and the shapefactor for the various possible species varied from 0.07 to64.98.I also test options for all cloud models to incorporatea fractional cloud coverage parameter. This assumes thatsome fraction f of the terminator has cloud, while the re-mainder (1- f ) is entirely cloud-free. P19 and MacDonald& Madhusudhan (2017) already include a fractional cloudparameter in their analysis; here, I investigate whether theinclusion of an extra free parameter is justified by the infor-mation content of current data. MNRAS000 , 1–15 (2019)
Joanna K. Barstow
Table 1.
Variable cloud parameters grouped by type, for each of the models considered. Quantities are as defined in the text. Priorranges as used in the retrieval are also included.Property B17 (updated) T18 F18 P19Opacity τ χ c κ a − , 10 − , 10 − , 10 − , 10 Scat index - γ -4 - a - γ
0, 14 Fixed 3, 7 0, 14Top pressure P top , all P top , grey None P top , grey − , 1.0 10 − , 1.0 N/A 10 − , 1.0Base pressure P base None None None P top , 1.0 N/A N/A N/AParticle size None r r NoneN/A 10 − , 10 − , 10 N/AShape factor None 50 Q NoneN/A Fixed 0.1, 65 N/A
Figure 1.
This provides a visual indication of how cloud structure is parameterised in each of the models discussed. The parametershighlighted in red are those that are allowed to vary in the retrieval.
Here I present and compare the results for each model, forHD 189733b in Section 3.1 and for HD 209458b in Sec-tion 3.2. With one exception, all four cloud parameteriza-tions can produce a fit to the data for all scenarios, andthe results are generally consistent. I find that there is goodevidence in favour of including a terminator cloud fractionparameter for HD 209458b, but a more confused picture forHD 189733b.A selection of retrievals are shown within the paper, and the full results can be found in an online repository .All corner plots are generated using the corner.py routine(Foreman-Mackey 2016).Both B17 and P19 found that the H O abundances forthese hot Jupiters were generally below that expected for asolar composition gas under the conditions of HD 189733band HD 209458b. The disequilibrium chemistry, solar com-position models of Moses et al. (2011) predict log(H O) val-ues of -3.42 and -3.45 for HD 189733 and HD 209458b re- https://tinyurl.com/qvpazxw MNRAS , 1–15 (2019) loud model choices and retrievals spectively. I compare these with the values retrieved in thisstudy below. The spectrum of HD 189733b displays a very steep slopeat visible wavelengths, and a muted but still present H Ofeature in the near infrared. The steep visible slope has gen-erally been attributed to the presence of high altitude clouds,and the results presented here corroborate this assumption.The extreme steepness of the visible slope presents achallenge to a completely agnostic approach to retrieval,which aims to place minimal prior constraints on the solu-tion. A very steep slope of this kind can be produced in twoways; either by modelling the cloud such that the extinctionefficiency drops off very rapidly with wavelength, or by in-creasing the temperature and thence the scale height of theatmosphere, which increases the amplitude of all features.For HD 189733b, if the temperature is allowed to varyfreely up to a threshold of 3000 K, then the retrieved termi-nator temperature is around 2000 K for all models. However,we actually have some information about the range of val-ues the terminator temperature can take, since we know themaximum equilibrium temperature of HD 189733b (assum-ing zero albedo) is around 1200 K (Sing et al. 2016; Barstowet al. 2017). In transit, for low spectral resolution, the regionof the atmosphere that is probed at visible to IR wavelengthstypically covers the range 10 − — 0.1 bar (Welbanks & Mad-husudhan 2019), which we expect to be broadly isothermal.Therefore, the temperature in the regions of the terminatorto which we are sensitive should be ≤ Hubble /WFC3 and the model wasnot therefore required to simulate a steep optical scatteringslope. However, this demonstrates that fixing the values ofmodel parameters should be undertaken with caution.In the restricted temperature case, temperatures at thehigh end of the prior range are still favoured. This begs thequestion as to why the retrieval converges preferentially onsolutions involving high temperatures rather than cases witha steep scattering slope produced by the cloud. I suggest thatthe reason for this is that the high temperature solution forthe large optical slope relies only on a single model parame-ter, whereas the solutions involving cloud (regardless of thecloud model chosen) rely on a specific combination of at leastthree parameters. Occam’s razor would therefore favour thetemperature solution, and this is only rejected in the lightof the prior knowledge we have of the planet’s temperature. I adopt the restricted temperature case as the moreplausible scenario, and key retrieved results are shown inTable 2.Other points of interest from these results are the con-sistency of the H O abundance results, and the retrievedscattering index, across B17, F18 and P19 for the homoge-neous cloud case. All models indicate an H O abundance ofapproximately 1/30 × solar, consistent with previous studiesB17 and P19, and cloud with a scattering index of around6.4. Rayleigh scattering corresponds to an index of 4, so thecloud particles present display super-Rayleigh behaviour andare likely to be small. This is discussed further in Section 4.2.The retrieved cloud top pressures must be comparedwith more caution, as this parameter represents differentthings across the four models. In B17 it represents the topof the entire cloud deck, whereas in T18 and P19 it is thetop of the grey, opaque cloud. Therefore, the retrieval of avery low top pressure in B17 compared with the others isconsistent.In summary, all models for HD 189733b indicate thepresence of small particle cloud high in the atmosphere, andabsence of an opaque grey cloud, and sub-solar H O.After performing the initial retrievals including homo-geneous cloud, I extend the analysis to also include a cloudfraction parameter for all cloud models. Key results fromthis analysis are presented in Table 2, along with the log ofthe Bayesian evidence for each model relative to the modelwith the highest evidence. This is also known as the Bayesfactor. Generally, if the Bayes factor difference is > > Model fits to HD 209458b proved more straightforward thanfor HD 189733b due to the lack of a substantial scatteringslope in the optical part of the spectrum. The retrieved tem-perature is in the range expected for the planet, even whena less restrictive temperature prior is used. The T18 modelis also able to produce a good fit to the data in this case.The key results for HD 209458b are presented in Ta-ble 3.2. As for HD 189733b, the best fitting model is high-lighted in green, and tests are run for both homogeneousand heterogeneous cloud.Again, the retrieved H O abundance is consistently subsolar across all models tested, at around 10 ppmv. Unlike theHD 189733b case, however, the evidence is more strongly infavour of a grey cloud deck, as both T18 and P19 modelshave top pressures that would make the grey cloud deck visi-
MNRAS , 1–15 (2019)
Joanna K. Barstow
Figure 2.
This figure compares posterior corner plots for the B17 cloud model with (bottom left) and without (top right) the informativetemperature prior. The secondary solution in the bottom left plot corresponds to the solution presented in B17, arrived at using a morerestricted version of the model. MNRAS , 1–15 (2019) loud model choices and retrievals Figure 3.
As Figure 2 but for the T18 modelMNRAS000
As Figure 2 but for the T18 modelMNRAS000 , 1–15 (2019)
Joanna K. Barstow
Figure 4.
Spectral fits to the HD 189733b data (red) with theT18 model with a restricted temperature prior (dark red) and abroad prior (purple). Spectra are generated using the median val-ues of the posterior distributions. The fit is generally poor for therestricted prior case as the relatively inflexible cloud parameteri-sation does not allow the short wavelength slope to be matched.Shaded regions on the spectrum plot indicate the 2- σ envelopefor each case Figure 5.
Spectral fits to the HD 189733b data for all models,for the restricted temperature case without fractional cloud. Themodel spectra are calculated using the median parameter values.
Figure 6.
As Figure 5 for HD 209458b. ble (see for example Figure 7). By contrast, the top pressureof the B17 model is higher than for HD 189733b, indicatingthat the upper part of the atmosphere is more likely to becloud-free.All four models agree well on cloud fraction for the het-erogeneous cloud case, and all four also show a substantial
Table 2.
Key retrieval results for HD 189733b. The differencein ln(Bayesian Evidence), ∆ ln(BE),is also presented, relative tothe best fitting model of the set. The top row of values for eachcase is assuming a homogeneous terminator cloud, with the lowerrow also retrieving a cloud fraction. The case with the highestBayesian evidence is highlighted in green. The T18 model is high-lighted in red due to the poor quality of the spectral fit.Property B17 T18 F18 P19Log(H O) -4.94 + . − . -5.56 + . − . -5.02 + . − . -4.98 + . − . -4.72 + . − . -7.50 + . − . -4.94 + . − . -3.58 + . − . Scat index 6.34 + . − . N/A 6.37 + . − . + . − . + . − . N/A 6.62 + . − . + . − . Top pressure -6.56 + . − . + . − . N/A 0.34 + . − . -7.2 + . − . -3.45 + . − . N/A -2.29 + . − . Cloud fraction N/A N/A N/A N/A0.68 + . − . + . − . + . − . + . − . ∆ ln(BE) -6.8 -41.9 -6.1 -6.3-6.4 -188.9 -8 0 Table 3.
Key retrieval results for HD 209458b, as Table 2.Property B17 T18 F18 P19Log(H O) -4.89 + . − . -5.02 + . − . -5.11 + . − . -4.95 + . − . -5.15 + . − . -5.19 + . − . -5.18 + . − . -5.06 + . − . Scat index 3.69 + . − . N/A 5.08 + . − . + . − . + . − . N/A 4.84 + . − . + . − . Top pressure -0.65 + . − . -0.79 + . − . N/A -0.61 + . − . -2.89 + . − . -3.64 + . − . N/A -5.6 + . − . Cloud fraction N/A N/A N/A N/A0.37 + . − . + . − . + . − . + . − . ∆ ln(BE) -10.9 -9.4 -7.3 -9-5.3 -1.3 -4 0 improvement in the goodness of fit when the cloud fractionparameter is included, with Bayesian evidence indicatingthat fractional cloud is strongly favoured. With fractionalcloud included, the cloud top pressure for the T18 and P19models is reduced (compare e.g. Figure 8 with Figure 7), aspartial cloud high up has a similar observational signatureto opaque cloud lower down.F18 and P19 have higher scattering indices than B17,but in the case of P19 there is relatively low opacity for theupper cloud/haze above the opaque grey cloud, as in thiscase the grey cloud dominates the signal.For the case where fractional cloud is included, the T18 MNRAS , 1–15 (2019) loud model choices and retrievals Figure 7.
Posterior probability distributions (excluding alkali metals) and best-fit spectrum for HD 209458b using the P19 cloud modelwith 100 % terminator cloud cover. The spectrum is generated using the median values of the posterior distributions. Shaded regions onthe spectrum plot indicate the 1 (darker) and 2 (paler) σ envelopes. and P19 models have similar evidence, whereas B17 andF18 are strongly and moderately disfavoured with respectto P19. T18 and P19 both include an opaque cloud deck,reinforcing the fact that HD 209458b is likely to have greycloud. Despite the fact that the four models include different waysof parameterising the cloud, for the most part a coherent pic-ture emerges. These are visually summarised for each planetin Figure 10.For both cases, the P19 model emerges as the one thatprovides the best fit to the observed spectra. The combina-
MNRAS000
MNRAS000 , 1–15 (2019) Joanna K. Barstow
Figure 8.
As Figure 7, but including a fractional cloud parameter. tion of grey cloud and overlying haze with a tunable scat-tering index seems to provide the most flexibility.Below, I discuss the major findings about H O abun-dance, cloud scattering properties and cloud location foreach planet. I contrast the findings with those presented inthe P19 paper. A similar comparison for the results fromF18 and T18 is not relevant, as these two papers dealt withonly a limited subset of the data considered here. O abundance
For both planets, the retrieved H O abundance is generallyvery robust to different assumptions about cloud, echoingsimilar findings by Mai & Line (2019). The exception to thisis the heterogeneous cloud P19 model for HD 189733b, whichhas a retrieved log(H O) abundance of -3.58, meaning theH O volume mixing ratio is an order of magnitude higherthan for all other models. This is also substantially higherthan the retrieved value from the P19 paper itself (-5.04).The trade-off that causes this can be traced to differencesin the retrieved cloud properties between the homogeneous
MNRAS , 1–15 (2019) loud model choices and retrievals Figure 9.
Posterior probability distributions (excluding alkali metals) and best-fit spectrum for HD 209458b using the F18 cloud modelwith fractional terminator cloud cover. The spectrum is generated using the median values of the posterior distributions. Shaded regionson the spectrum plot indicate the 1 (darker) and 2 (paler) σ envelopes. model in this study, the results of P19 and the heterogeneouscase here. The P19 study has a retrieved scattering index of ∼ , compared with 6.47 for the P19 parameterization inthis work for the homogeneous case and 10.29 for the het-erogeneous case. The higher the scattering index, the morerapidly the haze extinction efficiency drops off as a functionof wavelength. To compensate for the lack of cloud opacityat longer wavelengths when the scattering index is high, thegrey cloud top pressure for the heterogeneous P19 model inthis work is reduced to -2.29 from 0.34 in the homogeneous model. The grey cloud moving higher up means that a higherH O abundance is required to fit the 1.4 micron feature inthe WFC3 bandpass.Another key difference between this work and P19 isthat this work also includes the bandpass integrated pointsfrom the
HST /NICMOS instrument, as published by Pontet al. (2013) after Gibson et al. (2012). P19 did not includethese points, and if they are removed and the fits using theP19 model are repeated, the solution is closer to the P19
MNRAS000
MNRAS000 , 1–15 (2019) Joanna K. Barstow
Figure 10.
I show here a visual summary of the retrieved cloud structure and scattering behaviour. The cloud structure as shownemerges from the combined information of all four models, which indicates that HD 189733b has small-particle aerosols that cover atleast 60 % of the terminator, reaching to low pressures, but no grey cloud; and HD 209458b has opaque grey cloud deep in the atmospherecovering around 40 % of the terminator. The top panels show the wavelength dependence of the cloud extinction cross section for eachmodel case. The cross sections are normalised to the values at 1 µ m for comparison. The solid lines indicate the result for homogeneouscloud and the dotted lines for fractional cloud cover. published result, especially given the large error bars on theH O abundance (Table 4). All H O abundance values for HD 209458b are consis-tent within 1 σ for homogeneous and heterogeneous cloudmodels, and within 2 σ between homogeneous and hetero-geneous models. The H O abundance is also consistentlysubsolar, as for HD 189733b, compared with the value pre-dicted by Moses et al. (2011) of -3.45. The retrieved H O The fractional cloud case without NICMOS points for the P19model has bimodal probability distributions for the radius andcloud top pressures (figures in online material). The result of thisis that the spectrum generated using the median values actuallydoesn’t provide a good fit, as the median falls in between theprobability maxima. The maximum likelihood solution producesa spectrum that provides a much better fit to the observation. abundance from the P19 paper is -4.66 + . − . , slightly higherthan the values presented here, but the probability distribu-tions overlap at the 1 σ level.A likely explanation for any remaining differences be-tween this work and the P19 paper is our different ap-proaches to temperature profile parameterisation. In thiswork, I use a simple model assuming an isotherm plus anadiabat, whereas the more complex 6-parameter model ofP19 allows greater freedom in the structure of the T-p pro-file. Whilst in the region of greatest sensitivity the T-p pro-files are consistent with each other, variation in the loweratmosphere (where there is little constraint) could alter thedeep atmosphere scale height, which could in turn affect theother retrieved properties. MNRAS , 1–15 (2019) loud model choices and retrievals Table 4.
Comparison of retrieved results without the HD 189733bNICMOS data points from this work and P19. Error bars are notquoted for P19 values where they are not specified in detail inthe paper. Values without error bars are as read from probabilitydistribution histograms.Property 100% cloud Fractional cloud P19Log(H O) -5.21 + . − . -4.23 + . − . -5.04 + . − . Scat index 6.49 + . − . + . − . + . − . -1.36 + . − . + . − . Pinhas & Madhusudhan (2017) investigate the most repre-sentative scattering slope index for a range of possible cloudcompositions with different particle sizes, and find that val-ues of around 6 can be achieved with modal sizes of 0.01 µ mfor Na S and ZnS, and 0.1 µ m for MnS. These slopes areonly relevant across relatively narrow spectral ranges in theoptical, and a full study of the effective cross sections forthese species in Pinhas & Madhusudhan (2017) shows thatthe curve somewhat flattens out in the near infrared. In re-ality, the variation of extinction cross section as a functionof wavelength is much more complex than a simple powerlaw relationship.However, it is clear that an extremely steep scatteringslope is required to fit the optical spectrum for HD 189733b,that extends throughout the optical region. This slope doesnot appear to be achievable within any single species testedby Pinhas & Madhusudhan (2017). A secondary solution,indicated by the homogeneous B17 model (see Figure 2),is that a steep slope is created by a detached haze layerhigh in the atmosphere that becomes optically thin at longerwavelengths. Other, more complex solutions could includelayered clouds of different species.Previous explanations for the steep slope have also in-cluded unocculted starspots (McCullough et al. 2014), al-though the data used in this work were in fact already cor-rected for unocculted starspots according to the method out-lined by Pont et al. (2013). However, this method was onlyable to account for the variable level of starspots, whilst thebaseline spot coverage remains unknown and could still havean effect (Rackham et al. 2018).Our best hope for further understanding the cloud prop-erties of HD 189733b is that JWST will uncover spectral sig-natures of a specific condensate in the infrared (e.g. Wake-ford & Sing 2015); at this point, it is unclear which effect ofmany is responsible for the steep optical slope.No such explanations need to be invoked for HD209458b, as relatively flat spectra like this can be pro-duced by any condensate with a large spread of particlesizes. Whilst retrieved scattering indices are high for the P19model, this refers only to the upper layer of cloud which hasa low optical depth - it is the grey cloud deck for this and forthe T18 model that provides the majority of the cloud opac-ity. The scattering index for B17 and F18 is much lower forHD 209458b than for HD 189733b. I therefore conclude that the cloud on HD 209458b is consistent with cloud containinglarge aerosol particles. Again, the picture is more complicated for HD 189733b thanfor HD 209458b. Whilst all cases are consistent with a lowtop pressure for a small particle haze (from the B17 model),and a relatively high top pressure for any grey cloud thatmight be present (from T18 and P19), there is ambiguityabout the spatial location of the cloud. Whilst with the bestfitting P19 model there is a clear improvement when a frac-tional cloud parameter is included, this is not the case forthe B17 or F18 models.The retrieved cloud fraction from the P19 model is con-sistent with the value from the P19 paper itself, althoughthe probability distribution from the P19 paper is double-peaked and also has a secondary maximum close to 1.0. Thisreflects the spread of cloud fraction values from the othermodels tested in this work. There is also substantial degen-eracy between the grey cloud top pressure and cloud fractionfor the P19 model, as when a fractional cloud parameter isincluded the top pressure decreases from greater than 1 bar to less than 10 mbar (Table 2). The lower value here is notconsistent with the result from P19, that has the grey cloudtop pressure at around 1 bar, but removing the NICMOSpoints reduces the retrieved pressure in this work to a some-what closer value (Table 4). In any event, the degeneracyfirst pointed out in Line & Parmentier (2016) between globaland patchy cloud is evident here.On the other hand, there is strong and consistent ev-idence that HD 209458b has a terminator with only ∼ % cloud coverage. This is also in reasonable agreement withthe value from the P19 paper, which has a cloud fraction of51 % , with the probability distribution overlapping with theone in this work.The location of the cloud top for HD 209458b is a littlemore ambiguous than for HD 189733b, because once againthere is substantial degeneracy between the cloud fractionand the cloud top pressure. Lower cloud top pressures arepermitted for fractional terminator cloud coverage, as theirrelative effects on the cloud opacity cancel out. Indeed, forthe P19 model the top of the grey cloud deck could be veryhigh in the atmosphere, which may be somewhat implausibleas large particles would be unlikely to be lofted to pressuresof order a few µ bar. However, it should be noted that theconstraint for the P19 model is only an upper limit on thecloud top pressure.HD 189733b and HD 209458b have fairly similar equi-librium temperatures; HD 209458b is around 300 K warmer.They key difference relevant here is probably that HD209458b is a highly inflated planet, whereas HD 189733bhas a relatively high gravity for a hot Jupiter. This mightexplain why the cloud on HD 189733b appears to be com-posed of much smaller particles than on HD 209458b, aslarge particles would be more likely to rain out in a highergravity environment.Assuming that the cloud on HD 189733b and HD NEMESIS uses units of atmospheres for pressure. 1 atm =1.01325 bar.MNRAS000
Comparison of retrieved results without the HD 189733bNICMOS data points from this work and P19. Error bars are notquoted for P19 values where they are not specified in detail inthe paper. Values without error bars are as read from probabilitydistribution histograms.Property 100% cloud Fractional cloud P19Log(H O) -5.21 + . − . -4.23 + . − . -5.04 + . − . Scat index 6.49 + . − . + . − . + . − . -1.36 + . − . + . − . Pinhas & Madhusudhan (2017) investigate the most repre-sentative scattering slope index for a range of possible cloudcompositions with different particle sizes, and find that val-ues of around 6 can be achieved with modal sizes of 0.01 µ mfor Na S and ZnS, and 0.1 µ m for MnS. These slopes areonly relevant across relatively narrow spectral ranges in theoptical, and a full study of the effective cross sections forthese species in Pinhas & Madhusudhan (2017) shows thatthe curve somewhat flattens out in the near infrared. In re-ality, the variation of extinction cross section as a functionof wavelength is much more complex than a simple powerlaw relationship.However, it is clear that an extremely steep scatteringslope is required to fit the optical spectrum for HD 189733b,that extends throughout the optical region. This slope doesnot appear to be achievable within any single species testedby Pinhas & Madhusudhan (2017). A secondary solution,indicated by the homogeneous B17 model (see Figure 2),is that a steep slope is created by a detached haze layerhigh in the atmosphere that becomes optically thin at longerwavelengths. Other, more complex solutions could includelayered clouds of different species.Previous explanations for the steep slope have also in-cluded unocculted starspots (McCullough et al. 2014), al-though the data used in this work were in fact already cor-rected for unocculted starspots according to the method out-lined by Pont et al. (2013). However, this method was onlyable to account for the variable level of starspots, whilst thebaseline spot coverage remains unknown and could still havean effect (Rackham et al. 2018).Our best hope for further understanding the cloud prop-erties of HD 189733b is that JWST will uncover spectral sig-natures of a specific condensate in the infrared (e.g. Wake-ford & Sing 2015); at this point, it is unclear which effect ofmany is responsible for the steep optical slope.No such explanations need to be invoked for HD209458b, as relatively flat spectra like this can be pro-duced by any condensate with a large spread of particlesizes. Whilst retrieved scattering indices are high for the P19model, this refers only to the upper layer of cloud which hasa low optical depth - it is the grey cloud deck for this and forthe T18 model that provides the majority of the cloud opac-ity. The scattering index for B17 and F18 is much lower forHD 209458b than for HD 189733b. I therefore conclude that the cloud on HD 209458b is consistent with cloud containinglarge aerosol particles. Again, the picture is more complicated for HD 189733b thanfor HD 209458b. Whilst all cases are consistent with a lowtop pressure for a small particle haze (from the B17 model),and a relatively high top pressure for any grey cloud thatmight be present (from T18 and P19), there is ambiguityabout the spatial location of the cloud. Whilst with the bestfitting P19 model there is a clear improvement when a frac-tional cloud parameter is included, this is not the case forthe B17 or F18 models.The retrieved cloud fraction from the P19 model is con-sistent with the value from the P19 paper itself, althoughthe probability distribution from the P19 paper is double-peaked and also has a secondary maximum close to 1.0. Thisreflects the spread of cloud fraction values from the othermodels tested in this work. There is also substantial degen-eracy between the grey cloud top pressure and cloud fractionfor the P19 model, as when a fractional cloud parameter isincluded the top pressure decreases from greater than 1 bar to less than 10 mbar (Table 2). The lower value here is notconsistent with the result from P19, that has the grey cloudtop pressure at around 1 bar, but removing the NICMOSpoints reduces the retrieved pressure in this work to a some-what closer value (Table 4). In any event, the degeneracyfirst pointed out in Line & Parmentier (2016) between globaland patchy cloud is evident here.On the other hand, there is strong and consistent ev-idence that HD 209458b has a terminator with only ∼ % cloud coverage. This is also in reasonable agreement withthe value from the P19 paper, which has a cloud fraction of51 % , with the probability distribution overlapping with theone in this work.The location of the cloud top for HD 209458b is a littlemore ambiguous than for HD 189733b, because once againthere is substantial degeneracy between the cloud fractionand the cloud top pressure. Lower cloud top pressures arepermitted for fractional terminator cloud coverage, as theirrelative effects on the cloud opacity cancel out. Indeed, forthe P19 model the top of the grey cloud deck could be veryhigh in the atmosphere, which may be somewhat implausibleas large particles would be unlikely to be lofted to pressuresof order a few µ bar. However, it should be noted that theconstraint for the P19 model is only an upper limit on thecloud top pressure.HD 189733b and HD 209458b have fairly similar equi-librium temperatures; HD 209458b is around 300 K warmer.They key difference relevant here is probably that HD209458b is a highly inflated planet, whereas HD 189733bhas a relatively high gravity for a hot Jupiter. This mightexplain why the cloud on HD 189733b appears to be com-posed of much smaller particles than on HD 209458b, aslarge particles would be more likely to rain out in a highergravity environment.Assuming that the cloud on HD 189733b and HD NEMESIS uses units of atmospheres for pressure. 1 atm =1.01325 bar.MNRAS000 , 1–15 (2019) Joanna K. Barstow
In this paper, I have focused on the retrievals of H O abun-dance and cloud properties rather than the alkali metals.Na is detected in both atmospheres, whilst only an upperlimit for K can be obtained for HD 189733b. In general,the retrieved abundances are consistent between cloud mod-els, although the presence of fractional cloud affects the Naabundance for HD 189733b.For HD 189733b (without fractional cloud), Na and Kabundances are approximately 600 and 0.03 ppmv respec-tively; when fractional cloud is included, the Na abundanceis reduced to ∼
10 ppm. For HD 209458b, Na and K abun-dances are approximately 200 and 4 ppmv. Whilst the Kabundances are reasonably aligned with expectations (solarabundance from Asplund et al. (2009) is ∼ ∼ This paper has explored a range of cloud parameterisationsthat exist in the literature and applied them to spectra ofHD 189733b and HD 209458b. I find that, whilst each modelhas different approaches to representing cloud, the retrievalresults taken together present a surprisingly coherent andholistic picture of each planet. HD 189733b most likely hascloud made of small particles spread throughout the regionof the atmosphere to which we are sensitive, whereas HD209458b displays a thicker cloud with larger particles thatis restricted both to lower regions of the atmosphere andto roughly 40 % of the terminator. Both planets have H Oabundances that are consistently retrieved to be subsolar,confirming the previous findings of Barstow et al. (2017)and Pinhas et al. (2019).Despite their relatively similar equilibrium tempera-tures, and apparently similar chemistry extrapolated fromtheir H O abundances, the cloud properties of the two plan-ets indicate very different regimes. The lack of large particleson HD 189733b may be attributable to its higher gravity,whereas in the partially cloudy terminator on HD 209458b we may be seeing the effect of a slightly warmer eastern(evening) terminator.We now have access to an ever increasing number of hotJupiters with spectra covering the optical and infrared, sothis type of analysis can and should be expanded to covera broader range of targets. The James Webb Space Tele-scope will provide more insight by probing further into theinfrared, and could potentially reveal what these mysteriousclouds are made of.
ACKNOWLEDGEMENTS
I thank the anonymous referee for a very thorough reportand several helpful suggestions that substantially improvedthe clarity of the paper. I acknowledge the support of a RoyalAstronomical Society Research Fellowship. I thank Pat Irwinfor the use of NEMESIS and Dan Foreman-Mackey for theuse of the corner.py routine, which is available to downloadon GitHub: https://github.com/dfm/corner.py. I am grate-ful to the Center for Open Science for their provision of thefree Open Science Framework data hosting service.
DATA AVAILABILITY
The data for HD 189733b and HD 209458b used in thispaper are as published by Sing et al. (2016). The spectraare available to download here: https://pages.jh.edu/~dsing3/spectra/HD189733b_Sing_2015_Nature.csvhttps://pages.jh.edu/~dsing3/spectra/HD209458b_Sing_2015_Nature.csv
REFERENCES
Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A,47, 481Barber R. J., Tennyson J., Harris G. J., Tolchenov R. N., 2006,MNRAS, 368, 1087Barstow J. K., Heng K., 2020, arXiv e-prints, p. arXiv:2003.14311Barstow J. K., Aigrain S., Irwin P. G. J., Sing D. K., 2017, ApJ,834, 50Benneke B., 2015, arXiv e-prints, p. arXiv:1504.07655Blecic J., Dobbs-Dixon I., Greene T., 2017, ApJ, 848, 127Borysow A., 2002, Astronomy and Astrophysics, 390, 779Borysow A., Frommhold L., 1989, The Astrophysical Journal,341, 549Borysow A., Frommhold L., 1990, The Astrophysical Journal Let-ters, 348, L41Borysow A., Frommhold L., Moraldi M., 1989, The AstrophysicalJournal, 336, 495Borysow A., Jorgensen U. G., Fu Y., 2001, J. Quant. Spectrosc.Radiative Transfer, 68, 235Bouchy F., et al., 2005, A&A, 444, L15Buchner J., et al., 2014, A&A, 564, A125Charbonneau D., Brown T. M., Latham D. W., Mayor M., 2000,ApJ, 529, L45Cubillos P. E., et al., 2017, ApJ, 849, 145Feroz F., Hobson M. P., 2008, MNRAS, 384, 449Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2013, preprint,( arXiv:1306.2144 )Fisher C., Heng K., 2018, MNRAS, 481, 4698MNRAS , 1–15 (2019) loud model choices and retrievals Foreman-Mackey D., 2016, The Journal of Open Source Software,1, 24Gibson N. P., Aigrain S., Roberts S., Evans T. M., Osborne M.,Pont F., 2012, MNRAS, 419, 2683Heiter U., et al., 2008, in Journal of Physics Conference Series. p.012011, doi:10.1088/1742-6596/130/1/012011Helling C., et al., 2016, MNRAS, 460, 855H¨orst S. M., 2017, Journal of Geophysical Research (Planets),122, 432Irwin P. G. J., et al., 2008, JQSRT, 109, 1136Kitzmann D., Heng K., 2018, MNRAS, 475, 94Kitzmann D., Heng K., Oreshenko M., Grimm S. L., Apai D.,Bowler B. P., Burgasser A. J., Marley M. S., 2019, arXiv e-prints, p. arXiv:1910.01070Krissansen-Totton J., Garland R., Irwin P., Catling D. C., 2018,AJ, 156, 114Lacis A. A., Oinas V., 1991, J. Geophys. Res., 96, 9027Lavie B., et al., 2017, AJ, 154, 91Lee J.-M., Fletcher L. N., Irwin P. G. J., 2012, MNRAS, 420, 170Lee J.-M., Heng K., Irwin P. G. J., 2013, ApJ, 778, 97Line M. R., Parmentier V., 2016, ApJ, 820, 78Line M. R., et al., 2013, ApJ, 775, 137MacDonald R. J., Madhusudhan N., 2017, MNRAS, 469, 1979Madhusudhan N., Seager S., 2009, The Astrophysical Journal,707, 24Mai C., Line M. R., 2019, ApJ, 883, 144McCullough P. R., Crouzet N., Deming D., Madhusudhan N.,2014, ApJ, 791, 55Molli`ere P., et al., 2020, arXiv e-prints, p. arXiv:2006.09394Moses J. I., et al., 2011, ApJ, 737, 15Pinhas A., Madhusudhan N., 2017, MNRAS, 471, 4355Pinhas A., Madhusudhan N., Gandhi S., MacDonald R., 2019,MNRAS, 482, 1485Pont F., Sing D. K., Gibson N. P., Aigrain S., Henry G., HusnooN., 2013, MNRAS, 432, 2917Rackham B. V., Apai D., Giampapa M. S., 2018, ApJ, 853, 122Rodgers C. D., 2000, Inverse Methods for Atmospheric Sounding.World ScientificSing D. K., et al., 2016, Nature, 529, 59Tsiaras A., et al., 2018, AJ, 155, 156Wakeford H. R., Sing D. K., 2015, A&A, 573, A122Wakeford H. R., Visscher C., Lewis N. K., Kataria T., MarleyM. S., Fortney J. J., Mand ell A. M., 2017, MNRAS, 464,4247Waldmann I. P., Tinetti G., Rocchetto M., Barton E. J.,Yurchenko S. N., Tennyson J., 2015, ApJ, 802, 107Welbanks L., Madhusudhan N., 2019, AJ, 157, 206Welbanks L., Madhusudhan N., Allard N. F., Hubeny I., Spiegel-man F., Leininger T., 2019, ApJ, 887, L20This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000