Atmospheric Circulation of GJ436b

Atmospheric Circulation of Eccentric Hot Neptune GJ436b

Nikole K. Lewis and Adam P. Showman Department of Planetary Sciences and Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721 Jonathan J. Fortney Department of Astronomy & Astrophysics, University of California, Santa Cruz, CA 95064 Mark S. Marley and Richard S. Freedman11affiliation: SETI Institute, 515 N. Wishman Road, Mountain View, CA 94043 NASA Ames Research Center 245-3, Moffett Field, CA 94035 Katharina Lodders Washington University, St. Louis, MO 63130

GJ436b is a unique member of the transiting extrasolar planet population being one of the smallest and least irradiated and possessing an eccentric orbit. Because of its size, mass and density, GJ436b could plausibly have an atmospheric metallicity similar to Neptune (20-60 times solar abundances), which makes it an ideal target to study the effects of atmospheric metallicity on dynamics and radiative transfer in an extrasolar planetary atmosphere. We present three-dimensional atmospheric circulation models that include realistic non-gray radiative transfer for 1, 3, 10, 30, and 50 times solar atmospheric metallicity cases of GJ436b. Low metallicity models (1 and 3 times solar) show little day/night temperature variation and strong high-latitude jets. In contrast, higher metallicity models (30 and 50 times solar) exhibit day/night temperature variations and a strong equatorial jet. Spectra and light curves produced from these simulations show strong orbital phase dependencies in the 50 times solar case and negligible variations with orbital phase in the 1 times solar case. Comparisons between the predicted planet/star flux ratio from these models and current secondary eclipse measurements support a high metallicity atmosphere (30-50 times solar abundances) with disequilibrium carbon chemistry at play for GJ436b. Regardless of the actual atmospheric composition of GJ436b, our models serve to illuminate how metallicity influences the atmospheric circulation for a broad range of warm extrasolar planets.

planets and satellites: general, planets and satellites: individual: GJ436b, methods: numerical, atmospheric effects
slugcomment: Accepted to the Astrophysical Journal

1 Introduction

The “hot Neptune” GJ436b was first discovered by Butler et al. (2004) and later determined to transit its host star as seen from earth by Gillon et al. (2007b). Since the discovery of its transit and subsequent secondary eclipse, GJ436b has become a popular target for Hubble (Bean et al., 2008; Pont et al., 2009) and Spitzer (Gillon et al., 2007a; Deming et al., 2007; Demory et al., 2007; Stevenson et al., 2010) observations as well as modeling efforts (Spiegel et al., 2010; Madhusudhan & Seager, 2010). Given GJ436b’s mass (=0.0729M) and radius (=0.3767R), its interior must contain significant quantities of heavy elements in addition to hydrogen and helium (Adams et al., 2008; Figueira et al., 2009; Rogers & Seager, 2010; Nettelmann et al., 2010). This raises the possibility that, like Uranus and Neptune, whose atmospheric C/H ratios lie between 20 and 60 times solar (Gautier et al., 1995), the atmosphere of GJ436b is highly enriched in heavy elements. This makes GJ436b an excellent case study for atmospheric chemistry, radiative transfer, and global circulation that should differ significantly from the well studied “hot Jupiters” HD209458b and HD189733b.

Observations of HD189733b using the Spitzer Space Telescope provided the first clear evidence for atmospheric circulation on an extrasolar planet (Knutson et al., 2007, 2009). Most efforts to model atmospheric circulation for extrasolar planets have focused on hot Jupiters, specifically HD189733b and HD209458b (Showman & Guillot, 2002; Showman et al., 2008, 2009; Cho et al., 2003, 2008; Cooper & Showman, 2005, 2006; Dobbs-Dixon & Lin, 2008; Menou & Rauscher, 2009; Rauscher & Menou, 2010). Only a handful of studies have specifically investigated the effects of non-synchronous rotation (Cho et al., 2008; Showman et al., 2009), non-zero obliquity (Langton & Laughlin, 2007), and non-zero eccentricity (Langton & Laughlin, 2008). The possible effect of atmospheric composition, and hence opacity, on circulation patterns that may develop on extrasolar planets has been investigated to some extent by Dobbs-Dixon & Lin (2008) and Showman et al. (2009), but largely ignored in most of the current two-dimensional and three-dimensional atmospheric models. Atmospheric composition is key in determining opacity and radiative timescales that play a crucial role in the development of circulation on these planets.

Here we present three-dimensional atmospheric models for GJ436b that incorporate both equilibrium chemistry and realistic non-gray radiative transfer. Although the actual composition of GJ436b’s atmosphere is likely to deviate from an equilibrium chemistry solution as shown from the secondary eclipse observations of Stevenson et al. (2010), our investigation still serves to explore the effect metallicity can play in controlling the atmospheric circulation not only on GJ436b but for a broad range of gaseous extrasolar planets in a similar temperature range. Our models are not constrained to match measured chemical abundances and temperatures, but instead provide a systematic look at how changes in atmospheric metallicity over the range from 1 times to 50 times solar values affects the basic thermal and dynamical structure of the planet’s atmosphere. We do not expect that spectra and light curves from our model will provide a match to observational data, but instead illuminate some of the underlying atmospheric physics responsible for current observations and suggest areas of focus for future observations. Section 2 gives an overview of the three-dimensional coupled radiative transfer and atmospheric dynamics model used in this study. Section 3 presents the global thermal structures and winds that develop in each of our models along with predicted light curves and emission spectra. Sections 4 and 5 provide a brief discussion of the results and final conclusions.

2 Model

