Methane emission from a cool brown dwarf

Sample

JWST cycle 1 GO programme 2124 (PI Faherty) is obtaining G395H NIRSpec spectra for 12 cold brown dwarfs that span late-type T and Y dwarf classes. Sources were selected by their position in the mid-infrared (Spitzer) colour–magnitude diagram (M_[4.5] versus [3.6]–[4.5] where [3.6] and [4.5] are two Spitzer bands). Four colour bins were roughly defined to represent changing temperatures (from 800 to less than 350 K). Three sources were chosen in each colour bin such that one source was close to the median property line for the population⁷. The other two were, respectively, the brightest and faintest in M_[4.5] for that colour bin. Two of the sources chosen were CWISEP J193518.59-154620.3 (W1935; ref. ³¹) and WISE J222055.31-362817.4 (W2220; ref. ³²), which are the subject of this paper. At the time of the project design, we speculated that the [3.6]–[4.5] colour would define the temperature bin and we did not expect W1935 and W2220 to be so comparable. The results of this paper are strong evidence that M_[4.5] is a better temperature binning parameter. Extended Data Fig. 1 shows the Spitzer colour–magnitude diagram for cold brown dwarfs and highlights the sample for programme 2124 with the positions of W1935 and W2220 emphasized.

CWISEP J193518.59-154620.3

CWISEP J193518.59-154620.3 (W1935 for short) was first reported in ref. ³¹ after a concerted effort was applied to find cold compact sources within newly analysed Wide-field Infrared Survey Explorer (WISE) catalogue data^33,34. The object’s discovery was the result of a collaboration between the CatWISE team and the citizen science project Backyard Worlds: Planet 9 (ref. ³⁵). A Spitzer follow-up resulted in a Spitzer [3.6]–[4.5] colour of 3.24 ± 0.31 mag, making it one of the reddest mid-infrared sources known to date. Additional follow-up observations were done to obtain Spitzer data with a higher signal-to-noise ratio (S/N) in both channels, and the colour was refined to a Spitzer [3.6]–[4.5] colour of 2.984 ± 0.034 mag in ref. ³⁶. Reference ⁷ reported that the parallax for this source was 69.3 ± 3.8 mas and that the total proper motion was 293.4 ± 16.3 mas yr⁻¹. As noted by ref. ³⁶, this source has a particularly low tangential velocity (vtan) compared to other Y dwarfs analysed. The estimated spectral type from its photometry and parallax was > Y1, and the temperature was predicted to be 367 ± 79 K (ref. ⁷). Reference ³⁷ obtained Gemini NIRI imaging for the source and reported Mauna Kea Observatories (MKO) J = 23.93 ± 0.33 mag.

WISE J222055.31-362817.4

WISE J222055.31-362817.4 (W2220 for short) was first reported in ref. ³² after a search of the WISE catalogue³⁸ for cold compact objects. Reference ³² followed this up with observations by the Spitzer Space Telescope, Keck/NIRSpec-N3 and Keck/NIRSpec-N5 as well as AAT/IRIS2 and SOAR/OSIRIS. The source was confirmed as a cold brown dwarf with WISE W1, W2, Spitzer [3.6], [4.5], and MKO JH photometry as well as near-infrared spectra. References ^39,40 obtained a Hubble Space Telescope (HST) grism spectrum for the source, which confirmed the Y0 classification with higher S/N data. Follow-up photometric and astrometric imaging was done using both space-based (for example, Spitzer) and ground-based (for example, the FourStar Infrared Camera installed on Magellan) instruments and the most recent and highest S/N trigonometric parallax is 95.5 ± 2.1 mas (see refs. ^{7,40,41,42,43} for the astrometric history of the source). Reference ⁷ estimated a temperature of 452 ± 88 K for this source based on its position on the colour–magnitude diagram.

The data

