Spectrum of hot methane in astronomical objects using a comprehensive computed line list

Methane has five atoms and hence nine vibrational modes, more than are routinely treated using variational nuclear motion programs. However, it is a highly symmetric molecule belonging to the ${\mathcal{T}}_{\text{d}}$ point group which allows some simplification. The rovibrational energies, associated wave functions, and Einstein A coefficients, which are required to generate spectra, were computed using a specifically adapted version of the nuclear motion program TROVE (34) in conjunction with a new spectroscopically obtained potential energy surface (PES) and a previously calculated, ab initio dipole moment surface (DMS) (35). The energies were computed variationally, i.e., by diagonalizing a large set of (200) Hamiltonian matrices constructed using a suitable basis set in the symmetry adapted representation. The largest matrix considered (for rotational state J = 39) had 163,034 × 163,034 elements and was diagonalized using the (nonsparse) scalable linear algebra package (ScaLAPACK) direct eigensolver PDSYEVR. A direct diagonalization of large nonsparse matrices is the bottleneck in variational computation of comprehensive line lists for polyatomic molecules. These calculations were performed at the UK National Cosmology Supercomputer COSMOS and required about 1.5M CPU hours. Calculating the line intensities is also computationally expensive but more straightforward; these calculations are independent and thus can be very efficiently parallelized. This part of the project was performed at The Cambridge High Performance Computing Cluster Darwin (3.0M CPU hours). Full details of the computational method used to construct the line list will be published elsewhere (36).

To reduce computational requirements, we imposed an upper energy limit of E = hc 18,000 cm⁻¹. With this threshold for all rotational states up to J = 39, we obtained 6,603,166 energy levels. Einstein A coefficients were computed for all transitions involving rotational excitations J = 0 … 39 with the lower and upper state energies ranging up to hc 8,000 cm⁻¹ and hc 18,000 cm⁻¹, respectively, and covering the wavenumber range up to 12,000 cm⁻¹ (0.83 μm). A total of 9,819,605,160 transitions were computed. The energy values, transition wavenumbers, Einstein A coefficients, and degeneracy factors compose our 10to10 line list, which can be obtained from www.exomol.com.

The large basis set used in the diagonalization plus the combination of an accurate variational model, an empirical PES, and ab initio DMSs obtained at the high level of ab initio theory (37) all contribute to the high quality of the 10to10 line list. Further technical details are given in Supporting Information. This line list forms a key part of the ExoMol project (26), which aims to provide such lists for all molecules likely to be detected in exoplanetary atmospheres.

Typical agreement with the experimental spectra is shown in Fig. 1, where the absorption intensities for the ν₃ (stretching asymmetric) band at T = 296 K are shown as a stick diagram. Not only do the theoretical line positions agree well with the measurements, but the absolute intensities reproduce the experimental data with an accuracy comparable to that of the experiment, where the latter is available. Validation by experiment is important because it gives us confidence in our predictions of the methane opacity. Indeed, our intensities are based purely on the quality of the ab initio DMS, without any empirical adjustments.

The 10to10 line list offers at least 200 times more lines than any other previous compilation, including the most complete: the empirically constructed methane calculated spectroscopic database (MeCaSDa) (29). Most of these extra transitions come from spectral lines which involve high rotational levels and/or transitions between vibrationally excited states (hot bands). These transitions are generally very weak in room-temperature spectra but become important when the temperature increases. To demonstrate the importance of these transitions we have counted the number of active transitions as a function of temperature. At room temperature, only 1% of the 10to10 lines have intensities stronger than 10⁻³⁰ cm/molecule, which roughly corresponds to the scope of CH₄ in the HITRAN 2012 database (27). This shows that most 10to10 lines are unimportant for room temperature applications. However, at T = 1,500 K, 98% of 10to10 lines have absorption intensities stronger than 10⁻³² cm/molecule, 80% are stronger than 10⁻³⁰, and 40% are stronger than 10⁻²⁹ cm/molecule. In addition, the density of these lines is extremely high, up to 66,000 lines per cm⁻¹ at T = 1,500 K for the standard HITRAN cutoff of 10⁻²⁹ cm/molecule. This means that billions of transitions are required to model the high temperature spectrum of methane.