The atmospheric model used in this study is a three-dimensional (3D) coupled radiative transfer and dynamics model that was specifically developed with the study of extrasolar planetary atmospheres in mind. The Substellar and Planetary Atmospheric Radiation and Circulation (SPARC) model is described in detail in Showman et al. (2009) as applied to HD189733b and HD209458b. A basic overview of the SPARC model along with the specific changes made to the model setup for GJ436b are presented here for completeness. The SPARC model employs the MITgcm (Adcroft et al., 2004) to treat the atmospheric dynamics using the primitive equations, which are valid in stably stratified atmospheres where the horizontal dimensions of the flow greatly exceed the vertical dimension. For GJ436b, the horizontal length scale of the flow is m while the vertical scale height of the atmosphere is km. The simulations presented here take advantage of the cubed-sphere grid (Adcroft et al., 2004) at a resolution of C32 (roughly in latitude and longitude) to solve the relevant dynamic and energy equations. The vertical dimension in these simulations spans the pressure () range from 200 bar to 20 bar with 47 vertical levels, evenly spaced in . The boundary conditions in our simulations are an impermeable surface at the bottom and a zero pressure surface at the top both of which are free slip in horizontal velocity.

We have coupled the MITgcm to the non-gray radiative transfer model of Marley & McKay (1999) to realistically determine the magnitude of heating/cooling at each grid point. The radiative transfer model, a two-stream version of the Marley & McKay (1999) plane-parallel code, assumes local thermodynamic equilibrium and includes intensities over the wavelength range from 0.26 to 300 m. The opacity at each pressure-temperature-wavelength grid point is tabulated using the correlated-k method (Goody et al., 1989). Our extensive opacity database is described in Freedman et al. (2008). The chemical mixing ratios, which are computed assuming thermochemical equilibrium, are calculated as in Lodders & Fegley (2002, 2006). Calculated opacities assume a gaseous composition without particulate matter and account for the possibility of chemical rainout. Because GJ436b plausibly has an atmospheric chemistry that is enhanced in heavy elements, we developed opacity tables for 3 times (3), 10 times (10), 30 times (30), and 50 times (50) solar metallicity in addition to the 1 times (1) solar metallicity opacity table. In the enhanced metallicity opacity tables, all elements other than hydrogen and helium are assumed to be enhanced by the same factor over current solar values. The opacity databases of Freedman et al. (2008) were updated to include the opacity effects of CO, which is an important carbon bearing species at higher metallicities. The full opacity tables are divided into 30 wavelength bins as outlined in Showman et al. (2009). This binning of opacities allows for greater computational efficiency while only introducing small () deviations from the net radiative flux calculated with higher resolution opacity tables.

For each model, the winds are assumed to initially be zero everywhere and each column of the grid is assigned the same pressure-temperature profile. This initial pressure-temperature profile is derived from one-dimensional radiative-equilibrium calculations performed using the radiative transfer code in the absence of dynamics. Figure 1 shows the pressure-temperature profiles derived for each metallicity case of GJ436b. These pressure-temperature profiles were derived using the methodology presented in Fortney et al. (2008, 2005). The physical properties assumed for GJ436b and its host star (GJ436A) are presented in Table 1. Using these planetary and stellar parameters, the effective temperature () of GJ436b is calculated to be 649 K, assuming planet-wide redistribution of the incoming stellar flux. This corresponds to a mean photospheric level111Defined in this context as the atmospheric pressure where the local temperature equals the effective temperature. of 1 to 100 mbar depending on the assumed metallicity of the atmosphere (Figure 1).

Because GJ436b is known to have an eccentric orbit, we incorporated the effects of non-synchronous rotation and time-varying distance from the host star into the SPARC model. The most probable rotation rate for GJ436b was determined using the following pseudo-synchronous rotation relationship presented in Hut (1981):


where is the planetary rotation rate, is the orbital period of the planet, and is the eccentricity of the planetary orbit. In all cases considered here the obliquity of the planet is assumed to be zero. The time-varying distance of the planet with respect to its host star, , is determined using Kepler’s equation (Murray & Dermott, 1999) and used to update the incident flux on the planet at each radiative timestep. A diagram of GJ43b’s orbit is presented in Figure 2. To test the impact of pseudo-synchronous rotation and time-varying stellar insolation, additional simulations for the 1 and 30 solar metallicity cases were performed assuming synchronous rotation and zero eccentricity.

In our models, for computational efficiency, the radiative timestep used to update the radiative fluxes is longer than the timestep used to update the dynamics. Generally, as we increased the metallicity of the atmosphere, progressively shorter radiative and dynamical timesteps were needed to maintain stability. For the 1 and 3 solar metallicity cases a dynamic timestep of 25 s and a radiative timestep of 200 s were used. The 10 and 30 solar metallicity cases required a dynamic timestep of 20 s and a radiative timestep of 100 s while the 50 solar case required a dynamic timestep of 15 s and a radiative timestep of 60 s. Timestepping in our simulations is accomplished through a third-order Adams-Bashforth scheme (Durran, 1991). We applied a fourth-order Shapiro filter in the horizontal direction to both velocity components and the potential temperature over a timescale equivalent to twice the dynamical timestep in order to reduce small scale grid noise while minimally affecting the physical structure of the wind and temperature fields at the large scale.

We integrated each of our models until the velocities reached a stable configuration. Figure 3 show the root mean square (RMS) velocity as a function of pressure and simulated time, calculated according to:


where the integral is a global (horizontal) integral over the globe, is the horizontal area of the globe, is the east-west wind speed, and is the north-south wind speed. The high-frequency variations in the RMS velocity seen in the upper levels of both the 1 and 50 solar cases are largely due to variation in the incident stellar flux associated with the eccentric orbit of GJ436b. Notice that, in the observable atmosphere (pressures less than 100 mbar), the orbit-averaged winds become essentially steady within 2500 Earth days for solar metallicity and 1000 Earth days for 50 solar metallicity. RMS wind speeds typically reach 1 km s at photosphere levels. Any further increases in wind speeds will be small and confined to pressure well below the mean photosphere so as not to affect any synthetic observations derived from our simulations. As outlined in Showman et al. (2009) the energy available for the production of winds is limited largely by the global available potential energy within the atmosphere and to some extent energy losses due to the Shapiro filter which acts as a hyperviscosity. A full discussion of the energetics of our simulated GJ436b-like atmosphere is left for a future paper.

3 Results

