DETECTION OF H₂O AND EVIDENCE FOR TiO/VO IN AN ULTRA-HOT EXOPLANET ATMOSPHERE

Observations have revealed a diversity of atmospheres across the population of transiting gas giant exoplanets (Sing et al. 2016; Stevenson 2016). Transmission spectra at optical wavelengths show evidence for Rayleigh scattering by molecular hydrogen and high-altitude aerosols, and absorption by alkali metals (e.g., Charbonneau et al. 2002; Pont et al. 2008; Sing et al. 2013, 2015; Nikolov et al. 2014, 2015). In the near-infrared, water absorption has now been robustly measured for a number of planets, both in transmission and emission (e.g., Deming et al. 2013; Wakeford et al. 2013; McCullough et al. 2014; Stevenson et al. 2014; Kreidberg et al. 2015).

At temperatures above 2000 K, models predict that gaseous titanium oxide (TiO) and vanadium oxide (VO) are important absorbers, especially at optical wavelengths (e.g., Hubeny et al. 2003; Sharp & Burrows 2007; Fortney et al. 2008). Indeed, TiO/VO absorption is prominent in late M and early L dwarf atmospheres at these temperatures (Burrows & Sharp 1999; Kirkpatrick et al. 1999; Burrows et al. 2001). If TiO and VO are present in the upper atmosphere of an irradiated gas giant, it would likely generate a thermal inversion, as incoming stellar radiation is absorbed at low pressures (Hubeny et al. 2003; Fortney et al. 2008). Although thermal inversions have previously been claimed for a number of exoplanets, most of these results were based on Spitzer Space Telescope secondary eclipse measurements since called into question (Diamond-Lowe et al. 2014; Hansen et al. 2014; Evans et al. 2015). Furthermore, while Désert et al. (2008) have suggested a very low abundance of TiO/VO to explain absorption measured in the optical transmission spectrum of HD209458b (T_eq ~ 1500 K), searches in the much hotter WASP-12b (T_eq ~ 2600 K) and WASP-19b (T_eq ~ 2100 K) have failed to detect either species in transmission (Huitson et al. 2013; Sing et al. 2013).

One possibility is that TiO/VO is depleted from the upper atmosphere by cold trapping, either deeper in the dayside atmosphere or on the cooler nightside (Hubeny et al. 2003; Spiegel et al. 2009). The hottest planets, however, might avoid cold traps. Haynes et al. (2015) have recently published near-infrared secondary eclipse observations for the gas giant WASP-33b, which has an equilibrium temperature of T_eq ∼ 2700 K.⁷ These data reveal excess emission around 1.2 μみゅーm, which is consistent with TiO and would imply a thermal inversion. To confirm this hypothesis, further observations are required to spectrally resolve the TiO bandheads.

In this Letter, we present near-infrared transmission spectroscopy observations for another ultra-hot gas giant, WASP-121b (Delrez et al. 2016). This extreme planet is in a 1.3 day polar orbit around an F6V host star, with a semimajor axis of 0.025 au. This puts WASP-121b just beyond the Roche limit, where it is subject to intense tidal forces. With an equilibrium temperature of T_eq ∼ 2400 K, the atmosphere may be hot enough for gaseous TiO/VO to be abundant. WASP-121b is also one of the most inflated planets known, with a mass of 1.2 M_J and a radius of 1.8 R_J. These properties, combined with the brightness of the host star (J = 9.6 mag), make WASP-121b an excellent target for atmospheric characterization.

A single transit of WASP-121b was observed on 2016 February 6 UT using the Hubble Space Telescope (HST) Wide Field Camera 3 (WFC3) for Program 14468 (PI: Evans). Spectroscopic mode was used with grism G141 and forward scanning to allow longer exposures without saturating the detector (McCullough & MacKenty 2012). To reduce overheads, a 256 × 256 subarray containing the target spectrum was sampled with 15 non-destructive reads per exposure. Integration times were 103 s, using the SPARS10 readout mode.⁸ A scan rate of 0.12 arcsec s⁻¹ was adopted, giving scans along 100 pixel rows of the cross-dispersion axis for each exposure. Typical count levels remained below 2.2 × 10⁴ electrons, well within the linear regime of the detector.