JWST programme 2124 obtained both NIRSpec G395H spectra and MIRI F1000W, F1280W, and F1800W photometry to fill out the peak of the SED and the Rayleigh–Jeans tail of the SED for 12 brown dwarfs. NIRSpec data were obtained using the F290LP filter, the G395H grating, the S200A1 aperture and the SUB2048 subarray. The resultant wavelength coverage ranged from 2.87 to 5.14 μm with a resolving power of approximately 2,700. Acquisition images were first obtained for each target using the WATA method, the CLEAR filter and the NRSRAPID readout pattern. W2220 was observed with NIRSpec on 4 November 2022 with 28 groups per integration, ten integrations per exposure and three dithers for a summation of 30 total integrations in 1,356 s of exposure time. The recorded time including the overhead for the W2220 NIRSpec observation was 1.42 h. W1935 was observed with NIRSpec on 17 October 2022 with 46 groups per integration, 11 integrations per exposure and three dithers for a summation of 33 total integrations in 2,417 s of exposure time. The recorded time including the overhead for the W1935 NIRSpec observation was 1.76 h.

MIRI photometry was obtained with the F1000W, F1280W, and F1800W filters. For each filter, the FASTR1 readout pattern was chosen with a two-point dither pattern. W1935 was observed with MIRI on 20 September 2022 using 15 groups per integration for F1000W, 13 groups per integration for F1280W and 11 groups per integration for F1800W. The total exposure time plus the overhead for the MIRI observing of W1935 was 1.03 h. W2220 was observed with MIRI on 18 October 2022 using seven groups per integration for F1000W, seven groups per integration for F1280W, and ten groups per integration for F1800W. The total exposure time plus the overhead for the MIRI observation of W2220 was 0.54 h.

As noted in ‘Sample’ above, both W2220 and W1935 have previously published photometry and W2220 has previously published spectra. Extended Data Table 1 lists all observables, both previous and new in this paper, included in the analysis that follows.

Data reduction

We used the official JWST science calibration pipeline (v.1.10.0) to reduce the NIRSpec G395H spectra starting from uncalibrated data downloaded from the Mikulski Archive for Space Telescopes (MAST). The pipeline comprises three separate stages. Stage 1 performs detector-level corrections (for example, bias subtraction, dark subtraction and cosmic-ray detection) and ramp fitting to generate a count rate image for the individual uncalibrated image of each exposure. The resulting count rate images were calibrated by applying instrument-level and observing-mode corrections in Stage 2. Stage 3 combines several calibrated exposures and extracts the spectrum. We optimized the aperture extraction location by using the relative slit position of the target to account for inaccuracies in the object coordinates and the celestial World Coordinate System.

The flux uncertainties automatically propagated through the pipeline were all null for the extracted spectrum due to the most recent reference flat files for NIRSpec having NaN values. We recalculated the flux errors for the extracted spectrum by combining in quadrature the Poisson variance (FLUX_VAR_POISSON) and read noise variance (FLUX_VAR_RNOISE) provided in the extracted spectrum file.

Looking at the W1935 3.326 μm spectral feature alone, we carefully examined each dither position to confirm that the feature was present. Although this is, in general, an area of the spectrum with a low S/N (average of 5–10 across the entire feature), we found the methane emission persisted in individual single exposure dithers, thus confirming its presence.

For both W2220 and W1935, we used the MAST-produced pipeline reductions of MIRI photometry, choosing the aper total vegamag column as our preferred magnitude.

Radial velocity

Given the high resolution of G395H data, we were able to compute radial velocities for both targets. All values reported are from a correlation with the models of ref. ⁴⁴. Our values for W1935 and W2220, which are listed in Extended Data Table 1, are −36.9 ± 5.1 and −53.2 ± 2.8 km s⁻¹, respectively. We used these values in the banyan sigma kinematic code⁴⁵ to check whether our targets belonged to any known moving group. For both objects, the full kinematics yielded a 99.9% field population probability.

We also examined the total space velocity for each object and computed values of 42.02 ± 5.33 and 55.33 ± 2.82 km s⁻¹ for W1935 and W2220, respectively. These values are consistent with the field brown dwarf population (see, for instance, ref. ⁴⁶). Consequently, we have no evidence for youth or significant old age in the observables. Therefore, we choose a broad age range of 4.5 ± 4.0 Gyr to estimate fundamental parameters in the SED approach (see below).