The following sections overview the key results from the study of GJ436b’s atmospheric circulation at 1, 3, 10, 30, and 50 solar metallicity. Both the thermal structure and winds in these simulations have a strong dependence on the assumed composition of the atmosphere for GJ436b. Additionally, theoretical light curves and spectra are produced from our 3D model atmospheres and compared with available data.

3.1 Thermal Structure and Winds: Dependence on Metallicity

Figure 4 presents snapshots of the temperature and wind fields at three pressure levels in the atmosphere for the 1 and 30 solar cases near secondary eclipse when the full day-side of the planet faces Earth (Figure 2). Overall, the 30 solar case is significantly ( K) warmer than the 1 solar case at each pressure. The increased atmospheric opacity that comes with metallicity enhancements leads to an upward shift in the pressure-temperature profiles. This effect is self-consistently generated in the three-dimensional model integrations but can also be seen in the one-dimensional radiative-equilibrium solutions shown in Figure 1. Overall, the day/night temperature contrast in the upper layers of the atmosphere (1 mbar) and the equator/pole temperature contrast deeper in the atmosphere (30 mbar) increase with atmospheric metallicity. However, because the pressure at a given optical depth is smaller at high metallicity than low metallicity, the regions that develop significant day/night temperature contrasts shifts to higher altitude as metallicity increases. At the 1 bar level, the equator/pole temperature contrast in the 30 solar case is smaller than that in the 1 solar case. This lack of a strong temperature contrast at 1 bar in the 30 solar case occurs because this pressure is at a greater optical depth in the 30 solar case than in the 1 solar case, and thus occurs below the levels with the strongest heating/cooling.

It is also informative to compare the flow patterns indicated by the arrows in Figure 4 between the 1 and 30 solar cases. The overriding feature in all simulations is the development of a prograde (eastward) flow at low pressure. The flow patterns in the 30 solar case exhibit clear wavelike structures outside of the equatorial region at the 1 bar, 30 mbar, and 1 mbar levels. At the 1 bar level, the flow is predominately westward in the 30 solar case and predominately eastward in the 1 solar case. In both the 1 and 30 solar metallicity cases the strength of the winds indicated by the length of the wind vectors in Figure 4 is a stronger function of latitude than longitude, except at the highest levels (1 mbar) in the 30 solar case, which shows a significant day/night temperature contrast similar to what was seen in the simulations of the tidally locked hot Jupiters HD189733b and HD209458b from Showman et al. (2009).

We find that the atmospheric metallicity plays a key role in determining the jet structure for a planet with temperatures similar to those expected on GJ436b. This is demonstrated clearly in Figure 5, which shows the zonal-mean zonal wind222That is, the longitudinally averaged east-west wind, where eastward is defined positive and westward negative. See Holton (2004). versus latitude and pressure for each of the five atmospheric metallicities for GJ436b considered in this study. In the 1 solar case, strong high-latitude jets develop in the atmosphere with a weaker equatorial jet. Increasing the metallicity of the planet to 3 solar causes a strengthening of the equatorial jet and a weakening of the high-latitude jets. Once the metallicity of the atmosphere is increased to 10 solar or more, the high-latitude jets disappear and the equatorial jet becomes dominant. Overall, the maximum zonal wind speed increases with metallicity from roughly 1300 m s in the 1 solar case to over 2000 m s in the 50 solar case. The flow is subsonic everywhere in our 1, 3, and 10 solar metallicity cases. In our 30 and 50 solar metallicity cases, the winds are subsonic throughout most of the domain, but they become marginally supersonic at the very top of the domain (at pressures 1 mbar) in the equatorial region. Hydraulic jumps similar to those seen in the HD189733b and HD209458b cases presented in Showman et al. (2009) are present in the supersonic regions of the atmosphere in the 30 and 50 solar metallicity cases (see 1 mbar level of the 30 solar case in Figure 4). It is important to note that the flow in all of the metallicity cases are predominately eastward at pressure less than 1 bar and predominately westward at pressures greater than 1 bar. Momentum conservation requires that the eastward momentum in the upper atmosphere of these simulations comes from the deeper atmospheric layer, which requires the development of the mean westward flow at depth. The detailed mechanisms responsible for this momentum and the jet pumping mechanisms themselves will be discussed in a future paper.

In rapidly rotating atmospheres, atmospheric temperature gradients are linked to winds by dynamical balances, so it is interesting to next examine the atmospheric temperature structure in our simulations. Figure 6 shows the zonal-mean atmospheric temperature versus latitude and pressure for the 1 and 30 solar cases. In both cases, the deepest isotherms (for temperatures exceeding 1200 K) are flat, but isotherms between 600 and 1100 K are bowed upward, indicating a warm equator and cool poles. The relationship between this structure and the winds can be understood with the thermal-wind equation, which relates the latitudinal temperature gradients to the zonal wind and its derivative with pressure:


where , , , , , and are the zonal wind speed, latitude, planetary radius, planetary rotation rate, pressure, specific gas constant, and temperature, respectively. The latitudinal temperature gradient on the right-hand side is evaluated at constant pressure. This relationship derives from taking a vertical (pressure) derivative of the meridional momentum equation for a flow where the predominant zonal-mean meridional momentum balance is between the Coriolis, pressure-gradient, and curvature terms (called “gradient-wind” balance; see Holton 2004, p. 65-68). From Equation (3), one expects that, away from the equator, regions exhibiting vertical shear of the zonal wind must also exhibit latitudinal gradients of temperature. Comparing Figures 5 and 6 confirms that this is indeed the case: in the mid-latitudes, the 1 solar case exhibits the greatest vertical shear of the zonal wind in the pressure range of 0.1 to 3 bars (Figure 5), and this is the same pressure range over which the greatest latitudinal temperature gradients occur (Figure 6). For the 30 solar case, in the mid-latitudes, the regions exhibiting significant wind shear are shifted upward, occurring from 0.01 to less than 1 bar, and likewise this is the pressure range where significant latitudinal temperature gradients exist. The upward shift in the temperature gradients (and hence winds) at greater metallicity is the direct result of enhanced atmospheric opacities, which lead to shallower atmospheric heating (Fortney et al., 2008; Dobbs-Dixon & Lin, 2008).

