# 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

Simon Petrus, Niall Whiteford, Polychronis Patapis, Beth A. Biller, Andrew Skemer, Sasha Hinkley, Genaro Suárez, Paulina Palma-Bifani, Caroline V. Morley, Pascal Tremblin, et al.

► To cite this version:

Simon Petrus, Niall Whiteford, Polychronis Patapis, Beth A. Biller, Andrew Skemer, et al. 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 Astrophysical Journal Letters, 2024, 966, 10.3847/2041-8213/ad3e7c. insu-04726505

HAL Id: insu-04726505

https://insu.hal.science/insu-04726505v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# 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

Simon Petrus<sup>1,2,3</sup>, Niall Whiteford<sup>4</sup>, Polychronis Patapis<sup>5</sup>, Beth A. Biller<sup>6,7</sup>, Andrew Skemer<sup>8</sup>, Sasha Hinkley<sup>9</sup>, Genaro Suárez<sup>4</sup>, Paulina Palma-Bifani<sup>10,11</sup>, Caroline V. Morley<sup>12</sup>, Pascal Tremblin<sup>13</sup>, Benjamin Charnay<sup>11</sup>, Johanna M. Vos<sup>4,14</sup>, Jason J. Wang<sup>15,16</sup>, Jordan M. Stone<sup>17</sup>, Mickaël Bonnefoy<sup>18</sup>, Gaël Chauvin<sup>10</sup>, Brittany E. Miles<sup>8</sup>, Aarynn L. Carter<sup>8</sup>, Anna Lueber<sup>19,20</sup>, Christiane Helling<sup>21,22,23</sup>, Ben J. Sutlieff<sup>6,7</sup>, Markus Janson<sup>24</sup>, Eileen C. Gonzales<sup>25,26,70</sup>, Kielan K. W. Hoch<sup>27</sup>, Olivier Absil<sup>28</sup>, William O. Balmer<sup>27,29</sup>, Anthony Boccaletti<sup>11</sup>, Mariangela Bonavita<sup>30</sup>, Mark Booth<sup>31</sup>, Brendan P. Bowler<sup>12</sup>, Zackery W. Briesemeister<sup>32</sup>, Marta L. Bryan<sup>33</sup>, Per Calissendorff<sup>34</sup>, Faustine Cantalloube<sup>35</sup>, Christine H. Chen<sup>27,29</sup>, Elodie Choquet<sup>35</sup>, Valentin Christiaens<sup>28</sup>, Gabriele Cugno<sup>34</sup>, Thayne Currie<sup>36,37</sup>, Camilla Daniëlski<sup>38</sup>, Matthew De Furio<sup>34</sup>, Trent J. Dupuy<sup>6</sup>, Samuel M. Factor<sup>12</sup>, Jacqueline K. Faherty<sup>4</sup>, Michael P. Fitzgerald<sup>39</sup>, Jonathan J. Fortney<sup>8</sup>, Kyle Franson<sup>12</sup>, Julien H. Girard<sup>27</sup>, Carol A. Grady<sup>40</sup>, Thomas Henning<sup>41</sup>, Dean C. Hines<sup>27</sup>, Callie E. Hood<sup>8</sup>, Alex R. Howe<sup>32</sup>, Paul Kalas<sup>15,16</sup>, Jens Kammerer<sup>27,42</sup>, Grant M. Kennedy<sup>43</sup>, Matthew A. Kenworthy<sup>44</sup>, Pierre Kervella<sup>11</sup>, Minjae Kim<sup>43</sup>, Daniel Kitzmann<sup>20</sup>, Adam L. Kraus<sup>12</sup>, Masayuki Kuzuhara<sup>45</sup>, Pierre-Olivier Lagage<sup>11</sup>, Anne-Marie Lagrange<sup>11</sup>, Kellen Lawson<sup>32</sup>, Cecilia Lazzoni<sup>9</sup>, Jarron M. Leisenring<sup>46</sup>, Ben W. P. Lew<sup>47</sup>, Michael C. Liu<sup>48</sup>, Pengyu Liu<sup>6,7</sup>, Jorge Llop-Sayson<sup>16</sup>, James P. Lloyd<sup>26</sup>, Bruce Macintosh<sup>49</sup>, Mathilde Mâlin<sup>11</sup>, Elena Manjavacas<sup>50</sup>, Sebastián Marino<sup>9</sup>, Mark S. Marley<sup>51</sup>, Christian Marois<sup>52</sup>, Raquel A. Martinez<sup>53</sup>, Elisabeth C. Matthews<sup>41</sup>, Brenda C. Matthews<sup>52</sup>, Dimitri Mawet<sup>16,54</sup>, Johan Mazoyer<sup>11</sup>, Michael W. McElwain<sup>32</sup>, Stanimir Metchev<sup>55</sup>, Michael R. Meyer<sup>34</sup>, Maxwell A. Millar-Blanchaer<sup>56</sup>, Paul Mollière<sup>41</sup>, Sarah E. Moran<sup>51</sup>, Sagnick Mukherjee<sup>8</sup>, Eric Pantin<sup>11</sup>, Marshall D. Perrin<sup>27</sup>, Laurent Pueyo<sup>27</sup>, Sascha P. Quanz<sup>5</sup>, Andreas Quirrenbach<sup>57</sup>, Shrishmoy Ray<sup>9</sup>, Isabel Rebollido<sup>58</sup>, Jea Adams Redai<sup>59</sup>, Bin B. Ren<sup>10</sup>, Emily Rickman<sup>60</sup>, Steph Sallum<sup>53</sup>, Matthias Samland<sup>41</sup>, Benjamin Sargent<sup>27,29</sup>, Joshua E. Schlieder<sup>32</sup>, Karl R. Stapelfeldt<sup>55</sup>, Motohide Tamura<sup>61</sup>, Xianyu Tan<sup>62</sup>, Christopher A. Theissen<sup>63</sup>, Taichi Uyama<sup>64</sup>, Malavika Vasist<sup>28</sup>, Arthur Vigan<sup>35</sup>, Kevin Wagner<sup>46</sup>, Kimberly Ward-Duong<sup>65</sup>, Schuyler G. Wolff<sup>46</sup>, Kadin Worthen<sup>29</sup>, Mark C. Wyatt<sup>66</sup>, Marie Ygouf<sup>55</sup>, Alice Zurlo<sup>2,67,68</sup>, Xi Zhang<sup>8</sup>, Keming Zhang<sup>33</sup>, Zhoujian Zhang<sup>8</sup>, and Yifan Zhou<sup>69</sup>

<sup>1</sup>Núcleo Milenio Formación Planetaria—NPF, Universidad de Valparaíso, Av. Gran Bretaña 1111, Valparaíso, Chile; simon.petrus.pro@gmail.com

<sup>2</sup>Millennium Nucleus on Young Exoplanets and their Moons (YEMS), Chile

<sup>3</sup>Departamento de Física, Universidad de Santiago de Chile, Av. Victor Jara 3659, Santiago, Chile

<sup>4</sup>Department of Astrophysics, American Museum of Natural History, Central Park West at 79th Street, New York, NY 10034, USA

<sup>5</sup>Institute of Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland

<sup>6</sup>Scottish Universities Physics Alliance, Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK

<sup>7</sup>Centre for Exoplanet Science, University of Edinburgh, Edinburgh EH9 3HJ, UK

<sup>8</sup>University of California Santa Cruz, Santa Cruz, CA, USA

<sup>9</sup>University of Exeter, Astrophysics Group, Physics Building, Stocker Road, Exeter, EX4 4QL, UK

<sup>10</sup>Laboratoire Lagrange, Université Cote d'Azur, CNRS, Observatoire de la Cote d'Azur, 06304 Nice, France

<sup>11</sup>LESIA, Observatoire de Paris, Univ PSL, CNRS, Sorbonne Univ, Univ de Paris, 5 place Jules Janssen, 92195 Meudon, France

<sup>12</sup>University of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA

<sup>13</sup>Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

<sup>14</sup>School of Physics, Trinity College Dublin, The University of Dublin, Dublin 2, Ireland

<sup>15</sup>Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA

<sup>16</sup>Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

<sup>17</sup>Naval Research Laboratory, Remote Sensing Division, 4555 Overlook Avenue SW, Washington, DC 20375 USA

<sup>18</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>19</sup>LMU Munich, Faculty of Physics, University Observatory, Scheinerstr. 1, Munich D-81679, Germany

<sup>20</sup>University of Bern, Center for Space and Habitability, Gesellschaftsstrasse 6, CH-3012, Bern, Switzerland

<sup>21</sup>Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042, Graz, Austria

<sup>22</sup>Centre for Exoplanet Science, School of Physics & Astronomy, University of St Andrews, North Haugh, St Andrews, KY169SS, UK

<sup>23</sup>Institute for Theoretical Physics and Computational Physics, Graz University of Technology, Petersgasse 16/II, 8010, Graz, Austria

<sup>24</sup>Department of Astronomy, Stockholm University, AlbaNova University Center, SE-10691 Stockholm, Sweden

<sup>25</sup>Department of Physics and Astronomy, San Francisco State University, 1600 Holloway Avenue, San Francisco, CA 94132, USA

<sup>26</sup>Department of Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive, Ithaca, NY 14853, USA

<sup>27</sup>Space Telescope Science Institute, Baltimore, MD 21218, USA

<sup>28</sup>STAR Institute, Université de Liège, Allée du Six Août 19c, 4000 Liège, Belgium

<sup>29</sup>Department of Physics & Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA

<sup>30</sup>School of Physical Sciences, Faculty of Science, Technology, Engineering and Mathematics, The Open University, Walton Hall, Milton Keynes, MK7 6AA, UK
<sup>31</sup>UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK

<sup>32</sup>NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>33</sup>Department of Astronomy, 501 Campbell Hall, University of California Berkeley, Berkeley, CA 94720-3411, USA

<sup>34</sup>Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48103, USA

<sup>35</sup>Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>36</sup>Department of Physics and Astronomy, University of Texas-San Antonio, 1 UTSA Circle, San Antonio, TX, USA
37 Subaru Telescope, National Astronomical Observatory of Japan, 650 North A'ohōkū Place, Hilo, HI 96720, USA

38 Instituto de Astrofísica de Andalucía, CSIC, Glorieta de la Astronomía s/n, 18008, Granada, Spain