SED construction and results

Key to the analysis of this work was generating distance-calibrated SEDs inclusive of the newly obtained JWST data and all previously collected observables. To construct the SEDs, we used the open-source package SEDkit (ref. ¹⁴), which was first published and used by ref. ⁴⁷ to analyse the fundamental parameters of brown dwarfs. As described in ref. ⁴⁷, we first combined the parallax with spectra and photometry. For W1935, the only spectrum available was the G395H data. For W2220, an HST grism spectrum was available from ref. ³⁹ as well as the newly acquired JWST data. All photometry used in the SED for both objects is listed in Extended Data Table 1.

Using SEDkit and the input data, we constructed the SEDs, integrated under the data as described in ref. ⁴⁷ and calculated the L_bol value for each object. To examine the similarities and differences for the two objects, Fig. 1 shows the output SEDs for both W2220 and W1935 scaled to 10 pc. As described in the main text, the age estimate was paired with the evolutionary models of ref. ⁴⁸ to acquire a radius range, and we then semi-empirically calculated T_eff, log g and mass. All values are listed in Extended Data Table 1.

Figure 1 shows the G395H portion of the SEDs overplotted. The SEDs are both excellent fits for each other but show important differences. The 3.326 μm CH₄ band is striking given that it is clearly in the emission from W1935 but absent for W2220. There is also a difference in the 3.8–4.3 μm region, where W2220 appears to have the same molecular features (see Figs. 2 and 3 and Extended Data Figs. 4 and 5), but the gas is warmer. Hence, the flux is higher for W2220. which drives its bluer mid-infrared colour.

Retrieval analysis

We carried out retrieval analysis of the NIRSpec spectroscopy of our two targets using the Brewster package. Brewster is publicly available and was developed to model substellar atmospheres. It has been successfully applied to a range of brown dwarf atmospheres from cool T dwarfs to hot L-type subdwarfs, including cloudy and planetary mass objects^{15,16,17,18,19,49,50}. This is its first application in the Y dwarf regime.

Retrieval method

Brewster consists of a forward model and sampler. The forward model produces a synthetic spectrum for a given set of atmospheric parameters. This is coupled to a Bayesian sampler, which explores the parameter space and estimates the posterior distribution for the forward model parameters given the data.

The forward model is a one-dimensional, radiative-transfer, layered-atmosphere, model consisting of 64 layers distributed uniformly in log pressure space (in bars) between 2.3 and −4.0. The model evaluates the emergent flux using the two-stream source function technique of ref. ⁵¹, including scattering⁵². This requires each layer in the atmosphere to be specified in terms of the wavelength-dependent optical depth, single scattering albedo and asymmetry parameter, as well as the temperature of the gas in the layer. The resultant spectrum is then processed to allow direct comparison to the data and the calculation of a likelihood for the fit.

The parameters for the forward model fall into the following groups: gas-phase opacity, cloud opacity, temperature structure and global properties of the target. The gas-phase opacity is set by the choice of absorbing gases, their concentration and distribution in the model atmosphere. We included the following gases in this study: H₂O, CH₄, CO, CO₂, NH₃, H₂S and PH₃. In this case, we parameterized the concentrations of the absorbing gases as vertically constant volume mixing ratios as in previous studies with this code. This approach provides substantial flexibility in arriving at possible solutions while still retaining computational simplicity. We neglected cloud opacity in this study and do not discuss this aspect of Brewster any further here but leave that work for a future study.

We modelled the temperature structure using the same method as refs. ^53,54. Briefly, we parameterized the temperature at 13 evenly spaced points in log pressure and evaluated the temperature in our atmospheric layers using spline interpolation from these. We avoided implausible discontinuities and so-called ringing by penalizing the second derivative of the thermal profile in the likelihood by using the following log prior on the temperature:

This sums the discrete second derivative of the temperature T at each level i, and weights it by γ. The parameter γ sets the degree to which the likelihood will be penalized by wiggles in the thermal profile and is included in our retrieval as a free parameter. By including this free parameter, the data can set the degree of smoothing imposed on the profile. We follow ref. ⁵³ in adopting an inverse Γ distribution as the hyperprior for γ, with properties specified in Extended Data Table 3.

We parameterized the global properties of the target as follows. The radius of the source is encoded in the scale factor that is applied to the top-of-atmosphere flux produced by our radiative-transfer code to allow comparison to the flux received by NIRSpec. This scale factor is equal to R²/D², where R is the radius of the source and D is its heliocentric distance. Our forward model uses R²/D², and hence, this is the quantity that is estimated directly by the retrieval. We used our knowledge of the target’s distance to estimate the radius post hoc and incorporated the uncertainties in both the distance and the absolute flux calibration of the spectrum in the error estimate for the radius. The surface gravity for the source is parameterized as log g, where g is the gravitational acceleration at the altitude of our model atmosphere in cm s⁻². We combined this parameter with our post hoc estimate for the radius to estimate the mass of the source. Finally, we included parameters to apply rotational broadening and radial velocity shifts to our model spectrum. We used the rotational broadening code provided by ref. ⁵⁵ to achieve the former.

In addition to the free parameters described above, we hardcoded the following into our model atmosphere. In addition to the absorbing gases, we assumed that the atmosphere is dominated by H₂ gas, with He present at a solar H/He ratio. The total abundance of H₂/He was set by the remaining fraction after the absorbing gas fractions have been accounted for.

We calculated layer optical depths due to the absorbing gases using opacities sampled at a resolving power R = 30,000 taken from the compendium of refs. ^56,57, with updated opacities described in ref. ¹⁵ using the same method as ref. ¹⁶. We generated the latest version of CH₄ line list⁵⁸ with broadening coefficients relevant to H₂/He atmospheres using the computational methodology of ref. ⁵⁹ and incorporated PH₃ opacities from refs. ^48,60. We also included continuum opacities for H₂–H₂ and H₂–He collisionally induced absorption using the cross sections from refs. ^61,62 and Rayleigh scattering due to H₂, He and CH₄ but neglected the remaining gases. We also included continuum opacities due to bound–free and free–free absorption by H⁻ (refs. ^63,64) and free–free absorption by H₂⁻ (refs. ⁶⁵).

The emergent spectrum was convolved with a Gaussian kernel of width 2.2 pixels to simulate the instrumental profile of the NIRSpec G395H.

In this work, we coupled the forward model to the emcee sampler⁶⁶, which is the same method used in refs. ^{15,16,17,18,49}. We used a log-likelihood function to assess the fit of the data to the model:

where the index i refers to the ith of n spectral flux points y_i with errors s_i, which are compared to the forward model fluxes F_i for the current parameter set x. To (at least partially) account for model and other unaccounted sources of uncertainties, we inflated our errors using a tolerance parameter, such that our data error s_i is given by:

where σ_i is the measured error for the ith flux measurement and b is our tolerance parameter, which is retrieved^66,67,68. The full set of parameters used and their priors are listed in Extended Data Table 3.

We used 16 walkers per dimension. Following 10,000 iterations of burn-in, we ran our chains for blocks of 30,000 iterations, checking each time for convergence. In all cases convergence appeared to be achieved after the first 30,000 block, but we ran an additional 30,000 in each case nonetheless.

Retrieval results

We ran three models each for W1935 and W2220: (1) our baseline model as described earlier, (2) the baseline model adjusted with an additional prior such that T_i < T_i+1, which excludes the possibility of a temperature inversion and (3) the baseline model adjusted to include H₃⁺ as an absorbing gas at pressures less than 1 mbar. To distinguish the preferred model given our data, we calculated the Bayesian information criterion (BIC) for each case. For a set of models, the one with the smallest (typically most negative) BIC will be the preferred model with the strength of the preference dependent on the value of the difference in the BICs, ΔBIC. By convention, the ‘winning’ model is defined as having ΔBIC = 0 and lower-ranked models have ΔBIC quoted relative to that. Reference ⁶⁹ provided the following intervals for selecting between two models under the BIC:

• 0 < ΔBIC < 2: no preference worth mentioning;

• 2 < ΔBIC < 6: positive;

• 6 < ΔBIC < 10: strong;

• 10 < ΔBIC: very strong.

We found that for W1935, the baseline model is very strongly preferred over the model that includes H₃⁺. Both of these models are very strongly preferred over the model that rules out a temperature inversion. By contrast, for W2220 we found that there is no preference worth mentioning between models with and without a temperature inversion, which are both strongly preferred over the model that includes H₃⁺. This difference arises because W1935 shows emission in the CH₄ q-branch at 3.326 μm whereas W2220 does not. The ΔBIC values for each model are given in Extended Data Table 4.

Note that the inclusion of H₃⁺ in the model had no impact on the quality of the fit, the maximum likelihood or the values of the other parameters. The abundance of H₃⁺ was an upper limit of \(\log \,{f}_{{{\rm{H}}}_{3}^{+}}\approx -6\). We, thus, conclude that H₃⁺ is not detected in either source.

Parameter estimates (not related to the temperature structure) for W2220 are indistinguishable between the models that allow for a temperature inversion and those that do not. A comparison of the posterior distributions for models with and without a temperature inversion is shown in a corner plot in Extended Data Fig. 2.

The thermal profiles are also extremely similar, with only nonsignificant differences between the two. The two retrieved thermal profiles are compared in Extended Data Fig. 3. The two profiles are identical at pressures deeper than \(\log [P\,({\rm{bar}})]=-2.0\) and do not deviate significantly from one another at shallower pressures. Extended Data Fig. 4 shows that the retrieved model spectra are also indistinguishable and a similarly good fit to the data as the baseline model for W1935 (Fig. 2).

The posteriors for the non-temperature-related parameters for W1935 are also very similar between the two models, as shown in Extended Data Fig. 5. Although there are no significant differences between the compositions, there is a marginal trend across all the absorbing gases towards higher abundances in the baseline model (which allows for a temperature inversion).

Figures 2 and 3 demonstrate clearly that the model that allows a temperature inversion is the only one that is able to fit the CH₄ emission feature at 3.326 μm, which is absent in the spectrum of W2220.

Our retrieval results for both objects show some differences from those inferred from our SED-based luminosity and evolutionary model predictions (Extended Data Figs. 5 and 2 and Extended Data Table 1). In both cases, the masses implied by our retrievals are slightly higher than the 1σ upper limit suggested by evolutionary models. The difference for W1935 is 1.3σ, and for W2220, it is 1.1σ. It is not unusual for retrieval analyses to disagree with evolutionary models’ predictions of log g, radius and mass, particularly in the absence of any strong prior evidence to provide empirical constraints^15,16,50,70. In addition, these retrievals cover a relatively narrow wavelength range that lacks recognized gravity-sensitive spectral features. Establishing the presence and nature of any possible biases to higher gravity and mass estimates is beyond the scope of this work and does not impact our central results or conclusions.

Extended Data Fig. 6 compares the contribution functions for maximum likelihood retrieval models for W1935 and W2220. The contribution function in an atmospheric layer lying between pressures P₁ and P₂ is defined as:

C(λ, P) effectively maps the depth in the atmosphere from which the flux observed at a given wavelength originates. Extended Data Fig. 6 demonstrates that the CH₄ emission seen in W1935 originates from pressures shallower than 0.01 bar. Our model finds gas some 300 K hotter than is retrieved when no inversion is permitted.

For simplicity, we used only cloud-free retrieval models. It has been well documented that the omission of clouds in a retrieval of an atmosphere that contains clouds can bias the thermal profile to a more isothermal gradient, as the retrieval mimics the effect of clouds in screening hotter deeper layers in the opacity windows between molecular absorption features^15,16,70. Hence, the kinks seen in the retrieved thermal profiles may be suggestive that clouds should be considered for future retrieval studies in this temperature regime.