It is interesting to characterize the variations in wind speed that occur throughout the eccentric orbit due to the time-variable incident stellar flux. Figure 7 shows the RMS wind speed as a function of pressure and simulated time with outputs every two hours for a ten Earth day period after the simulations had reached equilibrium. Overall, wind speeds in these simulations of GJ436b vary with a frequency roughly equal to the orbital period at pressures above the mean photospheric level333In Figure 3, the high-frequency fluctuations appear to occur on periods of tens of Earth days, but this is an artifact that results from aliasing of the sampling frequency of s with the orbital period. (roughly 100 mbar in the 1 solar case and 10 mbar in the 50 solar case). Wind speeds are fairly constant at pressures below the mean photosphere for each metallicity case of GJ436b. The vertical lines in Figure 7 indicate the time of periapse passage, which are followed by a peak in the RMS wind speeds. The variation in the wind speeds between periapse and apoapse is several times larger in the 50 solar case compared with the 1 solar case. The higher atmospheric metallicity cases show greater variability in wind speeds as function of orbital phase, which could affect light curves especially if hotter or more eccentric systems are considered.

3.2 Effect of eccentricity and rotation rate

The models presented in Section 3.1 assume an eccentric orbit (hence time-variable stellar irradiation) and adopt the pseudo-synchronous rotation rate given in Table 1. Here, we explore the effect of eccentricity and rotation rate on the circulation by comparing our standard cases (Section 3.1) to cases with synchronous (rather than pseudo-synchronous) rotation rates and zero eccentricity. Figure 8 presents the thermal structure and winds at the 30 mbar level for the 1 and 30 solar cases assuming synchronous rotation (). The top panels of Figure 8 assume the nominal eccentric orbit of GJ436b while the bottom panels assume a circular orbit with the nominal semimajor axis of GJ436b (Table 1). The assumption of synchronous rotation and/or a circular orbit has little effect on the overall thermal structure and wind patterns that develop in these simulations. This is presumably because the eccentricity of GJ436b’s orbit is modest () and changes in the average stellar flux and pseudo-synchronous rotation rate from a circularized and tidally locked orbit are small.

Figure 9 presents the zonal-mean zonal winds for these same four cases (1 and 30 solar metallicity, with eccentricity of 0 or 0.15, all using the synchronous rotation period). Overall, the jet structures differ little from the nominal cases. It is interesting to note that the eastward equatorial jet in the 30 solar synchronous rotation cases has a maximum wind speed that is m s faster than what is seen in the non-synchronous case. Because GJ436b has a relatively small eccentricity orbit the effects of non-synchronous rotation and time-variable heating are only small perturbations on the synchronous rotation and circular orbit cases. Planets with higher eccentricities are likely to show a larger variation in the circulation patterns that develop compared with circularized and synchronous cases.

3.3 Light Curves and Spectra

The SPARC model is uniquely equipped to produce both theoretical light curves and spectra that account not only for radiative effects, but also dynamic movement in the atmosphere. Once each of the GJ436b models reached an equilibrium state, pressure and temperature profiles were recorded along each grid column at many points along the planet’s orbit (Figure 2). These pressure-temperature profiles were then used in high resolution spectral calculations to determine the emergent flux from each point on the planet and which portion of that emergent flux would be directed toward an earth observer including limb darkening/brightening effects. Spectra and light curve generation methods are fully described in Fortney et al. (2006). Figure 10 shows the theoretical light curves, expressed as the planet/star flux ratio, for the 1 and 50 solar cases as a function of orbital position. The 1 solar light curves are relatively flat due to the lack of a day/night temperature contrast as seen in Figure 4. The day/night temperature contrast is more prominent in the 50 solar case, which results in an increase in the planet/star flux ratio at secondary eclipse. The open and filled circles in Figure 10 represent output from three consecutive orbits separated by 100 simulated days to test for temporal variability in the light curves. Little variation is seen in the predicted light curves from orbit to orbit or over longer timescales.

Computed emission spectra from the points along the orbit shown in Figure 2 are presented in Figure 11 for both the 1 and 50 solar cases. As with all highly irradiated planets, the expectation is that the infrared spectra are carved predominantly by the absorption bands of HO vapor. These bands are most prominent in the near infrared between the J-, H-, and K-band flux peaks (it is at the flux peaks that the HO opacity is low). At wavelengths where HO opacity is low, the deeper, hotter layers of the atmosphere can be seen and the emitted flux is generally higher. However, other molecules can imprint absorption features on these emission peaks, and lessen them. This is true in the 4-5 m range where CO absorbs over the redder half of this range (4.5-5 m) along with CO (4.3 m), which is prominent at high metallicity. Given the assumption of chemical equilibrium, CH is abundant in all of the metallicity cases, which leads to a strong absorption band centered on 3.3 m, as well as a broad band from 7-8.5 m. However, as the atmospheric metallicity is increased, the CO abundance increases linearly, and the CO abundance increases quadratically (Lodders & Fegley, 2002; Zahnle et al., 2009b), which leads to a weakening of the CH bands and a strengthening of the CO and CO bands. In the 50 solar case a strong CO band at 4.4 m is clearly visible that does not appear in the 1 solar case. At high metallicity, this CO absorption band forces flux out a longer and shorter wavelengths, away from the 4-5 m flux peak.

In addition to Figure 11 showing emission spectra from 1 to 30 m, emitted flux distributions as a function of wavelength, with orbital phase, can also be analyzed. Figure 12 presents the planetary flux per unit wavelength for the 1 and 50 solar cases. In both the 1 and 50 solar models, the peak energy output of the planet occurs between 4 and 5 m. For the 50 model, the strong CO and CO bands just shortward of 5 m dramatically decrease the energy output there, forcing much of the energy out at both longer and shorter wavelengths. The dashed line in Figure 12 shows the integrated flux as a function of wavelength for secondary eclipse (blue dashed line, Figure 12). The flux from the planet in the 1-15 m range accounts for 95% of the planet’s emergent energy, which makes this an especially important wavelength range for determining the atmospheric properties of GJ436b.

