COMPOSITIONAL DIVERSITY IN THE ATMOSPHERES OF HOT NEPTUNES, WITH APPLICATION TO GJ 436b

Two chemical models are used in this study. The first is a chemical-equilibrium model using the NASA CEA code of Gordon & McBride (1994), and the second is a one-dimensional (1D) thermochemical and photochemical kinetics and transport model based on the Caltech/JPL KINETICS code of Allen et al. (1981). The kinetics/transport model is described more fully in Moses et al. (2011), Visscher & Moses (2011), and Moses et al. (2013). For the equilibrium calculations, we consider ∼500 gas-phase species and condensates containing the elements H, He, C, N, O, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, K, Ca, Ti, Cr, Mn, Fe, and Ni. For the kinetics/transport calculations, we solve the continuity equations for 92 atmospheric species via ∼1600 forward and reverse chemical-reaction pairs. Only species containing the elements H, He, C, N, and O are considered in the kinetics/transport model due to a lack of key rate-coefficient data for species containing the other elements. Note that the included elements are the dominant ones that will not be sequestered within condensates in the atmosphere of GJ 436b, and thus the GJ 436b results for the major gas-phase species are not expected to change significantly with the inclusion of additional elements.

The thermodynamic principle of microscopic reversibility, which is expected to be accurate even for complex multiple-potential-well chemical reactions (Miller et al. 2009), is assumed in the kinetics/transport model. Chemical equilibrium is achieved kinetically in the deepest, hottest regions of these models, but chemical reactions tend to be slower in the upper, cooler regions. When transport time scales drop below kinetic conversion time scales—at a pressure called the "quench pressure" or "quench level"—transport begins to dominate over chemical kinetics in controlling the species' vertical profiles (e.g., Prinn & Barshay 1977; Lewis & Fegley 1984; Fegley & Lodders 1994; Moses et al. 2010; Visscher & Moses 2011). When this situation occurs, the species can be "quenched" at mole fractions that remain constant with altitude above the quench level (and thus diverge from chemical-equilibrium predictions), as long as transport times scales remain shorter than chemical-kinetics time scales. Different assumptions about the rate coefficients for some key reactions can lead to differences in the predicted abundances of the quenched species (Visscher & Moses 2011; Venot et al. 2012); some of the current assumptions in published exoplanet models are reviewed by Moses (2013). The abundances of the carbon- and nitrogen-bearing species are particularly affected when kinetic interconversion between the dominant carbon constituents, CO and CH₄, and dominant nitrogen constituents, N₂ and NH₃, ceases to be effective.

Changes to the Moses et al. (2011) chemical mechanism, other than what is discussed in Moses et al. (2013), include the addition of O₃ and related reactions, which could potentially become important as the metallicity is increased. The rate coefficients and cross sections for these reactions are adopted from terrestrial atmospheric-chemistry compilations (e.g., Atkinson et al. 2004, 2006; Sander et al. 2011) and Venus photochemistry studies (e.g., Mills 1998). For the GJ 436b models discussed below, ozone does not build up to observable levels, and the added reactions do not measurably influence the observable species. Ions and ion chemistry are not included in the models, and our solutions at the highest altitudes in the model will be unrealistic due to our neglect of ion chemistry, high-temperature thermospheres, and possible hydrodynamic escape (e.g., García Muñoz 2007). The models contain 198 vertical levels separated uniformly in log pressure, with the hydrostatic equilibrium equation being used to solve for the background atmospheric parameters along the vertical grid. The chemical-equilibrium abundance profiles from the CEA code are adopted as initial conditions in the kinetics and transport models, with zero flux boundary conditions being assumed for all species at the top and bottom of the model. Our assumed solar composition is taken from the protosolar abundances listed in Table 10 of Lodders et al. (2009). Multiple Rayleigh scattering of incoming stellar radiation by gases is considered in the kinetics/transport models, but we assume that aerosols do not contribute to the atmospheric extinction.

We assume that vertical transport occurs through molecular and "eddy" diffusion, with the eddy diffusion coefficients K_zz being free parameters in the model. In the deep, convective portion of the atmosphere, free-convection and mixing-length theories (e.g., Stone 1976) predict relatively large eddy diffusion coefficients and short mixing time scales (e.g., K_zz ≈ 10⁹ cm²  s⁻¹ for the convective regions in GJ 436b), but K_zz values tend to be much smaller in the radiative regions in the upper troposphere and lower stratosphere. Higher up in the stratosphere, turbulence due to atmospheric tides and upward-propagating gravity waves (non-breaking as well as breaking waves) is expected to cause effective K_zz values to increase roughly with the inverse square root of atmospheric pressure (e.g., Lindzen 1981, and references therein)—a scaling that appears consistent with inferred vertical mixing in exoplanet general circulation models (GCMs; Parmentier et al. 2013; cf. also the K_zz profiles inferred from Showman et al. 2009, as shown in Moses et al. 2011). In our kinetics and transport models, the K_zz profile influences the quench behavior of molecules like CO, CH₄, and NH₃, and the observations themselves may ultimately provide the best means for defining both the K_zz values at the quench levels and the pressures at which those quench levels occur (e.g., Fegley & Lodders 1994; Bézard et al. 2002; Visscher et al. 2010b; Moses et al. 2010; Visscher & Moses 2011). However, the broadband photometric eclipse observations obtained to date for GJ 436b are not sufficient to constrain either the composition or thermal profile accurately enough to derive K_zz values in such a manner at this time. We therefore treat K_zz as a free parameter and will explicitly specify our assumptions for each model presented in Section 3.3.

Both the chemistry and the predicted spectrum of our hot-Neptune exoplanets will depend strongly on the adopted thermal structure. We do not self-consistently calculate temperatures within the kinetics/transport models. For our generic hot Neptunes, we explore a wide range of temperature–pressure conditions. For GJ 436b, we consider a variety of theoretically derived thermal profiles, including the dayside-average profiles at conditions of secondary eclipse from the 1 times and 50 times solar-metallicity GCMs of Lewis et al. (2010), as shown in Figure 1. For higher-metallicity GJ 436b scenarios, we use the PHOENIX atmospheric model (Hauschildt et al. 1999; Allard et al. 2001) in the presence of an external radiation field, as described in Barman et al. (2001) and Barman et al. (2005), to compute 1D temperature–pressure profiles under the assumption of inefficient day–night heat redistribution and efficient gravitational settling of condensates. The resulting profile for the case of 1000 times solar metallicity is shown in Figure 1. For all the non-solar-metallicity thermal models, solar ratios of elements other than hydrogen and helium are scaled by a constant factor. In certain instances, we also adopt thermal profiles from Madhusudhan & Seager (2011) that provide a good fit to the Stevenson et al. (2010) GJ 436b secondary-eclipse data, or we adopt profiles derived from new retrievals such as are described in Line et al. (2012, 2013; see also Section 3.2). We typically extend the thermal profiles upward in altitude nearly isothermally to a pressure of 10⁻¹¹ bar, where all major UV absorbers are optically thin, and downward assuming an adiabat to a lower boundary that reaches at least 2400 K, to ensure that the N₂–NH₃ quench levels are contained within the models (i.e., to encompass the pressure at which kinetic interconversion between N₂ and NH₃ slows down enough to prevent equilibrium from being maintained, see Moses et al. 2010, 2011). The adopted thermal profiles will be clearly described when we discuss the different GJ 436b models.

Zoom In Zoom Out Reset image size

