Clouds on the hot Jupiter HD189733b: constraints from the reflection spectrum
The hot Jupiter HD 189733b is probably the best studied of the known extrasolar planets, with published transit and eclipse spectra covering the near UV to mid-IR range. Recent work on the transmission spectrum has shown clear evidence for the presence of clouds in its atmosphere, which significantly increases the model atmosphere parameter space that must be explored in order to fully characterise this planet. In this work, we apply the NEMESIS atmospheric retrieval code to the recently published HST/STIS reflection spectrum, and also to the dayside thermal emission spectrum in the light of new Spitzer/IRAC measurements, as well as our own re-analysis of the HST/NICMOS data. We first use the STIS data to place some constraints on the nature of the cloud on HD 189733b, and explore solution degeneracy between different cloud properties and the abundance of Na in the atmosphere; as already noted in previous work, absorption due to Na plays a significant role in determining the shape of the reflection spectrum. We then perform a new retrieval of the temperature profile and abundances of HO, CO, CO and CH from the dayside thermal emission spectrum. Finally, we investigate the effect of including cloud in the model on this retrieval process. We find that the current quality of data does not warrant the extra complexity introduced by including cloud in the model; however, future data are likely to be of sufficient resolution and signal-to-noise that a more complete model, including scattering particles, will be required.
Subject headings:Methods: data analysis – planets and satellites: atmospheres – radiative transfer
Since its discovery in 2005 (Bouchy et al., 2005), the hot Jupiter HD 189733b has been repeatedly observed as it transits and is eclipsed by its parent star, leading to excellent coverage in both its transmission and eclipse spectra from the visible to mid-infrared. This has resulted in HD 189733b being probably the best characterised of all the known exoplanets; it is known from the observed slope of the reflectance spectrum that this unresolved planet would appear a deep shade of blue (Evans et al., 2013), a fact that testifies to the power of the transit spectroscopy technique. The transmission spectrum presented by Pont et al. (2013) shows clear evidence for haze or cloud in the atmosphere of this planet, making it the first transiting planet outside the solar system that is known to be cloudy. Lee et al. (2012) and Line et al. (2012) analyse the dayside spectrum from secondary eclipse observations, and retrieve the temperature-pressure profile and abundances of HO, CO, CO and CH. Knutson et al. (2012) use Spitzer/IRAC phase curves to investigate the longitudinal temperature variability, and de Wit et al. (2012) combine the phase curves with an analysis of the ingress/egress shape in eclipse to place further constraint on spatial variability.
The recent albedo spectrum from HST/STIS obtained by Evans et al. (2013) provides an opportunity to investigate the cloud structure on the dayside. Unlike the transmission spectrum investigated by Pont et al. (2013) and Lee et al. (2013b) which probes the limb of the exoplanet amosphere, the dayside reflection spectrum has near-nadir geometry, and so it is sensitive to deeper regions of the atmosphere. If the cloud is similar at the terminators and on the dayside, we can use this to further constrain its properties. Alternatively, we may see an entirely different cloud layer on the hotter dayside from that observed at the terminator. The albedo spectrum is useful for placing constraint on the cloud, as scattering particles have a significant effect on the optical reflectivity of an atmosphere. We expect there to be fewer gaseous absorbing species in the visible part of the spectrum than in the infrared; distinct absorption features of Na and K are seen in the transmission spectrum, but otherwise cloud appears to be the dominant opacity source in this region (Pont et al., 2013). We also expect little, if any, thermal contribution from the planet itself at these wavelengths due to its temperature.
In this work, we use the Non-linear optimal Estimator for MultivariateE spectral analySIS (NEMESIS) software (Irwin et al., 2008) to calculate synthetic spectra using a simple cloudy model atmosphere for HD 189733b, including multiple scattering. We vary cloud parameters and Na volume mixing ratio (VMR) and compare the resultant spectra to the STIS measurement; calculating the goodness-of-fit parameter allows us to determine the region of model parameter space that provides the best match. We then use a subset of our best-fitting models to examine the effect of clouds on our ability to accurately retrieve temperature and molecular abundances from the infrared dayside emission spectrum, following the work of Lee et al. (2012).
We use the NEMESIS spectral retrieval tool to produce forward models (predicted spectra for a range of model atmospheres) for comparison with the HST/STIS spectra, and also to retrieve the atmospheric state from the thermal emission spectrum as in Lee et al. (2012). NEMESIS was developed by Irwin et al. (2008) for atmospheric retrieval of solar system planets, and has since been extended to enable the same analysis for observations of transiting and directly imaged extrasolar planets (Lee et al., 2012, 2013a; Barstow et al., 2013b). NEMESIS is an Optimal Estimation retrieval model (Rodgers, 2000) and uses a correlated-k radiative transfer model (Lacis and Oinas, 1991; Goody and Yung, 1989). NEMESIS is not a radiative equilibrium model; instead, it simply uses the atmospheric model provided to compute the incoming and outgoing radiative flux. In an irradiated case, it will compute the incoming and scattered/reflected flux from the star, but it will not take into account the heating effect of the incoming stellar flux on the atmosphere. For the multiple scattering runs, NEMESIS uses beams over a user-specified number of zenith angles; azimuthal dependence is accounted for using Fourier decomposition. For more details of the scattering calculation in this work, see Section 1.1.1.
We use the same line and collision-induced absorption data as Lee et al. (2012),
Lee et al. (2013b) and Barstow
et al. (2013a). A list of sources for
absorption line data is given in Table 1. H-He
collision-induced absorption data are taken from the models of
Borysow and Frommhold (1989); Borysow et al. (1989); Borysow and Frommhold (1990); Borysow et al. (1997) and Borysow (2002). The reference
stellar spectrum is taken from the model set made available by
|HO||HITEMP2010 (Rothman et al., 2010)|
|CO||CDSD-1000 (Tashkun et al., 2003)|
|CO||HITRAN1995 (Rothman et al., 1995)|
|CH||STDS (Wenger and Champion, 1998)|
|Na||VALD (Heiter et al., 2008)|
|K||VALD (Heiter et al., 2008)|
Multiple scattering calculations
In order to reproduce an accurate reflection spectrum in the presence of optically thick cloud, it is necessary to include multiple scattering as many scattering events are likely. NEMESIS uses the matrix operator algorithm of Plass et al. (1973) for multiple scattering calculations, where the zenith angle integration is achieved using a 5-point Gaussian-Lobatto quadrature scheme and azimuth angle integration is achieved through Fourier decomposition, with the necessary number of Fourier components determined by the stellar and emission zenith angles. The analytical disc-averaged integration scheme used in Lee et al. (2012) for eclipse spectra is not used here as it is not applicable to scattering situations; instead, we represent the disc average for the synthetic STIS spectra by running multiple scattering calculations with the stellar zenith angles set to each of the 5 Gaussian-Lobatto quadrature angles and the azimuth angle set for back-scattering, since during secondary transit the observer is located in the same direction as the star and the stellar zenith angle is equal to the emission angle. The disc-average is then determined using a weighted average of these 5 different calculations, assuming that the atmospheric conditions are the same at all points on the disc. Retrievals including multiple scattering are computationally expensive. Therefore, we anticipate that the majority of the thermal emission retrieval calculations in the future will still be performed with an extinction-only approximation using the disc integration described by Lee et al. (2012). However, we also tested the effect of including multiple scattering, to test the sensitivity of the retrieval to differences in our modelling approach. The 5-angle approach used to calculate the STIS synthetic spectra is still too time-consuming for the emission spectrum retrieval, so in this case we approximate further by calculating the spectrum for a stellar zenith angle of 45 only, which represents the average angle of weighting function used to compute the disc-averaged spectrum since .
The scattering phase function of the particles is calculated using Mie theory and approximated by the Henyey-Greenstein parameterisation (Henyey and Greenstein, 1941), as in previous work on planetary clouds (e.g. Irwin et al. 2009, Barstow et al. 2012); expected deviations from the true phase function are small compared with the errors on the observed spectrum, so we consider this approximation to be valid. We use the double-peaked version of the phase function:
which represents the phase function as the sum of forward and backward scattering peaks. Here, is the cosine of the scattering angle, and are scattering asymmetry parameters for the forward and backward peaks respectively, and is the fractional contribution of the forward peak to the total phase function. For smaller particles approaching the Rayleigh scattering limit, the parameter is close to 0.5 and there are approximately equal contributions from the forward and backward scattering peaks; as the particle size increases relative to the wavelength of light the value of increases and the scattering becomes more asymmetric. However, in no case is the scattering completely isotropic.
We use the enstatite refractive index values of Scott and Duley (1996) and the MnS values from Huffman and Wild (1967) in our calculations of the scattering parameters. It is worth noting that, as we model realistic particles, they absorb as well as scatter incident radiation. The fraction of light that is absorbed rather than scattered is dependent on the particle size and also on the composition of the particles. In general, more of these realistic particles would be required than idealised, perfectly scattering particles to produce an equivalent planetary albedo in an atmospheric model, since for realistic particles some of the incident radiation is absorbed by the particles rather than scattered back to space. The single scattering albedo is also wavelength dependent, and this can therefore affect the shape of modelled planetary albedo and thermal emission spectra.
The HST/STIS data used in this work are presented in Evans et al. (2013). We use the tabulated eclipse depths binned in 6 channels, spanning the wavelength range 290—570 nm, to place constraint on the cloud properties for HD 189733b.
We also investigate the impact of including clouds on retrievals of temperature structure and atmospheric composition from the available thermal emission spectra. We use the same data as Lee et al. (2012), with the exception that we include the new values for the Spitzer/IRAC 3.6 and 4.5 m points presented in Knutson et al. (2012). Data used as in Lee et al. (2012) are the Spitzer/IRS spectrum from Grillmair et al. (2008), the Spitzer/IRS broadband point originally measured by Deming et al. (2006) and the broadband Spitzer MIPS and Spitzer/IRAC 5.8/8.0 m points from Charbonneau et al. (2008).
We also include a reanalysis of the HST/NICMOS observations of Swain et al. (2009). Gibson et al. (2011) re-analysed a number of exoplanet transmission spectra obtained with NICMOS, including that of XO-1b, originally published by Tinetti et al. (2010). This re-analysis showed that the shape of the NICMOS spectra is highly sensitive to the assumptions made about the instrumental systematics, and that the errors may have been under-estimated in the original publications. This was recently confirmed by new observations of XO-1b in transit with HST/WFC3 (Deming, 2013), which are less affected by systematics, and showed significant discrepancies with the original NICMOS dataset. We expect that the NICMOS emission spectrum of HD189733b may suffer from the same problems as the aforementioned transmission spectra, but it is the only dataset published to date in the crucial near-infrared wavelength range. We therefore re-analysed it using a method proposed by Gibson et al. (2012), which makes minimal assumptions about the nature of the systematics. This re-analysis, which we describe briefly in Section 2.1, leads to larger, but somewhat more robust error bars, and thus more conservative results for the retrieval analysis.
We omit the ground-based K- and L-band data of Waldmann et al. (2012), as these data contain an extremely strong emission feature at 3.3 m which is attributed to non-Local Thermodynamic Equilibrium (non-LTE) processes; modelling this feature is therefore beyond the limits of our code. The attribution of this feature as resulting from the atmosphere of HD 189733b has also been disputed by Mandell et al. (2011) as these authors did not detect it; they suggest it may be a result of telluric absorption. This conflict is so far unresolved and as a result we do not include these data in our analysis.
2.1. Reanalysis of NICMOS emission spectrum
We re-analysed the NICMOS emission spectrum presented in Swain et al. (2009) using the method presented by Gibson et al. (2012). We briefly outline the process here, but refer the interested reader to Gibson et al. (2012) for further information on the data reduction and the details of the systematics modelling method. We model the secondary eclipse and systematics simultaneously in each of 18 wavelength channels, using Gaussian Processes (GP) to model the dependence of the systematics on the instrumental parameters which are known to affect them: the orbital phase of HST and its square ; the and -position of the spectral trace on the detector, its width and its angle relative to the x-axis. The ephemeris, planet-to-star radius ratio, inclination and system scale () were fixed to the values obtained by Winn et al. (2006), These are also consistent with the more recent determination of these parameters by Southworth (2008). The eccentricity was assumed to be zero, as there is no evidence in the literature for a non-zero eccentricity. We used three Markov-Chain Monte Carlo (MCMC) chains of 40,000 steps each to marginalise over all the parameters except the eclipse depth in each channel.
Figure 1 shows the spectrum obtained after re-analysis, with the original from Swain et al. (2009) for comparison. The two versions of the spectrum are largely similar in shape, but our error bars are about a factor of 3 larger, which is consistent with the findings of Gibson et al. (2012). The wavelength range probed by NICMOS is particularly sensitive to the gaseous abundances on the dayside of HD 189733b (Lee et al., 2012), so the increased error bars will impact the constraints we can extract from the thermal emission spectrum. The spectrum and error bars are provided in numeric form in Table 2.
We note that the choice of which instrumental parameters to include in the systematics model for NICMOS affects both the shape of the resultant spectrum and the errors, whether using traditional linear de-correlation methods (Gibson et al., 2011) or our GP approach. However, in a GP model, the length scale associated with each parameter gives an indication of the relevance of that parameter (see Gibson et al. 2012 for details). In the present case, the best-fit length scales for all the instrument parameters we tried including in the model were found to be of the same order of magnitude, which is why we included them all in the final model. Nonetheless, this highlights the fact that the NICMOS results should be viewed with caution, as the spectrum is extremely sensitive to the details of the systematics removal.
The reanalysed spectrum and 1 errors are provided in Table 2.
|Wavelength||Flux ratio (10)||Error|
3. Cloud models
We use the HST/STIS data to investigate the range of cloud models that provide a reasonable match with HD 189733b’s reflection spectrum. We use the best-fit atmospheric temperature profile retrieved by Lee (2012); the other parameters with a significant effect in the STIS wavelength range are all related to cloud, except for the abundance of sodium. Sodium has been detected in transmission spectra of HD 189733b (Redfield et al., 2008; Huitson et al., 2012), but it is not possible to obtain an absolute constraint on the sodium abundance from these measurements. Atmospheric models by Sudarsky et al. (2000) and Heng and Demory (2013) suggest that Na should have a measurable effect on the reflection spectrum. The abundance is left as a free parameter in the model.
Clouds are extremely complex phenomena in radiative transfer. The resultant spectrum for a cloudy atmosphere depends on the scattering properties of the cloud, which in turn are affected by the composition and size of the cloud particles. The altitude and optical depth of the cloud affect how far light from the star will penetrate the atmosphere, which affects the fraction of the light that is absorbed before it can be reflected back. In order to investigate the potential cloud structure in the simplest possible way, the cloud has been modelled for a fixed set of particle sizes, optical depths and cloud altitudes; the free parameters in our model for the STIS spectrum are listed in Table 3, with the range of values taken by each. The scattering efficiency and phase function for the cloud are calculated using Mie theory with a double-peaked Henyey-Greenstein approximation to the phase function; the scattering properties are therefore a consequence of the size and composition of the cloud particles.
All other parameters, including the atmospheric temperature and HO, CO, CO and CH VMRs, are fixed during this analysis at the best-fit values of Lee (2012); none of these gases are expected to have strong absorption features between 0.3 and 0.6, so their only effect on the spectrum is a small effect on the mean molecular weight. We also include Rayleigh scattering from H and He, and the abundances of these gases are also fixed.
|Base pressure (mbar)||1000, 100, 10, 1, uniform|
|Particle size (m)||0.01, 0.03, 0.1, 0.3, 1, 3, 10|
|Optical depth at 0.25 m||0.1, 0.2, 0.4, 0.6, 0.8, 1, 10|
|Na VMR (ppmv)||0.5, 5, 50, 500|
For most of the tests described below, we assume that the cloud particles are made of enstatite, MgSiO (Lecavelier Des Etangs et al. 2008; refractive indices used are taken from Scott and Duley 1996). This is likely to be condensed at the stratospheric temperature on HD 189733b (Fortney et al., 2010), but it is not the only chemical species that may be relevant; Morley et al. (2012) suggest that constituents such as MnS may also form clouds on hot Jupiters, and we consider the impact of changing the composition of the cloud particles in Section 3.2. Conversely, we do not expect to see absorption due to gaseous metal oxides such as TiO in the atmosphere as Ti is likely to be in condensate form at the expected photospheric temperature of HD 189733b (Fortney et al., 2010). This may mean that any condensates present also contain Ti oxides, as suggested by Helling and Woitke (2006).
To test the effect of particle size, we use a range of monodisperse (single-size) cloud populations from 10 nm to 10 m in radius.For the STIS wavelength range, the Rayleigh scattering limit is approached for a particle size of 10 nm, so the scattering behaviour as a function of wavelength will be similar for any particles smaller than 10 nm; the only difference would be that a larger number of smaller particles would be required to achieve the same optical depth. The nadir optical depths (referenced at 0.25 m) are varied between 0.1 and 10, and the cloud is located in a layer 1 decade thick in pressure coordinates (i.e. between 1000 and 100 mbar, or 100 and 10 mbar, etc.). We locate the cloud base at 1000, 100, 10 and 1 mbar, and also test the effect of distributing the cloud particles uniformly throughout the atmosphere. The cloud particle number density scale height is assumed to be the same as the pressure scale height in all cases.
Synthetic reflection spectra are generated for all combinations of these three cloud parameters and the Na abundance parameter, using NEMESIS. We list the values used for each parameter in Table 3.
We also generate a series of cloud-free models to investigate the impact of changing the sodium abundance by itself. It is clear that for a range of sodium abundances the STIS spectrum can be reproduced without the need for clouds (Figure 2; the effect of sodium on the spectrum is as predicted by Sudarsky et al. 2000) but if the atmosphere in fact contains more/less sodium additional scatterers/absorbers would be required (Heng and Demory, 2013).
For each of the cloudy models, the goodness-of-fit parameter is computed to indicate the quality of the fit to the 6 data points in the STIS spectrum. We use this as a means of comparing models across our parameter space; our 4 variables are Na VMR, cloud particle radius, cloud altitude and cloud optical depth. For a good fit, therefore, the should approach 2 (). Contour plots of the are shown for each parameter combination in Figure 3, as well as the number of models for each value of each single parameter with a of less than 10 and less than 5. The contour plots show two-dimensional cuts through the four-dimensional parameter space, with two parameters varying in each plot and the other two held fixed. The values for each fixed parameter are indicated by the dashed lines in the histogram plots. This indicates that the best-fit region of parameter space occurs mostly for either small particles ( m in size) or 10 m particles, 10 solar Na (50 ppmv) and optical depths of 1 or less. The preference for small particles (0.01—0.1 m) is consistent with the findings of Lecavelier Des Etangs et al. (2008).
A subsolar sodium abundance (0.1, or 0.5 ppmv) is insufficient to reproduce the spectrum for the cases we test here. Examples of models with different cloud properties and Na VMRs are shown in Figure 4; models with 5—500 ppmv Na could all produce an adequate fit to the spectrum depending on the assumed cloud properties, but a model with 0.5 ppmv Na does not produce a sufficiently large absorption feature at wavelengths 0.5 m. In general, the cloud altitude is the model parameter with the smallest effect, except for the case where the sodium abundance is 100 solar (Figure 5). The degeneracy between sodium abundance and the cloud properties becomes complex when the sodium abundance is high, because the cloud must make the atmosphere optically thick at precisely the right altitude in order for the model to fit. For 50 ppmv sodium a cloud-free model can fit the spectrum with a of less than 5 (Figure 2), but the increase in sodium abundance means that an additional reflective component is required in the model atmosphere for 500 ppmv.
The data do not allow us to distinguish between the finite deck-type models, where the cloud is confined to a certain pressure range, and models where the cloud particles are uniformly distributed throughout the atmosphere as for many cases we can obtain a good fit using either (Figure 3). Based on the results of Parmentier et al. (2013), a uniform cloud of particles up to 0.1 m in size is plausible; we may expect a cloud formed of larger particles to be located lower in the atmosphere due to sedimentation and rain out. For particle sizes of up to 0.1 m, the uniform cloud model contains fewer assumptions (as we do not define a cloud base or cloud top pressure) and Occam’s razor suggests that this should be preferred for any given values of Na VMR, particle radius and optical depth, where the quality of fit is equal.
Heng and Demory (2013) also compare a set of cloud models to the HST/STIS spectrum, and find similar degeneracy between Na abundance and cloud. They find that the abundance of Na relative to the abundance of cloud particles decreases by a factor 10 if the cloud particle size decreases by a factor 10 from 100 nm to 10 nm. For sub-m particles, we find a similar result. The best fit Na VMR, and therefore the number of Na atoms, is roughly constant as a function of particle size if optical depth and cloud altitude are fixed (Figure 3), and we require a factor 10 increase in the number of cloud particles as the size decreases from 100 nm to 10 nm to maintain a fixed optical depth. As in the work of Heng and Demory (2013), we find that the shape of the reflection spectrum is largely dictated by the particle size and the Na VMR.
3.1. K, TiO and VO
So far, we have only examined the effect of Na on the visible spectrum, and have ignored other species. We include K in the model atmosphere (0.1 ppmv), and Figure 6 shows the effect of varying the abundance over 4 orders of magnitude for a cloud-free model. It can be seen that K has a much smaller effect on the shape of the spectrum than Na (Figure 2), and the magnitude of the change in flux ratio is smaller than the error bars on the STIS spectrum, justifying our decision not to vary K in our main analysis.
HD 189733b is not expected to have a hot enough photosphere for absorption due to gaseous TiO and VO to be present in spectra (Fortney et al., 2010), and no features due to these gases have been observed in the transmission spectrum; however, it is still worthwhile testing the effect of these gases on the spectrum should they unexpectedly be present. We show a series of cloud-free synthetic spectra including absorption due to TiO and VO (using absorption data from Freedman 2011) in Figure 7. The effect is much more significant than the effect of varying the K abundance, and including some TiO and VO as well as Na (10 ppb) improves the spectral fit from the case without TiO and VO. This should not be taken as evidence that TiO and VO are present on HD 189733b, as the problem is degenerate and we still expect this planet to be too cold for these species to be present in the gas phase; in addition, TiO and VO have not yet been detected in the atmospheres of hotter planets (e.g. Sing et al. 2013), so the presence of these molecules in planetary atmospheres is questionable.
3.2. Cloud composition
We have so far assumed that any cloud on HD 189733b is made of enstatite. Morley et al. (2012) propose a range of other minerals that are likely to form clouds on brown dwarfs, and MnS is a possible candidate for brown dwarfs with similar temperatures to HD 189733b. We repeat the analysis for a sodium VMR of 5 ppmv, with clouds composed of MnS. The results are presented in Figure 8.
It is clear that the particle size/optical depth parameter space still encompasses the most important spectral variability for the cloud, as discussed in the Appendix of Lee et al. (2013a). However, choosing MnS instead of enstatite results in a greater range of cloud models that produce a fit to the data with 10, but no models that provide a very good fit to the data (5). This is due to the fact that MnS has a different refractive index to enstatite, which affects the scattering efficiency and phase function of the particles, and serves to illustrate the complexity of the cloud parameter space.
3.3. Comparison with terminator cloud models
The results presented here illustrate the highly degenerate nature of any cloudy atmosphere solutions for HD 189733b. If we are to assume that the dayside atmosphere is similar to the terminator atmosphere, then we can place some further constraint on the cloud by comparing with results from the transmission spectrum. The analysis of Lee et al. (2013b) shows that the best-fit models for the terminator atmosphere, for a uniformly-distributed cloud, have particles 0.1 m in size, with an optical depth of 0.01 at 1.0 m. This corresponds to our set of solutions for uniformly-distributed 0.1 m particles with an optical depth of 0.5 at 0.25 m, as the extinction cross-section of 0.1 m enstatite particles is a factor 50 smaller at 1.0 m than at 0.25 m. A of 6 is achieved for these particles if the Na VMR is 5 ppmv, and 3 for a VMR of 50 ppmv. If we can expect the cloud structure to be the same at the terminator and on the dayside, we could therefore place some limits on the Na VMR.
Whether or not we expect the same cloud to form is a question for general circulation models of HD 189733b to solve. So far, few have dealt with cloudy atmospheres, with some exceptions. Parmentier et al. (2013) present a 3D circulation model for the somewhat hotter planet HD 209458b, including passive tracer particles of various sizes to examine the relative effects of circulation and sedimentation on aerosols and condensates; they find that, whereas larger particles sediment out and show significant longitudinal variability in abundance, the number of 0.1 m-sized tracer particles remains relatively uniform with longitude. This indicates that a cloud composed of 0.1 m particles could easily be approximately uniform throughout the atmosphere of a hot Jupiter, and therefore we may be observing the same cloud deck on the dayside as Lee et al. (2013b) observe at the terminator, although it is likely that different altitudes are probed.
4. Thermal emission retrieval
We now take a subset of the cloud models used in Section 3 and examine their effect on retrievals of temperature and molecular abundances from the infrared secondary eclipse spectrum (Figure 9). The retrieval is performed as in Lee et al. (2012) using hemispheric integration for non-scattering runs and a 45 emission angle approximation for scattering runs. We also use our reanalysed NICMOS spectrum and the warm Spitzer 3.6, 4.5 and 8 m points presented in Knutson et al. (2012). We include the cloud model in the model atmosphere, but do not make any attempt to retrieve it; instead, we examine its effect on the other parameters in the retrieval.
Firstly, we repeat the cloud-free analysis of Lee et al. (2012) with the addition of the reanalysed NICMOS data points; Lee et al. (2012) state that the NICMOS spectrum provides information about the temperature profile in the deep atmosphere, and also the abundances of trace gases, and as these data have changed we expect to retrieve a different temperature profile and abundances. We perform the retrieval of temperature as a function of pressure, over 50 atmospheric levels with a correlation length of 1.5 in ln(pressure), and altitude-independent abundances of HO, CO, CO and CH using a range of different temperature and gas abundance priors; since the problem is underconstrained, it is necessary to test the influence of the prior on the retrieved property. We present the results of this analysis for temperature in Figure 10 and Table 4. The best-fit retrieved temperature profile (Figure 10 top left) and gas abundances are averaged over all runs. The temperature sensitivity weighting function averaged over all wavelengths is shown qualitatively in Figure 10 by a shaded bar, with the peak (lightest colour) occuring between 1 and 100 mbar; this corresponds to the region of the atmosphere probed by the majority of the data, and therefore the temperature is least dependent on the temperature prior in this region; this is shown in Figure 10 top right, in which the retrieved profiles for different priors converge between 1 and 100 mbar.
The temperature structure is most sensitive to the Knutson et al. (2012) Spitzer/IRAC 3.6 and 4.5 m points as these data have by far the smallest fractional error, meaning that they drive the temperature retrieval. This accounts for the differences in the shape of the retrieved T-p profile compared with that of Lee et al. (2012). The 3.6 m weighting function peaks at around 1 mbar and constrains the stratospheric temperature, whereas the 4.5 m weighting function peaks at the tropopause, at around 100 mbar.
|Mean||Min||Max||Lee et al.|
In Figure 10 it can be seen that the greatest degeneracy between temperature and gas abundance occurs for HO and CO, with the precise shape of the retrieved T-p profile between 1 and 100 mbar varying as a function of the HO and CO abundance priors. This is because the peak of the temperature weighting function and therefore the pressure probed varies with atmospheric opacity, and these two gases are the most spectrally active. The increase in the error bars for the NICMOS spectrum means that the gas abundance retrievals are more sensitive to the prior than found by Lee et al. (2012), so for a small HO or CO prior abundance the retrieved abundance is smaller, the atmospheric opacity is less and the pressure probed is higher. The reverse is true for high prior abundances of HO or CO. The effects of changing the CO and CH abundances are smaller, but similar.
The average temperature profile retrieved is largely similar to that of Line et al. (2013), except that we see a slight decrease in temperature above the 1 mbar level compared with Line et al. (2013). This is probably due to the fact that we retrieve temperature for each atmospheric level individually, with an assumed correlation length (1.5 in log-pressure coordinates), whereas Line et al. (2013) parameterise the temperature profile. Therefore, the upper atmosphere temperature retrieved by Line et al. (2013) cannot vary freely with respect to the temperature between 1 and 100 mbar, where the sensitivity is greatest, whilst we retrieve a continuous profile.
The mean, standard deviation and range of retrieved gas VMRs are given in Table 4. The standard deviations are highest for CO and CO, implying a high dependence on the prior. The retrieved values for HO and CO are similar to that of Lee et al. (2012), but the value for CO is smaller, although due to the dependence on the prior this result must be viewed with caution. However, we retrieve a much higher CH abundance than Lee et al. (2012), roughly a factor of 10 higher than the Spitzer-derived upper limit of Madhusudhan and Seager (2009). Our higher value for CH is similar to the result of Line et al. (2013), and so it is likely to be a result of the revised 3.6 and 4.5 m points which are also used by these authors.
We now seek to investigate the effect of cloud on the retrieved temperature profile and gas abundances. We use the average temperature prior and a priori gas VMRs of 10 for all retrieved species. The parameters for the set of models we use are shown in Table 5. These correspond to some of the best-fitting solutions from the previous analysis that have very different infrared optical depths, so are likely to affect the thermal emission retrieval differently.
|Model||Base p||0.25 m O. D.||Particle radius|
|2||1000 mbar||1||0.01 m|
|3||1000 mbar||1||3 m|
|4||10 mbar||10||0.03 m|
We show the retrieved temperature profiles and gas VMRs for each of the cloud models, as well as the gas abundances for the cloud-free case (Figure 11, Table 6). We perform each of the retrievals twice, first using the fast, extinction-only disc integration technique as in Lee et al. (2012), and secondly using multiple scattering as in Section 3, except we approximate the disc with a single 45 emission/incidence angle as the 5-angle calculation is extremely time-consuming for a full retrieval.
|Model||HO VMR (10)||CO VMR (10)||CO VMR (10)||CH VMR (10)|
|4||3.82.0/6.83.7||3117/137||1.31.3/1.31.3||0.31 0.09/0.31 0.11|
It is clear from Figure 11 and Table 6 that the thermal emission retrieval is largely insensitive to the inclusion of clouds corresponding to our best-fitting models. This is due to the fact that the majority of best-fitting models with optical depths greater than 1 at 0.25 m have small particle sizes, and therefore the extinction drops rapidly as a function of wavelength. These particles have little effect at wavelengths longer than 1 m because the particle size is small with respect to the wavelength of the light, and so the optical depth in the infrared is negligible even though it is significant at shorter wavelengths. This means that even higher optical depths of 0.1 m-sized particles would be unlikely to change the temperature retrieval.
The exception is the retrieved HO and CO abundances and temperature profile for cloud model 4, using an extinction-only assumption; including a high (10—1 mbar) cloud of 0.03 m-sized particles results in a larger retrieved CO abundance, which in turn affects the shape of the temperature profile in this pressure range (as in Figure 10). Where solutions are degenerate, as in this case, it is possible for small changes in one model variable (here, the cloud optical depth) to result in different best-fit solutions being obtained. In Figure 12, we show the synthetic spectrum for the model 4 retrieved case, and the spectrum for a forward model with retrieved temperature profile and gas abundances as for a cloud-free model but including a model 4 cloud. The spectra are very similar, with of 138 and 142 respectively, demonstrating that the difference in retrieved CO abundance for this case should not be treated as significant. Additionally, the multiple scattering retrieval for model 4 yields a lower abundance of CO and a temperature profile more similar to those retrieved using other cloud models. For the case of HD 189733b, therefore, it seems that an accurate retrieval of temperature and atmospheric composition from the thermal emission spectrum may not be critically dependent on an understanding of the cloud properties.
We have used a very simple set of cloud models to investigate the acceptable parameter space for the HD 189733b STIS spectrum. The quality of the data at this stage do not warrant more detailed modelling, but it is instructive to consider the most important steps that could be taken to improve this model in anticipation of future results.
5.1. Monodisperse particles
We have assumed that all the particles in the cloud are of a single size. This is of course not likely to be the case in practice; cloud models based on microphysical processes indicate that particle sizes will change as a function of altitude, and at a given altitude particle size distributions may be relatively broad, for example as discussed by Helling et al. (2008). Coalescence processes, condensational growth and sedimentation will result in the presence of larger grains towards the bottom of the cloud deck.
Broadening the size distribution will make the cloud optical depth more uniform as a function of wavelength; extinction due to cloud particles is most efficient when the particle size is comparable to the wavelength of light. This may result in the cloud having a larger effect on the thermal emission spectrum than we find in this work; although in most cases we find that particle sizes of larger than 0.1 m do not provide a good fit to the data at short wavelengths, this does not preclude larger particles being present as they would scatter the short wavelength light less strongly than smaller particles.
When the data are available that would allow us to test a more complex model, we consider that investigating different particle size distributions would be the most important step to take, as particle size is the property that has the largest effect on the spectrum.
5.2. Particle composition
The composition of clouds on cool brown dwarfs and hot exoplanets is treated somewhat differently in different models. Morley et al. (2012) consider a range of possible cloud compositions for specific temperatures, under the assumption that if clouds of more than one species are present they will form at different temperatures, and therefore altitudes. Helling et al. (2008) have a very different model, in which multiple species condense onto the same TiO ‘seed particle’ (Helling and Woitke, 2006), thus forming grains of mixed composition. This significantly complicates the calculation of spectral properties for these clouds; however, Helling et al. (2008) find that, for a low-gravity atmosphere with T1300 K, the enstatite is the constituent with the largest effect on the spectrum, meaning that enstatite probably represents the best approximation to the multi-component grains.
We tested the effect of including MnS, as proposed by Morley et al. (2012) for cool brown dwarfs, and find that the different spectral properties of the material do affect the spectrum, but to a lesser extent than the size of the particles. Therefore, particle composition is an interesting route for further investigation, but is of secondary importance to developing more detailed representations of particle size.
5.3. Effect of temperature profile
We have modelled the STIS reflection spectrum using the best-fit temperature profile from Lee (2012), but have then retrieved a different profile using the thermal emission spectrum. We tested the effect on the best-fit cloud properties of including the altered temperature profile in our analysis of the cloud. We found that the changes to the spectral shape were negligible compared with the error bars on the STIS spectrum (Figure 13), and the values for different cloud models do not change greatly. Therefore small inaccuracies in the temperature profile we use would not significantly affect our conclusions; the same cloud models provide good/poor fits to the measured spectrum regardless of the temperature profile used.
The reason for the limited effect of temperature is likely to be that the absorption due to Na at the altitudes to which we are sensitive is dominated by pressure-broadened line wings. Temperature variation is more likely to affect the line core absorption at higher altitudes, but at the resolution of the STIS spectrum such variation would not be seen. The STIS measurements also lie in a wavelength region where the signal is dominated by reflected starlight, rather than light emitted from the planet.
We find that a large range of enstatite cloud models can fit the measured HST/STIS spectrum. Cloud-free model atmospheres are also acceptable solutions, although we consider this to be less plausible due to the clear evidence from the transmission spectrum that HD 189733b is cloudy (Pont et al., 2013; Lee et al., 2013b) and the likelihood that small particles are evenly distributed in hot Jupiter atmospheres (Parmentier et al., 2013). Small enstatite particles (0.1 m) and 50 ppmv of Na provide the best fit to the STIS data of the examples we test; we find an overlap with the models of Lee et al. (2013b) for a uniformly-distributed 0.1 m cloud with an optical depth of 0.5 at 0.25 m. However, the problem is extremely degenerate and we cannot exclude solutions with larger cloud particles.
The retrieval of temperature and atmospheric composition from the thermal emission spectrum is relatively insensitive to the inclusion of cloud in the model atmosphere, for our best-fitting models. This suggests that for the case of HD 189733b accurate retrieval of temperature and gaseous abundances from the thermal emission spectrum is possible, even without detailed knowledge of the cloud properties. Solution degeneracy prevents firm conclusions from being drawn about the nature of the cloud on HD 189733b; the current quality and coverage of spectroscopic data for most exoplanets is therefore insufficient to simultaneously constrain temperature structure, gaseous abundances and multiple cloud properties. Given the additional complexity (and therefore number of degenerate solutions) introduced to the retrieval problem when clouds are included, this implies that the best approach with the currently available secondary eclipse data is to use cloud-free model atmospheres for temperature retrieval. As data quality improves, alternative, more detailed cloud models must be explored and their effect on the emission spectrum reassessed.
JKB acknowledges the support of the John Fell Oxford University Press (OUP) Research Fund for this research and LNF is supported by a Royal Society Research Fellowship. We thank the anonymous reviewer for his/her report, and Tom Evans, Frédéric Pont, David Sing and Kevin Heng for helpful discussions about the visible albedo observations.
- F. Bouchy, S. Udry, M. Mayor, C. Moutou, F. Pont, N. Iribarne, R. Da Silva, S. Ilovaisky, D. Queloz, N. C. Santos, et al., A&A 444, L15 (2005), URL http://dx.doi.org/10.1051/0004-6361:200500201.
- T. M. Evans, F. Pont, D. K. Sing, S. Aigrain, J. K. Barstow, J.-M. Désert, N. Gibson, K. Heng, H. A. Knutson, and A. Lecavelier des Etangs, The Astrophysical Journal Letters 772, L16 (2013), 1307.3239.
- F. Pont, D. K. Sing, N. P. Gibson, S. Aigrain, G. Henry, and N. Husnoo, Monthly Notices of the Royal Astronomical Society 432, 2917 (2013), 1210.4163.
- J.-M. Lee, L. N. Fletcher, and P. G. J. Irwin, Monthly Notices of the Royal Astronomical Society 420, 170 (2012), 1110.2934.
- M. R. Line, X. Zhang, G. Vasisht, V. Natraj, P. Chen, and Y. L. Yung, The Astrophysical Journal 749, 93 (2012), 1111.2612.
- H. A. Knutson, N. Lewis, J. J. Fortney, A. Burrows, A. P. Showman, N. B. Cowan, E. Agol, S. Aigrain, D. Charbonneau, D. Deming, et al., The Astrophysical Journal 754, 22 (2012), 1206.6887.
- J. de Wit, M. Gillon, B.-O. Demory, and S. Seager, A&A 548, A128 (2012), 1202.3829.
- J.-M. Lee, P. G. J. Irwin, Leigh, N. Fletcher, K. Heng, and J. K. Barstow, ArXiv e-prints (2013b), 1310.5868.
- P. G. J. Irwin, N. A. Teanby, R. de Kok, L. N. Fletcher, C. J. A. Howett, C. C. C. Tsang, C. F. Wilson, S. B. Calcutt, C. A. Nixon, and P. D. Parrish, Journal of Quantitative Spectorscopy and Radiative Transfer 109, 1136 (2008).
- J.-M. Lee, K. Heng, and P. G. J. Irwin, The Astrophysical Journal 778, 97 (2013a), 1307.1404.
- J. K. Barstow, S. Aigrain, P. G. J. Irwin, L. N. Fletcher, and J.-M. Lee, MNRAS (2013b), 1306.6567.
- C. D. Rodgers, Inverse Methods for Atmospheric Sounding (World Scientific, 2000).
- A. A. Lacis and V. Oinas, J. Geophys. Res. 96, 9027 (1991).
- R. M. Goody and Y. L. Yung, Atmospheric radiation : theoretical basis (1989).
- J. K. Barstow, S. Aigrain, P. G. J. Irwin, N. Bowles, L. N. Fletcher, and J.-M. Lee, MNRAS (2013a).
- A. Borysow and L. Frommhold, The Astrophysical Journal 341, 549 (1989).
- A. Borysow, L. Frommhold, and M. Moraldi, The Astrophysical Journal 336, 495 (1989).
- A. Borysow and L. Frommhold, The Astrophysical Journal Letters 348, L41 (1990).
- A. Borysow, U. G. Jorgensen, and C. Zheng, Astronomy and Astrophysics 324, 185 (1997).
- A. Borysow, Astronomy and Astrophysics 390, 779 (2002).
- E. K. Baines, H. A. McAlister, T. A. ten Brummelaar, N. H. Turner, J. Sturmann, L. Sturmann, P. J. Goldfinger, and S. T. Ridgway, The Astrophysical Journal 680, 728 (2008), URL http://stacks.iop.org/0004-637X/680/i=1/a=728.
- L. S. Rothman, I. E. Gordon, R. J. Barber, H. Dothe, R. R. Gamache, A. Goldman, V. I. Perevalov, S. A. Tashkun, and J. Tennyson, Journal of Quantitative Spectroscopy and Radiative Transfer 111, 2139 (2010).
- S. A. Tashkun, V. I. Perevalov, J.-L. Teffo, A. D. Bykov, and N. N. Lavrentieva, Journal of Quantitative Spectroscopy and Radiative Transfer 82, 165 (2003).
- L. S. Rothman, R. B. Wattson, R. Gamache, J. W. Schroeder, and A. McCann, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, edited by J. C. Dainty (1995), vol. 2471 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, pp. 105–111.
- C. Wenger and J. P. Champion, Journal of Quantitative Spectroscopy and Radiative Transfer 59, 471 (1998).
- U. Heiter, P. Barklem, L. Fossati, R. Kildiyarova, O. Kochukhov, F. Kupka, M. Obbrugger, N. Piskunov, B. Plez, T. Ryabchikova, et al., Journal of Physics Conference Series 130, 012011 (2008).
- G. N. Plass, G. W. Kattawar, and F. E. Catchings, Applied Optics 12, 314 (1973).
- L. G. Henyey and J. L. Greenstein, Astrophysical Journal 93, 70 (1941).
- P. G. J. Irwin, N. A. Teanby, and G. R. Davis, Icarus 203, 287 (2009).
- J. K. Barstow, C. C. C. Tsang, C. F. Wilson, P. G. J. Irwin, F. W. Taylor, K. McGouldrick, P. Drossart, G. Piccioni, and S. Tellmann, Icarus 217, 542 (2012).
- A. Scott and W. W. Duley, Astrophysical Journal Supplement 105, 401 (1996).
- D. R. Huffman and R. L. Wild, Physical Review 156, 989 (1967), URL http://link.aps.org/doi/10.1103/PhysRev.156.989.
- C. J. Grillmair, A. Burrows, D. Charbonneau, L. Armus, J. Stauffer, V. Meadows, J. van Cleve, K. von Braun, and D. Levine, Nature 456, 767 (2008), 0901.4774.
- D. Deming, J. Harrington, S. Seager, and L. J. Richardson, The Astrophysical Journal 644, 560 (2006), arXiv:astro-ph/0602443.
- D. Charbonneau, H. A. Knutson, T. Barman, L. E. Allen, M. Mayor, S. T. Megeath, D. Queloz, and S. Udry, The Astrophysical Journal 686, 1341 (2008), 0802.0845.
- M. R. Swain, G. Vasisht, G. Tinetti, J. Bouwman, P. Chen, Y. Yung, D. Deming, and P. Deroo, The Astrophysical Journal Letters 690, L114 (2009).
- N. P. Gibson, F. Pont, and S. Aigrain, Monthly Notices of the Royal Astronomical Society 411, 2199 (2011), 1010.1753.
- G. Tinetti, P. Deroo, M. R. Swain, C. A. Griffith, G. Vasisht, L. R. Brown, C. Burke, and P. McCullough, The Astrophysical Journal Letters 712, L139 (2010), 1002.2434.
- D. e. a. Deming, The Astrophysical Journal 774, 95 (2013), 1302.1141.
- N. P. Gibson, S. Aigrain, S. Roberts, T. M. Evans, M. Osborne, and F. Pont, Monthly Notices of the Royal Astronomical Society 419, 2683 (2012), 1109.3251.
- I. P. Waldmann, G. Tinetti, P. Drossart, M. R. Swain, P. Deroo, and C. A. Griffith, The Astrophysical Journal 744, 35 (2012), 1104.0570.
- A. M. Mandell, L. Drake Deming, G. A. Blake, H. A. Knutson, M. J. Mumma, G. L. Villanueva, and C. Salyk, The Astrophysical Journal 728, 18 (2011), 1011.5507.
- J. N. Winn, J. A. Johnson, G. W. Marcy, R. P. Butler, S. S. Vogt, G. W. Henry, A. Roussanova, M. J. Holman, K. Enya, N. Narita, et al., The Astrophysical Journal Letters 653, L69 (2006), astro-ph/0609506.
- J. Southworth, MNRAS 386, 1644 (2008), 0802.3764.
- J.-M. Lee, Ph.D. thesis, University of Oxford (2012).
- S. Redfield, M. Endl, W. D. Cochran, and L. Koesterke, The Astrophysical Journal Letters 673, L87 (2008), 0712.0761.
- C. M. Huitson, D. K. Sing, A. Vidal-Madjar, G. E. Ballester, A. Lecavelier des Etangs, J.-M. Désert, and F. Pont, Monthly Notices of the Royal Astronomical Society 422, 2477 (2012), 1202.4721.
- D. Sudarsky, A. Burrows, and P. Pinto, The Astrophysical Journal 538, 885 (2000), astro-ph/9910504.
- K. Heng and B.-O. Demory, The Astrophysical Journal 777, 100 (2013), 1309.5956.
- A. Lecavelier Des Etangs, F. Pont, A. Vidal-Madjar, and D. Sing, Astronomy and Astrophysics 481, L83 (2008), 0802.3228.
- J. J. Fortney, M. Shabram, A. P. Showman, Y. Lian, R. S. Freedman, M. S. Marley, and N. K. Lewis, The Astrophysical Journal 709, 1396 (2010), 0912.2350.
- C. V. Morley, J. J. Fortney, M. S. Marley, C. Visscher, D. Saumon, and S. K. Leggett, The Astrophysical Journal 756, 172 (2012), 1206.4313.
- C. Helling and P. Woitke, Astronomy and Astrophysics 455, 325 (2006).
- V. Parmentier, A. P. Showman, and Y. Lian, ArXiv e-prints (2013), 1301.4522.
- R. Freedman, Private Communication (2011).
- D. K. Sing, A. Lecavelier des Etangs, J. J. Fortney, A. S. Burrows, F. Pont, H. R. Wakeford, G. E. Ballester, N. Nikolov, G. W. Henry, S. Aigrain, et al., Monthly Notices of the Royal Astronomical Society 436, 2956 (2013), 1309.5261.
- M. R. Line, H. Knutson, A. Wolf, and Y. Yung, ArXiv e-prints (2013), 1309.6663.
- N. Madhusudhan and S. Seager, The Astrophysical Journal 707, 24 (2009), 0910.1347.
- C. Helling, P. Woitke, and W.-F. Thi, Astronomy and Astrophysics 485, 547 (2008), 0803.4315.