It is interesting to note the increased variability in the emitted flux from the planet as a function of orbital position in the 50 solar case compared with the 1 solar case in Figures 11 and 12. The flux emitted from the 50 solar case at secondary eclipse (blue line, Figure 11) is lacking in many of the predominant spectral features seen at other orbital phases, due to a shallower day-side temperature gradient. The absorption features due to CH are much weaker at secondary eclipse indicating a reduction in CH abundance on the day-side compared to the night-side which is seen during transit (orange line, Figure 11). Observing the flux emitted from GJ436b as a function of wavelength at several different points along its orbit could reveal a great deal about its overall chemical composition.

4 Discussion

The atmospheric models presented here are not only useful for exploring circulation regimes, chemistry, and radiative transfer, but can also provide insight into current observations and help guide future observations of GJ436b. Figure 10 also includes the available Spitzer secondary eclipse measurements from Stevenson et al. (2010). In all metallicity cases, our predicted planet/star flux ratio falls short of the measured 3.6, 5.8, 8.0, 16.0, and 24.0 m values and is higher than the observed 4.5 m value. However, it is useful to note that 50 solar model planet/star flux ratio comes much closer to matching the observations than the 1 solar model, which could hint at a high metallicity as already suggested by Stevenson et al. (2010). High metallicity solutions (30 solar) are also favored by Spiegel et al. (2010) to match the 8 m observations of GJ436b from Deming et al. (2007). The predicted planet/star flux ratios in the 50 solar case are within two sigma of the values measured by Stevenson et al. (2010) at all bandpasses except 3.6 m.

In comparing our predicted planet/star flux ratios with those observed by Stevenson et al. (2010) it is important to remember that we have not altered the temperature or chemistry in our models in an attempt to match observations. The interplay between the equilibrium chemistry mixing ratios, the absorption and emission of flux, and the atmospheric dynamics dictates the temperature structure and emergent spectrum. As pointed out in Stevenson et al. (2010) it is likely that disequilibrium chemistry plays a strong role in the atmosphere of GJ436b. They suggest a CO/CH ratio that is many orders of magnitude larger than what one would predict from an equilibrium chemistry model. Enhancement of CO at the expense of CH is well known in the atmospheres of the solar system’s giant planets and brown dwarfs, but not to such an extreme degree. Additionally, they suggest that the photochemical destruction of CH is important to further lessen this molecule’s importance in the planet’s atmosphere. The carbon chemistry of an atmosphere will strongly affect flux measurements in the 3.6, 4.5, 5.8, and 8.0 m bandpasses. As shown in Figure 11 and discussed in Section 3.3, lowering the amount of CH in the atmosphere will increase the emergent flux from the planet in the 3.6 and 8.0 m bandpasses. Higher order hydrocarbons produced by mixing and photochemistry (Zahnle et al., 2009a), such as CH, CH, and CH, are also strong absorbers throughout the near- and mid-infrared. Increasing the amount of CO, along with CO, in the atmosphere will decrease the emergent flux from the planet in the 4.5 m bandpass while at the same time increasing the emergent flux at wavelength on either side of this bandpass, including the 3.6 and 5.8 m bandpasses.

Although our models do not include disequilibrium chemistry, they can provide an important constraint for disequilibrium chemistry models—namely, estimates of dynamical timescales and vertical mixing rates. Disequilibrium chemistry is expected in all regions of the atmosphere where dynamic timescales, , are shorter than chemical timescales, . The dynamic timescale is given simply by


where is the relevant length scale in the horizontal or vertical direction and is horizontal or vertical wind speed. For one-dimensional photochemical models, the timescales considered are only those in the vertical direction in which case , where is the vertical velocity as a function of height or pressure in the atmosphere. Given the values of from our models, the vertical velocity can be estimated as


where is the scale height of the atmosphere, which is a function of the pressure at given atmospheric level. Figure 13 shows the RMS values for calculated as global, day-side, and night-side averages for both the 1 and 50 solar cases at secondary eclipse. The RMS values for were calculated as , where the integral is a horizontal integral (at constant pressure) over the day-side, night-side, or entire globe as appropriate. In both the 1 and 50 solar cases the vertical wind speeds increase monotonically from pressures around 1 bar to 0.1 mbar with peak speeds around . The vertical eddy diffusion coefficient, , can be calculated from these vertical wind speed profiles simply by multiplying by the relevant vertical length scale, . As shown in Smith (1998), in order to calculate the correct value of at each pressure level both dynamic and chemical timescales must be considered, but to first order especially near the quench level where . Since km for GJ436b, one can expect values from 10 cm s near 100 bar to 10 cm s near 0.1 mbar. These values of are in line with those favored by Madhusudhan & Seager (2010) to explain the disequilibrium CO/CH ratio needed to fit the Stevenson et al. (2010) observations.

Regardless of the detailed chemistry, our models suggest that the basic circulation regime of planets in the temperature range of GJ436b depend strongly on overall metallicity. High metallicity cases (30–50 solar) produce a dominant eastward equatorial jet that peaks on or near the equator (Figure 5). This pattern resembles that previously obtained for more strongly irradiated hot Jupiters (Cooper & Showman, 2005; Showman et al., 2008, 2009; Rauscher & Menou, 2010). At low metallicity, on the other hand, the equatorial jet weakens and fast eastward jets develop in midlatitudes—a jet pattern distinct from those previously reported in the hot-Jupiter modeling literature. These changes suggest that the atmosphere experiences a regime shift where different dynamical mechanisms control the jet structure at low versus high metallicity. It may be that each regime is relevant to a range of planets of differing effective temperatures and compositions, providing a motivation to better understand their underlying dynamics. We will discuss the specific mechanisms responsible for this regime shift in a future paper.