Exposures were taken contiguously over five HST orbits, with 17 exposures obtained per orbit. The first orbit allowed the telescope to settle into its new pointing position and was discarded during the light curve analysis due to the large amplitude instrumental systematics it exhibited, as is standard practice for exoplanet transit observations (e.g., Deming et al. 2013; Huitson et al. 2013; Stevenson et al. 2014; Haynes et al. 2015). Of the remaining four HST orbits, the first and fourth provided an out-of-transit baseline, while the second and third occurred during the transit.

Our analysis commenced with the IMA data files produced by the CALWF3 pipeline (v3.1.6), which already have basic calibrations such as flat fielding and bias subtraction applied. We extracted flux for WASP-121 from each exposure by taking the difference between successive non-destructive reads. For each read difference, we removed the background by taking the median flux in a box of pixels well away from the stellar spectra. Typical background levels were 115 electrons. We then determined the flux-weighted center of the WASP-121 scan and set to zero all pixel values located more than 35 pixels above and below along the cross-dispersion axis. Application of this tophat filter had the effect of masking the flux contributions from nearby contaminants, including a faint star separated by 7 arcsec on the sky. It had the additional advantage of eliminating many of the pixels affected by cosmic rays. The final reconstructed images were produced by adding together the read differences for each exposure. We extracted spectra from the reconstructed images by summing the flux within rectangular apertures centered on the scanned spectra and spanning the full dispersion axis. We experimented with aperture widths along the cross-dispersion axis ranging from 100 to 200 pixels in increments of 10 pixels, and found that an aperture of 160 pixels minimized the residuals in the final white light curve fit (Section 4). The wavelength solution was determined by cross-correlating the first spectrum, with dispersion 4.64 nm/pixel, against an ATLAS stellar model (Kurucz 1979; Castelli & Kurucz 2004) with properties similar to WASP-121 (T_eff = 6500 K, log g = 4.0 cgs, v_turb = 2 km s⁻¹) modulated by the throughput of the G141 grism. Prior to cross-correlation, both the measured spectrum and ATLAS spectrum were smoothed using a Gaussian filter with a full width half maximum of 20 nm. As a result, the cross-correlation was most sensitive to the steep edges of the G141 response curve rather than individual stellar lines. We repeated this process for all remaining spectra to determine shifts along the dispersion axis over the course of the observations, which were found to be within 0.23 pixels.

White light curves were generated for each trial aperture by summing the flux for each spectrum along the dispersion axis. We discarded the first exposure from each orbit, as these had significantly lower counts than subsequent exposures. Two additional frames were flagged as outliers, with closer inspection revealing that one was affected by a cosmic ray and the other exhibited an anomalous scan. The resulting light curve obtained for the 160 pixel aperture is shown in Figure 1. Instrumental systematics that correlate with the HST orbital phase are evident, caused by the varying thermal environment experienced by the spacecraft during its orbit (e.g., Sing et al. 2013; Wakeford et al. 2016).

Zoom In Zoom Out Reset image size

To model the white light curves, we adopted a Mandel & Agol (2002) analytic function for the planet signal. For stellar limb darkening, we used the four-parameter nonlinear limb law of Claret (2000) with coefficients c₁, c₂, c₃, and c₄ listed in Table 1. The latter were obtained by fitting to the limb-darkened intensities of the ATLAS stellar model described in Section 3 multiplied by the G141 throughput profile.

Table 1.  White Light Curve Fit Results and Adopted Parameters

Parameter	Value
Rp/R⋆	${0.12109}_{-0.00032}^{+0.00031}$
Tmid (HJDUTC)	${2457424.88307}_{-0.00011}^{+0.00010}$
a/R⋆	3.754
b	0.160
i (°)	87.557
c1	0.582
c2	0.151
c3	−0.435
c4	0.199

Download table as:  ASCIITypeset image