39 University of California, Los Angeles, 430 Portola Plaza, Box 951547, Los Angeles, CA 90095-1547, USA

40 Eureka Scientific, 2452 Delmer Street, Suite 1, Oakland, CA 96402, USA

41 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

42 European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748, Garching, Germany

43 Department of Physics, University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK

44 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

45 Astrobiology Center of NINS, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan

46 Steward Observatory and the Department of Astronomy, The University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

47 Bay Area Environmental Research Institute and NASA Ames Research Center, Moffett Field, CA 94035, USA

48 Institute for Astronomy, University of Hawai'i, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

49 Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA

50 AURA for the European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

51 Department of Planetary Sciences, Lunar & Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA

52 Herzberg Astronomy & Astrophysics Research Centre, National Research Council of Canada, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada

53 Department of Physics and Astronomy, 4129 Frederick Reines Hall, University of California, Irvine, CA 92697, USA

54 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

55 Western University, Department of Physics & Astronomy and Institute for Earth and Space Exploration, 1151 Richmond Street, London, ON N6A 3K7, Canada

56 Department of Physics, University of California, Santa Barbara, CA 93106, USA

57 Landesernwarte, Zentrum für Astronomie der Universität Heidelberg, Königstuhl 12, 69117 Heidelberg, Germany

58 Centro de Astrobiología (CAB CSIC-INTA) ESAC Campus Camino Bajo del Castillo, s/n, Villanueva de la Cañada, 28692, Madrid, Spain

59 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

60 European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

61 The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

62 Tsung-Dao Lee Institute, Shanghai Jiao Tong University, 520 Shengrong Road, Shanghai, People's Republic of China

63 Department of Astronomy & Astrophysics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0424, USA

64 Infrared Processing and Analysis Center, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA

65 Department of Astronomy, Smith College, Northampton, MA 01063, USA

66 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

67 Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile

68 Escuela de Ingeniería Industrial, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile

69 University of Virginia, Department of Astronomy, 530 McCormick Road, Charlottesville, VA 22904, USA

Received 2023 December 7; revised 2024 January 15; accepted 2024 January 23; published 2024 April 25

## Abstract

The unprecedented medium-resolution ($R_\lambda \sim 1500–3500$) near- and mid-infrared (1–18 µm) spectrum provided by JWST for the young (140 ± 20 Myr) low-mass (12–20 $M_{Jup}$) L–T transition (L7) companion VHS 1256 b gives access to a catalog of molecular absorptions. In this study, we present a comprehensive analysis of this data set utilizing a forward-modeling approach applying our Bayesian framework, ForMoSA. We explore five distinct atmospheric models to assess their performance in estimating key atmospheric parameters: $T_{eff}$, log(g), [M/H], C/O, $\gamma$, $f_{sed}$, and R. Our findings reveal that each parameter’s estimate is significantly influenced by factors such as the wavelength range considered and the model chosen for the fit. This is attributed to systematic errors in the models and their challenges in accurately replicating the complex atmospheric structure of VHS 1256 b, notably the complexity of its clouds and dust distribution. To propagate the impact of these systematic uncertainties on our atmospheric property estimates, we introduce innovative fitting methodologies based on independent fits performed on different spectral windows. We finally derived a $T_{eff}$ consistent with the spectral type of the target, considering its young age, which is confirmed by our estimate of log(g). Despite the exceptional data quality, attaining robust estimates for chemical abundances [M/H] and C/O, often employed as indicators of formation history, remains challenging. Nevertheless, the pioneering case of JWST’s data for VHS 1256 b has paved the way for future acquisitions of substellar spectra that will be systematically analyzed to directly compare the properties of these objects and correct the systematics in the models.

Unified Astronomy Thesaurus concepts: Exoplanet atmospheres (487); Infrared spectroscopy (2285); Direct imaging (387); L dwarfs (894); Nested sampling (1894); James Webb Space Telescope (2291); Astronomy data modeling (1859)

## 1. Introduction

The diversity of exoplanets that have been discovered to date has reshaped our understanding of planetary systems and has raised new questions regarding their formation and evolution.

Although different formation scenarios, such as core accretion (Pollack et al. 1996), gravitational instability (Boss 2000), and gravoturbulent fragmentation of molecular clouds (Padoan et al. 2005), have been proposed, constraining the evolutionary history of observed planets remains challenging. Indeed, the initial physical and chemical conditions of the circumstellar disk in which they formed are no longer directly observable, and the dynamical evolution of the system (migration, ejection, or planetary capture) is unknown. However, formation models have identified key parameters in the atmosphere of young planets that can potentially serve as formation tracers, such as

<figure><img src="image_3.png" /></figure>

<footnote n="70">70 51 Pegasi b Fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
the metallicity [M/H] (Ormel et al. 2021), the carbon-to-oxygen ratio C/O (Öberg et al. 2011), and, more recently, the isotopic ratio of ¹²CO/¹³CO (Zhang et al. 2021). By detecting specific molecules in the atmospheres of transiting planets, it has become possible to estimate these observables (Hoeij-makers et al. 2019; Spake et al. 2021). However, this technique is limited to highly irradiated planets close to their host stars. Direct imaging allows for the characterization of the atmospheres of planets with larger separations, but due to the high contrast with their host star, only sufficiently bright companions can be accessed. Presently, there are fewer than 30 low-mass companions for which spectroscopic data have been obtained. Only a few of these observations have reached sufficient spectral resolution ($R_\lambda > 1000$) to attempt estimating chemical abundances in their atmospheres, including β Pic b (Snellen et al. 2014), HD 106906 b (Daemgen et al. 2017), 1RXS J1609 b (Lafrenière et al. 2010), Delorme 1(AB) b (Eriksson et al. 2020; Betti et al. 2022; Ringqvist et al. 2023), HR 8799 b, c (Konopacky et al. 2013; Ruffio et al. 2019; Wang et al. 2021, 2023), $\kappa$ And b (Wilcomb et al. 2020), TYC 8998 b (Zhang et al. 2021), HIP 65426 b (Petrus et al. 2021), AB Pic b (Palma-Bifani et al. 2023), and, lastly, VHS 1256 b (Hoch et al. 2022; Petrus et al. 2023), which is the subject of this work.

The system VHS J125601.92–125723.9 (hereafter VHS 1256) has been identified as a hierarchical system with the detection of the companion VHS 1256 b, orbiting at ~8''06 ± 0''03 (projected physical separation of 179 ± 9 au) from the tight (0.1 as) M7.5+M7.5 binary VHS 1256 AB (Gauza et al. 2015; Stone et al. 2016). Its distance and age have been estimated to 21.14 ± 0.22 pc (Brown et al. 2021) and 140 ± 20 Myr (Dupuy et al. 2023), respectively.

Bowler et al. (2020) and Zhou et al. (2020) monitored the companion with the Hubble Space Telescope and Spitzer, respectively, and revealed significant wavelength-dependent variability, corresponding to an estimated rotation period ranging between 21 and 24 hr. Zhou et al. (2022) confirmed this hourly period variability and further identified a higher temporal baseline variability with multiepoch tracking over a 2 yr period. Both publications emphasized VHS 1256 b as one of the most variable substellar objects discovered thus far (over 20% between 1.1 and 1.7 µm) and interpreted this variability as the signature of a complex and heterogeneous atmospheric structure (inhomogeneous cloud cover, heterogeneous temperature, etc.). These conclusions align with the results provided by Gauza et al. (2015), who identified VHS 1256 b as an L7-type object at the boundary of the L–T transition (between 1000 and 1400 K). Through this temperature range, planetary atmospheres are expected to undergo significant modifications to their cloud structure in terms of vertical distribution (condensation and/or sedimentation of FeH, TiO, and VO; Cushing et al. 2008) and/or horizontal distribution (appearance of holes in the cloud cover; Burgasser et al. 2002; Marley et al. 2010). Furthermore, by employing the methodology detailed in Allers & Liu (2013), Gauza et al. (2015) noted indicators of low surface gravity within the spectrum of this object, consistent with its young age. Notably, they observed a characteristic triangular H-band shape. This low gravity implied a reddening effect (see Figure 1 in Miles et al. 2023) that provides additional evidence supporting the presence of cloudy layers.

Several studies have analyzed medium- and high-resolution spectra of VHS 1256 b to better understand this complexity. Using Keck Near-Infrared Spectrograph (NIRSpec) data in the H band ($R_\lambda \sim 25,000$), Bryan et al. (2018) estimated a projected rotational velocity $v \sin(i)$ of $13.5 \pm 4.1$ km s⁻¹, which was consistent with the rotational period estimated by Zhou et al. (2020). This $v \sin(i)$ measurement also enabled Zhou et al. (2020) to estimate that VHS 1256 b is very likely to be observed equatorially. Keck/NIRSpec L-band spectroscopy ($\Delta\lambda = 2.9–4.4$ µm; $R_\lambda \sim 1300$) enabled the detection of methane, which revealed nonequilibrium chemistry between CO and CH₄ (Miles et al. 2018). Hoch et al. (2022) compared medium-resolution Keck/OSIRIS spectroscopy in the K band ($R_\lambda \sim 4000$) with a custom grid of precomputed synthetic spectra to derive a $T_{\text{eff}}$ range of 1200–1300 K, a log(g) range of 3.25–3.75 dex, and a C/O of 0.590 ± 0.354. These atmospheric parameter estimates were confirmed by Petrus et al. (2023), who fit a medium-resolution spectrum from the X-Shooter ($R_\lambda \sim 8000$) spanning 1.10–2.48 µm simultaneously with the ATMO grid (Tremblin et al. 2015) to estimate $T_{\text{eff}} = 1380 \pm 54$ K, log(g) = 3.97 ± 0.48 dex, [M/H] = 0.21 ± 0.29, and C/O > 0.63. They concluded that while the estimates of [M/H] and C/O suggested supersolar values, indicating significant enrichment in solids during its formation, the hypothesis of a formation by fragmentation of the initial molecular clouds could not be dismissed due to the characteristics of the system’s architecture (hierarchical system) and the lack of sufficient precision in estimates of these formation tracers.