The atmospheric circulation models presented here represent a first step in better understanding the dynamical, radiative, and chemical mechanisms that shape the atmosphere of GJ436b. Although chemical disequilibrium is likely to play an important role, it is important to first consider —as we have done here— a baseline atmospheric model where the effects of chemical equilibrium opacity are tested before more complex chemistry, clouds, and other factors are added and complicate the interpretation of the results. It is unlikely that the presence of disequilibrium chemistry in the atmosphere of GJ436b would strongly affect the global circulation patterns seen in this study. Madhusudhan & Seager (2010) suggest that a high metallicity (30 solar) atmosphere, with chemical disequilibrium induced both by thermal quenching and photochemistry, can best explain the observations of Stevenson et al. (2010). If this is the case, then one would expect global circulation patterns on GJ436b akin to our 30 and 50 solar models with a strong equatorial jet and non-negligible day/night temperature contrasts. As shown in Figures 10 and 11 one would expect to see variations in the planet/star flux ratio as a function of orbital phase for a high metallicity GJ436b-like atmosphere. This variation in the planet/star flux ratio could reveal a great deal about circulation patterns on GJ436b along with day/night chemistry differences. For this reason we recommend that full-orbit light curves of the planet in the 3.6 m bandpass be obtained during the Warm Spitzer mission. Looking further into the future, with the James Webb Space Telescope (JWST) a combination of NIRSpec or NIRCam data from 3-5 m, as well as MIRI data shortward of 15 m, will sample the vast majority of the planet’s emitted energy. Observations as a function of orbital phase with these instruments could tightly constrain the energy budget of the planet. Detection of the CO band depth at 4.4 m with NIRSpec would be also an important metallicity indicator and probe of the carbon chemistry of GJ436b.

5 Conclusions

In this study we have investigated atmospheric circulation for the hot Neptune GJ436b for various assumed atmospheric metallicities using a three-dimensional coupled radiative transfer and general circulation model. We have found that assumed atmospheric composition for a planet like GJ436b can have strong effects on both day/night recirculation and jet patterns within the atmosphere. Light curves and spectra produced from our models show that enhancements in atmospheric metallicity over solar produce an increase as well as orbital phase variations in the planet/star flux ratio. Given the possible strong dependence of emitted flux with orbital phase, GJ436b could provide an important probe of atmospheric chemistry outside of our solar system. Moreover, we showed that for warm gaseous extrasolar planets in the effective temperature range of 600-900 K, there exists a regime shift as metallicity is increased from a circulation dominated by mid-latitude jets and minimal longitudinal temperature differences to one dominated by an equatorial jet and large day-night temperature differences. Although these models of GJ436b do not include all the processes at play in any given atmosphere, the basic trends for atmospheric circulation as a function of metallicity will be useful in better understanding GJ436b and many other extrasolar planetary atmospheres.