For the white light curve fitting we treated the data as a Gaussian process (GP), using the approach described in Gibson et al. (2012). For the mean function we adopted a Mandel & Agol (2002) transit model multiplied by a linear time trend and for the covariance matrix we used a Matérn νにゅー = 3/2 kernel (see Gibson et al. 2013) with HST orbital phase and dispersion shift as the input variables. We allowed the linear trend and covariance parameters to vary in the fitting, along with the planet-to-star radius ratio R_p/R_⋆ and transit mid-time T_mid. The orbital period P, normalized semimajor axis a/R_⋆, orbital inclination i, and eccentricity e were fixed to values reported in Delrez et al. (2016) and listed in Table 1. Uniform priors were adopted for R_p/R_⋆, T_mid, and the linear trend parameters. Gamma priors of the form Gam(αあるふぁ = 1, βべーた = 100) and Gam(αあるふぁ = 1, βべーた = 1) were adopted for the covariance amplitude and inverse correlation length scales, respectively, giving preference to simpler systematics models, i.e., smaller covariance amplitudes and longer correlation length scales. Maximum likelihood solutions were located using nonlinear optimization, as implemented by the fmin routine of the scipy.optimize Python software package.⁹ The likelihood distribution was then marginalized using affine-invariant Markov chain Monte Carlo (MCMC) as implemented by the emcee software package (Foreman-Mackey et al. 2013). This was done by initializing two groups of walkers in close proximity to the maximum likelihood solution and allowing them to run for 300 steps. For both walker groups, correlation length scales were calculated for each parameter and a burn-in phase of three times the longest correlation length scale was discarded from all walker chains before combining them into a single chain. The resulting posterior samples displayed good mixing and convergence. Best-fit model residuals for the light curve generated with the 160 pixel aperture (Section 3) had the lowest scatter, within 5% of the photon noise floor. We adopt this reduction for all subsequent analysis and report the results in Table 1 with best-fit model shown in Figure 1.

After fitting the white light curve, spectroscopic light curves were produced using a similar approach to Deming et al. (2013). First, a reference spectrum was produced by taking the average of the out-of-transit spectra. Each individual spectrum was then shifted laterally in wavelength and stretched vertically in flux to match the reference spectrum, using linear least squares. The residuals of these fits were binned into 28 spectroscopic channels across the 1.12–1.64 μみゅーm wavelength range, each spanning four columns of the dispersion axis. The spectroscopic residuals were then added to a transit signal with R_p/R_⋆ and T_mid set to the white light curve best-fit values (Table 1) and limb darkening appropriate to the wavelength channel, giving the final spectroscopic light curves shown in Figure 2. This process reduces systematics that are common-mode in wavelength as well as those arising due to the spectra drifting along the dispersion axis. It retains, however, channel-to-channel differences in the flux level, including possible transit depth variations caused by the wavelength-dependent opacity of the planetary atmosphere.

Zoom In Zoom Out Reset image size

We fit the spectroscopic light curves using the same GP treatment described above for the white light curve. For these fits, R_p/R_⋆ was allowed to vary separately for each channel, while T_mid was held fixed to the best-fit white light curve value. Again, a nonlinear limb-darkening law was used with coefficients fixed to the values reported in Table 2, which were determined using the same method described in Section 3 for the white light curve. We report the results for R_p/R_⋆ in Table 2 and show best-fit models in Figure 2.

Table 2.  Inferred Transmission Spectrum and Adopted Nonlinear Limb Darkening Coefficients