It is common astronomical practice to estimate the integrated flux within so-called J, H, K, L, M … colors, or bands, which correspond to the transparency windows of the Earth atmosphere in the infrared. Here we present integrated opacities of CH₄ at different temperatures for the three main methane spectroscopic windows: the J, H, and K bands, which are defined as 1.1–1.4 μm, 1.5–1.8 μm, and 2.0–2.4 μm, respectively. Fig. 2 compares 10to10 integrated intensities with estimates obtained from the HITRAN and MeCaSDa databases. Agreement is good at room temperature for the H and K, but all previous compilations show a rapid drop in the methane opacity with temperature, whereas 10to10 suggests that it is approximately flat. The 10to10 line list gives enhanced absorption in the J band for all temperatures, which is to be expected because this region is poorly sampled experimentally (20). At high temperatures (T > 1,500 K), however, the 10to10 intensities deviate from a flat line, showing incompleteness of the 10to10 line list for such temperatures. To improve the temperature coverage to, say, 2,000 K, the lower energy threshold has to be also extended at least up to 10,000 cm⁻¹, and the rotational excitations have to be extended to about J = 45.

In Fig. 3 we present the results of modeling a spectrum of the bright T4.5 dwarf 2MASS 0559-14 taken from the Infrared Telescope Facility spectral library (39, 40). At this spectral type, methane bands at around 1.6 and 2.2 μm are becoming prominent features in the spectrum. Previous analyses of the spectrum of this object have derived effective temperatures in the range 1,000–1,200 K and gravities of log g (in cgs units) = 4.5–5.0 (41, 42). The model spectrum was calculated using the versatile software for transfer of atmospheric radiation (VSTAR) and the methods described by Bailey and Kedziora-Chudczer (22). The model is based on a pressure–temperature structure for effective temperature 1,200 K and log g = 5 (43) and assumes equilibrium chemistry. The new model also uses updated absorption coefficients for the H₂–H₂ collision induced absorption (CIA) (44). The blue spectrum in Fig. 3 corresponds to the previous model which used the spherical top data system (STDS) methane line list (45), a precursor to the recent MeCaSDa line list (29). The red spectrum is the new model using the 10to10 line list: The comparison of the new VSTAR model using 10to10 with the old VSTAR model using the STDS list provides the direct comparison of the effect of just changing the opacities. Also included on the plot are results obtained using the BT-Settl model (38) for T_eff = 1,200 K and log g = 5, which is also based on the old line lists and similarly fails to match the observations in the regions of strong methane absorption.

It can be seen that the new model using the 10to10 list fits the data much better, particularly in the 1.6–1.8 μm region. In this region the STDS based lists used previously include no hot bands and clearly fail to properly represent the absorption in this region. Fig. 4 shows two expanded views of the spectral regions sensitive to methane and shows that much of the detailed line structure is reproduced by the model. The VSTAR model used here, which is a cloud-free model, fits the data for this object significantly better than the BT-Settl model, which includes clouds.

A line list as large as 10to10 presents something of a challenge for a line-by-line modeling code such as VSTAR. Initially, we tried to model the spectrum using a line intensity cutoff of 10⁻²⁷ cm/molecule at 1,500 K to reduce the number of lines used in the line-by-line code. However, this procedure was found to produce significantly less absorption than the use of the full line list. Even though the lines below this cutoff are very weak, the very large number of them means they still contribute significant total absorption. To speed up the line-by-line calculations, we have divided the lines into strong and weak lines. Only the strong methane lines are modeled with a full line shape calculation. For the many weak lines we simply add the total line absorption into the single wavelength bin in which the line center lies. Given that there are many more lines than there are wavelength bins in our spectrum, these lines are clearly never going to be individually resolved. In the models presented here the boundary between strong and weak lines was set at 10⁻²⁶ cm/molecule at the temperature of the atmospheric layer being modeled.

It remains very challenging to obtain even relatively crude spectra of exoplanets. The most productive technique has been monitoring the variation in observed starlight for those planets that transit their host star as viewed from Earth, for both the primary and secondary eclipse. Interpreting these observations requires the construction of detailed radiative transfer models. tau (46, 47) is a one-dimensional radiative transfer code for transmission spectroscopy of exoplanets, especially designed for hot Jupiters with near-stellar orbits. It uses a line-by-line integration scheme to model transmission of the radiation from the parent star through the atmosphere of the orbiting planet, equating physically to observations of the radius ratio (transit depth) as a function of wavelength in the primary transit geometry. Abundances of absorbing molecules in the atmosphere can hence be estimated by hypothesizing compositions and comparing to any available observations. Specifically, the algorithm calculates the optical depth of the atmosphere (and hence effective radius of the planet) at a particular wavelength, given model trace molecular abundances and the atmospheric structure and absorbing behavior (in the form of line lists) of those molecules. A transit depth can hence be calculated as the ratio of the squared radii of the planet and the star, and a spectrum created of absorption as a function of wavelength.