More recently, VHS 1256 b has been observed with JWST, which combined medium-resolution spectroscopy from NIR-Spec and the Mid-Infrared Instrument (MIRI) (Miles et al. 2023; see Section 2.1). The spectra were integrated to obtain a robust estimate of the bolometric luminosity log($L_{bol}/L_\odot$) = −4.550 ± 0.009. Dupuy et al. (2023) combined this luminosity with their estimate of the age to derive the mass of VHS 1256 b using the hybrid cloudy evolutionary model of Saumon & Marley (2008). They identified a bimodal solution with a lower-mass regime at 12 ± 0.1 $M_{\text{Jup}}$ and a higher-mass regime at 16 ± 1 $M_{\text{Jup}}$. Given this bimodality, we can confidently place the mass between 12 and 20 $M_{\text{Jup}}$, but, as we do not know whether this object is currently burning deuterium, we do not know which of the bimodal peaks it truly falls in. Dupuy et al. (2023) also estimated a radius R of 1.30 $R_{\text{Jup}}$ and 1.22 $R_{\text{Jup}}$, a $T_{\text{eff}}$ of 1153 ± 5 K and 1194 ± 9 K, and a log(g) of 4.268 ± 0.006 dex and 4.45 ± 0.03 dex according to the lower-mass and higher-mass solution, respectively. Miles et al. (2023) also compared this spectrum with a grid of synthetic spectra from Mukherjee et al. (2023) to estimate a $T_{\text{eff}}$ of 1100 K, a log(g) of 4.5 dex, and a radius of 1.27 $R_{\text{Jup}}$ and to highlight a complex atmosphere with a different cloud deck and varying regimes of vertical mixing. Additionally, the MIRI data enabled the detection of a strong silicate absorption at around 10 µm, supporting the presence of silicate clouds in its atmosphere.

These previous studies indicate the usefulness of employing grids of precomputed synthetic spectra for characterizing the atmosphere of VHS 1256 b. These grids are generated from atmospheric models that incorporate our current understanding of the physics and chemistry involved. However, it has been observed that these models sometimes encounter challenges in accurately reproducing certain spectral features and extended
wavelength ranges. These systematic errors can be attributed to
the use of strong assumptions when constraining the physical
and chemical processes occurring in planetary atmospheres,
such as cloud formation and evolution, chemical disequili-
brium, vertical mixing, and dust sedimentation (see
Section 2.2). The imposition of these assumptions leads to a
limited number of free parameters (<6), and the quality of the
fit heavily relies on the quality of the models. The impacts of
these systematic errors on the estimation of atmospheric
properties have already been recognized for L and M spectral
types (Cushing et al. 2008; Stephens et al. 2009; Bonnefoy
et al. 2014; Manjavacas et al. 2014; Lachapelle et al. 2015;
Bayo et al. 2017; Petrus et al. 2020, 2023; Suárez et al. 2021;
Lueber et al. 2023; Palma-Bifani et al. 2023). These studies
have interpreted these discrepancies as indications of a
deficiency in the modeling of dust and haze within the
synthetic atmosphere models.

With the wavelength coverage and data quality of the spectra
now provided by JWST, it is imperative to establish a
comprehensive understanding of the limitations of self-
consistent models of atmospheres and to develop analysis
methods that take appropriate account of their systematic
errors. This is crucial to ensure the production of reliable and
robust estimates of the atmospheric properties of planetary-
mass objects and to avoid overinterpretation of these properties,
especially when dealing with spectra covering a narrower
spectral interval. In this paper, we exploit the unprecedented
data provided by JWST to propose an original method that
aims to identify and propagate the systematic errors of the
models directly into the error bars of the atmospheric
parameters estimated from the forward-modeling approach.
Five different grids of precomputed synthetic spectra are used,
generated from five different self-consistent models, to
compare the results they yield and the interpretation that can
be made to constrain the formation pathway of VHS 1256 b.
Section 2 describes the data set and models used in this study.
Section 3 exploits the broad wavelength coverage provided by
JWST to study the dispersion of parameter estimates along the
spectral energy distribution (SED). Section 4 focuses on
medium-resolution spectral features to estimate chemical
abundances. Finally, Section 5 is devoted to a discussion of
the results, and a conclusion is given in Section 6.

## 2. Description of the Data and the Models

2.1. Data

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

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Spectral resolution as a function of the wavelength coverage allowed by the NIRSpec + MIRI data. The observed data are represented by the black lines. The five models are depicted by the colored lines: ATMO in green, Exo-REM in blue, Sonora in yellow, BT-Settl in red, and DRIFT-PHOENIX in purple. For the models, the spectral resolution has been calculated assuming a Nyquist sampling rate.</figcaption></figure>

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).

## 2.2. Models

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
<table><caption>**Table 1**
Parameter Spaces Explored by the Models</caption><thead><tr><th rowspan="3"></th><th colspan="2">T<sub>eff</sub></th><th colspan="2">log(g)</th><th colspan="2">[M/H]</th><th colspan="2">C/O</th><th colspan="2">γ</th><th colspan="2">f<sub>sed</sub></th></tr><tr><th colspan="2">(K)</th><th colspan="2">(dex)</th><th colspan="2">(M/H)</th><th colspan="2">(M/H)</th><th colspan="2">(M/H)</th><th colspan="2">(M/H)</th></tr><tr><th>Range</th><th>Step</th><th>Range</th><th>Step</th><th>Range</th><th>Step</th><th>Range</th><th>Step</th><th>Range</th><th>Step</th><th>Range</th><th>Step</th></tr></thead><tbody><tr><td>ATMO</td><td>[800, 3000]</td><td>100</td><td>[2.5, 5.5]</td><td>0.5</td><td>[-0.6, 0.6]</td><td>0.3</td><td>[0.3, 0.7]</td><td>0.25</td><td>[1.01, 1.05]</td><td>0.02</td><td>...</td><td>...</td></tr><tr><td>Exo-REM</td><td>[400, 2000]</td><td>100</td><td>[3.0, 5.0]</td><td>0.5</td><td>[-0.5, 0.5]</td><td>0.5</td><td>[0.1, 0.8]</td><td>0.1</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Sonora</td><td>[900, 2400]</td><td>100</td><td>[3.5, 5.5]</td><td>0.5</td><td>[-0.5, 0.5]</td><td>0.5</td><td>...</td><td>...</td><td>...</td><td>...</td><td>[1.0, 8.0]</td><td>1.0</td></tr><tr><td>BT-Settl</td><td>[1200, 3000]</td><td>100</td><td>[2.5, 5.5]</td><td>0.5</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>DRIFT-PHOENIX</td><td>[1000, 3000]</td><td>100</td><td>[3.0, 6.0]</td><td>0.5</td><td>[-0.6, 0.3]</td><td>0.3</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** Effective temperature ($T_{\text{eff}}$), surface gravity (log(g)), metallicity ([M/H]), carbon-oxygen ratio (C/O), adiabatic index ($\gamma$), and sedimentation factor ($f_{\text{sed}}$).