λらむだ ( μみゅーm)	Rp/R⋆	c1	c2	c3	c4
0.392–0.481	${0.12375}_{-0.00069}^{+0.00061}$	0.337	0.653	−0.098	−0.041
0.555–0.670	${0.12521}_{-0.00069}^{+0.00065}$	0.465	0.553	−0.408	0.104
0.836–0.943	${0.12298}_{-0.00117}^{+0.00114}$	0.492	0.309	−0.311	0.083
1.121–1.139	${0.12033}_{-0.00038}^{+0.00038}$	0.507	0.201	−0.268	0.086
1.139–1.158	${0.12107}_{-0.00036}^{+0.00035}$	0.501	0.233	−0.318	0.109
1.158–1.177	${0.12122}_{-0.00033}^{+0.00035}$	0.502	0.239	−0.346	0.124
1.177–1.195	${0.12149}_{-0.00034}^{+0.00035}$	0.500	0.255	−0.377	0.139
1.195–1.214	${0.12116}_{-0.00034}^{+0.00034}$	0.497	0.270	−0.404	0.151
1.214–1.232	${0.12113}_{-0.00036}^{+0.00033}$	0.491	0.319	−0.482	0.186
1.232–1.251	${0.12143}_{-0.00035}^{+0.00034}$	0.487	0.338	−0.519	0.205
1.251–1.269	${0.12070}_{-0.00036}^{+0.00037}$	0.496	0.325	−0.518	0.205
1.269–1.288	${0.12036}_{-0.00035}^{+0.00036}$	0.543	0.236	−0.486	0.194
1.288–1.307	${0.12033}_{-0.00037}^{+0.00036}$	0.500	0.343	−0.565	0.227
1.307–1.325	${0.12035}_{-0.00033}^{+0.00033}$	0.499	0.364	−0.620	0.259
1.325–1.344	${0.12118}_{-0.00034}^{+0.00034}$	0.507	0.356	−0.628	0.266
1.344–1.362	${0.12201}_{-0.00034}^{+0.00035}$	0.516	0.342	−0.630	0.271
1.362–1.381	${0.12193}_{-0.00035}^{+0.00034}$	0.528	0.314	−0.609	0.264
1.381–1.399	${0.12138}_{-0.00036}^{+0.00035}$	0.534	0.327	−0.659	0.293
1.399–1.418	${0.12118}_{-0.00035}^{+0.00037}$	0.549	0.293	−0.634	0.286
1.418–1.437	${0.12156}_{-0.00036}^{+0.00035}$	0.567	0.248	−0.601	0.278
1.437–1.455	${0.12164}_{-0.00036}^{+0.00035}$	0.590	0.180	−0.534	0.252
1.455–1.474	${0.12167}_{-0.00036}^{+0.00035}$	0.608	0.154	−0.521	0.249
1.474–1.492	${0.12088}_{-0.00040}^{+0.00039}$	0.622	0.107	−0.474	0.232
1.492–1.511	${0.12160}_{-0.00037}^{+0.00037}$	0.634	0.082	−0.459	0.229
1.511–1.529	${0.12101}_{-0.00039}^{+0.00039}$	0.667	−0.002	−0.388	0.209
1.529–1.548	${0.12053}_{-0.00039}^{+0.00042}$	0.690	−0.058	−0.352	0.202
1.548–1.567	${0.12123}_{-0.00042}^{+0.00044}$	0.715	−0.098	−0.334	0.199
1.567–1.585	${0.12057}_{-0.00041}^{+0.00041}$	0.723	−0.147	−0.276	0.176
1.585–1.604	${0.12030}_{-0.00041}^{+0.00042}$	0.723	−0.170	−0.241	0.161
1.604–1.622	${0.12071}_{-0.00041}^{+0.00042}$	0.751	−0.227	−0.201	0.150
1.622–1.641	${0.11954}_{-0.00043}^{+0.00042}$	0.754	−0.214	−0.214	0.154

Download table as:  ASCIITypeset image

Following Wakeford et al. (2016), we also performed an independent analysis of the white and spectroscopic light curves using the method of Gibson (2014) in which marginalization was performed over a grid of parametric systematics models. This gave results in good agreement with the GP analysis.

In addition to the spectroscopic WFC3 data, we fit the ground-based photometric transit data presented in Delrez et al. (2016), comprised of three B band light curves, two r' light curves, and four z' light curves. For all ground-based light curves we used a GP systematics model with a transit signal multiplied by a linear time trend as the mean function. For the covariance matrix we adopted a Matérn νにゅー = 3/2 kernel with multiple input variables such as the xy coordinates of the star on the detector, seeing, and airmass. Light curves in the same bandpass were fit simultaneously, with R_p/R_⋆ shared but T_mid, linear trends, covariance parameters, and white noise levels allowed to vary separately. Nonlinear limb-darkening laws were adopted with coefficients fixed to the values reported in Table 2, which were taken from the catalogs of Claret (2000, 2003). Remaining transit parameters were fixed to the values adopted for the WFC3 analysis, namely those reported by Delrez et al. (2016) and listed in Table 1. Fitting was performed by locating the maximum likelihood and marginalizing over the parameter space with affine-invariant MCMC, as described above. Results are reported in Table 2.

The measured transmission spectrum for WASP-121b is plotted in Figure 3. The WFC3 data show structure, with two broad features centered around 1.2 and 1.4 μみゅーm. The scale of these variations is on the order of a few atmospheric scale heights, consistent with expectations for molecular absorption. Additionally, the effective planetary radius is found to be significantly larger at optical wavelengths relative to near-infrared wavelengths.

Zoom In Zoom Out Reset image size

