The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems. V. Do Self-consistent Atmospheric Models Represent JWST Spectra? A Showcase with VHS 1256–1257 b

The data set analyzed in this study was first presented in Miles et al. (2023) as the first spectroscopic data set of a substellar object obtained through direct imaging with JWST (ERS #1386; Hinkley et al. 2022). The observations were conducted using two instruments: NIRSpec (Böker et al. 2022; Jakobsen et al. 2022) in integral field unit mode and MIRI (Wright et al. 2023) in medium-resolution spectrometer mode (Argyriou et al. 2023). This data set offers the widest wavelength coverage to date for this kind of object, spanning from ∼1 to 18 μm. It is important to note that the signal-to-noise ratio (S/N) dropped significantly in channels 4A, 4B, and 4C of MIRI, which covers wavelengths above 18 μm. As a result, this channel was not included in this analysis. The spectral resolution of the data ranges from ∼1500 to 3500 (see Figure 1 for a visualization of the resolution). The data set is composed of 12 individual channels, with three acquired by NIRSpec and nine obtained by MIRI. These channels were sequentially observed over a time span from UT 09:43:42 to 13:56:05 on 2022 July 5, corresponding to ∼20% of the rotation period of the object. The data were reduced using adapted version 1.7.2 of the standard JWST pipeline for NIRSpec (Bushouse et al. 2022a) and adapted version 1.8.1 for MIRI (Bushouse et al. 2022b). The different steps of the reduction are detailed in Miles et al. (2023).

Zoom In Zoom Out Reset image size

Since the initial application of atmospheric models to imaged planetary-mass objects in the late 1990s (Allard et al. 1996; Marley et al. 1996; Tsuji et al. 1996), various families of 1D models have been developed to replicate the spectral characteristics of progressively colder atmospheres. To address current challenges in atmospheric modeling (cloud formation and evolution, nonequilibrium chemistry, etc.), these models employ different strategies, varying in the complexity and inclusion of physical processes, as well as the range of parameters they explore. Additionally, these models can cover different wavelength ranges and operate at various spectral resolutions. In this study, we considered five distinct grids of synthetic spectra, generated by four cloudy and one cloud-free model. The rest of this section is dedicated to their introduction. The parameter spaces they investigate are summarized in Table 1. Their spectral resolutions have been estimated from the sampling of the synthetic spectra they provided, assuming a Nyquist sampling rate. They are illustrated in Figure 1.

• 1.  
  ATMO (Tremblin et al. 2015). This is a cloudless model that focuses on two nonequilibrium chemical reactions: CO/CH₄ and N₂/NH₃. These instabilities can induce diabatic convection (fingering convection) due to the difference in mean molecular weights between CO and CH₄ at the L–T transition and N₂ and NH₃ at the T/Y transition. That affects the temperature–pressure (T-P) profile in the atmosphere, reducing the temperature gradient and providing an alternative explanation for observed reddening that does not invoke clouds. The adiabatic index γ drives this temperature gradient, while log(g) determines the eddy coefficient responsible for vertical mixing (and thus nonequilibrium chemistry). This approach is debated (Leconte 2018), and the lack of clouds is now challenged by the clear detection of silicate features. But the output spectra reproduce the observations very well, in particular at the L–T transition (e.g., Petrus et al. 2023), and it remains to be understood how much this proposed ingredient could play a role in changing the T-P profile and producing redder near-infrared (NIR) and dimmer spectra, as produced by ATMO. This model incorporates 277 species, assuming elemental abundances from Caffau et al. (2011) and including those related to nonequilibrium chemistry as defined by the model of Venot et al. (2012) originally developed for hot Jupiters, while the other species precipitate under the photosphere. Opacity sources encompass collision-induced absorptions of H₂–H₂ and H₂–He, 13 molecules (H₂O, CO₂, CO, CH₄, NH₃, TiO, VO, FeH, PH₃, H₂S, HCN, C₂H₂, and SO₂), seven atoms (Na, K, Li, Rb, Cs, Fe, and H⁻), and the Rayleigh scattering opacities for H₂, He, CO, N₂, CH₄, NH₃, H₂O, CO₂, H₂S, and SO₂. This model explores five parameters: T_eff from 800 to 3000 K, log(g) from 2.5 to 5.5 dex, [M/H] from −0.6 to 0.6, C/O from 0.3 to 0.7, and γ from 1.01 to 1.05.
• 2.  
  Exo-REM (Charnay et al. 2018). This model adopts a simplified approach to microphysics to facilitate constraining and interpreting atmospheric parameters. It calculates the flux iteratively through 65 sublayers of the photosphere under the assumption of radiative–convective equilibrium. Nonequilibrium chemistry is considered, following Zahnle & Marley (2014), for a limited number of molecules (CO, CH₄, CO₂, and NH₃), and the cloud model incorporates the formation of iron, Na₂S, KCl, silicates, and water. Opacity sources include collision-induced absorptions of H₂–H₂ and H₂–He, rovibrational bands from nine molecules (H₂O, CH₄, CO, CO₂, NH₃, PH₃, TiO, VO, and FeH), and resonant lines from Na and K. Vertical mixing is parameterized using an eddy-mixing coefficient derived from cloud-free simulations. The chemical abundances of each element are defined according to Lodders (2010). By simulating the dominant physical and chemical processes and incorporating the feedback effect of clouds on the thermal structure of the atmosphere, Exo-REM reproduces the spectra of objects at the L–T transition. The parameter space explored includes T_eff ranging from 400 to 2000 K, log(g) from 3.0 to 5.0 dex, [M/H] from −0.5 to 0.5, and C/O from 0.1 to 0.8.
• 3.  
  Sonora (C. Morley et al. 2024, in preparation). Sonora Diamondback is a new grid of models within the Sonora family of models (Sonora-Bobcat: Marley et al. 2021; Sonora-Cholla: Karalidi et al. 2021). It is based on the radiative–convective equilibrium model described in Marley & McKay (1999) and used to model brown dwarfs and exoplanets (e.g., Fortney et al. 2008; Saumon & Marley 2008; Morley et al. 2012). Clouds are parameterized following the approach in Ackerman & Marley (2001). Opacities are included for 15 molecules and atoms, as well as the collision-induced opacity of hydrogen and helium, and the solar abundances from Lodders (2010) are considered. Vertical mixing in the cloud model is calculated using the mixing-length theory. Models assume chemical equilibrium throughout the atmosphere. The parameter space considered includes temperatures from 900 to 2400 K, log(g) from 3.5 to 5.5 dex, [M/H] from −0.5 to +0.5, and the cloud sedimentation efficiency parameter f_sed from 1 (fully cloudy) to 8 (cloudless).
• 4.  
  BT-Settl (Allard et al. 2012). This model estimates the abundance and size distributions of dust grains by comparing timescales of condensation, coalescence, mixing, and gravitational settling for 55 types of solids in different atmospheric layers. The radiative transfer is then calculated using the PHOENIX code. Nonequilibrium chemistry is permitted for several molecules, including CO, CH₄, CO₂, N₂, and NH₃. Vertical mixing is accounted for using the mixing-length theory under hydrostatic and chemical equilibrium. We used the CIFIST version of BT-Settl, which incorporates solar chemical abundances defined by Caffau et al. (2011). The explored parameter space includes T_eff ranging from 1200 to 7000 K (with the high end limited to 3000 K for computational efficiency) and log(g) ranging from 2.5 to 5.5 dex.
• 5.  
  DRIFT-PHOENIX (Helling et al. 2008; Witte et al. 2009, 2011). This model combines the stationary nonequilibrium cloud model called DRIFT (Woitke & Helling 2003, 2004; Helling & Woitke 2006) that incorporates growth, evaporation, and gravitational settling of seven solids (TiO₂, Al₂O₃, Fe, SiO₂, MgO, MgSiO₃, and Mg₂SiO₄) with the code PHOENIX (Hauschildt et al. 1997; Allard et al. 2001) that calculates radiative transfer and hydrostatic and chemical equilibrium and employs mixing-length theory in 256 atmospheric layers. To replenish the upper layers depleted of grains due to condensation, vertical mixing by convection is included. The grain coagulation is not taken into account. The parameter space covered by this model ranges from T_eff of 1000–3000 K, log(g) from 3.0 to 6.0 dex, and metallicity [M/H] from −3.0 to 6.0. For [M/H] = 0.0, the solar element abundances from Grevesse et al. (1992) are used, implying a fixed C/O = 0.55 in the grid.

To compare these models with the data, we used the forward-modeling code ForMoSA. ⁷¹ Initially introduced in Petrus et al. (2020), this code has been updated for the analysis of the X-Shooter data related to VHS 1256 b (Petrus et al. 2023). This updated version is the one considered in this work. This publicly available tool utilizes the Bayesian inversion technique known as "nested sampling" (Skilling 2004) to efficiently explore the complex parameter space provided by grids of precomputed synthetic spectra. By doing so, it proposes an estimation of the atmospheric properties of brown dwarfs and directly imaged planets but also other parameters such as the radius, the interstellar extinction, the radial velocity, and the projected rotational velocity. Regarding the spectral resolution of the data analyzed in this study, only the radius will be constrained here, in addition to the parameters explored by the grids. In order to optimize the analysis process, the model grids presented in this study have been converted to the xarray ⁷² format. This format allows for efficient manipulation, including interpolation and labeling of the various dimensions allowed by the grids. The spectral resolutions of both the data and synthetic spectra are compared at each wavelength to align them with optimized values. When the model has a lower spectral resolution compared to the data (as observed in certain ATMO and Exo-REM wavelengths; see Figure 1), the data's resolution is adjusted accordingly. This adjustment involves convolving both the observed and synthetic spectra with a Gaussian distribution, using the calculated FWHM to achieve the desired resolution. Subsequently, the Python module spectres ⁷³ is employed to resample the spectra onto a wavelength grid, ensuring a defined Nyquist sampling and achieving the desired spectral resolution with wavelengths. This final step optimizes computation time by limiting the number of data points considered during Bayesian inversion, all while preserving the information. The nested sampling inversion process is then carried out using the Python module nestle. ⁷⁴ For each fit, we considered priors uniformly distributed through the parameter space defined by the grids, a uniform prior between 0 and 10 R_Jup for the radius, and a Gaussian prior for the distance with μ = 21.15 pc and σ = 0.22 pc, corresponding to the Gaia DR3 measurement and error, respectively.

To assess the influence of systematic effects on the estimation of parameters, we conducted a comparative analysis between the results obtained from fitting the full wavelength coverage permitted by JWST and the results obtained from fitting different spectral windows defined along the SED. We applied each model introduced in Section 2.2 for this purpose.

In the initial step, we combined the NIRSpec and MIRI spectra. For regions where the wavelength coverage overlapped between the channels, we selected the data with the highest spectral resolution. The fitting procedure was then carried out for each model, considering a wavelength range from 0.97 to 18.02 μm. The corresponding fits are represented by unicolor solid lines in Figure 2, and the corresponding sets of parameters are summarized in Tables 4 and 5.

Zoom In Zoom Out Reset image size

Next, we divided the spectrum into 15 distinct spectral windows and conducted independent fits for each window. To account for potential calibration discrepancies between channels, we utilized the wavelength coverage of each channel to define the specific range for each window. For the NIRSpec data, we further divided them into subwindows to examine the model performance across different infrared bands, while in the case of MIRI, each channel was treated as an individual spectral window. The complete list of spectral windows can be found in Tables 2 and 3, and the corresponding best fits and parameter estimates are illustrated in Figure 2 using bicolor lines and in Tables 4 and 5, respectively.

One of the advantages of this approach is the ability to assess the contribution of each spectral window to the likelihood calculation, denoted as C_win. This allows us to identify the spectral windows that have the most significant impact on the inversion results when considering the full SED. The calculation of C_win is performed using the following formula:

$$\begin{eqnarray}&&{C}_{\mathrm{win}}\propto \displaystyle \sum _{\lambda }{\rm{S}}/{{\rm{N}}}_{\lambda }^{2},\end{eqnarray} \tag{ 1 }$$

where S/N_λ is the S/N associated with each wavelength, assuming completely uncorrelated errors. The sum of S/N² is considered to align with the definition of the likelihood function utilized during the fits, which is directly proportional to a χ² value. Consequently, C_win is maximized for spectral windows that have a greater number of data points and higher S/N. This information is illustrated in the top right panel of Figure 2. In the case of JWST's data, the fit is primarily driven by the NIRSpec G395HF/290LP channel, which contributes approximately 60% to the overall fit. On the other hand, all MIRI data collectively contribute around 0.26% to the overall fit.

In this approach, we also gain the ability to estimate the sensitivity of the models to each parameter depending on the spectral window used. For a given spectral window Δλ, model m, and parameter p, we calculated this sensitivity S_Δλ,m,p using the following relation:

$$\begin{eqnarray}&&{S}_{{\rm{\Delta }}\lambda ,m,p}=\displaystyle \sum _{\lambda }\left(\displaystyle \frac{\mathrm{abs}({F}_{m,\lambda }({p}_{\max })-{F}_{m,\lambda }({p}_{\min }))}{{F}_{m,\lambda }({p}_{\max })}\right)\times \displaystyle \frac{1}{{N}_{{\rm{\Delta }}\lambda }},\end{eqnarray} \tag{ 2 }$$

where ${F}_{m,\lambda }({p}_{\min })$ and ${F}_{m,\lambda }({p}_{\max })$ are the synthetic flux generated by the model m at the wavelength λ with the minimum and maximum values of the parameter p explored by m, respectively, and N_Δλ is the number of data points contained in Δλ. The remaining parameters are estimated using the posteriors of the fits corresponding to the specific spectral window and model considered. Therefore, these index values are only relevant for objects with atmospheric properties similar to VHS 1256 b. If, for a given model m, the spectral information F_m,λ (p) within the considered window is not affected significantly by the parameter p, the difference ${F}_{m,\lambda }({p}_{\max })$ – ${F}_{m,\lambda }({p}_{\min })$ (and consequently S_Δλ,m,p ) will tend toward 0. Conversely, a substantial S_Δλ,m,p value will result if there is significant spectral variation within the window corresponding to variations in the parameter p. The sensitivity is visualized in Figures 3–5 through the size of the data points and will be further discussed in the subsequent sections.

Zoom In Zoom Out Reset image size

All five models exhibit deviations from the observed data when the entire wavelength range is considered for the fit, particularly at longer wavelengths (λ > 7 μm), which carry less weight in the likelihood calculation. To assess the goodness of fit, we calculated the reduced χ² values for the fits using the full SED, denoted as ${\chi }_{r,\mathrm{full}}^{2}$, as well as for the fits using the spectral windows, denoted as ${\chi }_{r,\mathrm{win}}^{2}$ (see Figure 2). For the latter case, we constructed a synthetic SED by combining the best fit from each spectral window, following the merging procedure described earlier, enabling a direct comparison with the fit of the full SED. Among the models tested, the Exo-REM model demonstrates the best performance, yielding ${\chi }_{r,\mathrm{full}}^{2}$ = 246 and ${\chi }_{r,\mathrm{win}}^{2}$ = 57. The utilization of smaller wavelength ranges in the independent fits (spectral windows) results in a decrease in the overall ${\chi }_{r,\mathrm{win}}^{2}$ from a factor of 1.5 for the DRIFT-PHOENIX model to 6 for the BT-Settl model. This improvement is not surprising, as this procedure is similar to a fit in which the number of free parameters is multiplied by the number of spectral windows. Consequently, there is a variation in the estimated parameters depending on the specific spectral window considered.

As illustrated in Figures 3–5, the high S/N provided by JWST allows for precise estimation of the atmospheric properties of VHS 1256 b in each fit, resulting in very small error bars. However, the significant dispersion observed in these estimates across the different spectral windows indicates that these errors cannot be representative of this dispersion. Furthermore, as depicted in these figures, the sensitivity of the models to each parameter varies depending on the spectral window considered. We propose here a procedure that provides a robust estimate of each parameter Θ_win,f, where Θ represents the parameter of interest. This procedure considers both the dispersion and the models' sensitivity of each parameter to redefine a more robust error bar.

For each parameter estimated with each model, we defined a sample composed of subsamples of points randomly drawn from the posteriors obtained with each spectral window. We excluded those affected by the silicate absorption between 7.5 and 10.5 μm, which are known to not be reproduced by the current generation of models. The number of points in each subsample is determined proportionally based on the calculated sensitivity S_Δλ,m,p using Equation (2). The final value and error of each parameter are defined as the mean and standard deviation, respectively, of this sample and are presented in Figures 3–5 as the solid gray line and gray area, respectively. They are also reported in Tables 4 and 5.

3.2.1. The Effective Temperature T_eff

3.2.2. The Surface Gravity log(g)

3.2.3. The Metallicity [M/H]

The exploration of [M/H] is performed using the grids Exo-REM, Sonora, DRIFT-PHOENIX, and ATMO. When the full wavelength range is considered, the ATMO model indicates a solar [M/H], while the Sonora model suggests a supersolar metallicity. On the other hand, both the Exo-REM and DRIFT-PHOENIX models converge toward a subsolar metallicity (Figure 4). These results vary when considering the spectral windows. In this case, the Exo-REM model estimates a solar [M/H], while the Sonora model suggests a supersolar [M/H]. Conversely, both the DRIFT-PHOENIX and ATMO models converge toward a subsolar [M/H]. The three identified regimes of metallicity are further confirmed by the estimates of [M/H]_win,f. Similar to the log(g) parameter, the [M/H] estimation is primarily constrained by the NIR data (λ < 5 μm).

Zoom In Zoom Out Reset image size

3.2.4. The Carbon–Oxygen Ratio C/O

3.2.5. The Adiabatic Index γ and the Sedimentation Factor f_sed

The parameters γ and f_sed are specific to the ATMO and Sonora models, respectively (see Figure 5). When considering the NIRSpec data, their estimates exhibit distinct trends, with a medium/high value for γ and a low value for f_sed. However, these parameters display similar variations when analyzed with the MIRI data. Specifically, there is an initial increase in both parameters up to λ ∼ 9 μm, followed by a sharp decline in the MIR range. Similarly to C/O, the new errors for γ_win,f and f_sed,win,f both fail to represent the dispersion adequately. The f_sed values are particularly intriguing, indicating a minimal impact of clouds at λ < 3 μm (f_sed ∼ 2), a significant impact at λ = [3–10] μm (f_sed > 5), and no discernible impact at λ > 10 μm (f_sed ∼ 1). This may reveal the complexity of the cloud clover in the atmosphere of VHS 1256 b, which stands at the L–T transition. Additionally, the estimation of γ (∼1.03) and f_sed (∼1) using the full SED aligns coherently with the results observed for late L objects (Stephens et al. 2009; Tremblin et al. 2017).

Zoom In Zoom Out Reset image size

3.2.6. The Radius R

The radius is not a parameter explored by the model grids. Instead, it is estimated using the dilution factor C_K = (R/d)², which takes into account the flux dilution between the synthetic spectra generated at the outer boundary of the atmosphere and the observed spectrum. Here, R represents the planet's radius and d represents the distance to the planet. To propagate the error from the distance to the radius, we define a Gaussian prior for the distance centered at 21.14 pc with a standard deviation of 0.22 pc, according to the estimate from Brown et al. (2021). We impose a flat prior to the radius to leave it unconstrained. The results are illustrated in Figure 6. It is observed that the radius estimated by each model is consistently underestimated, except for the Exo-REM grid, which provides a radius in good agreement with the predictions of evolutionary models. When the full SED is used for the fit, the radius estimates from the Sonora and ATMO grids also align with the results from Dupuy et al. (2023) and Miles et al. (2023). However, a significant dispersion remains, resulting in substantial error bars for the final values (see Table 5).

Zoom In Zoom Out Reset image size

In one of our previous works, Petrus et al. (2021) demonstrated the challenges associated with fitting spectral features detected at medium resolution using precomputed grids of synthetic spectra. They highlighted that even small inhomogeneities in the grids could remain and propagate through interpolation and introduce biases in the resulting posteriors. Reproducing the depth of the absorption features and accounting for systematics in the continuum shape were also identified as potential sources of bias. To mitigate these issues, we have developed a two-step procedure that analyzes each spectral feature independently.

We initially conducted a fit on the wavelength ranges, Δλ_continuum, that avoided the considered spectral features (red spectral windows in Figure 7). From this fit, we extracted the T_eff, which was then used to define a local pseudocontinuum at the positions of the spectral features. This approach serves as an alternative to the definition of the pseudocontinuum as a low-resolution version of the spectrum (e.g., Petrus et al. 2023), which may be suitable for higher spectral resolution data but could be unreliable for the resolution of JWST (R_λ < 3000). Additionally, we extracted the dilution factor C_K to ensure that the depth of the absorption is solely related to the atmospheric properties explored by the grid. In the second fit, we only considered the wavelength ranges, Δλ_fit, which focused on the spectral features (green spectral windows in Figure 7). We fixed the values of T_eff and C_K to the estimates obtained from the first fit. All spectral windows considered are summarized in Table 6.

Zoom In Zoom Out Reset image size

For each spectral feature, we illustrate in Figure 7 the synthetic spectrum interpolated from the model grids that maximizes the likelihood when compared to the observed data. The reduced ${\chi }_{{\rm{r}}}^{2}$, calculated using Δλ_fit, is provided to quantify and compare the quality of each fit. Its value greater than 1 suggests that the models may have limitations in replicating these spectral signatures and/or that there may have been an underestimation of the error bars during data extraction. The corresponding set of parameters is represented in Figure 8. Similar to the method described in Section 3, we have observed variations in the estimates of atmospheric parameters depending on the spectral feature considered for the fit, as well as the model used. Globally, the models are capable of reproducing this information at medium resolution. However, certain spectral features appear to be more challenging to reproduce, such as the potassium line at 1.243 μm. Now, let us delve into the detailed performance of each individual model.

Zoom In Zoom Out Reset image size

4.2.1.  ATMO

4.2.2.  Exo-REM

4.2.3.  Sonora

4.2.4.  BT-Settl

4.2.5.  DRIFT-PHOENIX

This JWST spectrum covers 98% of the bolometric flux of VHS 1256 b. By filling the remaining 2% with a model and integrating this full SED, Miles et al. (2023) calculated a log(L_bol/L_⊙) of −4.550 ± 0.009, which was consistent with previous estimates provided by Hoch et al. (2022) and Petrus et al. (2023). In our study, we derived the bolometric luminosity directly from our estimates of T_eff and the radius R. For each spectral window defined in Section 3.1, we calculated L_bol using the Stefan–Boltzmann law: ${L}_{\mathrm{bol}}=4\pi \sigma \times {R}^{2}\times {T}_{\mathrm{eff}}^{4}$, where σ is the Stefan–Boltzmann constant. These results are visualized in Figure 9, and as expected, we observe dispersion across the different spectral windows, influenced by the dispersion observed for the T_eff and the radius estimates. The fits using the full SED yielded bolometric luminosity in good agreement with the estimates of Miles et al. (2023) for each model but with a very small relative error (less than 0.03%). Using the same method described in Section 3.2, we calculated robust values of L_bol for each model and reported them in Table 5. Exo-REM and ATMO are the two models that provide the most stable L_bol over wavelength coverage.

Zoom In Zoom Out Reset image size

We also calculated the mass from our estimates of the log(g) and the radius, applying the gravitational law:

$$\begin{eqnarray}&&m\,g=\displaystyle \frac{{\rm{G}}\,M\,m}{{R}^{2}}\longleftrightarrow M=\displaystyle \frac{g\,{R}^{2}}{{\rm{G}}},\end{eqnarray} \tag{ 3 }$$

where G is the gravitational constant. As with the bolometric luminosity, we also observe dispersion in the mass estimates (Figure 9). We take this dispersion into account to provide the final values reported in Table 5. Interestingly, for the grids BT-Settl, Exo-REM, and ATMO, the mass values obtained from each spectral window appear to be too low and unphysical (<5 M_Jup) when compared to the previous estimate made by Miles et al. (2023). They considered the age of the system, the measured L_bol, and the hybrid cloudy evolutionary model of Saumon & Marley (2008) to estimate the mass of VHS 1256 b to be lower than 20 M_Jup with two potential solutions: one at 12 ± 0.5 M_Jup and another at 16 ± 2 M_Jup. Nonetheless, this comparison must be approached with caution due to the inherent uncertainties within the evolutionary models, encompassing factors such as initial entropy variations and the uncertainty in age determinations, including the potential for age discrepancies between the star and its companion. The trend of log(g) as shown in Figure 3 seems to impact the Sonora grid, indicating a higher mass when using MIRI data and a lower mass when using NIRSpec data. Lastly, the substantial dispersion in log(g) obtained with the DRIFT-PHOENIX grid results in a wide range of mass estimates.

The results presented in Sections 3 and 4 have revealed that the estimates of atmospheric properties for VHS 1256 b, derived using the forward-modeling approach, exhibit a significant dependence on several factors. These include the wavelength coverage, the specific model used, and the type of information considered for the fit (large spectral windows or specific spectral features). The dispersion observed in Figures 3–6 can be attributed to three main sources: the object's characteristics, the data reduction, and the models' performances and properties (chemical abundances, physics considered).

First, VHS 1256 b has been identified as the most variable substellar object by Bowler et al. (2020) and Zhou et al. (2022), showing a luminosity variation of over 20% between 1.1 and 1.7 μm over an 8 hr period and ∼5.8% at 4.5 μm (Zhou et al. 2020). This variability has been attributed to an inhomogeneous cloud cover. A spectroscopic time follow-up of this object is underway to analyze its cloud and atmospheric structure, and the results will be presented in dedicated papers. Using the three-sinusoid model proposed by Zhou et al. (2022), Miles et al. (2023) conducted a simulation to estimate the expected variability during the 4 hr time acquisition of the JWST data used in this study. They estimated a maximum variability of less than 5% for 1 μm < λ < 3 μm, about 1.5% for 3 μm < λ < 5 μm, and negligible variability for λ > 5 μm. However, Figure 9 reveals an unexpected increase in the best-fit bolometric luminosity as a function of wavelength for fitted spectral channels between 1 and 7 μm, with approximate percentage changes of 330%, 390%, 400%, 40%, and 210% for the grids BT-Settl, Sonora, DRIFT-PHOENIX, Exo-REM, and ATMO, respectively, which cannot be attributed to the variability of VHS 1256 b.

An alternative explanation for this luminosity dispersion might be attributed to the extraction methodology. The data set employed in this analysis comprises distinct channels observed independently (three for NIRSpec and nine for MIRI). As the spectral coverage provided by these channels is employed to define the spectral windows used in the fitting process in Section 3, inaccuracies in flux calibration during extraction could potentially influence the inferred luminosity. However, it appears improbable, given that this discrepancy is discernible between the sub-windows of each NIRSpec channels, despite a similar extraction process, and therefore a similar flux calibration (see Table 2). Furthermore, the magnitude of the dispersion varies across models, indicating a model-dependent interpretation.

Previous studies have invoked systematic errors within models of atmospheres to explain biases in estimating atmospheric parameters for early-type and hotter objects. Petrus et al. (2020) attributed them to deficiencies in dust grain modeling, which were partially mitigated by introducing interstellar extinction (A_V ) as an additional free parameter in the fitting process. This approach was further employed by Hurt et al. (2023) in the characterization of 90 late M and L dwarfs. They confirmed that the fit quality significantly improved when incorporating an extra source of extinction. Notably, no interstellar extinction was detected for these observed targets. While both Petrus et al. (2020) and Hurt et al. (2023) focused on data within the wavelength range λ < 2.5 μm, we extended this concept to a broader wavelength range. We incorporated the interstellar extinction law as defined by Draine (2003) as a free parameter in the fitting process, using our five different models.

The results, as presented in Appendix B, demonstrate that exploring the impact of an interstellar extinction leads to an enhancement in the quality of the fits for the models BT-Settl, Exo-REM, and ATMO, evident through reduced ${\chi }_{{\rm{r}}}^{2}$ values. This improvement is not observed in the case of DRIFT-PHOENIX and Sonora, where the calculated ${\chi }_{{\rm{r}}}^{2}$ values across the entire SED are higher when considering A_V . This phenomenon can be explained by the significant heterogeneity inherent in the data set concerning spectral resolution and S/N, as illustrated in the top right panel of Figure 2. Notably, focusing exclusively on the channel G395HF/290LP of NIRSpec (2.87–5.27 μm), which holds the greatest sensitivity in the likelihood computation, yields a reduction in ${\chi }_{{\rm{r}}}^{2}$ values (${\chi }_{{\rm{r}}}^{2}$ = 126 and ${\chi }_{{\rm{r}},{A}_{{V}}}^{2}$ = 116 for DRIFT-PHOENIX and ${\chi }_{{\rm{r}}}^{2}$ = 297 and ${\chi }_{{\rm{r}},{A}_{{V}}}^{2}$ = 148 for Sonora). With each model, substantial deviations persist when considering A_V , particularly at λ > 5 μm, underscoring that the considered extinction law, derived from a population of dust grains within the interstellar medium, fails to accurately capture the dust and cloud characteristics that shape the atmosphere of VHS 1256 b. We have also observed a negative value of A_V when the model Exo-REM is considered, indicating that this model is excessively red. This underscores the critical significance of advancing the development of extinction laws specifically tailored for characterizing atmospheric dust and haze.

Moreover, the luminosity depicted in Figure 9 is calculated based on the estimated values of T_eff and radius considering the Stefan–Boltzmann law. Consequently, the dispersion of the luminosity within the 1–7 μm range directly correlates with variations in these two parameters. Upon comparing the dispersions of T_eff and radius illustrated in Figures 3 and 6, respectively, it becomes apparent that the T_eff estimates remain relatively consistent across all models, except for Sonora. Conversely, the radius estimates show a general increase across all models, except for Sonora. One plausible explanation for this behavior is the potential variance in cloud cover between the Sonora model and the remaining models, attributed to the parameter f_sed, which governs the sedimentation process of the diverse clouds generated (Luna & Morley 2021). Indeed, it is expected that different cloud coverage would strongly impact the shape of the SED at low spectral resolution (N. Whiteford et al. 2024, in preparation) but also should be critical in the reproduction of the spectral information at higher resolution (see Section 4) and consequently in the robust estimate of chemical abundances. Moreover, the cloudless model ATMO, grounded in chemical disequilibrium physics, offers good performances in reproducing the data. We suggest that this cloudless approach should be integrated with cloud models, as a synergistic interplay between the two factors is plausible and can yield combined effects. Indeed, in the case of VHS 1256 b, the presence of clouds in the photosphere is directly detected from the spectra, as evidenced by the silicate absorption. This implies that a fully cloud-free model is evidently not self-consistent in reproducing the atmosphere of this kind of object.

Lastly, it is important not to disregard the potential influence of the spectral covariance and the choice of Bayesian priors, which are known to impact the estimation of atmospheric parameters (Greco & Brandt 2016). In this work, we have opted for a conservative approach by defining flat priors constrained by the size of the grids considered. Some inversion codes, like STARFISH (Czekala et al. 2015), have integrated the impact of an error covariance matrix into likelihood calculations. This code also quantifies the influence of systematic errors in models using a Markov Chain Monte Carlo framework. In our future work, we plan to implement and test these functionalities while considering nested sampling in our code ForMoSA (M. Ravet et al. 2024, in preparation).

JWST's data provided the first detection of silicate absorption in the atmosphere of a low-mass companion (Miles et al. 2023), directly revealing the presence of clouds in the atmosphere of VHS 1256 b. This feature had already been detected in the atmosphere of isolated objects, in particular at the L–T transition (e.g., Cushing et al. 2006; Looper et al. 2008; Suárez & Metchev 2022). At these temperatures (∼1200–1700 K), clouds of forsterite (Mg₂SiO₄), enstatite (MgSiO₃), and quartz (SiO₂) can condense following nonequilibrium chemistry (Lodders 2002; Helling & Woitke 2006). This has been detected empirically using atmospheric retrieval methods (Burningham et al. 2021; Vos et al. 2023) and will be explored in the retrieval analysis of the data set of VHS 1256 b used in this paper (N. Whiteford et al. 2024, in preparation). In their work, Suárez & Metchev (2022) used a library of low-resolution spectra (R_λ ∼ 90), obtained with Spitzer, to define a spectral index in order to investigate the depth of the silicate absorption as a function of the spectral type. Their library was composed of ∼110 field and young isolated planetary-mass objects covering a wide range of spectral types from M5 to T9. Their method was based on multiple linear interpolations to estimate the pseudocontinuum and the absorption depth at the position of the feature (at 9 μm). The spectral index was defined as the ratio of both fluxes. Thanks to the diversity of their spectral library, they identified positive correlations between their silicate index and the NIR color excess and between their silicate index and the photometric variability, proving that silicate clouds had a great impact on the atmospheric properties of planetary-mass objects. The silicate index definition was improved in Suárez & Metchev (2023) considering the silicate absorption extends up to 13 μm in relatively young objects, like VHS 1256 b, and a better approach to estimate the continuum. This led to the main conclusion that the silicate absorption is sensitive to surface gravity. It is redder, broader, and more asymmetric in low surface gravity dwarfs. This work was followed by Suárez et al. (2023), who identified, in the same Spitzer library, a dependence between the depth of the silicate absorption and the viewing geometry of the target, confirming that equatorial latitudes are cloudier. In the case of VHS 1256 b, which is seen equator-on (90°${}_{-28}^{+0}$; Zhou et al. 2020), they find a strong silicate absorption and reddening. However, the inability to date to observe the silicate absorption between 8.0 and 12.0 μm at medium resolution (R_λ > 100) has restricted the comparison of this molecular absorption's properties (depth, shape, etc.) with atmospheric model predictions. This limitation has slowed down the development of these models, as they currently cannot replicate this feature. Consequently, our comprehension of the physical and chemical mechanisms underlying its appearance, notably its chemical composition, remains hindered. In this section, we propose a new method based on the forward-modeling analysis to define a new silicate index.

5.3.1. Method

To avoid the region where the models fail to reproduce the silicate absorption, we defined two windows on the left side (6.0–8.0 μm) and the right side (12.0–14.0 μm) of the feature. The data were fitted using ForMoSA, considering only the wavelengths covered by these two windows. The parameters estimated from this fit were then used to estimate the pseudocontinuum between 8 and 12 μm. We performed this fit for each model defined in Section 2.2, as well as for each isolated object in the Spitzer library and for the JWST data of VHS 1256 b. We limited the Spitzer library to objects with spectral types between M6 and T8. The final sample consisted of 12 young objects (age <400 Myr) and 93 field objects. The index I_Si was calculated as the equivalent width of the silicate absorption using the formula

$$\begin{eqnarray}&&{I}_{\mathrm{Si}}=\mathrm{EW}(\mathrm{Si})={\int }_{\lambda }\left(1-\displaystyle \frac{{F}_{\lambda }}{{C}_{\lambda }}\right)d\lambda ,\end{eqnarray} \tag{ 4 }$$

where C_λ and F_λ are the fluxes of the estimated pseudocontinuum and the data, respectively, at the wavelength λ, with λ ∈ [8.0, 12.0] μm. The error bar on I_Si is estimated by the calculation of ${I}_{\mathrm{Si},\min }$ and ${I}_{\mathrm{Si},\max }$, considering F − F_err and F + F_err, respectively, where F_err is the error of the data. This method is illustrated in Figure 10.

Zoom In Zoom Out Reset image size

5.3.2. Results

In their work, Suárez & Metchev (2022) found that the silicate index was higher on average between the spectral types L4 and L6 compared to colder and warmer objects in their library. They interpreted this phenomenon as the result of Si being unable to condense into clouds for T_eff > 2000 K and the sedimentation of clouds below the photosphere after the L–T transition. Figure 11 displays I_Si as a function of spectral type for each object in the Spitzer library, including VHS 1256 b. The top right panel reproduces the silicate index as defined by Suárez & Metchev (2023), enabling a direct comparison between the two approaches. Similar results were obtained for each model used. An I_Si ∼ 0 is calculated for M-type objects, indicating a lack of significant silicate absorption in this spectral range. Additionally, the detection of CH₄ after the L–T transition, as well as the potential contribution of the NH₃ absorption at 10.5 μm, which is present in the spectra of T dwarfs, may lead to negative values for the estimated equivalent widths. However, it should be noted that the validity of the results from the five models is limited to the spectral types they cover, as indicated by the colored area. Furthermore, we observe that the silicate index is generally higher for young objects compared to field objects, as recently highlighted by Suárez & Metchev (2023). This is likely caused by the higher surface gravity of older objects, resulting in their contraction, which increases the sedimentation rate of the silicate clouds below the photosphere.

Zoom In Zoom Out Reset image size

The silicate index calculated for VHS 1256 b, following the equivalent-width procedure, aligns with the observed trends. The L7 spectral type and young age of this object favor the presence of a silicate absorption feature, resulting in a high value of I_Si for each model. As observed in Figure 10, the absorption feature of VHS 1256 b seems to be redder and broader compared to the absorption detected in the spectrum of other objects (e.g., 2MASS 2148+4003). The reason for this distinct shape needs to be investigated, but it could be attributed to complex silicate clouds with varying abundances of different silicate elements. Indeed, Cushing et al. (2006) showed that forsterite (Mg₂SiO₄) was redder than enstatite (MgSiO₃). The in-depth analysis of the silicate feature is a relatively new research direction that is currently in progress (Burningham et al. 2021; Suárez & Metchev 2023; Suárez et al. 2023). With this analysis, we have demonstrated that the current generation of atmospheric models can be utilized to accurately estimate the pseudocontinuum at the position of this feature, taking advantage of their inability to reproduce the silicate absorption feature. This allows for robust quantification of the equivalent width and facilitates discussions on the spectral type and age of observed objects.

This work highlighted the current limitations of forward modeling in characterizing the atmospheres of low-mass objects. One approach to overcome these limitations is to improve the quality of the existing generation of atmospheric models. As demonstrated in this study, it is crucial to compare the outcomes obtained by considering various model families. Therefore, simultaneous improvements are necessary for each model. Various enhancement points have been identified.

The physics governing the clouds (formation and evolution) requires improvement. Currently, the particle size distribution is defined by a single broad lognormal distribution, which may need adaptation because it suggests smaller particles at higher atmospheric levels. Incorporating microphysical models that present smaller particles in a "nucleation mode" and larger particles in a "growth/settling mode" could be an interesting solution (see Helling & Woitke 2006). Having a cloud model with smaller particles at the right heights might aid in reproducing the ∼10 μm silicate feature and potentially explain some of the observed interstellar reddening from optical to NIR wavelengths.

As demonstrated by ATMO, diabatic convection can replicate cloud effects in the atmosphere, substantially altering temperature structures, particularly at the L–T transition. However, the presence of clouds is now confirmed through the detection of silicate absorption. Simultaneous integration of disequilibrium and clouds marks a promising advancement, slated for implementation in Sonora and refinement in Exo-REM. Additionally, introducing physically motivated adjustments to temperature structures, including effects like mean molecular weight gradients not presently considered, could impact the T-P profile in the mid- and upper atmosphere, influencing cloud physics and estimates of chemical abundances.

All the model families considered in this study are 1D, combining simulations of the vertical atmosphere structure with radiative transfer codes to generate synthetic spectra. VHS 1256 b exhibits significant variability, suggesting a complex atmosphere potentially marked by an inhomogeneous cloud cover. It is evident that 1D models alone will not suffice to capture the complexity of these objects existing in three dimensions. Currently, significant efforts are underway to interpret such variability by developing a new generation of 3D models of atmospheres (Showman et al. 2020; Tan & Showman 2021a, 2021b; Plummer & Wang 2022) that aim to map the surfaces of exoplanets observed by JWST and, in the near future, by the Extremely Large Telescope (ELT).

Finally, the recent identification of isotopologue ratios as tracers of formation motivates the development of models that can explore their potential different values. Recently, isotopologues (¹³C, ¹⁸O, and ¹⁷O) were identified in the atmosphere of VHS 1256 b through comparisons of JWST data with synthetic spectra generated within a retrieval framework (Gandhi et al. 2023). It is important to investigate the capability of self-consistent models in detecting these chemical elements, estimating their abundance, and studying the impact of systematic errors in the models on these values. Expansion of grids is currently underway to address this new avenue of research.