Although the thermal profiles derived from the three-dimensional (3D) GCM and the 1D PHOENIX models have some notable differences for any given metallicity (not shown in Figure 1), both models predict a general upward vertical shift in the thermal profile as the metallicity is increased above solar. An increased metallicity results in increased mole fractions of opacity sources like water and other key molecular absorbers at lower pressures, which moves the optical-depth unity level (and "photosphere" in general) upward to higher altitudes. However, because the overall atmospheric number density decreases with increasing altitude, there is a limit to how high the photosphere can be shifted upward in altitude from this increased metallicity (see also Lewis et al. 2010; Miller-Ricci & Fortney 2010) because optical thickness cannot be approached at very high altitudes. As an example, the PHOENIX-based models show virtually no differences in the thermal profiles between 1000 times and 10,000 times solar metallicities. The chemical equilibrium CH₄-CO equal-abundance boundary also shifts downward in altitude with increasing metallicity (see Figure 1).

These two metallicity-dependent effects, first noted by Lodders & Fegley (2002) and discussed in relation to GJ 436b by Lewis et al. (2010), motivate our investigation and demonstrate why higher-metallicity models are more likely to explain the inferred large CO/CH₄ ratio needed to reproduce the secondary-eclipse data from GJ 436b (Stevenson et al. 2010; Madhusudhan & Seager 2011). For instance, Figure 1 illustrates that the thermal profile for the 1 times solar-metallicity model lies completely within the CH₄-dominated regime, so that no matter what the rate of vertical mixing or at what pressure the CO quenches, the carbon monoxide mole fraction will never exceed the CH₄ mole fraction for this solar-metallicity thermal profile. Because the dayside atmosphere is expected to be hotter than the terminators, the likelihood of having CO dominate at conditions relevant to the transit is even smaller if the atmosphere has a solar-like metallicity. In contrast, the thermal structure of the 1000 times solar model shown in Figure 1 lies completely within the CO-dominated regime, down to at least the 30 bar level. Since the CO–CH₄ quench point (i.e., where transport processes dominate over the chemical-kinetic interconversion of CO and CH₄) is likely to be within the ∼0.1–30 bar region on GJ 436b (Moses et al. 2011; Visscher & Moses 2011; Line et al. 2011), the quenching will occur in the CO-dominated regime for the 1000 times solar-metallicity model, and carbon monoxide will dominate over methane.

In general, the hotter the photosphere of GJ 436b, the more abundant CO is likely to be. Aside from the effect of increased metallicities, hotter photospheres can result from the addition of tidal heating or a large residual internal heat source. GJ 436 is a relatively quiet star and slow rotator (Saffe et al. 2005; Demory et al. 2007; Jenkins et al. 2009; Sanz-Forcada et al. 2010; Knutson et al. 2010, 2011), suggesting that the system is relatively old (e.g., 6$^{+4} _{-5}$ Gyr, according to Torres et al. 2008). For an older planet of GJ 436b's mass, the interior would be expected to have cooled significantly such that T_int ≈ 60 K for GJ 436b (comparable to that of Neptune with its T_int ≈ 50 K), defined in terms of an intrinsic internal heat flux of $\sigma T_{{\rm int}}^4$ (e.g., Fortney et al. 2007; Marley et al. 2007; Baraffe et al. 2008; Rogers & Seager 2010a). However, the orbit of GJ 436b has a significant eccentricity of 0.146 (von Braun et al. 2012), suggesting that the atmosphere is being tidally heated. Since orbital circulation times due to tidal dissipation are of order 30 Myr for GJ 436b (Deming et al. 2007), the eccentricity is likely being continually forced by one or more additional planets in the system (e.g., Deming et al. 2007; Demory et al. 2007; Stevenson et al. 2012). Depending on where the tidal energy is dissipated within the planet, the additional tidal heating could increase temperatures at the CO–CH₄ quench point, pushing the thermal profile into the CO stability field. Such a situation is shown in Figure 2 for a deep-seated intrinsic heat source on a 50 times solar-metallicity GJ 436b.

Zoom In Zoom Out Reset image size

Note, however, that even the highest T_int profile shown in Figure 2 still resides close enough to the CO = CH₄ curve in the region from a few tens of bars to a few tenths of a bar that the quenched CH₄ abundance will be relatively large, even if CO dominates. For example, for the most favorable case of T_int = 300 K with vigorous mixing (e.g., K_zz = 10¹¹ cm² s⁻¹), which would result a quench point at relatively high temperatures and pressures, the quenched CH₄ mole fraction is 4 × 10⁻⁵, which is a factor or 40 larger than the derived 1 ppm upper limit from Spitzer eclipse observations Madhusudhan & Seager (2011). A higher metallicity and a large internal heat source would make higher CO/CH₄ ratios more likely.