We used tau (46, 47) to generate models to simulate the atmospheric transmission of the warm Neptune GJ 436b (48) and of the hot Jupiter HD 189733b (49) (Figs. 5 and 6, respectively). Both these exoplanets have been extensively observed (6, 12, 15, 50–63). Assuming for GJ 436b an equilibrium temperature of about 700 K, the atmosphere of this planet is expected to contain significant quantities of methane as a trace absorber in chemical equilibrium (64–66) and some water. In models for GJ 436b containing only methane, the differences in the line lists result in an observed theoretical radius ratio increase Δ(R_p/R_∗) of up to 0.3% at 5 μm and below 2 μm (after smoothing). Similar models on the well-studied exoplanet HD 189733b, a hot Jupiter with an effective temperature 1,200 K, found an increase of up to 0.6% (Fig. 6). If due solely to methane opacity, this would require an increase in modeled abundance of ∼5 times using the HITRAN opacities to match the level predicted by 10to10. A large portion of these regions is, however, screened by water. This can be also seen in Figs. 5 and 6, where the results of the CH₄ + H₂O models are presented. In particular, the 5-μm regions are completely dominated by the water absorption, making the effect from the new methane data at this wavelength negligible. However, the windows around 1.6 μm and 2.2 μm show a significant difference with an increase of the radius ratio of up to 0.1% for the CH₄ (10to10) + H₂O model of the hot Jupiter HD 189733b in Fig. 6. The corresponding mixing ratios of CH₄ and H₂O were obtained by fitting to the observations available in the literature. In addition to the 10to10 and HITRAN 2012 line lists for methane, the BT2 line list (67) was used to model absorption by water. It should be noted that the large uncertainties of the existing observation data points, some of which appear to contradict each other (e.g., Fig. 5), lead to a very large number of possible solutions. Figs. 5 and 6 show only one selection of models as an example illustrating the significance of the potential errors involved in using line lists that are known to have missing opacities. These errors may not have been obvious in previous studies because the methane discrepancies are often masked to a large extent by the inclusion of realistic quantities of water at around the 10⁻⁵ abundance level.

We present a new comprehensive line list for methane suitable for modeling absorptions up to 1,500 K. An important consideration is that use of 10to10 increases opacities generated by methane, which changes the level of the atmosphere at which methane is absorbed or emitted. Future exoplanet models will have to account for the new opacities to retain self-consistency. The increased methane opacities are expected to have some effect on the pressure–temperature profile and thermal evolution of brown dwarfs and exoplanets. However, this is beyond the scope of the models considered here. Such effects can be properly evaluated when the new opacities are included in self-consistent structure models for the atmospheres and interiors. However, we would expect the changes to be relatively small as the extra opacity seen with the new line list is a relatively small fraction of the total opacity due to all sources, such as water, methane, and CIA.

Despite the size of the 10to10 line list it is still not complete. We can estimate the effect of the incompleteness of the 10to10 line list by comparison with the high-temperature partition function (68). The critical parameter here is the constraint introduced by our use of a lower energy threshold at hc 8,000 cm⁻¹ on the states that can be thermally occupied. The partition function of CH₄ computed using 10to10 levels lying below this threshold suggests that at 1,500 K these sample about 85% of the total contribution. This means that 10to10 remains incomplete and that for higher temperatures an even larger number of lines will be required for accurate modeling of the opacity of methane in hot media such as exoplanets and (cool) stars.

This work made extensive use of the high performance computing (HPC) systems Darwin and COSMOS as part of the DiRAC (Distributed Research Utilising Advanced Computing) facility. DiRAC is the UK HPC facility for particle physics, astrophysics, and cosmology and is supported by the Science and Technology Research Council (STFC) and the Department for Business Innovation and Skills (BIS) through the BIS National E-infrastructure capital Grant ST/J005673/1 and STFC Grants ST/H008586/1 and ST/K00333X/1. We thank University College London for use of the Legion High Performance Computer for performing the electronic structure calculations. This work was supported by the UK STFC, European Research Council Advanced Investigator Project 267219, and the Australian Research Council through Discovery Grant DP110103167.