Before exploring implications for the planetary atmosphere, we considered the possibility of stellar activity causing the different radii measured at optical and near-infrared wavelengths. Making the conservative assumption that star spots are non-luminous, we calculate that the unocculted spot coverage would need to decrease by an amount equivalent to 7% of the visible stellar disk between the optical and near-infrared observation epochs. This seems highly unlikely, especially given that WASP-121 appears photometrically stable at the millimag level (Delrez et al. 2016). Also, analyses of the individual light curves for each bandpass produced consistent results at the 1σしぐま level (top panel of Figure 3), suggesting that stellar activity does not affect the inferred radius significantly from epoch to epoch. We therefore conclude that variations in the inferred radius shown in Figure 3 are due to the planetary atmosphere.

We used the one-dimensional radiative transfer code ATMO to generate model spectra and investigate atmospheric opacity sources (Amundsen et al. 2014; Tremblin et al. 2015, 2016). Isothermal temperature-pressure profiles were assumed, with the temperature allowed to vary in the fitting. We produced models with and without TiO, VO, and FeH. The latter species are observed in M/L dwarfs (Burrows et al. 2001) and models predict they may also be present in hot gas giant planets (e.g., Sharp & Burrows 2007). For models including TiO/VO/FeH, we varied the relative abundances of TiO/H₂O, VO/H₂O, and Fe/H₂O during fitting. Solar abundances under thermochemical equilibrium were adopted for other gas species, including H₂O. We note that for transmission spectra the absolute abundances of individual species are often unconstrained owing to a well-known degeneracy with the absolute pressure level, but the relative abundances between species can be well-measured (Lecavelier Des Etangs et al. 2008). In our fitting, TiO/VO/FeH abundances were adjusted using the analytical relation of Lecavelier Des Etangs et al. (2008) for the wavelength-dependent transit radius of an isothermal atmosphere. We also considered a simple haze-scattering model in the form of a Rayleigh profile with collisional cross-section area allowed to vary as a free parameter as described in Sing et al. (2016).

Best fits for a number of illustrative models are shown in Figure 3, that were obtained using the nonlinear least squares IDL routine mpfit (Markwardt 2009). We find that a model excluding TiO, VO, FeH, and haze is unable to account for the data (the purple line). In particular, it fails to reproduce the optical photometry and short wavelength WFC3 data. A good fit is achieved, however, when TiO, VO, and FeH opacity is included (the red line). For this model, we find a temperature of 1500 ± 230 K and abundances relative to thermochemical equilibrium of ∼7× solar for TiO/H₂O, ∼5× solar for VO/H₂O, and ∼0.2× solar for FeH/H₂O. The inferred temperature is lower than the dayside equilibrium temperature of T_eq ∼ 2400 K, which may indicate a significantly cooler nightside hemisphere. The strong absorption of TiO and VO in the optical accounts for the larger radii measured at these wavelengths relative to the near-infrared. To illustrate the effect of FeH, we also generated a model with it removed, but TiO and VO retained (the orange line). Such a model is unable to reproduce the WFC3 transmission spectrum at wavelengths near 1.3 μみゅーm, where FeH has prominent absorption. Lastly, a model excluding TiO, VO, and FeH, but including haze, gives a relatively poor fit to the data (the green line). Most significantly, haze scattering underpredicts the r' opacity by 2.6σしぐま. To assess the relative quality of fit between the haze and TiO/VO/FeH models, we computed the Bayesian information criterion (BIC) for each, obtaining BIC = 50.4 for the haze model and BIC = 42.6 for the TiO/VO/FeH model. Following Kass & Raftery (1995) we approximate the Bayes factor as $\mathrm{exp}(-{\rm{\Delta }}{\rm{BIC}}/2)\;=\;0.02$ in favor of the TiO/VO/FeH model over the haze model. This strongly supports the hypothesis that TiO and VO rather than haze are responsible for the enhanced opacity measured at optical wavelengths.

All models favor H₂O absorption, which matches the measured spectrum well between 1.3 and 1.64 μみゅーm. To quantify the significance of this detection, we compared the quality of fit obtained for the best-fit TiO/VO/FeH model described above with that obtained for the same model but H₂O removed. We determine that H₂O is detected at a confidence level of 5.4σしぐま.