This work was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program (Grant NNX08AX02H) and Origins Program (Grant NNX08AF27G). Work by KL was supported by NSF grant AST 0707377 and by the National Science Foundation, while working at the Foundation. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The authors would like to thank Yuan Lian for his helpful insights into working with the MITgcm and the anonymous referee for their helpful comments and suggestions.


  • Adams et al. (2008) Adams, E. R., Seager, S., & Elkins-Tanton, L. 2008, ApJ, 673, 1160
  • Adcroft et al. (2004) Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. 2004, Monthly Weather Review, 132, 2845
  • Bean et al. (2008) Bean, J. L., Benedict, G. F., Charbonneau, D., Homeier, D., Taylor, D. C., McArthur, B., Seifahrt, A., Dreizler, S., & Reiners, A. 2008, A&A, 486, 1039
  • Butler et al. (2004) Butler, R. P., Vogt, S. S., Marcy, G. W., Fischer, D. A., Wright, J. T., Henry, G. W., Laughlin, G., & Lissauer, J. J. 2004, ApJ, 617, 580
  • Cho et al. (2003) Cho, J. Y.-K., Menou, K., Hansen, B. M. S., & Seager, S. 2003, ApJ, 587, L117
  • Cho et al. (2008) —. 2008, ApJ, 675, 817
  • Cooper & Showman (2005) Cooper, C. S. & Showman, A. P. 2005, ApJ, 629, L45
  • Cooper & Showman (2006) —. 2006, ApJ, 649, 1048
  • Deming et al. (2007) Deming, D., Harrington, J., Laughlin, G., Seager, S., Navarro, S. B., Bowman, W. C., & Horning, K. 2007, ApJ, 667, L199
  • Demory et al. (2007) Demory, B.-O., Gillon, M., Barman, T., Bonfils, X., Mayor, M., Mazeh, T., Queloz, D., Udry, S., Bouchy, F., Delfosse, X., Forveille, T., Mallmann, F., Pepe, F., & Perrier, C. 2007, A&A, 475, 1125
  • Dobbs-Dixon & Lin (2008) Dobbs-Dixon, I. & Lin, D. N. C. 2008, ApJ, 673, 513
  • Durran (1991) Durran, D. R. 1991, Monthly Weather Review, 119, 702
  • Figueira et al. (2009) Figueira, P., Pont, F., Mordasini, C., Alibert, Y., Georgy, C., & Benz, W. 2009, A&A, 493, 671
  • Fortney et al. (2006) Fortney, J. J., Cooper, C. S., Showman, A. P., Marley, M. S., & Freedman, R. S. 2006, ApJ, 652, 746
  • Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  • Fortney et al. (2005) Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69
  • Freedman et al. (2008) Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504
  • Gautier et al. (1995) Gautier, D., Conrath, B. J., Owen, T., de Pater, I., & Atreya, S. K. 1995, in Neptune and Triton, ed. D. P. Cruikshank, M. S. Matthews, & A. M. Schumann, 547–611
  • Gillon et al. (2007a) Gillon, M., Demory, B.-O., Barman, T., Bonfils, X., Mazeh, T., Pont, F., Udry, S., Mayor, M., & Queloz, D. 2007a, A&A, 471, L51
  • Gillon et al. (2007b) Gillon, M., Pont, F., Demory, B.-O., Mallmann, F., Mayor, M., Mazeh, T., Queloz, D., Shporer, A., Udry, S., & Vuissoz, C. 2007b, A&A, 472, L13
  • Goody et al. (1989) Goody, R., West, R., Chen, L., & Crisp, D. 1989, Journal of Quantitative Spectroscopy and Radiative Transfer, 42, 539
  • Holton (2004) Holton, J. R. 2004, An Introduction to Dynamic Meterology, 4th edn. (Elsevier Academic Press)
  • Hut (1981) Hut, P. 1981, A&A, 99, 126
  • Knutson et al. (2007) Knutson, H. A., Charbonneau, D., Allen, L. E., Fortney, J. J., Agol, E., Cowan, N. B., Showman, A. P., Cooper, C. S., & Megeath, S. T. 2007, Nature, 447, 183
  • Knutson et al. (2009) Knutson, H. A., Charbonneau, D., Cowan, N. B., Fortney, J. J., Showman, A. P., Agol, E., Henry, G. W., Everett, M. E., & Allen, L. E. 2009, ApJ, 690, 822
  • Langton & Laughlin (2007) Langton, J. & Laughlin, G. 2007, ApJ, 657, L113
  • Langton & Laughlin (2008) —. 2008, ApJ, 674, 1106
  • Lodders & Fegley (2002) Lodders, K. & Fegley, B. 2002, Icarus, 155, 393
  • Lodders & Fegley (2006) Lodders, K. & Fegley, Jr., B. Chemistry of Low Mass Substellar Objects, ed. J. W. Mason (Springer Verlag), 1–+
  • Madhusudhan & Seager (2010) Madhusudhan, N. & Seager, S. 2010, ArXiv e-prints
  • Marley & McKay (1999) Marley, M. S. & McKay, C. P. 1999, Icarus, 138, 268
  • Menou & Rauscher (2009) Menou, K. & Rauscher, E. 2009, ApJ, 700, 887
  • Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics (Cambridge University Press)
  • Nettelmann et al. (2010) Nettelmann, N., Kramm, U., Redmer, R., & Neuhaeuser, R. 2010, ArXiv e-prints
  • Pont et al. (2009) Pont, F., Gilliland, R. L., Knutson, H., Holman, M., & Charbonneau, D. 2009, MNRAS, 393, L6
  • Rauscher & Menou (2010) Rauscher, E. & Menou, K. 2010, ApJ, 714, 1334
  • Rogers & Seager (2010) Rogers, L. A. & Seager, S. 2010, ApJ, 712, 974
  • Showman et al. (2008) Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559
  • Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., Marley, M. S., Freedman, R. S., Knutson, H. A., & Charbonneau, D. 2009, ApJ, 699, 564
  • Showman & Guillot (2002) Showman, A. P. & Guillot, T. 2002, A&A, 385, 166
  • Smith (1998) Smith, M. D. 1998, Icarus, 132, 176
  • Stevenson et al. (2010) Stevenson, K. B., Harrington, J., Nymeyer, S., Madhusudhan, N., Seager, S., Bowman, W. C., Hardy, R. A., Deming, D., Rauscher, E., & Lust, N. B. 2010, Nature, 464, 1161
  • Torres et al. (2008) Torres, G., Winn, J. N., & Holman, M. J. 2008, ApJ, 677, 1324
  • Zahnle et al. (2009a) Zahnle, K., Marley, M. S., & Fortney, J. J. 2009a, ArXiv e-prints
  • Zahnle et al. (2009b) Zahnle, K., Marley, M. S., Freedman, R. S., Lodders, K., & Fortney, J. J. 2009b, ApJ, 701, L20
  • Spiegel et al. (2010) Spiegel, D. S., Burrows, A., Ibgui, L., Hubeny,I.,& Milsom, J. A. 2010, ApJ, 709, 149