Another input to the kinetics/transport models is the ultraviolet flux from the host star, which has been classified in the literature as an early M dwarf of spectral type M2.5V to M3.5V (Kirkpatrick et al. 1991; Hawley et al. 1996; Maness et al. 2007; Jenkins et al. 2009). Figure 3 illustrates our adopted flux for GJ 436 (normalized to 1 AU), as derived from a compilation of various observational and theoretical sources. At wavelengths from 0 through the Lyman beta line at 1025.7 Å, we use the GJ 436 synthetic spectrum from the X-exoplanets archive at the Centro de Astrobiología (Sanz-Forcada et al. 2011), assuming the units in the downloadable file are photons cm⁻² s⁻¹ per wavelength bin, which then reproduces the correct integrated X-ray and EUV luminosity from Table 6 of Sanz-Forcada et al. (2011). At wavelengths between Lyman beta and the Lyman alpha line at 1215.7 Å, we use the solar spectrum of Woods & Rottman (2002) at solar-cycle minimum, scaled downward by a factor of 6.3 to transition smoothly to the longer-wavelength flux. For the Lyman alpha line itself, we use the reconstructed GJ 436 flux from Ehrenreich et al. (2011). For wavelengths between 1215.7 Å and 1800 Å, we use International Ultraviolet Explorer (IUE) data for GL 15B (a M3.5V star) from the MAST archive (http://archive.stsci.edu), and for wavelengths longer than 1800 Å, we use the spectrum of GL 15B from the Next Generation Spectral Library (Heap & Lindler 2010). For our dayside atmospheric models, we scale this normalized 1 AU spectrum to the 0.027 AU orbital distance relevant to the GJ 436b secondary eclipse. The key spectral region as far as the neutral atmospheric chemistry is concerned is the wavelength region from Lyman alpha out to ∼2400 Å. We use a fixed solar zenith angle θ = 48° to simulate the eclipse conditions in the calculations (see Moses et al. 2011). A directly measured ultraviolet spectrum for GJ 436 has recently been made available by France et al. (2013), and we emphasize that such studies are of great utility to theoretical photochemical models for exoplanets.

Zoom In Zoom Out Reset image size

To calculate the emergent planetary spectrum for GJ 436b from the assumed thermal structure and derived abundance profiles from the chemical models, we use a plane-parallel radiative-transfer code as described in Line et al. (2013). Opacity from H₂, He, H₂O, CO, CO₂, CH₄, and NH₃ is included in the modeling, although NH₃ itself is not included in the retrievals. The source line parameters are described in Line et al. (2013). Collision-induced absorption from H₂–H₂ and H₂–He is included in the models, but collision-induced absorption from other relevant collisional pairs such as H₂–CO, CO–CO, H₂O–CO is not. As atmospheric metallicity is increased, molecules like CO₂, CO, and H₂O increase in importance relative to H₂, and a future study of the full spectral effects of the changes in collision-induced absorption and overall continuum as the composition shifts would be interesting. Local thermodynamic equilibrium is assumed throughout, and scattering is ignored. In a manner similar to Sharp & Burrows (2007), absorption cross sections are precomputed using a line-by-line code at high spectral resolution and tabulated on a grid with 1 cm⁻¹ wavenumber resolution and 20 evenly spaced temperature points between 500–3000 K and log(pressure) points between 50–10⁻⁶ bar (see http://www.atm.ox.ac.uk/RFM/). The cross sections from this pre-tabulated grid can then be interpolated for the pressure–temperature-abundance conditions relevant to the model atmosphere, which has 90 grid layers equally spaced between 50–10⁻⁶ bar. As is typical, we plot the model emission spectrum in terms of the flux of the planet divided by the flux of the star. The flux emergent from the entire disk of the planet is calculated using a four-point gaussian quadrature. The stellar flux is derived from a PHOENIX model (Hauschildt et al. 1999; Allard et al. 2001) assuming T_eff = 3350 K (Maness et al. 2007), although we should note that the latest determinations of the stellar effective temperature from von Braun et al. (2012) are somewhat hotter at T_eff = 3416 ± 54 K (see also Bean et al. 2006; Nettelmann et al. 2010; Southworth 2010).

With the core-accretion model of giant-planet formation (Mizuno et al. 1978; Bodenheimer & Pollack 1986; Pollack et al. 1996), a rocky or rock-ice protoplanetary core initially forms and grows from the accretion of solid planetesimals within a protoplanetary disk, with gaseous envelopes forming around these emerging protoplanets through core outgassing, direct accretion of nebular gas, and ablation of incoming solid planetesimals within the gaseous envelope. Rapid accretion of the surrounding largely H₂ and He nebular gas occurs when the protoplanet reaches a certain critical mass (∼10 M_⊕ in these traditional core-accretion models). The standard explanation for the "ice giants" Uranus and Neptune being so much more enriched in heavy elements than Jupiter and Saturn is that the accretion rate of solids was slow enough for the proto-Uranus and -Neptune that they did not completely reach the runaway gas-accretion-phase before the nebular gas was dispersed from the disk (e.g., Lissauer & Stevenson 2007). Elements heavier than hydrogen and helium therefore make up 80%–85% of Uranus and Neptune by mass (Podolak et al. 1995; Hubbard et al. 1995; Fortney & Nettelmann 2010). Depending on the degree of mixing between the core and atmosphere, the outermost gaseous envelope could have a heavy-element enrichment by mass (Z_env) different from that of the bulk planet (e.g., Nettelmann et al. 2013); however, observations of CH₄ on the giant planets support the picture of a greater atmospheric Z_env for the ice giants Uranus and Neptune in comparison with the gas giants Jupiter and Saturn (e.g., Moses et al. 2004; Fouchet et al. 2009; Fegley et al. 1991; Gautier et al. 1995; Karkoschka & Tomasko 2009, 2011).

A metal-enriched atmospheric envelope might also be true for hot-Neptune exoplanets, despite potentially widely different evolutionary and migration histories. Although the hot Neptunes for which we have good constraints on both mass and radius show a large scatter in the mass-versus-radius relation, indicating likely different compositions, Weiss et al. (2013) demonstrate that there is a general trend toward increasing density with decreasing mass for planetary masses below ∼150 M_⊕ (see also Miller & Fortney 2011). This trend suggests that smaller planets become progressively more enriched in heavy elements, similar to the situation in our own solar system. Planet-formation and population-synthesis models (e.g., Alibert et al. 2005; Figueira et al. 2009; Mordasini et al. 2012a, 2012b) also reflect this trend, with Fortney et al. (2013) predicting a significant increase in Z_env of Neptune-mass planets in comparison with planets more massive than ∼100 M_⊕. There is also the possibility of in situ formation and capture of gas at small radial distances (e.g., Hansen & Murray 2012), which can lead to a variety of bulk gas fractions depending on planetary size and orbital distance.

To predict the atmospheric composition of any particular hot Neptune, we would need to know the properties of the protoplanetary disk in which the planet formed, the formation location within the disk, the planet's migration and impact history, and the details of the subsequent atmospheric evolution (e.g., interior outgassing, atmospheric escape, impact delivery/erosion, irradiation history, tidal heating, climate evolution, magnetospheric interactions, disequilibrium chemistry, etc.). Given the stochastic nature of some of these evolutionary processes and a lack of information about others, the task of predicting any particular atmospheric composition is exceedingly difficult (although the attempt can still be made, e.g., Alibert et al. 2006; Mousis et al. 2011; Madhusudhan et al. 2011b; moreover, population-synthesis models along the lines of those mentioned above can provide valuable insights into atmospheric properties of the ensemble). Instead of pursuing these types of models, we go through the simple exercise of investigating the expected equilibrium composition of Neptune-class exoplanets as a function of temperature, bulk metallicity, and bulk elemental ratios in the atmosphere. The increase in metallicity in these calculations then becomes a proxy for the evolution of smaller Neptune-mass planets, which are more likely to have high Z_env due to either inefficient gas accretion or efficient hydrogen escape.

Figure 4 illustrates how the atmospheric composition is expected to change as a function of either the atmospheric metallicity or the atmospheric C/O ratio for a few different temperatures at a pressure of 100 mbar. The 100 mbar pressure was selected here because it represents a typical infrared photospheric emission region in transiting-planet atmospheres. For these calculations, we use the protosolar abundances of Lodders et al. (2009) and assume an increase in metallicity occurs uniformly for all species except H, He, and Ne; moreover, we define the oxygen abundance through the C/O ratio, so the "metallicity" in this context refers to the C/H ratio (or X/H ratio, where X is any species except H, He, Ne, or O). This "C/H metallicity" should not be confused with bulk atmospheric metallicity, as the overall heavy-element enrichment will change with the C/O ratio—at very low C/O ratios, for example, the overall atmospheric metallicity, as defined by the abundance of C + O + all heavy elements in relation to that of the Sun, can be quite a bit higher than the C/H metallicity due to the exaggerated importance of the enhanced O. Condensates are not assumed to rain out in this equilibrium model. Given such simplifications and given that solar ratios of elements may not be preserved as a planetary atmosphere forms and evolves, the relevance of these calculations to real planets is questionable, but the exercise does effectively demonstrate that the atmospheric composition of hot Neptunes could be highly variable with atmospheric properties.

Zoom In Zoom Out Reset image size

Exo-Neptunes with solar-like elemental ratios and moderately low metallicities could have hydrogen-dominated atmospheres very reminiscent of our own Neptune, but water and methane will make up an increasing fraction of the atmosphere with increasing metallicity, until H₂ itself becomes less important. At high-enough metallicities, exo-Neptunes can have CO₂-dominated atmospheres, qualitatively reminiscent of Venus. Carbon monoxide becomes an increasingly important constituent at higher temperatures and higher metallicities, and even O₂ can become dominant in equilibrium at very high metallicities and very low C/O ratios, further emphasizing that O₂ by itself is not necessarily a good indicator of biological activity on exoplanets (e.g., see Selsis et al. 2002; Segura et al. 2007; Schaefer et al. 2012; Hu et al. 2012).

These scenarios are borne out more clearly in Figure 5, which shows how the equilibrium abundances of several atmospheric species vary as a function of temperature and metallicity for an assumed pressure of 100 mbar and an assumed C/O ratio that is maintained at the protosolar value (C/O = 0.46) with an increase in metallicity (i.e., all elements except H, He, and Ne are increased by a constant factor). Because we do not have a complete set of condensates in these equilibrium calculations (and in particular, we are missing several Ca–Ti–Al silicates, see Lodders 2010), these gas-phase mixing ratios are not completely accurate, but the general trends with temperature and metallicity hold true. For example, Figure 5 demonstrates that H₂ becomes relatively less abundant at all temperatures as the metallicity is increased, whereas the water (H₂O) mole fraction increases roughly linearly with increasing metallicity at most temperatures, until H₂O starts to decrease at very high metallicities (e.g., ≳1000 times solar) due to the overall decrease in the bulk H fraction. Methane (CH₄) is more stable at low temperatures, where its abundance also increases with increasing metallicity until the bulk hydrogen abundance drops and makes less H available to form methane. Carbon monoxide (CO) becomes more stable at high temperatures and high metallicities, where it is a significant atmospheric component under these conditions. Carbon dioxide (CO₂) is not very abundant at low metallicities, especially at low temperatures, but CO₂ increases significantly with increasing metallicity, becoming the dominant constituent at very high metallicities for all temperatures. Carbonyl sulfide (OCS) has a behavior similar to CO₂, although it is never as abundant. Solid graphite is stable in the lower-right corner of these plots in Figure 5 (high metallicities, low temperatures) and sequesters a notable fraction of the available carbon, thereby lowering the gas-phase C/O ratio in this region of temperature–metallicity space. Ammonia (NH₃) is abundant in equilibrium at the lowest temperatures considered, but N₂ becomes the most abundant nitrogen component at all temperatures as the metallicity is increased, although it never dominates the overall atmospheric composition at the conditions considered, assuming no further atmospheric evolution. The abundance of hydrogen cyanide (HCN; not shown in Figure 5) never rivals that of N₂ and NH₃, although its abundance increases with increasing temperature at moderately high metallicities. Hydrogen sulfide (H₂S) is the dominant sulfur constituent under the considered conditions, increasingly roughly linearly with metallicity until hydrogen becomes scarce.

Zoom In Zoom Out Reset image size

Note that atmospheric metallicities of 10, 100, 1000, and 10,000 times solar correspond to Z_env = 0.13, 0.61, 0.94, and 0.99, respectively, for solar ratios of the full suite of elements considered in the models (see Section 2.1). This range of Z_env variation is expected for Neptune-sized planets in the population-synthesis formation and evolution models of Fortney et al. (2013), with Z_env = 0.6–0.9 being the most common for such planets. We might therefore expect the whole range of possible phase space in Figure 5 to be represented in the hot-Neptune exoplanet population. Although Uranus and Neptune are good examples of the cool, moderate-metallicity, H₂-dominated end members, many hot-Neptune exoplanets would have atmospheric compositions that are not found in our own solar system, with various H₂O-, CO-, or CO₂-rich possibilities being particularly worth mentioning. We should also note that metallicities ≳500–600× solar are not possible if the "metals" are brought in predominantly via water or very water-rich volatiles (and/or through H-rich species such as CH₄, NH₃, and H₂S)—assuming no further atmospheric evolution due to hydrogen escape or other fractionation processes—because such volatiles would also deliver large amounts of hydrogen to the atmosphere (Nettelmann et al. 2011).

Aside from bulk metallicity, the atmospheric composition is also dependent on the C/O ratio (Lodders & Fegley 1997; Seager et al. 2005; Kuchner & Seager 2005; Lodders 2010; Line et al. 2010; Madhusudhan et al. 2011a; Madhusudhan 2012; Kopparapu et al. 2012; Moses et al. 2013). It is unclear what C/O ratio to expect for hot Neptunes, as that will depend greatly on the planet's original formation location (especially in relation to specific condensation fronts in the protoplanetary disk), as well as to its evolutionary history (e.g., Moses et al. 2013, and references therein). Atmospheric C/O ratios <0.1 could occur if the disk or feeding zones were oxygen-rich or if planetesimals composed largely of water ice or clathrate hydrates dominated the delivery of heavy elements to the protoplanetary envelope. Similarly, the bulk C/O ratio could conceivably be greater than unity on giant planets (e.g., Lodders 2004; Madhusudhan et al. 2011a, 2011b; Madhusudhan 2012; Öberg et al. 2011; Mousis et al. 2012; Kopparapu et al. 2012; Moses et al. 2013) if the planet accreted carbon-rich solids inward of the water-ice line in the disk; if the planet accreted CO-rich, H₂O-poor gas from a region between the H₂O and CO ice lines, from a region in the innermost disk, or from within a heterogeneous disk or water-poor feeding zones; or if the disk were carbon-rich in the first place. However, graphite condensation followed by gravitational settling of the condensates (i.e., precipitation or "rain out") is expected to deplete the gaseous carbon and maintain the gas-phase C/O ratio near 1 above any graphite clouds in solar-metallicity exoplanet atmospheres (Lodders & Fegley 1997; Seager et al. 2005; Moses et al. 2013); therefore, we do not consider C/O ratios >1 in our current models. If the planet were so hot that graphite never becomes stable, if the graphite clouds were located within the photosphere, or if the graphite clouds were kinetically inhibited from forming even under conditions of thermodynamic-equilibrium stability, then photospheric gas-phase C/O ratios >1 would be possible (see Section 5 for more discussion).

In Figure 6, we show how the equilibrium abundances of several atmospheric constituents vary as a function of temperature and C/O ratio for an assumed 100 mbar pressure and a moderately high assumed C/H metallicity of 300 times solar (i.e., where all elements other than H, He, Ne, and O are assumed to be 300 times the protosolar abundances of Lodders et al. 2009). Recall here that we define the oxygen abundance through the C/O ratio, such that the very low C/O ratios at the left edge of these plots also correspond to large overall atmospheric metallicities, in addition to O being a dominant element. For example, our nominal solar C/H ratio is 2.78 × 10⁻⁴ and O/H ratio is 6.06 × 10⁻⁴. For a uniform enrichment of all elements of 300 times solar, the corresponding C/H and O/H ratios would be 8.33 × 10⁻² and 1.82 × 10⁻¹, respectively. For a C/H ratio of 300 times solar but oxygen defined via the C/O ratio from O/H = (C/H)/(C/O), a C/O ratio of 0.08 would correspond to an O/H ratio of 1.04, which is ∼1700 times the solar O/H ratio. Under such very low C/O ratio conditions (i.e., C/O ≲ 0.08) in this otherwise 300 times solar case, hydrogen is no longer the dominant element, and the H₂ abundance drops precipitously, as seen in Figure 6. Molecular oxygen then becomes the dominant gas at C/O ≲ 0.04 in this scenario—if such conditions have any relevance to real atmospheres—but H₂O quickly takes over with increasing C/O ratio to dominate at 0.05 ≲ C/O ≲ 0.19 (or at even greater C/O ratios for lower temperatures), whereas H₂ dominates at moderate-to-high C/O ratios (C/O ≳ 0.2). Carbon monoxide becomes an important constituent at moderate-to-high temperatures and moderate-to-high C/O ratios, CH₄ and NH₃ are important at low temperatures for all but the lowest C/O ratios, and CO₂ is important under most conditions except for the lowest temperatures considered at high C/O ratios. HCN increases in significance at high temperatures and high C/O ratios. Molecular nitrogen (not shown in Figure 6) is relatively unaffected by the atmospheric C/O ratio. Condensed graphite is stable above ∼750 K for the higher C/O ratios.

Zoom In Zoom Out Reset image size

Note the major and sudden shift from a more oxidized to a more reduced atmosphere at C/O ≈ 0.7 at high temperatures in Figure 6 for this 300 times solar-metallicity model. Graphite becomes stable at C/O ≳ 0.7 as this shift occurs, methane and HCN become more prevalent, and water and CO₂ become less prevalent. This transition is equivalent to the one that occurs at C/O ratios near 1.0 for a solar-composition, solar-metallicity gas, as discussed extensively by Lodders & Fegley (1997), Madhusudhan (2012), and Moses et al. (2013). For any given pressure in the planet's photospheric region, this key transition occurs at lower C/O ratios for higher metallicities. Thus the transition to the Madhusudhan (2012) "carbon-rich" C2 class would occur at lower C/O ratios for higher-metallicity planets, and metallicity can be an additional complication for any composition-based classification scheme.

The phase boundaries and stability regimes for the dominant carbon-bearing gases under the conditions plotted in Figures 5 and 6 are further illustrated in Figure 7. The conditions under which graphite is stable are also shown (see Section 5 for additional discussion). Note that carbon dioxide dominates at very high metallicities and/or very low C/O ratios. At more moderate metallicities and C/O ratios, CH₄ dominates at low temperatures and CO at high temperatures.

Zoom In Zoom Out Reset image size

The chemical-equilibrium results are also sensitive to pressure, as is shown in Lodders & Fegley (2002) and Visscher et al. (2006), for example. However, given that transport-induced quenching will affect the predicted abundances as a function of pressure in real atmospheres, we have simply chosen a single representative photospheric pressure for the above figures. The quench level may occur at higher pressures (lower altitudes) in hot-Neptune atmospheres, thereby affecting the predicted compositions, and the infrared photosphere of high-metallicity exoplanets may reside at lower pressures than our nominal choice of 100 mbar. In fact, we emphasize that these simple equilibrium models are largely phenomenological and are not designed to represent all the complex processes that have gone into shaping the atmospheric composition of actual exoplanet atmospheres. These equilibrium models do serve a useful purpose, though, in illustrating the possible diversity of hot Neptunes and in highlighting specific trends such as the increasing dominance of CO₂ at very high metallicities, the importance of H₂O at moderate-to-high metallicities for a variety of other conditions, and the change in the relative importance of methane and CO with increasing temperature. These general trends can be useful for considerations of the likely bulk atmospheric properties of specific hot Neptunes such as GJ 436b, based on the compositional clues provided by transit and eclipse data.

The secondary-eclipse data for GJ 436b (Stevenson et al. 2010; see also Beaulieu et al. 2011; Knutson et al. 2011; Stevenson et al. 2012) provide important constraints for compositional models. Madhusudhan & Seager (2011) use a Markov Chain Monte Carlo (MCMC) method to help identify the range of parameter space allowed for GJ 436b from comparisons of synthetic emission spectra with the Spitzer photometric data. Consistent with the earlier conclusions of Stevenson et al. (2010), Madhusudhan & Seager (2011) find that plausible GJ 436b models require a low methane abundance (e.g., mole fractions of 10⁻⁶–10⁻⁷) and a large a CO abundance (mole fraction ⩾10⁻³), along with a H₂O mole fraction ⩽10⁻⁴ and a CO₂ mole fraction in the range ∼10⁻⁶–10⁻⁴ in order to adequately reproduce the Stevenson et al. (2010) eclipse data. As is noted by Stevenson et al. (2010) and Madhusudhan & Seager (2011), such compositions are inconsistent with low-metallicity equilibrium models. Line et al. (2011) and Venot et al. (2013) further demonstrate that disequilibrium chemical processes like photochemistry and transport-induced quenching do not help resolve this problem.

The relatively sparse spectral coverage, low signal-to-noise ratio, and systematic uncertainties of the Spitzer secondary-eclipse data from GJ 436b and other exoplanets make retrieving atmospheric information difficult, as many solutions are statistically valid (Madhusudhan & Seager 2009, 2011; Lee et al. 2012; Line et al. 2012, 2013; Benneke & Seager 2012; Barstow et al. 2013). As a check on the derived best-fit abundances and thermal profile for GJ 436b, we apply the CHIMERA retrieval code of Line et al. (2013) to the Stevenson et al. (2010) secondary-eclipse data. Briefly, CHIMERA employs a suite of Bayesian retrieval algorithms—optimal estimation, bootstrap Monte Carlo, and differential-evolution MCMC—to determine the allowed range of temperatures and gas-phase mixing ratios for GJ 436b. Here we discuss the results of the differential-evolution MCMC approach. With this technique, the posterior probability distribution for each of the parameters that controls the temperature structure (based on the analytic parameterizations of Guillot 2010) and assumed constant-with-altitude mole fractions for H₂O, CH₄, CO, and CO₂ is characterized using a genetic algorithm that generates ∼10⁵ models (see Line et al. 2013, for further details). The statistics from the retrieval are shown in Figures 8–10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, Figure 8 shows both the retrieved thermal structure and the corresponding spectra for the ensemble of models for GJ 436b generated from the differential-evolution MCMC approach. Rather than plotting the many thousands of spectra and thermal profiles generated from this retrieval approach, we instead plot the median of the spectra (or temperatures in the plot insert) in black, along with the 1σ and 2σ spread in the spectra (or temperatures) in dark gray and light gray, respectively. In essence, these spreads reflect the fact that if we were to draw a random set of parameters from the posterior probability distributions, there would be a 95% chance that the flux at any one wavelength corresponding to the Spitzer bandpasses would fall within the 2σ spread, and so on. Note that although the synthetic spectra qualitatively reproduce the 3.6 to 4.5 μm flux ratio from the Spitzer photometric data, the non-detection of the eclipse in the 4.5 μm bandpass is particularly difficult to reproduce (see also Stevenson et al. 2010; Madhusudhan & Seager 2011).

Next, Figure 9 shows histograms of the H₂O, CH₄, CO, and CO₂ mole fractions for GJ 436b from the marginalized posterior probability distribution, as derived from the differential-evolution MCMC approach. These distributions show that H₂O and CO are relatively well constrained from the GJ 436b secondary-eclipse spectra, with the most probable H₂O mole fraction residing within the range of a few times 10⁻⁷ to a few times 10⁻⁴, and the most probable CO mole fraction restricted to ≳10⁻⁴, although the solutions for both H₂O and CO contain an extended, highly unconstrained tail at lower mixing ratios. The methane distribution also shows that the CH₄ mole fraction is restricted to values less than 10⁻⁶. These results are consistent with the analysis of Madhusudhan & Seager (2011). The posterior CO₂ distribution has an interesting double-peaked structure that makes the precise value for the CO₂ mole fraction less well constrained, but viable solutions are found for mole fractions between 10⁻⁸ and 10⁻². This result for CO₂ is also consistent with the χ²/N_obs ⩽ 2 results of Madhusudhan & Seager (2011), where N_obs is the number of available photometric data points. Our results confirm the picture of a GJ 436b atmosphere that possesses a large CO abundance, a very low methane abundance, and a moderately low water abundance.

Finally, Figure 10 shows the correlations in the retrieved abundances amongst the different gases. The CH₄ abundance remains low, regardless of the abundance of the other species, and thus the CH₄ mole fraction appears uncorrelated with either CO, CO₂, or H₂O. An apparent correlation between CO and CO₂ results from the very low flux (non-detection) in the 4.5 μm Spitzer channel where both CO and CO₂ absorb: if the CO mole fraction is large enough, the CO₂ mole fraction is relatively unconstrained and can be low, but if the CO mole fraction is small, the CO₂ mole fraction must be large in order to explain the very low 4.5 μm flux. Note also that solutions with low CO and high CO₂ do not fit the data as well as models with relatively high CO, perhaps because the high-CO₂ solutions end up with too much absorption in the 16 μm Spitzer channel. Figure 10 also illustrates a correlation between CO₂ and H₂O. If water is fairly abundant, the CO₂ abundance is not well constrained, but if the water abundance is low, the CO₂ abundance is more tightly constrained to fall within a mole-fraction range of ∼10⁻⁸–10⁻⁶. At the highest water abundances, there is also a positive correlation between H₂O and CO₂ in that very high water abundances require very high CO₂ abundances. This may reflect a strong correlation between H₂O (as the gas that has the largest overall influence on the spectrum) and temperatures, in that a very large water abundance can only be accommodated through high temperatures, at which point more CO₂ is needed to keep the low flux at 4.5 μm. Given that more solutions are found for low-to-moderate water abundances, the low-abundance peak for CO₂ in Figure 9 dominates over the secondary high-CO₂-abundance peak.

These new retrievals do not help resolve the apparent puzzles with respect to atmospheric composition on GJ 436b. Recalling Figures 5 and 6, a very low methane abundance can occur in equilibrium at very high temperatures at low metallicities (i.e., T ≳ 1200 K in the photosphere, hotter than is likely for GJ 436b) or at more moderate temperatures when the C/O ratio is low or when the atmospheric metallicity is very high. However, under those conditions, the H₂O abundance is likely to be significantly larger than is indicated by all the retrievals (above, and Madhusudhan & Seager 2011). Similarly, the preference for high CO mole fractions in the retrievals suggests a high metallicity for GJ 436b, which again implies unacceptably large water (and potentially CO₂) abundances. In fact, looking at Figures 5 and 6, we need to emphasize that there are no equilibrium chemistry solutions—or even, as we will show below, no disequilibrium chemistry solutions—that fall within the favored temperature-abundance ranges indicated by the retrievals. This result illustrates one potential drawback of retrievals from sparse, noisy, systematics-prone data: the solutions may reproduce the data well but not make any physical sense. Going to higher metallicities does help in general in that the photosphere is both hotter and more likely to reside in the CO-dominated regime (see Figure 1), but a high water abundance is a necessary consequence of high metallicity, unless the metallicity is so large (≳10,000 times solar) that the atmosphere contains very little H. That solution, which corresponds to Z_env = 0.994, (i.e., 0.6% H/He by mass), is unfortunately ruled out by GJ 436b's mass–radius combination, according to most interior models (cf. Adams et al. 2008; Baraffe et al. 2008; Figueira et al. 2009; Rogers & Seager 2010a; Nettelmann et al. 2010; Miller & Fortney 2011).

Figure 11 illustrates the possible C/O ratio versus metallicity phase space that can accommodate a low methane abundance on GJ 436b for the assumption of a moderately hot, high-altitude photosphere (e.g., assuming 900 K at 10 mbar), which may be relevant for a high-metallicity GJ 436b atmosphere. The shaded regions in this plot are consistent with the retrieved CH₄ mole fractions of ≲ 1 ppm, which can occur at low atmospheric C/O ratios and/or high metallicities. Again, there are no equilibrium solutions within the shaded regions of Figure 11 that are entirely consistent with the retrievals, with the main problem being too much H₂O and too much CO₂. However, the atmosphere is not likely to be in equilibrium (Line et al. 2010; Moses 2013), and we use these constraints as a guide for our disequilibrium kinetics and transport models, in an attempt to find forward models with more physically meaningful atmospheric properties that can provide a reasonable fit to the GJ 436b Spitzer secondary-eclipse data.

Zoom In Zoom Out Reset image size

Photochemistry and transport-induced quenching can affect the predicted photospheric abundances on extrasolar giant planets, resulting in mixing ratios that are often orders of magnitude different from chemical-equilibrium expectations (Liang et al. 2003, 2004; Zahnle et al. 2009; Line et al. 2010, 2011; Moses et al. 2011, 2013; Miller-Ricci Kempton et al. 2012; Kopparapu et al. 2012; Hu et al. 2012; Venot et al. 2012, 2013). Most disequilibrium models to date have focused either on giant planets with a near-solar-like complement of elements or on hydrogen-poor terrestrial exoplanets or super-Earths. We now test various scenarios with moderate and high metallicities (but still containing a non-negligible hydrogen mass fraction) for GJ 436b using disequilibrium kinetics/transport models to see how photochemistry and transport-induced quenching can alter the predicted transit and eclipse spectra. We start with the first-principles-based temperature profiles shown in Figure 1 and adopt a solar C/O ratio and a constant eddy diffusion coefficient of K_zz = 10⁹ cm²  s⁻¹ for these initial models. Figure 12 shows how the model results for CH₄, H₂O, CO, and CO₂ compare with the mixing ratios retrieved from the Spitzer secondary-eclipse data. This figure demonstrates that even when photochemistry and transport-induced quenching are considered, none of these models can reproduce the abundances derived from the retrievals.

Zoom In Zoom Out Reset image size

The mixing-ratio profiles for other potentially interesting species in these disequilibrium models are shown in Figure 13. The species profiles for the 1 times solar model are qualitatively similar to those of Line et al. (2011): behind H₂ and He, the next most abundant gases are H₂O, CH₄, and NH₃. Water survives at a near-equilibrium abundance until very high altitudes (∼1 μbar), at which point photolysis and other destruction mechanisms start to irreversibly convert the H₂O to CO and other oxygen-bearing species. Although photolysis and other chemical mechanisms operate to destroy H₂O at lower altitudes, water is efficiently recycled in the background H₂ atmosphere, and H₂O remains the dominant infrared opacity source in these 1 times solar models. Methane is even less stable than water at high altitudes, due primarily to the large H abundance released by H₂O photolysis and subsequent catalytic destruction of H₂ (e.g., Liang et al. 2003; Moses et al. 2011). The carbon that was in high-altitude CH₄ gets photochemically converted to CO, C₂H₂, HCN, and atomic carbon, primarily (see Moses et al. 2011 for details). Ammonia is photolyzed by longer-wavelength UV radiation that can penetrate a bit deeper in the atmosphere, and the nitrogen liberated by the NH₃ photodestruction largely ends up as HCN and atomic N. There are some quantitative differences between our 1 times solar model and that of Line et al. (2011) due to different adopted reaction rate coefficients, K_zz assumptions, stellar ultraviolet flux, and atmospheric temperatures, but the results from both Line et al. (2011) and our own models indicate that CH₄ survives photochemical destruction throughout the bulk of the ∼0.0001–1 bar photosphere on the 1 times solar GJ 436b. As a result, the 1 times solar model has more methane than is indicated by the retrieval analyses of the Spitzer secondary-eclipse data (see Figure 12, Section 3.2, and Madhusudhan & Seager 2011).

Zoom In Zoom Out Reset image size

The resulting spectrum for the 1 times solar model does not compare well with the Spitzer secondary-eclipse data (Figure 14). One obvious problem is a reversed 3.6 to 4.5 μm flux ratio in comparison with observations. The CO–CH₄ quench point in the 1 times solar model occurs within the CH₄-dominated regime, so that the quenched CO abundance is relatively small, and the methane abundance is large. Since the equilibrium abundance of CO is everywhere lower than that of CH₄ with the 1× solar thermal profile (see Figure 1), there is no eddy K_zz value we could adopt that would change this conclusion. Some CO is produced at high altitudes from CH₄ and H₂O photochemistry, but the CO column abundance is never large enough to influence the relative 3.6 to 4.5 μm flux ratio. The large methane abundance therefore results in a lower flux (more absorption) at 3.6 μm in comparison with 4.5 μm. Another obvious problem is the overall lower brightness temperature at most wavelengths in comparison with the data. The thermal profile from Lewis et al. (2010) adopted in the 1 times solar model is not significantly colder than the profiles derived from the thermal retrievals shown in Madhusudhan & Seager (2011) or in Section 3.2, but the 1 times solar model does have a water abundance at the upper end (or greater than) the retrievals indicate (see Figure 12). Water accounts for most of the absorption in the 1 times solar model, and the relatively large water abundance leads to more absorption at most wavelengths than can be accommodated by the observations—hence the requirement of a low water abundance from the retrievals in the first place. Since photochemistry does not permanently destroy the photospheric water and methane in our 1 times solar models, such models cannot explain the Spitzer secondary-eclipse data.

Zoom In Zoom Out Reset image size

The quench point for the 50 times solar model is closer to the equilibrium CO–CH₄ equal-abundance boundary due to an overall hotter atmosphere (see Figure 1), and both CO and CH₄ end up being major atmospheric constituents (see Figure 13). Carbon dioxide and molecular nitrogen also become more important constituents at the higher 50 times solar metallicity, while ammonia becomes relatively less important due to the N₂–NH₃ quench point falling within the N₂-dominated regime. Because NH₃ is more photochemically active than N₂, the photochemical production of HCN and complex nitriles does not increase significantly at high altitudes in the higher-metallicity model because the overall NH₃ mixing ratio has not changed much between 1 times and 50 times solar metallicities. However, the column abundance of HCN actually increases with the increased metallicity here because of thermochemical kinetics and the quenching of CH₄ and NH₃ at higher-than-equilibrium abundances (i.e., HCN in the photosphere maintains a pseudo-equilibrium with the quenched NH₃ and CH₄; see Moses et al. 2011, 2013). Carbon dioxide is produced effectively at high altitudes from the photochemistry of CO and H₂O, but equilibrium through the net reaction CO + H₂O $\leftrightarrows$ CO₂ + H₂ is maintained kinetically through much of the photosphere. Although the CO mole fraction in the 50× solar model now falls within the range required by the retrievals discussion in Section 3.2, both H₂O and CH₄ are much more abundant than is allowed by the retrievals, and the 50× solar model does not fit the Spitzer secondary-eclipse data well.

At metallicities of 1000 times solar, hydrogen and helium now make up only ∼6% of the atmosphere by mass. The methane abundance has dropped significantly in the 1000 times solar-metallicity model because the CH₄–CO quench point lies within the CO-dominated regime. The resulting CO/CH₄ ratio is finally in the right direction to explain the observed relative 3.6 to 4.5 μm flux ratio from the Spitzer eclipse data. Carbon monoxide begins to rival H₂ as the dominant constituent, and H₂O, CO₂, and N₂ are all very abundant. The atmosphere as a whole is more oxidized, and species like O₂ and NO that are produced from photochemistry at high altitudes are able to survive more readily than they do in a more reducing atmosphere. The fact that elements other than O, C, N, and H have been neglected from these disequilibrium models (particularly Cl and S) will have an impact on the resulting abundances in these higher-metallicity models due to omitted catalytic cycles (e.g., Yung & DeMore 1999), and atmospheric escape and other evolutionary processes will likely alter this simple picture of steadily increasing metallicity, but the general supplantation of H₂ and the increasing dominance of CO₂ with increasing metallicity is a robust conclusion. At atmospheric metallicities of 10,000× solar, for instance, hydrogen and helium make up only ∼0.6% of the atmosphere by mass, and Figure 13 shows that CO₂ has solidly replaced H₂ as the dominant constituent, and even CO, N₂, and H₂O are more abundant than H₂.

Figure 14 shows that the 1000 times solar-metallicity model fares better in reproducing the Spitzer secondary-eclipse data than the 1 times solar-metallicity model did. The flux in the 3.6 μm channel is predicted to be greater than that in the 4.5 μm model now with the 1000 times solar model, due to the CH₄ abundance being much less than the CO abundance. The 1000 times solar model still provides too much absorption at 3.6 μm to be consistent with the eclipse depth as determined by Stevenson et al. (2010), indicating too much methane in the model; however, we should note that Beaulieu et al. (2011) derive both a lower 3.6 μm flux and a larger uncertainty in the eclipse depth at this wavelength (and at 4.5 μm), and the 1000× solar model more readily satisfies those constraints. More intriguing is the excellent fit in the 5.8 and 8.0 μm channels where water in this model is providing the dominant opacity source. This good fit with a water abundance that exceeds a 10% mixing ratio—in contrast to the much lower abundances derived from the retrievals—reflects the influence of the higher-temperature photosphere that results from the increased atmospheric opacity at high metallicities. The water abundance and photospheric temperatures are closely linked, and a hotter photosphere requires more water to remain consistent with the 5.8 and 8.0 μm eclipse depths, whereas, conversely, a colder photosphere requires less water. Still, this 1000 times solar model provides a relatively poor fit to the eclipse data of Stevenson et al. (2010, 2012), providing a χ²/N_obs ∼ 11.5, with the main problem being a 3.6 μm flux that is much lower than observed, but also a 4.5 μm flux that is in excess of the observational upper limit and a 16 μm flux that is much lower than observed (due to CO₂ absorption).

Although the 1000 times solar-metallicity model itself does not provide a good fit to the Stevenson et al. (2010, 2012) secondary-eclipse data, and none of the models in Figure 12 fall within the constraints imposed by the retrievals for all the species, the model-data comparisons do suggest further avenues to explore with higher-metallicity models. The larger observed brightness temperature at 3.6 μm in relation to that at 5.8 and 8 μm in the Spitzer data suggests that the atmosphere has a sharper thermal gradient than the one derived from the 1D PHOENIX model, given that the contribution functions at the three wavelengths are not hugely separated in terms of the pressures at which they peak in the high-metallicity models. We note that both the dayside-average profiles from the 3D GCMs of Lewis et al. (2010) and the thermal profiles derived from the retrievals exhibit a larger thermal gradient in the relevant photospheric region than the 1000 times solar 1D PHOENIX models do. A larger thermal gradient could also result in a hotter atmosphere at the CH₄–CO quench point, if the quench point occurs at a high-enough altitude, pushing the atmosphere toward greater CO dominance at the expense of CH₄. That would improve the relative 3.6 to 4.5 μm flux ratio in comparison with observations. To investigate such a scenario, we simply scale one of the thermal profiles derived from the Madhusudhan & Seager (2011) MCMC analysis upward uniformly in altitude to lower pressures and run different disequilibrium chemistry models, assuming different bulk atmospheric properties. The results from one such model are shown in Figure 15.

Zoom In Zoom Out Reset image size

For the model presented in this figure, the atmospheric C/H and N/H metallicity are assumed to be 300 times solar, the C/O ratio above the silicate clouds is assumed to be 0.6 (which actually requires a slightly sub-solar initial bulk C/O ratio if ∼21% of the oxygen is sequestered in silicates and other condensates at deeper atmospheric levels, as expected by Lodders 2010 and Visscher et al. 2010a), and the eddy diffusion coefficient follows the relation K_zz = 10⁷ cm²  s⁻¹ at pressures P ⩾ 10⁻⁴ bar, with K_zz increasing with the inverse square root of atmospheric pressure at P < 10⁻⁴ bar. The thermal profile as shown in the insert of Figure 4 of Madhusudhan & Seager (2011) has been shifted upward uniformly by −1.2 in log₁₀(P) (i.e., pressures have been multiplied by 10^−1.2) to account for the higher-altitude photosphere that will result from the higher metallicity. The 300 times solar atmospheric metallicity in this model is consistent with the interior models for GJ 436b (Adams et al. 2008; Baraffe et al. 2008; Figueira et al. 2009; Rogers & Seager 2010a; Nettelmann et al. 2010; Miller & Fortney 2011), and the lower value of K_zz in the lower atmosphere allows CH₄ to remain in equilibrium until ∼0.1 bar, at which point the methane quenches at the desired <1 ppm abundance. Judging from the PHOENIX models, the thermal profile at photospheric pressures may be hotter in this model than is expected on average in the dayside hemisphere of GJ 436b at this metallicity, so the likelihood of this scenario will need to be investigated by further 3D circulation modeling, but given that the observed emission will derive preferentially from the hottest regions on the observed disk, the profile does not appear to be unreasonably hot.

Figure 15 shows that H₂ is still the dominant atmospheric constituent (by number) in the atmosphere in this 300 times solar-metallicity model, although both CO and H₂O at mole fractions of ∼10% each are encroaching on the H₂ dominance. Kinetic production of CO₂ from H₂O and CO occurs in the middle and upper photosphere in this model, allowing the CO₂ abundance to exceed that of water at high altitudes. Quenching of ammonia at non-negligible values occurs and is responsible for maintaining HCN at greater than ppm levels. The CH₄ abundance remains low, and the emission spectrum (Figure 16) is dominated by opacity from H₂O, CO₂, and CO. This model finally has a low-enough methane abundance that the 3.6 μm flux from the Spitzer secondary-eclipse analysis of Stevenson et al. (2010) is better reproduced. The resulting 3.6 to 4.5 μm flux ratio is also more in line with observations, although this model, like all others including the retrievals, still produces too much flux in the 4.5 μm band. This result is more a function of the atmosphere being too hot where the contribution function peaks at 4.5 μm, rather than being due to insufficient CO and CO₂, as can be readily seen in Figures 12 and 15. In fact, the very large CO₂ abundance in this 300 times solar-metallicity model produces too much absorption in the 16 μm Spitzer channel, where the contribution function peaks at cooler, higher altitudes due to the large CO₂ column abundance. The non-detection of the GJ 436b secondary eclipse at 4.5 μm is difficult to explain by any known physically reasonable scenario without adversely affecting the fit at other wavelengths. In terms of the statistical fit to all the Stevenson et al. (2010, 2012) secondary-eclipse depths, this model fares better, at χ²/N_obs = 2.8. The primary problems are insufficient flux at 3.6, 16, and 24 μm, and too much flux at 4.5 μm.

Zoom In Zoom Out Reset image size

This particular forward model is not unique in terms of fitting the data, and there are numerous models with atmospheric metallicities in the couple-hundred to couple-thousand times solar range that fit the Spitzer data to within χ²/N_obs ⩽ 3 (see, for example the 2000 times solar-metallicity model we presented in Richardson et al. 2013), but we have not found any disequilibrium chemistry models that fit to within χ²/N_obs ⩽ 2. All of these high-metallicity models produce too much CO₂ to be consistent with the observed flux in the 16-μm channel. Models with high metallicities in the 100–10,000 times solar range also end up with a large water abundance, which then requires a hot photosphere to explain the 5.8, 8.0, and 24 μm fluxes, at which point there is too much emission at 4.5 μm. Slightly cooler models than the one shown in Figures 15 and 16 tend to fit better at 4.5 μm, but then the fit tends to be worse at 3.6 and 24 μm. Hotter models fit the 24 μm flux better, but then the fit at 4.5, 5.8, and 8.0 μm is worse. Going to a higher C/O ratio will help remove excess CO₂ to improve the fit at 16 μm, but then CH₄ becomes too abundant to explain the 3.6 μm flux. Lower C/O ratios help keep the methane abundance low, but then the problems with excess H₂O and CO₂ are exacerbated. Therefore, while these high-metallicity models are promising in their ability to qualitatively reproduce several aspects of the observed secondary-eclipse observations, several quantitative problems remain. The solutions derived from the retrieval methods (Section 3.2 and Madhusudhan & Seager 2011) provide a better overall fit to the data, but problems with the fit still remain at 4.5 and/or 16 μm; more importantly, these solutions have problems with plausibility in that it is difficult to theoretically explain how such abundances (e.g., the inferred CO/H₂O ratio) could be obtained chemically or physically in an atmosphere with any reasonable bulk properties. Therefore, despite the less-than-perfect fit to the existing GJ 436b secondary-eclipse data, such high-metallicity models should not be dismissed out of hand.

The higher the atmospheric metallicity, the higher the mean molecular weight, and the smaller the atmospheric scale height for any given temperature and gravity profile. High-metallicity atmospheres will then have flatter transmission spectra (e.g., Miller-Ricci et al. 2009; Fortney et al. 2013). Examples of transmission spectra from some of our GJ 436b models are shown in Figure 17, in comparison with observations. The relatively flat spectrum observed in HST/NICMOS data (Pont et al. 2009; Gibson et al. 2011) seems more consistent with the high-metallicity models (see also H. A. Knutson et al. 2013, in preparation, for comparisons with recent HST/WFC3 data, which are also quite flat in the 1.1–1.6 μm region), and such models also seem to compare well to the results from the Knutson et al. (2011) Spitzer transit analysis (see the open squares at 3.6, 4.5, and 8 μm in Figure 17). However, the high-metallicity models do not compare well with the Beaulieu et al. (2011) Spitzer transit analysis (see the open diamonds in Figure 17) or the K_s-band transit data of Cáceres et al. (2009), and in fact none of the models provides a reasonable fit to all the available transit data for GJ 436b. Given the observed variability in the apparent transit depths with time from different observations (e.g., Knutson et al. 2011) and the lack of agreement between transit results analyzed with different procedures (e.g., Pont et al. 2009; Gibson et al. 2011; Knutson et al. 2011; Beaulieu et al. 2011), the poor fits of these models to the transit data are not particularly surprising. Obtaining simultaneous spectral observations at multiple wavelengths should help resolve some of these issues and help distinguish between competing models (see also Benneke & Seager 2013).

Zoom In Zoom Out Reset image size

Note that we have not adjusted for potential differences in temperatures or chemical abundances between the terminator and dayside atmospheres when performing these calculations. The dayside atmosphere observed at secondary eclipse is likely hotter than the limb atmosphere at the terminators probed in the planetary transit (see Lewis et al. 2010). The non-negligible eccentricity of GJ 436b will also lead to variations in temperature over the orbit, which is captured in the GCMs of Lewis et al. (2010). It is likely that horizontal quenching will help homogenize the atmospheric abundances with respect to longitude due to strong zonal winds or rapid thermal changes over the orbit (Cooper & Showman 2006; Iro & Deming 2010; Visscher 2012; Agúndez et al. 2012), but without more detailed multi-dimensional and temporal modeling, it is difficult to predict the exact terminator composition. Regardless of the exact molecular composition, however, the high-metallicity models will have a higher atmospheric mean molecular mass, and a flatter overall transit spectrum.