adiabatic index $\gamma$ 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$_2$–H$_2$ and H$_2$–He, 13 molecules (H$_2$O, CO$_2$, CO, CH$_4$, NH$_3$, TiO, VO, FeH, PH$_3$, H$_2$S, HCN, C$_2$H$_2$, and SO$_2$), seven atoms (Na, K, Li, Rb, Cs, Fe, and H$^-$), and the Rayleigh scattering opacities for H$_2$, He, CO, N$_2$, CH$_4$, NH$_3$, H$_2$O, CO$_2$, H$_2$S, and SO$_2$. This model explores five parameters: $T_{\text{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 $\gamma$ 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$_4$, CO$_2$, and NH$_3$), and the cloud model incorporates the formation of iron, Na$_2$S, KCl, silicates, and water. Opacity sources include collision-induced absorptions of H$_2$–H$_2$ and H$_2$–He, rovibrational bands from nine molecules (H$_2$O, CH$_4$, CO, CO$_2$, NH$_3$, PH$_3$, 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_{\text{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_{\text{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$_4$, CO$_2$, N$_2$, and NH$_3$. 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_{\text{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$_2$, Al$_2$O$_3$, Fe, SiO$_2$, MgO, MgSiO$_3$, and Mg$_2$SiO$_4$) 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.

## 2.3. Parameter Space Exploration with ForMoSA

To compare these models with the data, we used the forward-modeling code ForMoSA.<fnref n="71" /> 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`<fnref n="72" /> 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`<fnref n="73" /> 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`.<fnref n="74" /> 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 $\mu = 21.15$ pc and $\sigma = 0.22$ pc, corresponding to the Gaia DR3 measurement and error, respectively.

## 3. Fits as a Function of the Wavelength Range

The forward-modeling approach used in this study offers several advantages, including the ability to fit high-resolution data across a broad range of wavelengths within a relatively short computing time (less than 1 day). However, this approach has its challenges, as the quality of the fit is directly influenced

by the quality of the model, which can be difficult to fully quantify. Petrus et al. (2020) demonstrated that systematic errors present in the model BT-Settl greatly influenced the posteriors of the parameter estimates explored during the inversion of nine spectra of young M7–M9 dwarfs. The fits converged to different sets of parameters depending on the spectral band considered. These differences considerably dominated the errors on these parameters derived from the Bayesian algorithm. This work was further extended to the analysis of isolated brown dwarfs at the L–T transition observed with Spitzer (Suárez et al. 2021), L dwarfs using X-Shooter data of VHS 1256 b with the ATMO grid (Petrus et al. 2023), and Very Large Telescope/SINFONI data of AB Pic b with the BT-Settl and Exo-REM grids (Palma-Bifani et al. 2023). These subsequent studies confirmed that the systematic errors observed in the previous work were indeed a recurring issue in the forward-modeling approach. In this section, we present an extension of the method proposed in these three previous studies. We address the systematics in the models to provide a robust estimation of the atmospheric properties of VHS 1256 b, taking advantage of the extensive spectral coverage offered by JWST.

### 3.1. Description of the Method

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 NIRSpect 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.

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 NIRSpect 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:

$$C_{win} \propto \sum_{\lambda} S/N_{\lambda}^2 \quad (1)$$

where $S/N_{\lambda}$ is the S/N associated with each wavelength, assuming completely uncorrelated errors. The sum of $S/N^2$ is considered to align with the definition of the likelihood

<footnote n="71">⁷¹ https://formosa.readthedocs.io</footnote>

<footnote n="72">⁷² https://docs.xarray.dev</footnote>

<footnote n="73">⁷³ https://pypi.org/project/spectres/</footnote>

<footnote n="74">⁷⁴ http://kylebarbary.com/nestle/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison between the data and the interpolated synthetic spectra generated with the set of parameters that maximized the likelihood for each model. The best fit using the full SED is depicted with unicolor lines, while the best fits using different spectral windows are illustrated with bicolor lines. The data are represented by black lines. All spectra have been normalized with the bolometric luminosity of VHS 1256 b to simplify the figure, and an offset was applied. The use of spectral windows enhances the quality of the fit for each model. We note that a possible underestimation of the error bars of the data may play a role in the high value of $\chi_r^2$, which is provided as a comparative measure of the performance of the best-fit models, as well as the number of degrees of freedom (DoF) that have been used to calculate it. The top right panel shows the relative contribution of each spectral window to the fit using the full SED. It corresponds to $C_{win}/C_{win,ot}$ when $C_{win}$ is calculated from Equation (1) and $C_{win,ot}$ is the sum of $C_{win}$ over the 15 spectral windows. The NIRSpec channel G395HF/290LP, from 2.87 to 5.27 µm, with its higher spectral resolution and S/N, carries the greatest weight in the calculation of the likelihood. See Appendix A for a zoomed version of this plot.</figcaption></figure>

function utilized during the fits, which is directly proportional to a $\chi^2$ 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 $\Delta\lambda$, model
<table><caption>**Table 2**
Spectral Windows Defined for the NIRSpect Data</caption><thead><tr><th>NIRSpec</th><th></th><th>Window 1<br>(μm)</th><th>Window 2<br>(μm)</th></tr></thead><tbody><tr><td rowspan="3">Channels</td><td>G140HF/100LP</td><td>0.97–1.40</td><td>1.40–1.89</td></tr><tr><td>G235HF/170LP</td><td>1.89–2.45</td><td>2.45–3.17</td></tr><tr><td>G395HF/290LP</td><td>2.87–4.10</td><td>4.10–5.27</td></tr></tbody></table>

Note. Each NIRSpec channel has been split into two windows, corresponding to different NIR bands.

<table><caption>**Table 3**
Spectral Windows Defined for MIRI Data</caption><thead><tr><th rowspan="2">MIRI</th><th rowspan="2">Channels</th><th colspan="3">Grating Settings</th></tr><tr><th>A<br>(μm)</th><th>B<br>(μm)</th><th>C<br>(μm)</th></tr></thead><tbody><tr><td>1</td><td>4.90–5.74</td><td>5.66–6.63</td><td>6.53–7.65</td></tr><tr><td>2</td><td>7.51–8.77</td><td>8.67–10.13</td><td>10.02–11.70</td></tr><tr><td>3</td><td>11.55–13.47</td><td>13.34–15.57</td><td>15.41–17.98</td></tr></tbody></table>

Note. The nine windows correspond to the three channels used to acquire the data combined with the three grating settings.

m, and parameter p, we calculated this sensitivity $S_{\Delta\lambda,m,p}$ using the following relation:

$$S_{\Delta\lambda,m,p} = \sum_{\lambda} \left( \frac{\text{abs}(F_{m,\lambda}(p_{\text{max}}) - F_{m,\lambda}(p_{\text{min}}))}{F_{m,\lambda}(p_{\text{max}})} \right) \times \frac{1}{N_{\Delta\lambda}}, \quad (2)$$

where $F_{m,\lambda}(p_{\text{min}})$ and $F_{m,\lambda}(p_{\text{max}})$ are the synthetic flux generated by the model m at the wavelength $\lambda$ with the minimum and maximum values of the parameter p explored by m, respectively, and $N_{\Delta\lambda}$ is the number of data points contained in $\Delta\lambda$. 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,\lambda}(p)$ within the considered window is not affected significantly by the parameter p, the difference $F_{m,\lambda}(p_{\text{max}}) - F_{m,\lambda}(p_{\text{min}})$ (and consequently $S_{\Delta\lambda,m,p}$) will tend toward 0. Conversely, a substantial $S_{\Delta\lambda,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.

## 3.2. Results: Parameters as a Function of Wavelengths

All five models exhibit deviations from the observed data when the entire wavelength range is considered for the fit, particularly at longer wavelengths ($\lambda > 7$ µm), which carry less weight in the likelihood calculation. To assess the goodness of fit, we calculated the reduced $\chi^2$ values for the fits using the full SED, denoted as $\chi^2_{r,\text{full}}$, as well as for the fits using the spectral windows, denoted as $\chi^2_{r,\text{win}}$ (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^2_{r,\text{full}} = 246$ and $\chi^2_{r,\text{win}} = 57$. The utilization of smaller wavelength ranges in the independent fits (spectral windows) results in a decrease in the overall $\chi^2_{r,\text{win}}$ 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 $\Theta_{\text{win,f}}$, where $\Theta$ 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_{\Delta\lambda,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_{\text{eff}}$

Considering the young age of VHS 1256 b, the estimate of $T_{\text{eff,full}}$ considering the full wavelength range is consistent with its spectral type and with the prediction of the evolutionary models (Dupuy et al. 2023) for the grids Sonora, Exo-REM, and ATMO (Figure 3). The BT-Settl and DRIFT-PHOENIX grids appear to converge toward higher $T_{\text{eff}}$ values than those predicted by evolutionary models. However, it is noteworthy that for these two models, ForMoSA converges to $T_{\text{eff}}$ values close to their lower limit. The ATMO and Exo-REM grids demonstrate greater stability across the different spectral windows compared to BT-Settl, Sonora, and DRIFT-PHOENIX. This leads to smaller error bars on the final values $T_{\text{eff,win,f}}$ which remains consistent with $T_{\text{eff,full}}$ for each model. The large error bar at $\lambda > 13$ µm is attributed to the low S/N at these wavelengths. The sensitivity of $T_{\text{eff}}$ remains relatively constant across all spectral windows.

### 3.2.2. The Surface Gravity $\log(g)$

The results using the full wavelength range for the models Exo-REM, ATMO, Sonora, and BT-Settl converge toward low values of $\log(g)$ (<4.5 dex), whereas DRIFT-PHOENIX tends to converge toward high values (>4.5 dex; Figure 3). Significant dispersion was observed across the various spectral windows, covering the entire range of $\log(g)$ values explored,
<table><caption>**Table 4**
Summary of the Estimations of the Parameters Explored by the Grids</caption><thead><tr><th></th><th>Fit</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>log(g)<br>(dex)</th><th>(M/H)</th><th>C/O</th><th>γ</th><th><i>f</i><sub>sed</sub></th></tr></thead><tbody><tr><td rowspan="2">ATMO</td><td>Full</td><td>1203 ± 1</td><td>2.54 ± 0.01</td><td>-0.01 ± 0.01</td><td>&lt;0.30</td><td>1.030 ± 0.001</td><td>...</td></tr><tr><td>Windows</td><td>1318 ± 166</td><td>3.59 ± 0.53</td><td>-0.32 ± 0.28</td><td>0.46 ± 0.16</td><td>1.035 ± 0.011</td><td>...</td></tr><tr><td rowspan="2">Exo-REM</td><td>Full</td><td>1164 ± 1</td><td>&lt;3</td><td>-0.19 ± 0.01</td><td>0.50 ± 0.01</td><td>...</td><td>...</td></tr><tr><td>Windows</td><td>1153 ± 152</td><td>3.70 ± 0.34</td><td>0.03 ± 0.17</td><td>0.38 ± 0.20</td><td>...</td><td>...</td></tr><tr><td rowspan="2">Sonora</td><td>Full</td><td>1116 ± 1</td><td>&lt;3.50</td><td>&gt;0.50</td><td>...</td><td>...</td><td>1.01 ± 0.01</td></tr><tr><td>Windows</td><td>1349 ± 208</td><td>4.50 ± 0.63</td><td>0.41 ± 0.14</td><td>...</td><td>...</td><td>2.97 ± 2.09</td></tr><tr><td rowspan="2">BT-Settl</td><td>Full</td><td>1482 ± 1</td><td>3.50 ± 0.01</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Windows</td><td>1560 ± 203</td><td>3.62 ± 0.60</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">DRIFT-PHOENIX</td><td>Full</td><td>1451 ± 1</td><td>4.89 ± 0.01</td><td>-0.28 ± 0.01</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Windows</td><td>1540 ± 286</td><td>4.54 ± 1.13</td><td>-0.36 ± 0.26</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** Estimations of the parameters explored by the grids are provided for both the full wavelength range and the combined spectral windows. For the combined spectral windows, error bars have been calculated by considering the dispersion of the individual fits from each spectral window. The error bars of the fits on the full SED reflect the high S/N of the data.

<table><caption>**Table 5**
Summary of the Estimations of the Parameters Not Explored by the Grids</caption><thead><tr><th></th><th>Fit</th><th><i>R</i><br>(<i>R</i><sub>Jup</sub>)</th><th>log(<i>L</i><sub>bol</sub>/<i>L</i><sub>☉</sub>)</th><th><i>M</i><br>(<i>M</i><sub>Jup</sub>)</th></tr></thead><tbody><tr><td rowspan="2">ATMO</td><td>Full</td><td>1.19 ± 0.01</td><td>-4.57 ± 0.01</td><td>0.19 ± 0.01</td></tr><tr><td>Windows</td><td>1.09 ± 0.23</td><td>-4.56 ± 0.20</td><td>2.24 ± 1.96</td></tr><tr><td rowspan="2">Exo-REM</td><td>Full</td><td>1.29 ± 0.01</td><td>-4.56 ± 0.01</td><td>0.64 ± 0.01</td></tr><tr><td>Windows</td><td>1.39 ± 0.28</td><td>-4.55 ± 0.14</td><td>7.87 ± 6.33</td></tr><tr><td rowspan="2">Sonora</td><td>Full</td><td>1.34 ± 0.01</td><td>-4.59 ± 0.01</td><td>2.17 ± 0.08</td></tr><tr><td>Windows</td><td>1.12 ± 0.16</td><td>-4.48 ± 0.25</td><td>&lt;94.40</td></tr><tr><td rowspan="2">BT-Settl</td><td>Full</td><td>0.81 ± 0.01</td><td>-4.54 ± 0.01</td><td>0.80 ± 0.01</td></tr><tr><td>Windows</td><td>0.95 ± 0.16</td><td>-4.37 ± 0.29</td><td>&lt;20.90</td></tr><tr><td rowspan="2">DRIFT-PHOENIX</td><td>Full</td><td>0.79 ± 0.01</td><td>-4.60 ± 0.01</td><td>18.73 ± 0.69</td></tr><tr><td>Windows</td><td>0.93 ± 0.15</td><td>-4.43 ± 0.31</td><td>&lt;249.42</td></tr></tbody></table>

**Note.** Estimations of the parameters not explored by the grids are provided for both the full wavelength range and the combined spectral windows. The radius *R* is determined from the dilution factor *C*<sub>*K*</sub>, the luminosity *L* is calculated using the Stefan-Boltzmann law and the estimates of the radius and the *T*<sub>eff</sub>, and the mass *M* is calculated using the gravitational law and the estimates of the radius and the log(g).

particularly for the BT-Settl, Sonora, and DRIFT-PHOENIX grids. Smaller dispersion was also evident in the other grids, Exo-REM and ATMO. However, trends could be identified, with a preference for lower surface gravities along the SED. The grid Sonora provides low log(g) for wavelengths shorter than 4 µm but higher log(g) for longer wavelengths. The sensitivity analysis confirmed that the log(g) parameter is primarily constrained by the NIR data (λ < 5 µm). Due to the wide dispersion observed across the different spectral windows, it is challenging to provide a robust value for log(g) except for Exo-REM and ATMO, which converge to consistent values. Although these log(g)<sub>win,f</sub> point to low values (<4.12 dex), typical of young objects, they are significantly lower than the predictions of evolutionary models provided by Dupuy et al. (2023).

### 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]<sub>win,f</sub>. Similar to the log(g) parameter, the [M/H] estimation is primarily constrained by the NIR data (λ < 5 µm).

### 3.2.4. The Carbon-Oxygen Ratio C/O

Among the considered grids, only Exo-REM and ATMO enable the exploration of the C/O ratio (Figure 4). The Exo-REM model appears to exhibit a preference for low C/O ratios, although the wide dispersion of estimates makes it difficult to establish a definitive trend. The situation is even more complex with the ATMO model, as the estimated C/O ratios span the entire range of values explored by the grid. Consequently, this parameter cannot be reliably estimated using this method. Nonetheless, it appears that the NIR data exhibit the highest sensitivity to this parameter

### 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 γ<sub>win,f</sub> 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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Estimation of the $T_{\text{eff}}$ (left panel) and $\log(g)$ (right panel) using the full SED during the fit (dashed colored line) and using the spectral windows (colored dots). The size of the dots indicates the sensitivity of each spectral window to each parameter. If the sensitivity is too low to be represented by a dot, a colored cross is used. Each panel corresponds to a different model, with the parameter range indicated by the filled colored area. The final value and error extracted from these results are shown by the solid gray line and gray area, respectively. The dotted black lines represent the value estimated by Dupuy et al. (2023), who interpolated the evolutionary model from Saumon & Marley (2008) assuming the bolometric luminosity and age of VHS 1256 b. The two solutions they have found are given.</figcaption></figure>

aligns coherently with the results observed for late L objects (Stephens et al. 2009; Tremblin et al. 2017).

### 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)^2$, 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).

## 4. Focusing on Spectral Features

The medium spectral resolution provided by JWST has enabled the redetection of various atomic and molecular absorptions in the spectrum of VHS 1256 b. Notably, the two potassium (K I) doublets at 1.173 and 1.248 µm are present, previously identified as robust indicators of surface gravity and metallicity through empirical and synthetic analysis of spectra from planetary-mass objects (Allers & Liu 2013; Petrus et al. 2023). Additionally, we identify carbon monoxide (CO) overtones at approximately 2.3 µm, which can be used to measure the C/O ratio (Konopacky et al. 2013; Nowak et al. 2020; Petrus et al. 2021; Hoch et al. 2022; Petrus et al. 2023). The methane (CH₄) absorption between 2.8 and 3.8 µm is also detected, confirming the lower-resolution detection of these species in the atmosphere of VHS 1256 b by Miles et al. (2018), who suggested nonequilibrium chemistry between CO and CH₄ and a significant vertical mixing ($K_{zz} > 10^8$ cm² s⁻¹) in the atmosphere to explain the presence of this feature. Finally, this JWST spectrum enabled the detection of a highly resolved CO absorption forest centered at 4.6 µm, already detected in the low-resolution spectra of isolated brown dwarfs by Sorahana & Yamamura (2012). In this section, we conduct a comprehensive analysis of these prominent spectral features to
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Figure 3 but for the [M/H] (left) and the C/O (right). The dotted black lines represent the solar values 0.0 and 0.55 for [M/H] and C/O, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Same as Figure 3 but for the adiabatic index $\gamma$ (left) and the sedimentation factor $f_{sed}$ (right).</figcaption></figure>

determine the extent to which they can be utilized to constrain atmospheric properties through the forward-modeling approach. We assess the performance of each model described in Section 2.2.

## 4.1. *Description of the Method*

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, $\Delta\lambda_{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_{\lambda} < 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, $\Delta\lambda_{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.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Figure 3 but for the radius *R*. This parameter is estimated from the dilution factor $C_K$ considering the distance from Brown et al. (2021).</figcaption></figure>

<table><caption>**Table 6**
Spectral Windows Were Employed to Fit the Four Spectral Features: K I lines, CO Overtones, CH$_4$ Absorption, and CO Forest</caption><thead><tr><th>Chemical Element</th><th>Δλ<sub>continuum</sub><br>(μm)</th><th>Δλ<sub>fit</sub><br>(μm)</th></tr></thead><tbody><tr><td>K I lines</td><td>[1.160, 1.166] + [1.180, 1.241]<br>+ [1.255, 1.260]</td><td>[1.166, 1.180]<br>+ [1.241, 1.255]</td></tr><tr><td>CO overtones</td><td>[2.250, 2.290] + [2.305, 2.319]<br>+ [2.330, 2.400]</td><td>[2.290, 2.305]<br>+ [2.319, 2.330]</td></tr><tr><td>CH<sub>4</sub> absorption</td><td>[3.100, 3.300] + [3.375, 3.600]</td><td>[3.300, 3.375]</td></tr><tr><td>CO forest</td><td>[4.300, 5.000]</td><td>[4.300, 5.000]</td></tr></tbody></table>

**Note.** The Δλ<sub>continuum</sub> is used to estimate the T<sub>eff</sub> and the radius. These derived values serve as constraints within the Δλ<sub>fit</sub> for fitting these features.

## 4.2. Results

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 χ², calculated using Δλ<sub>fit</sub>, 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.

### 4.2.1. ATMO

The T<sub>eff</sub> derived from the various spectral windows excluding the spectral features is in agreement with the estimate obtained in Section 3, ranging between 1242 and 1520 K. This consistency is expected, as we employed a similar methodology by utilizing a fit based on Δλ<sub>continuum</sub>, which encompasses a broad range of wavelengths. The low values of log(g) (less than 4.37 dex) also align well with the age of the system. However, it is worth noting that a remarkably low and unphysical estimate (2.56 dex) is obtained when considering methane absorption. Subsolar metallicities are estimated for the companion, except for the potassium absorption lines, which yielded a supersolar value. The reported C/O ratios range from solar to supersolar, except for the CH<sub>4</sub> feature that provides a subsolar value. Furthermore, the adiabatic index (γ) shows variability across the entire grid size in the different fits, making it challenging to identify a clear trend. Among the models used, ATMO demonstrates the ability to reproduce this information at medium resolution with the lowest χ². However, it seems to inaccurately reproduce the CO overtone at 2.30 μm.

### 4.2.2. Exo-REM

Similarly to ATMO, Exo-REM allowed us to estimate consistent T<sub>eff</sub> values, ranging from 1233 to 1460 K, and low log(g) values, lower than 4.17 dex. We also encountered discrepancies with the very low T<sub>eff</sub> values provided by methane absorption and the very low log(g) values provided by both methane absorption and CO overtones. The estimated metallicity is mostly solar, except for the fit considering the CO absorption forest, which yielded a supersolar value at the grid's edge. However, caution should be exercised when interpreting this solar value, as it appears that the fits converged to a single node of the grid, possibly due to a bias in the interpolation step caused by small inhomogeneities in the grid. Two modes of C/O ratios are identified: one supersolar, when CO is considered in the fits, and one subsolar, when methane and potassium are used. Exo-REM demonstrated a successful reproduction of each spectral feature with a relatively low χ² compared to the other models, except for the potassium absorption lines, where the depth is not fully reached.

### 4.2.3. Sonora

The T<sub>eff</sub> estimates obtained from Sonora range from 1170 to 1401 K, with a deviant value of 1837 K when considering the CH<sub>4</sub> feature. These results are consistent with the two previous models and with Section 3. In contrast, we observed high surface gravities, with values exceeding 4.5 dex, which are inconsistent with the age of the system. Additionally, a supersolar metallicity is found, except when fitting using CH<sub>4</sub>. Similar to the adiabatic index γ explored by ATMO, the sedimentation factor f<sub>sed</sub> also varies across the size of the grid.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison between the data and the interpolated synthetic spectra generated with the set of parameters that maximized the likelihood for each model and for each spectral feature considered. The red area represents the wavelength coverage used to estimate $T_{\text{eff}}$ and the dilution factor $C_K$, while the green area represents the wavelengths used to fit the spectral features with $T_{\text{eff}}$ and $C_K$ fixed. The calculated $\chi_r^2$ values for each fit are provided, allowing for a comparison of each model's ability to reproduce the data.</figcaption></figure>

Sonora is capable of reproducing the different spectral features and generally exhibits the second-lowest $\chi_r^2$ compared to the other grids.

### 4.2.4. BT-Settl

This model estimates $T_{\text{eff}}$ to be between 1426 and 1718 K, with a low surface gravity of <4.5 dex, which is coherent with the spectral type and age of the object, respectively. Despite its ability to reproduce each spectral feature, the calculated $\chi_r^2$ is generally higher than that of the three previous models we discussed, indicating lower performance.

### 4.2.5. DRIFT-PHOENIX

As shown in Figure 7 and by its high $\chi_r^2$ values, this particular version of DRIFT-PHOENIX was unsuccessful in accurately reproducing the spectral information at medium

resolution. Consequently, making a robust estimate of atmospheric parameters becomes challenging.

## 5. Discussion

In Sections 3 and 4, we have presented two different approaches to leverage the spectral information observable from imaged exoplanets with JWST, aiming to estimate their atmospheric properties. The results of this study indicate a dispersion of these estimated parameters based on the wavelength coverage used for the fit, the considered spectral features, and the adopted model. We have proposed a method to incorporate this dispersion into a new error bar, attempting to provide a robust estimate for each parameter. This section is dedicated to discussing these results, investigating the origin of the dispersion, and analyzing the detected silicate absorption at 10 µm.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Parameter estimates provided by each model. The results obtained using each spectral feature are represented by different line formats. The colored areas represent the parameter spaces explored by each grid. The black lines indicate the estimates of T<sub>eff</sub> and log(g) from Dupuy et al. (2023) for the two mass scenarios, as well as the solar values of [M/H] and C/O.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9. Same as Figure 3 but for the bolometric luminosity $L_{bol}$ (left panel) and the mass $M$ (right panel). These parameters are estimated from the Stefan–Boltzmann law and the gravitational law, respectively.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

## 5.1. The Bolometric Luminosity and the Mass

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 \pm 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_{bol} = 4\pi\sigma \times R^2 \times T_{eff}^4$, where $\sigma$ 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.

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

$$m g = \frac{G M m}{R^2} \leftrightarrow M = \frac{g R^2}{G}, \qquad (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 \pm 0.5 M_{Jup}$ and another at $16 \pm 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.
## 5.2. The Origins of the Dispersion

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 $\lambda < 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_r^2$ values. This improvement is not observed in the case of DRIFT-PHOENIX and Sonora, where the calculated $\chi_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_r^2$ values ($\chi_r^2 = 126$ and $\chi_{r,A_V}^2 = 116$ for DRIFT-PHOENIX and $\chi_r^2 = 297$ and $\chi_{r,A_V}^2 = 148$ for Sonora). With each model, substantial deviations persist when considering $A_V$, particularly at $\lambda > 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).

### 5.3. Exploitation of the Silicate Absorption

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_λ \sim 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°+0°-28°; 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.

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Comparison between three examples of silicate absorption identified in the Spitzer library and the detection in the JWST data of VHS 1256 b. The spectral resolution of the JWST data has been reduced to match that of Spitzer for a fair comparison. The red area represents the spectral window used to estimate the pseudocontinuum, while the green area indicates the spectral window used to calculate the spectral index. In this figure, we specifically present the case of BT-Settl as an illustrative example.</figcaption></figure>

#### 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

$$I_{Si} = EW(Si) = \int_{\lambda} \left( 1 - \frac{F_{\lambda}}{C_{\lambda}} \right) d\lambda, \quad (4)$$

where $C_λ$ and $F_λ$ are the fluxes of the estimated pseudocontinuum and the data, respectively, at the wavelength $λ$, with $λ \in [8.0, 12.0]$ µm. The error bar on $I_{Si}$ is estimated by the calculation of $I_{Si,min}$ and $I_{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.

#### 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
<figure><img src="image_15.png" /><figcaption>**Figure 11.** The silicate index as a function of spectral type for the objects in the Spitzer library and VHS 1256 b. The top right panel displays the results obtained using the method defined by Suárez & Metchev (2023) for comparison. The other panels present the results obtained in this study, where the silicate index is calculated as an equivalent width. The open circles and filled circles represent young and field objects, respectively, and the star represents the estimate for VHS 1256 b. The colored areas indicate the spectral type ranges covered by the five different models used in the analysis.</figcaption></figure>

Suárez & Metchev (2023), enabling a direct comparison between the two approaches. Similar results were obtained for each model used. An $I_{Si} \sim 0$ is calculated for M-type objects, indicating a lack of significant silicate absorption in this spectral range. Additionally, the detection of CH$_4$ after the L-T transition, as well as the potential contribution of the NH$_3$ 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.

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$_2$SiO$_4$) was redder than enstatite (MgSiO$_3$). 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.

5.4. *Current and Future Improvements in Models*

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 temper-
ature 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, iso-
topologues (¹³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.

**6. Summary and Conclusion**

The main objective of this study was to conduct a
comprehensive analysis, employing the forward-modeling
approach, of the most detailed spectrum obtained to date for
an imaged planetary-mass companion. The aim was to assess
the capabilities and limitations of the current generation of five
different self-consistent atmospheric models and determine the
extent to which we can characterize the atmospheres of such
objects.

In Section 3, we presented an innovative method that took
advantage of the extended wavelength coverage provided by
JWST. By dividing the SED into multiple spectral windows,
we explored the dispersion of parameter estimates across these
windows and propagated this dispersion to determine new error
bars for each atmospheric property. Our analysis revealed that
relatively robust estimates could be obtained for T<sub>eff</sub>, log(g),
and radius. However, we observed that the log(g) and radius
tended to be underestimated compared to the predictions of
evolutionary models. This discrepancy resulted in under-
estimated values of the mass. On the other hand, the calculated
luminosity, derived from T<sub>eff</sub> and radius, was in good
agreement with evolutionary models. Interestingly, we
observed a curious trend of increasingly large fitted bolometric
luminosities as a function of wavelength region fit, which may
be attributed to a misrepresentation of dust, haze, and clouds in
the models of atmospheres. Estimating chemical abundances
([M/H] and C/O), as well as parameters such as γ and f<sub>sed</sub>,
proved to be more challenging due to the significant dispersion
observed along the SED. Furthermore, we found that the choice
of the model used for the fit had a significant impact on the final
parameter values, underscoring the importance of careful
model selection in interpreting atmospheric properties, con-
firming the results of Lueber et al. (2023) regarding the
estimation of log(g).

In Section 4, we employed a targeted approach to constrain
the [M/H] and C/O ratios by focusing on specific molecular
and atomic absorption features known to be sensitive to these
chemical abundances. By optimizing the fitting strategy for
each absorption feature, we obtained consistent estimates of
T<sub>eff</sub> and log(g) with those derived in Section 3. Regarding [M/
H], different scenarios emerged depending on the choice of the
atmospheric model. The Sonora model suggested a supersolar
metallicity, while the DRIFT-PHOENIX and ATMO models
indicated a subsolar metallicity. Lastly, the Exo-REM model
tended to converge toward solar values. Moreover, obtaining a
robust estimate of the C/O ratio proved to be challenging as
well. The ATMO model suggested supersolar values, similar to
the Exo-REM model when considering CO, and coherent with
previous studies (Hoch et al. 2022; Petrus et al. 2023).
Nevertheless, in both cases, the inversion reached the edge of
the grid.

After considering various sources of dispersion in the
estimation of atmospheric parameters, such as the high
variability of VHS 1256 b and potential biases during data
reduction, we turned our attention to systematic errors in the
five different models used. Indeed, there is a noticeable
discrepancy in accurately reproducing the complexity of the
atmosphere of VHS 1256 b, encompassing factors such as the
clouds' chemical composition and dispersion, the properties of
dust and haze, and various physical processes (vertical mixing, sedimentation, chemical disequilibrium, etc.). These systematics inherently contribute to biases in the determination of atmospheric parameters. Identifying and rectifying them will constitute one of the main challenges of exoplanetary atmospheric characterization in the coming decade.

To characterize the silicate absorption feature detected, we took advantage of the fact that its absorption band at ~10 µm is not reproduced by the models. This allowed us to estimate the pseudocontinuum between 8 and 12 µm, which enabled the calculation of the equivalent width of the silicate absorption. We compared it to those obtained from the silicate absorption detected in the spectra of isolated low-mass objects observed at low resolution with Spitzer. Through this analysis, we were able to confirm the correlation between the depth of the silicate absorption and the spectral type, as previously identified by Suárez & Metchev (2022). Additionally, we identified a connection between the silicate absorption feature and the age of the objects.

All of these results have demonstrated the capability of models to partially reproduce the observed data. The use of JWST data has allowed the development of optimized inversion procedures, which take into account the systematics present in the models. By performing a multimodel strategy and propagating these systematics into the error bars of the estimated parameters, we ensure a more accurate and realistic characterization of the atmospheric properties and avoid any potential overinterpretation of the results. However, constraining the formation history of VHS 1256 b using its atmospheric properties estimated from the forward-modeling approach remains challenging. A supersolar metallicity is usually used to quantify the solid material that is accreted during the formation, and the C/O is considered as a tracer of the birth location of the planet (see Nowak et al. 2020; Petrus et al. 2021; Hoch et al. 2022; Palma-Bifani et al. 2023), but with this study, we have shown that the estimates of these commonly used tracers of formation were strongly dependent on the model and the spectral information used for the fit. Moreover, Mollière et al. (2022) have demonstrated that the relationship between chemical abundances and the formation mechanisms described by the models (e.g., Öberg et al. 2011) is more complex than initially anticipated, and it is dependent on various intricate parameters, particularly in the case of formation within a disk. The chemical composition of the protoplanetary disk from which VHS 1256 b may have formed is unknown, but it needs to be estimated in order to establish reference values. An alternative approach is to estimate it from the host star. However, in our case, VHS 1256 (AB) is a short-period binary, which poses significant challenges in accurately determining its chemical abundances. Furthermore, Mollière et al. (2022) have demonstrated the significant influence of various factors on the evolution of the disk structure, such as gap opening, winds, and self-shadowing. They have also highlighted the impact of pebble evaporation inside ice lines on the accretion of chemical elements (in both gas and solid phases) during the planet's formation. These processes can potentially affect the overall metallicity and the carbon-to-oxygen ratio (C/O). Lastly, the authors emphasize that the chemical composition estimated solely from spectroscopic data may not necessarily reflect the bulk composition of the planet but rather describes the photospheric composition. In order to truly connect the measured atmospheric composition and the

original disk composition, it is crucial to gain a comprehensive understanding of the mechanisms that facilitate vertical chemical exchange between the upper and deeper layers of the planetary atmosphere.

Due to the presence of these unknown factors, it seems challenging to definitively determine the mode of formation of VHS 1256 b, and more generally of planetary-mass objects, by exclusively exploring its atmosphere. Nevertheless, the development and systematic application of spectral inversion methods, such as the one presented in this study, to the populations of imaged planetary-mass objects (both star companions and isolated objects) will enable a “comparative” analysis of their properties, evolving into a “statistical” study over the course of the JWST mission and once ELTs become operational around 2030.

## Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. We are truly grateful for the countless hours that thousands of people have devoted to the design, construction, and commissioning of JWST. This project was supported by a grant from STScI (JWST-ERS-01386) under NASA contract NAS5-03127. This work benefited from the 2022 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation. S.P. acknowledges the support of ANID, -Millennium Science Initiative Program-Center Code NCN19_171. This project received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (COBREX; grant agreement No. 885593) and from the ANR project FRAME (ANR-20-CE31-0012). J.M.V. acknowledges support from a Royal Society—Science Foundation Ireland University Research Fellowship (URF\1\221932). S.M. is supported by a Royal Society University Research Fellowship (URF-R1-221669). M.B. received funding from the European Union's Horizon 2020 research and innovation program (AtLAST; grant agreement No. 951815). R.A.M. is supported by the National Science Foundation MPS-Ascend Postdoctoral Research Fellowship under grant No. 2213312. E.G. acknowledges support from the Heising-Simons Foundation for this research. I.R. is supported by grant FJC2021-047860-I and PID2021-127289NB-I00 financed by MCIN/AEI /10.13039/501100011033 and the European Union NextGenerationEU/PRTR. A.Z. and S.P. acknowledge support from ANID—Millennium Science Initiative Program—Center Code NCN2021_080. All of the data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via 10.17909/ceq5-9g20.

## Appendix A
### Focus on the Fit Comparison

Figures 12–15 allow a comparison between the data and the best-fitting models with more details than in Figure 2.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Same as Figure 2. A zoom has been performed to illustrate the differences between the models and the data with more details. In this figure, the zoom was chosen to represent the data acquired with the channel G140HF/100LP of NIRSpec between 0.97 and 1.89 µm.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Same as Figure 12 but with the zoom chosen to represent the data acquired with the channel G235HF/170LP of NIRSpec between 1.89 and 3.17 µm.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Same as Figure 12 but with the zoom chosen to represent the data acquired with the channel G395HF/290LP of NIRSpec between 3.17 and 5.27 µm.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Same as Figure 12 but with the zoom chosen to represent the data acquired with MIRI between 4.90 and 17.98 µm.</figcaption></figure>
## Appendix B
Interstellar Extinction as a Free Parameter

Figure 16 illustrates the best-fitting model, achieved when integrating the interstellar extinction $A_V$ as defined by Draine (2003) into the list of parameters explored by ForMoSA,

following the concept introduced by Petrus et al. (2020). Notably, the $\chi^2$ value is substantially diminished in comparison to cases where $A_V$ is omitted (with the exception of DRIFT-PHOENIX and Sonora). However, despite these improvements, notable deviations persist in the reproduction of the observed data.

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Same figure as Figure 2 but with the exploration of the interstellar extinction $A_V$. The estimate of the other parameters, obtained with each model, are given in each panel, respectively.</figcaption></figure>
ORCID iDs

Simon Petrus https://orcid.org/0000-0003-0331-3654

Niall Whiteford https://orcid.org/0000-0001-8818-1544

Polychronis Patapis https://orcid.org/0000-0001-8718-3732

Beth A. Biller https://orcid.org/0000-0003-4614-7035

Andrew Skemer https://orcid.org/0000-0001-6098-3924

Sasha Hinkley https://orcid.org/0000-0001-8074-2562

Genaro Suárez https://orcid.org/0000-0002-2011-4924

Paulina Palma-Bifani https://orcid.org/0000-0002-6217-6867

Caroline V. Morley https://orcid.org/0000-0002-4404-0456

Pascal Tremblin https://orcid.org/0000-0001-6172-3403

Benjamin Charnay https://orcid.org/0000-0003-0977-6545

Johanna M. Vos https://orcid.org/0000-0003-0489-1528

Jason J. Wang https://orcid.org/0000-0003-0774-6502

Jordan M. Stone https://orcid.org/0000-0003-0454-3718

Mickaël Bonnefoy https://orcid.org/0000-0001-5579-5339

Gaël Chauvin https://orcid.org/0000-0003-4022-8598

Brittany E. Miles https://orcid.org/0000-0002-5500-4602

Aarynn L. Carter https://orcid.org/0000-0001-5365-4815

Anna Lueber https://orcid.org/0000-0001-6960-0256

Christiane Helling https://orcid.org/0000-0002-8275-1371

Ben J. Sutlieff https://orcid.org/0000-0002-9962-132X

Markus Janson https://orcid.org/0000-0001-8345-593X

Eileen C. Gonzales https://orcid.org/0000-0003-4636-6676

Kielan K. W. Hoch https://orcid.org/0000-0002-9803-8255

Olivier Absil https://orcid.org/0000-0002-4006-6237

William O. Balmer https://orcid.org/0000-0001-6396-8439

Anthony Boccaletti https://orcid.org/0000-0001-9353-2724

Mariangela Bonavita https://orcid.org/0000-0002-7520-8389

Mark Booth https://orcid.org/0000-0001-8568-6336

Brendan P. Bowler https://orcid.org/0000-0003-2649-2288

Zackery W. Briesemeister https://orcid.org/0000-0002-1764-2494

Marta L. Bryan https://orcid.org/0000-0002-6076-5967

Per Calissendorff https://orcid.org/0000-0002-5335-0616

Faustine Cantalloube https://orcid.org/0000-0002-3968-3780

Christine H. Chen https://orcid.org/0000-0002-8382-0447

Elodie Choquet https://orcid.org/0000-0002-9173-0740

Valentin Christiaens https://orcid.org/0000-0002-0101-8814

Gabriele Cugno https://orcid.org/0000-0001-7255-3251

Thayne Currie https://orcid.org/0000-0002-7405-3119

Camilla Danielski https://orcid.org/0000-0002-3729-2663

Matthew De Furio https://orcid.org/0000-0003-1863-4960

Trent J. Dupuy https://orcid.org/0000-0001-9823-1445

Samuel M. Factor https://orcid.org/0000-0002-8332-8516

Jacqueline K. Faherty https://orcid.org/0000-0001-6251-0573

Michael P. Fitzgerald https://orcid.org/0000-0002-0176-8973

Jonathan J. Fortney https://orcid.org/0000-0002-9843-4354

Kyle Franson https://orcid.org/0000-0003-4557-414X

Julien H. Girard https://orcid.org/0000-0001-8627-0404

Carol A. Grady https://orcid.org/0000-0001-5440-1879

Thomas Henning https://orcid.org/0000-0002-1493-300X

Dean C. Hines https://orcid.org/0000-0003-4653-6161

Callie E. Hood https://orcid.org/0000-0003-1150-7889

Alex R. Howe https://orcid.org/0000-0002-4884-7150

Paul Kalas https://orcid.org/0000-0002-6221-5360

Jens Kammerer https://orcid.org/0000-0003-2769-0438

Grant M. Kennedy https://orcid.org/0000-0001-6831-7547

Matthew A. Kenworthy https://orcid.org/0000-0002-7064-8270

Pierre Kervella https://orcid.org/0000-0003-0626-1749

Minjae Kim https://orcid.org/0000-0001-6218-2004

Daniel Kitzmann https://orcid.org/0000-0003-4269-3311

Adam L. Kraus https://orcid.org/0000-0001-9811-568X

Masayuki Kuzuhara https://orcid.org/0000-0002-4677-9182

Anne-Marie Lagrange https://orcid.org/0000-0002-2189-2365

Kellen Lawson https://orcid.org/0000-0002-6964-8732

Cecilia Lazzoni https://orcid.org/0000-0001-7819-9003

Jarron M. Leisenring https://orcid.org/0000-0002-0834-6140

Ben W. P. Lew https://orcid.org/0000-0003-1487-6452

Michael C. Liu https://orcid.org/0000-0003-2232-7664

Pengyu Liu https://orcid.org/0000-0001-7047-0874

Jorge Llop-Sayson https://orcid.org/0000-0002-3414-784X

Bruce Macintosh https://orcid.org/0000-0003-1212-7538

Mathilde Mâlin https://orcid.org/0000-0002-2918-8479

Elena Manjavacas https://orcid.org/0000-0003-0192-6887

Sebastián Marino https://orcid.org/0000-0002-5352-2924

Mark S. Marley https://orcid.org/0000-0002-5251-2943

Christian Marois https://orcid.org/0000-0002-4164-4182

Raquel A. Martinez https://orcid.org/0000-0001-6301-896X

Elisabeth C. Matthews https://orcid.org/0000-0003-0593-1560

Brenda C. Matthews https://orcid.org/0000-0003-3017-9577

Dimitri Mawet https://orcid.org/0000-0002-8895-4735

Johan Mazoyer https://orcid.org/0000-0002-9133-3091

Michael W. McElwain https://orcid.org/0000-0003-0241-8956

Stanimir Metchev https://orcid.org/0000-0003-3050-8203

Michael R. Meyer https://orcid.org/0000-0003-1227-3084

Maxwell A. Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

Paul Mollière https://orcid.org/0000-0003-4096-7067

Sarah E. Moran https://orcid.org/0000-0002-6721-3284

Sagnick Mukherjee https://orcid.org/0000-0003-1622-1302

Eric Pantin https://orcid.org/0000-0001-6472-2844

Marshall D. Perrin https://orcid.org/0000-0002-3191-8151

Laurent Pueyo https://orcid.org/0000-0003-3818-408X

Sascha P. Quanz https://orcid.org/0000-0003-3829-7412

Andreas Quirrenbach https://orcid.org/0000-0002-3302-1962

Shrishmoy Ray https://orcid.org/0000-0003-2259-3911

Isabel Rebollido https://orcid.org/0000-0002-4388-6417

Jea Adams Redai https://orcid.org/0000-0002-4489-3168

Bin B. Ren https://orcid.org/0000-0003-1698-9696

Emily Rickman https://orcid.org/0000-0003-4203-9715

Steph Sallum https://orcid.org/0000-0001-6871-6775

Matthias Samland https://orcid.org/0000-0001-9992-4067

Benjamin Sargent https://orcid.org/0000-0001-9855-8261

Joshua E. Schlieder https://orcid.org/0000-0001-5347-7062

Karl R. Stapelfeldt https://orcid.org/0000-0002-2805-7338

Motohide Tamura https://orcid.org/0000-0002-6510-0681

Xianyu Tan https://orcid.org/0000-0003-2278-6932

Christopher A. Theissen https://orcid.org/0000-0002-9807-5435

Taichi Uyama https://orcid.org/0000-0002-6879-3030
Malavika Vasist https://orcid.org/0000-0002-4511-3602

Arthur Vigan https://orcid.org/0000-0002-5902-7828

Kevin Wagner https://orcid.org/0000-0002-4309-6343

Kimberly Ward-Duong https://orcid.org/0000-0002-4479-8291

Schuyler G. Wolff https://orcid.org/0000-0002-9977-8255

Kadin Worthen https://orcid.org/0000-0002-5885-5779

Mark C. Wyatt https://orcid.org/0000-0001-9064-5598

Marie Ygouf https://orcid.org/0000-0001-7591-2731

Alice Zurlo https://orcid.org/0000-0002-5903-8316

Xi Zhang https://orcid.org/0000-0002-8706-6963

Keming Zhang https://orcid.org/0000-0002-9870-5695

Zhoujian Zhang https://orcid.org/0000-0002-3726-4881

Yifan Zhou https://orcid.org/0000-0003-2969-6040

## References

Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872

Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357

Allard, F., Hauschildt, P. H., Baraffe, I., & Chabrier, G. 1996, ApJL, 465, L123

Allard, F., Homeier, D., Freytag, B., & Sharp, C. M. 2012, in EAS Publications Series, ed. C. Reylé, C. Charbonnel, & M. Schultheis, Vol. 57 (Les Ulis: EAS Publication Series), 3

Allers, K. N., & Liu, M. C. 2013, ApJ, 772, 79

Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111

Bayo, A., Barrado, D., Allard, F., et al. 2017, MNRAS, 465, 760

Betti, S. K., Follette, K. B., Ward-Duong, K., et al. 2022, ApJL, 935, L18

Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82

Bonnefoy, M., Chauvin, G., Lagrange, A. M., et al. 2014, A&A, 562, A127

Boss, A. P. 2000, ApJ, 536, L101

Bowler, B. P., Zhou, Y., Morley, C. V., et al. 2020, ApJL, 893, L30

Brown, A. G. A., Gaia Collaboration, Gaia Collaboration, Vallenari, A., et al. 2021, A&A, 649, A1

Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, NatAs, 2, 138

Burgasser, A. J., Marley, M. S., Ackerman, A. S., et al. 2002, ApJL, 571, L151

Burningham, B., Faherty, J. K., Gonzales, E. C., et al. 2021, MNRAS, 506, 1944

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022a, JWST Calibration Pipeline, v1.7.0, Zenodo, doi:10.5281/zenodo.7071140

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022b, JWST Calibration Pipeline, v1.8.1, Zenodo, doi:10.5281/zenodo.7215690

Caffau, E., Ludwig, H. G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, SoPh, 268, 255

Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172

Cushing, M. C., Marley, M. S., Saumon, D., et al. 2008, ApJ, 678, 1372

Cushing, M. C., Roellig, T. L., Marley, M. S., et al. 2006, ApJ, 648, 614

Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2015, ApJ, 812, 128

Daemgen, S., Todorov, K., Quanz, S. P., et al. 2017, A&A, 608, A71

Draine, B. T. 2003, ARA&A, 41, 241

Dupuy, T. J., Liu, M. C., Evans, E. L., et al. 2023, MNRAS, 519, 1688

Eriksson, S. C., Asensio Torres, R., Janson, M., et al. 2020, A&A, 638, L6

Fortney, J. J., Marley, M. S., Saumon, D., & Lodders, K. 2008, ApJ, 683, 1104

Gandhi, S., de Regt, S., Snellen, I., et al. 2023, ApJL, 957, L36

Gauza, B., Béjar, V. J. S., Pérez-Garrido, A., et al. 2015, ApJ, 804, 96

Greco, J. P., & Brandt, T. D. 2016, ApJ, 833, 134

Grevesse, N., Noels, A., & Sauval, A. J. 1992, in ESA Special Publication, Vol. 348, Coronal Streamers, Coronal Loops, and Coronal and Solar Wind Composition, ed. C. Mattok (Paris: ESA), 305

Hauschildt, P. H., Baron, E., & Allard, F. 1997, ApJ, 483, 390

Helling, C., Dehn, M., Woitke, P., & Hauschildt, P. H. 2008, ApJL, 675, L105

Helling, C., & Woitke, P. 2006, A&A, 455, 325

Hinkley, S., Carter, A. L., Ray, S., et al. 2022, PASP, 134, 095003

Hoch, K. K. W., Konopacky, Q. M., Barman, T. S., et al. 2022, AJ, 164, 155

Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165

Hurt, S. A., Liu, M. C., Zhang, Z., et al. 2023, 961, 121

Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80

Karalidi, T., Marley, M., Fortney, J. J., et al. 2021, ApJ, 923, 269

Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Sci, 339, 1398

Lachapelle, F.-R., Lafrenière, D., Gagné, J., et al. 2015, ApJ, 802, 61

Lafrenière, D., Jayawardhana, R., & van Kerkwijk, M. H. 2010, ApJ, 719, 497

Leconte, J. 2018, ApJL, 853, L30

Lodders, K. 2002, ApJ, 577, 974

Lodders, K. 2010, Astrophysics and Space Science Proceedings, Principles and Perspectives in Cosmochemistry, Vol. 16 (Berlin: Springer), 379

Looper, D. L., Kirkpatrick, J. D., Cutri, R. M., et al. 2008, ApJ, 686, 528

Lueber, A., Kitzmann, D., Fisher, C. E., et al. 2023, ApJ, 954, 22

Luna, J. L., & Marley, C. V. 2021, ApJ, 920, 146

Manjavacas, E., Bonnefoy, M., Schlieder, J. E., et al. 2014, A&A, 564, A55

Marley, M. S., & McKay, C. P. 1999, Icar, 138, 268

Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJL, 723, L171

Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Sci, 272, 1919

Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85

Miles, B. E., Biller, B. A., Patapis, P., et al. 2023, ApJL, 946, L6

Miles, B. E., Skemer, A. J., Barman, T. S., Allers, K. N., & Stone, J. M. 2018, ApJ, 869, 18

Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74

Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172

Mukherjee, S., Batalha, N. E., Fortney, J. J., & Marley, M. S. 2023, ApJ, 942, 71

Nowak, M., GRAVITY Collaboration, GRAVITY Collaboration, Lacour, S., et al. 2020, A&A, 633, A110

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A&A, 647, A175

Padoan, P., Kritsuk, A., Norman, M. L., & Nordlund, Å. 2005, ApJL, 622, L61

Palma-Bifani, P., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, A90

Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2020, A&A, 633, A124

Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59

Petrus, S., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, L9

Plummer, M. K., & Wang, J. 2022, ApJ, 933, 163

Pollack, J. B., Murray-Clay, R., & Bodenheimer, P., et al. 1996, Icar, 124, 62

Ringqvist, S. C., Viswanath, G., Aoyama, Y., et al. 2023, A&A, 669, L12

Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200

Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327

Showman, A. P., Tan, X., & Parmentier, V. 2020, SSRv, 216, 139

Skilling, J. 2004, in AIP Conf. Ser. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint (Melville, NY: AIP), 395

Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Natur, 509, 63

Sorahana, S., & Yamamura, I. 2012, ApJ, 760, 151

Spake, J. J., Sing, D. K., Wakeford, H. R., et al. 2021, MNRAS, 500, 4042

Stephens, D. C., Leggett, S. K., Cushing, M. C., et al. 2009, ApJ, 702, 154

Stone, J. M., Skemer, A. J., Kratter, K. M., et al. 2016, ApJL, 818, L12

Suárez, G., & Metchev, S. 2022, MNRAS, 513, 5701

Suárez, G., & Metchev, S. 2023, MNRAS, 523, 4739

Suárez, G., Metchev, S., Leggett, S. K., Saumon, D., & Marley, M. S. 2021, ApJ, 920, 99

Suárez, G., Vos, J. M., Metchev, S., Faherty, J. K., & Cruz, K. 2023, ApJL, 954, L6

Tan, X., & Showman, A. P. 2021a, MNRAS, 502, 678

Tan, X., & Showman, A. P. 2021b, MNRAS, 502, 2198

Tremblin, P., Hubundsen, D. S., Mourier, P., et al. 2015, ApJL, 804, L17

Tremblin, P., Chabrier, G., Baraffe, I., et al. 2017, ApJ, 850, 46

Tsuji, T., Ohnaka, K., Aoki, W., & Nakajima, T. 1996, A&A, 308, L29

Venot, O., Hébrard, E., Agúndez, M., et al. 2012, A&A, 546, A43

Vos, J. M., Burningham, B., Faherty, J. K., et al. 2023, ApJ, 944, 138

Wang, J., Wang, J. J., Ruffio, J.-B., et al. 2023, AJ, 165, 4

Wang, J. J., Delorme, J.-R., Ruffio, J.-B., et al. 2021, Proc. SPIE, 11823, 1182302

Wilcomb, K. K., Konopacky, Q. M., Barman, T. S., et al. 2020, AJ, 160, 207

Witte, S., Helling, C., Barman, T., Heidrich, N., & Hauschildt, P. H. 2011, A&A, 529, A44

Witte, S., Helling, C., & Hauschildt, P. H. 2009, A&A, 506, 1367

Woitke, P., & Helling, C. 2003, A&A, 399, 297

Woitke, P., & Helling, C. 2004, A&A, 414, 335

Wright, G. S., Rieke, G. H., Glasse, A., et al. 2013, PASP, 135, 048003

Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41

Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021, Natur, 595, 370

Zhou, Y., Bowler, B. P., Apai, D., et al. 2022, AJ, 164, 239

Zhou, Y., Bowler, B. P., Morley, C. V., et al. 2020, AJ, 160, 77