Our observations reveal H₂O absorption in the atmosphere of WASP-121b with an amplitude of >2 gas pressure scale heights. This result is consistent with a clear atmosphere, assuming TiO/VO and FeH absorption is responsible for the additional opacity measured across the 1.12–1.3 μみゅーm wavelength range. In contrast, WASP-12b is observed to have a hazy atmosphere without TiO/VO (Sing et al. 2013), despite having similar properties to WASP-121b, reinforcing the emerging picture of hot gas giant diversity (Sing et al. 2016). We note, however, that the optical data cannot currently exclude a model including both TiO/VO absorption and haze scattering for WASP-121b. Such a model would be very similar to our best-fit model (the red line in Figure 3), except that the effective radius would continue increasing for wavelengths shortward of the B bandpass. To confidently discount such a scenario, it will be necessary to spectrally resolve the characteristic steep rise in opacity at ∼0.35 μみゅーm due to TiO/VO with further observations.

Until then, the evidence for TiO/VO absorption in the atmosphere of WASP-121b remains tentative, as it primarily hinges upon the relative radius measured from the ground in the r' bandpass. We emphasize, however, that this measurement is based on two light curves that produce consistent results when analyzed separately and simultaneously (the top panel of Figure 3). In addition, Delrez et al. (2016) report ground-based secondary eclipse measurements in the z' bandpass that imply a brightness temperature higher than the equilibrium temperature predicted for zero Bond albedo and instantaneous heat re-radiation. This could potentially be explained by a thermal inversion, with TiO being observed in emission in the z' bandpass. Again, follow-up observations that spectrally resolve the TiO/VO features will be necessary for confirmation.

To date, HST spectroscopy data have been published for three exoplanets with T_eq > 2000 K. Of these, TiO/VO absorption has been ruled out for WASP-12b (T_eq ∼ 2600 K; Sing et al. 2013) and WASP-19b (T_eq ∼ 2100 K; Huitson et al. 2013), while evidence for TiO emission has been reported for WASP-33b (T_eq ∼ 2700 K; Haynes et al. 2015). It is possible that WASP-33b is hot enough to maintain TiO/VO in the gas phase throughout its dayside atmosphere, whereas the day–night terminators of WASP-12b and WASP-19b are too cool, or cold trapping occurs elsewhere in their atmospheres. However, the fact that WASP-121b has a lower equilibrium temperature (T_eq ∼ 2400 K) than WASP-12b, and yet shows evidence for TiO/VO absorption at low pressures, indicates that additional factors are at play. For instance, models predict larger temperature contrasts between the dayside and nightside hemispheres with increasing stellar irradiation and planet metallicity (Dobbs-Dixon & Lin 2008; Fortney et al. 2008; Kataria et al. 2015). It is thus plausible that WASP-12b, with its higher irradiation and host star metallicity, has a cooler nightside than WASP-121b, which could result in a nightside cold trap for WASP-12b but not WASP-121b. Compositional differences could also have a significant influence on the temperature profiles for each planet (Parmentier & Guillot 2014; Tremblin et al. 2016), which in turn determines the susceptibility to cold trapping.

We presented a transmission spectrum for WASP-121b spanning the 1.12–1.64 μみゅーm wavelength range. Absorption by H₂O is detected at high confidence (5.4σしぐま), consistent with predictions of clear atmosphere models. Deeper transits measured at optical wavelengths relative to the near-infrared strongly favor models including TiO/VO absorption, but scattering by a high-altitude haze cannot yet be definitively excluded. We also find evidence for FeH absorption in the WFC3 bandpass. WASP-121b is one of the most favorable targets available for both transmission and emission spectroscopy and offers particular promise for exploring the link between strong optical absorbers, such as TiO/VO, and thermal inversions in hot gas giant atmospheres.

The authors would like to thank the referee for their prompt and thoughtful review. This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program GO-14468. The authors are grateful to the WASP-121 discovery team for generously providing the ground-based photometric light curves. The research leading to these results received funding from the European Research Council under the European Union Seventh Framework Program (FP7/2007-2013) ERC grant agreement no. 336792. H.R.W. acknowledges support by an appointment to the NASA Postdoctoral Program at Goddard Space Flight Center, administered by ORAU and USRA through a contract with NASA. N.P.G. gratefully acknowledges support from the Royal Society in the form of a University Research Fellowship.