Figure 1: Initial pressure temperature profiles assumed for each metallicity case of GJ436b. Lines of equal abundance for CH vs CO and NH vs N are shown to highlight the dominant carbon and nitrogen bearing species at each pressure level for each metallicity case. As the metallicity in the atmosphere is increased form 1 to 50 solar, the dominant carbon bearing species changes from CH to CO. Diamonds represent the level of the mean photosphere (), which decreases in pressure as the metallicity is increased. Profiles assume planet-wide redistribution of absorbed incident energy. (a color version of this figure is available in the online journal.)
Figure 2: Orbit of GJ436b. The true anomaly, , represents the angular distance of the planet from periapse. Assuming a longitude of pericenter, , of 343 from Deming et al. (2007), transit occurs at and secondary eclipse occurs at . Dots along the orbital path represent points where data was extracted to produce Figures 10, 11, and 12. Colored dots represent points near periapse (red), transit (orange), apoapse (green), and secondary eclipse (blue), which correspond to the colored spectra presented in Figures 11 and 12. Figure is to scale with the small purple dot after periapse representing the size of GJ436b in relation to its host star and orbit. (a color version of this figure is available in the online journal.)
Figure 3: RMS velocity (colorscale) as a function of pressure and simulated time for the 1 (top) and 50 (bottom) solar cases of GJ436b. The RMS velocity at each pressure level is calculated from the instantaneous wind speeds recorded every 510 s. Simulations are continued until the winds reach a relatively flat profile at all pressure levels. (a color version of this figure is available in the online journal.)
Figure 4: Temperature (colorscale) and winds (arrows) for the 1 (left) and 30 (right) solar metallicity cases of GJ436b. For both 1 and 30 solar cases the thermal structure and winds are shown at the 1 mbar (top), 30 mbar (middle), and 1 bar (bottom) levels of the simulation. The longitude of the substellar point is indicated by the solid vertical line in each panel. Each panel is a snap shot of the atmospheric state taken near secondary eclipse (, Figure 2). The horizontal resolution of these runs is C32 (roughly 12864 in longitude and latitude) with 47 vertical layers. (a color version of this figure is available in the online journal.)
Figure 5: Zonal-mean zonal winds for the five atmospheric metallicities considered in this study for GJ436b assuming pseudo-synchronous rotation. The wind speeds presented here represent 100 day averages of the zonal winds taken after each simulation was considered to have reached an equilibrium state. The colorbar shows the strength of the zonally averaged winds in m s. Contours are spaced by 100 m s. Positive wind speeds are eastward, while negative wind speeds are westward. Note the significant change in the jet structure as a function of atmospheric metallicity.(a color version of this figure is available in the online journal.)
Figure 6: Zonal-mean temperatures (colorscale) as a function of pressure and latitude for the pseudo-synchronous 1 (top) and 30 (bottom) solar cases of GJ436b. Contours represent isotherms and are spaced by 50 K. Gradients in temperature along an isobaric surface tend to drive zonal winds according to equation (3). The 30 solar case shows stronger thermal gradients at lower pressures that extend into lower latitudes when compared with the 1 solar case. This change in the thermal gradient structure can be directly related to the zonal wind profiles seen in Figure 5.(a color version of this figure is available in the online journal.)
Figure 7: RMS velocity (colorscale) as a function of pressure and simulated time for the 1 (top) and 50 (bottom) solar cases of GJ436b as calculated from high-cadence simulation outputs every 7200 s. Vertical lines indicate periapse passage. Overall, wind speeds in the atmosphere vary with a period equal to the orbital period. The peak in the wind speeds typically occurs 4-8 hours after periapse passage. (a color version of this figure is available in the online journal.)
Figure 8: Temperature (colorscale) and winds (arrows) for the 1 (left) and 30 (right) solar metallicity cases of GJ436b assuming synchronous rotation at the 30 mbar level. The top panels represent simulations where synchronous rotation was assumed, but the nominal eccentric orbit was maintained. The bottom panels represent simulations where synchronous rotation and a circular orbit with the nominal semimajor axis (Table 1) were assumed. In the eccentric cases (top) the panels represent a snap shot taken near secondary eclipse (, Figure 2). The solid vertical line in each panel represents the longitude of the substellar point.(a color version of this figure is available in the online journal.)
Figure 9: Zonal-mean zonal winds for the synchronous rotation () cases at 1 (left) and 30 (right) solar metallicity GJ436b atmospheres. The top panels represent the jet structure for a synchronously rotating GJ436b in an eccentric orbit, while the bottom panels assume a circular orbit with the same semimajor axis. The wind speeds presented here represent 100 day averages of the zonal winds taken after each simulation was considered to have reached an equilibrium state. Note that jet structures for the synchronous cases are very similar to the zonal-mean zonal wind plots presented in Figure 5 for the same atmospheric metallicity.(a color version of this figure is available in the online journal.)
Figure 10: Planet/Star flux ratio as a function of orbital phase in each of the Spitzer bandpasses for the 1 and 50 solar metallicity cases of GJ436b. The mean anomaly is an angle that increases linearly with time from periapse passage. For GJ436b, transit and secondary eclipse occur at a mean anomalies of 90 and 303 respectively. The filled and open circles represent data extracted with 100 simulated days of separation to test for temporal variability, which appears to be minimal. The secondary eclipse measurements from Stevenson et al. (2010) are shown for comparison. (a color version of this figure is available in the online journal.)
Figure 11: Flux per unit frequency, (erg s cm Hz), as a function of wavelength for the 1 (top) and 50 (bottom) solar metallicity cases of GJ436b. The black spectra presented for each case are taken from 32 locations along a single orbit as shown in Figure 2. The central wavelengths of J-, H-, K-, L-, and M-bandpasses are indicted by the corresponding letters at the top of the plot. Dotted lines at the bottom indicate the bandpasses of the four Spitzer IRAC bands from 3-9 m, IRS blue filter at 16 m, and the MIPS 24 m bandpass. Colored spectra are taken from secondary eclipse (blue), periapse (red), transit (orange), and apoapse (green). The 50 solar case shows a strong dependence of the emergent flux density on orbital position while the 1 solar case does not.(a color version of this figure is available in the online journal.)
Figure 12: Flux per unit wavelength, (erg s cm m), as a function of wavelength for the 1 (top) and 50 (bottom) solar metallicity cases of GJ436b. The spectra presented for each case are taken from several points along the orbit as shown in Figure 2, with secondary eclipse (blue), periapse (red), transit (orange), and apoapse (green) highlighted. The dashed line represents the integrated flux as a function of wavelength from the planet at secondary eclipse and is tied to the axes on the right with 1.0 indicating the total integrated flux over all wavelengths. Dotted lines at the bottom indicate the bandpasses of the four Spitzer IRAC bands from 3-9 m. Solid gray bars at the top indicate the wavelength coverage of 3 planned instruments for the JWST: NIRCam, NIRSpec, and MIRI. The vertical position of the bars is arbitrary and does not signify instrument sensitivity.(a color version of this figure is available in the online journal.)
Figure 13: RMS vertical velocity () as a function of pressure for both the 1 (left) and 50 solar metallicity cases of GJ436b at secondary eclipse. Global (green line), day-side (red line), and night-side (blue line) averages of the vertical velocity are presented. These profiles are useful in determining the vertical mixing rate in the atmosphere as a function of pressure for use in disequilibrium chemistry and photochemical models. Note that the decrease in near the top of the domain results from the existence of the model’s upper boundary where is forced to be zero. The vertical velocities would likely continue increasing with altitude if the top of the model had been placed at even lower pressures. (a color version of this figure is available in the online journal.)
Parameter Value
(R) 0.3767
(M) 0.0729
(m s) 12.79
(AU) 0.02872
(deg) 343
(days) 2.64385
(days) 2.32851
(R) 0.464
(M) 0.452
(K) 3350

Note. – Planetary and stellar parameters taken from Torres et al. (2008). Values for and were taken from Deming et al. (2007).

Table 1: GJ436A/b parameters.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description