© The Authors 2026

1 Introduction

Although only a few planetary-mass companions have been detected and characterized using direct imaging¹ (approximately 80), this sample includes diverse objects with effective temperatures (T_eff) ranging from 270 to 2700 K and orbital separations spanning from a few tens to thousands of au. The majority of these companions are young (<300 Myr) and the observed diversity of properties (e.g., log(LZL_⊙), T_eff, log(g), C/O, [M/H], e) is a product of each object’s formation and evolution history (e.g., Spiegel & Burrows 2012; Allers & Liu 2013; Ruffio et al. 2026).

The COCONUTS-2 (L 34-26 and WISEPA J075108.79-763449.6) system stands out as an outlier in the field of direct imaging. It consists of a cold ∼8 M_Jup (Zhang et al. 2025), T9 (Kirkpatrick et al. 2011) companion orbiting an M3 star with an estimated age of 414 ± 23 Myr (Kiman et al. 2026), thus providing a unique opportunity to study the atmosphere of a mature super-Jupiter. Originally thought to be two separate free-floating objects, Zhang et al. (2020, 2021a) demonstrated, through proper motion and parallax measurements, that COCONUTS-2 A and b are gravitationally bound. They estimated the semimajor axis of COCONUTS-2 b to be ${7506}_{-2060}^{+5205}$ au and its effective temperature to be T_eff = 434±9 K, making it the second-coldest exoplanet and the one with the largest directly imaged separation discovered to date. The large separation poses a challenge for planet formation models. COCONUTS-2b may have formed within the circumstellar disk of A, either through core accretion (Pollack et al. 1996; Ida & Lin 2008; Alibert et al. 2005, 2013; Benz et al. 2014) or disk instability (Kuiper 1951; Cameron 1978; Boss 1997), followed by outward migration due to dynamical interactions (Vorobyov 2013). Alternatively, COCONUTS-2b could have formed via gravitational instability (Padoan & Nordlund 2002) in a collapsing molecular cloud, similarly to a stellar binary (Zuckerman & Song 2009; Kirkpatrick et al. 2011). A less likely scenario is a capture during a flyby, as suggested by Marocco et al. (2024). The recent extensive study from Zhang et al. (2025) using the combination of WISE and Spitzer photometry with Gemini/FLAMINGO-2 spectroscopy (0.98-2.51 μm, R_λ ~ 2001200) alongside multiple atmospheric grids revealed that COCONUTS-2 b likely has a substellar or near-stellar C/O ratio and metallicity ([M/H]). This may suggest that it accreted oxygen-rich and carbon-depleted gas from within ice lines (Line et al. 2021; August et al. 2023), and was not influenced by core erosion processes or planetesimal bombardment (e.g., Alibert et al. 2013; Miller & Fortney 2011; Thorngren & Fortney 2019; Madhusudhan 2019; Schneider & Bitsch 2021; Zhang et al. 2023). This could indicate a planet-like formation followed by outward migration. However, the spread of values in their retrieved metallicities across models (ranging from -0.4 to 0.1) remains fairly consistent with the metallicity of COCONUTS-2 A (0.00 ± 0.08; Hojjatpanah et al. 2019), and therefore cannot be used to rule out a scenario of a binary-like formation, as the metallicities of the host star and the companion coincide.

We present in this study new observations from JWST/MIRI-LRS that we analyzed jointly with existing spectroscopic data to consolidate the estimate of the physical properties of COCONUTS-2 b. The combination of the very high S/N of these data and the low effective temperature of COCONUTS-2 b presents a challenge for classical forward modeling approaches (Petrus et al. 2024). These temperature and age regimes remain largely unexplored, and self-consistent models may suffer from systematic deviations, either due to missing or incomplete physics (e.g., disequilibrium chemistry, inhomogeneous atmospheres) or due to limitations in parameter space exploration. Our analysis method is presented in Sect. 2, and the results of our analysis are given in Sect. 3. We discuss the implications of our results in Sect. 4, and present our conclusions in Sect. 5. Additional material is provided in the Appendices.

2 Methods

2.1 Observation and data reduction

We used both archival and newly obtained spectro-photometric data on the target. This section describes each of them. Table 2 provides a concise summary of all of them.

Photometry. COCONUTS-2 b was observed with both Spitzer and WISE as part of a survey of 184 late T and Y dwarfs presented in Kirkpatrick et al. (2019). We used the Spitzer [3.6] and [4.5] points, as well as the WISE W1 (3.32 ± 0.33 μm) and W2 (4.56 ± 0.52 μm). The W3 and W4 bands were excluded as their bandpasses are not fully covered by the Sonora Elf Owl grid; W4 lies entirely outside the modeled wavelength range. We note that this has little impact on the retrieved parameters for the other grids.

Gemini/FLAMINGOS-2. FLAMINGOS-2 is a low- to medium-resolution spectrograph located on the 8.1 m Gemini-South telescope. We used archival FLAMINGOS-2 observations of COCONUTS-2 b obtained as part of programs GS-2021B-FT-111 and GS-2021B-FT-204 and first presented in Zhang et al. 2025. The instrument provides simultaneous coverage of the Y-, J-, H-, and K-bands (1.01-2.50 μm). COCONUTS-2 b was observed on 16 and 19 December 2021 with the 2 pixel-wide slit (0.36” × 263”) and an ABBA nodding pattern, resulting in an average spectral resolution of R_λ ~ 900 (Zhang et al. 2025). Because COCONUTS-2 A and b are widely separated, the adaptative optics (AO) system was not needed. To extract the final spectrum, a standard data reduction was applied using the Gemini-IRAF² package.

JWST/NIRSpec. As part of the “Explaining the Diversity of Cold Worlds” survey (Program ID: 2124, PI: J. Faherty, DOI), COCONUTS-2 b was observed with JWST/NIRSpec. Observations took place on 9 July 2023 for a total integration time of 2777.706 s and were reduced and presented in Kiman et al. 2026. They used the G395H grating with the F290LP filter element resulting in an average spectral resolution of R_λ ~ 3000 and a wavelength coverage of 2.88-5.14 μm. We convolved the spectra to a lower resolution of R_λ ~ 1500 in our analysis to match the sampling of the ATMO2020++ grids.

JWST/MIRI-LRS. MIRI-LRS on the JWST is a low-resolution spectrograph operating in the mid-infrared. It was used to observe COCONUTS-2 b on 29 March 2024 (Program ID: 3514, PI: M. Bonnefoy, DOI) for a total integration time of 3335.598 s. The observation strategy was based on using a two-dither pattern with two integrations of 300 groups per exposure. The data were reduced using a combination of the official JWST STScI pipeline³, supplemented with customized functions described in Voyer et al. (2025), Sect. 2.2. As noted in other studies (Xuan et al. 2024; Voyer et al. 2025), the retrievals showed a significant wavelength-dependent offset between the line positions of the data and the theoretical models. To correct these offsets we include a third-order polynomial for the wavelength correction in the retrieval. The coefficients of this correction are retrieved once and then fixed for subsequent analysis following the methods in Xuan et al. (2024); Voyer et al. (2025). The wavelength solution from Table 1⁴ of Voyer et al. (2025) was used for the pixel-to-wavelength calibration. Beyond 12.5 μm, the extracted spectrum is dominated by background emission. With only two integrations per exposure, the background cannot be reliably averaged out, and the resulting data quality is insufficient for scientific analysis. We therefore excluded these longer wavelengths from the analysis.

2.2 Framework

To infer the atmospheric properties of COCONUTS-2 b, we used a custom version of the Forward Modeling tool for Spectral Analysis, ForMoSA⁵ (Ravet et al. 2025). It relies on the nested-sampling (Skilling 2004) approach using the Python package PyMultinest⁶ (Buchner 2016). Each inversion utilizes 1000 live points, adaptive sampling efficiency mode, an evidence tolerance of 0.5; uniform priors were used for the atmospheric parameters. We also fit for the planet radius R constrained by the flux dilution factor ${\left(\frac{\text{R}}{d}\right)}^2$ using a fixed distance of 10.888 (±0.002) pc (see Table 1). Nested sampling allows the direct computation of the log-Bayesian evidence, ln Z, thus enabling a robust quantitative comparison between models via the log-Bayes factor, ln B_1,2 = ln Z₁ - lnZ₂. This metric inherently accounts for both model fit and complexity (i.e., number of free parameters). To facilitate an intuitive interpretation, the Bayes factor is often approximately translated into a sigma significance level, representing the strength of preference for one model over another (Sellke et al. 2001; Trotta 2008; Benneke & Seager 2013). However, as recently pointed out by Kipping & Benneke (2025), such a conversion may only provide an upper limit on the sigma values. The true values are generally lower, which can lead to an overestimation of the significance of model preferences. We therefore chose to rely on the log-Bayesian evidence as a more conservative proxy for model significance. Importantly, the trend we find is also consistent with the Akaike information criterion (AIC; Akaike 1974), the Bayesian information criterion (BIC; e.g., Mollière et al. 2025), and the simplified Bayesian predictive information criterion (BPICS; Ando 2011; Thorngren et al. 2026) across all tested models (see Table 3). For the rest of the study, we define these quantities relative to the favored model:

$$\begin{array}{rl}& \ln \mathcal{B}=\ln {\mathcal{Z}}_{max}-\ln \mathcal{Z},\\ & \mathrm{\Delta }\text{AIC}=\text{AIC}-{\text{AIC}}_{min},\\ & \mathrm{\Delta }\text{BIC}=\text{BIC}-{\text{BIC}}_{min},\\ & \mathrm{\Delta }\text{BPICS}=\text{BPICS}-{\text{BPICS}}_{min}.\end{array}$$(1)

The subscripts max and min denote the maximum and minimum values across all models. With this convention, the preferred model has ln B = 0, while all other models have ln B > 0.

2.3 Likelihood mapping

Since using only the injected error bars coming from the data-reduction pipelines likely underestimates the total uncertainty, we introduce in ForMoSA the following modified log-likelihood function:

$$\begin{array}{rl}& \ln {\mathcal{L}}^{\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{c}}=\sum _i\ln {\mathcal{L}}_i^{\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{c}}\\ & \ln {\mathcal{L}}_i^{\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{c}}=-\frac{N_i}{2}\ln \left(\frac{{\chi }_i^2}{N_i}\right)-\frac{1}{\text{2}}\ln (|{\sigma }_{\mathbf{i}}|)-\frac{N_i}{2}\ln (2\pi )-\frac{N_i}{2}\\ & \text{with}{\chi }_i^2=\sum {\left(\frac{{\mathbf{d}}_{\mathbf{i}}-{\mathbf{m}}_{\mathbf{i}}}{{\sigma }_{\mathbf{i}}}\right)}^2.\end{array}$$(2)

Here the index i takes values from 1 to 4 and corresponds to the three spectroscopic observations and four photometric points; d_i, σ_i, and m_i refer to the data, error, and model vectors, respectively, each of length N_i. This likelihood formulation proposed by Ruffio et al. (2019) is obtained by marginalizing the standard loglikelihood over a noise scaling factor σ_s = sσ. The maximum likelihood estimator for this scaling is

$$\widehat{s_i}=\frac{{\chi }_i^2}{N_i}.$$(3)

This approach partially compensates for model-data mismatches by inflating the error bars, thus propagating them into the posterior distributions. Other approaches exist in the literature, including alternative multiplicative or additive prescriptions (e.g., Piette & Madhusudhan 2020; Zhang et al. 2025). Here we adopted a marginalization over a simple multiplicative factor as a computationally efficient choice that does not assume a particular functional form of the inflation. In a forward modeling framework, the quality of the fit is primarily limited by the fidelity of the atmospheric models rather than by the formal observational uncertainties (Petrus et al. 2025), making the inclusion of these noise scaling parameters a pragmatic way to account for residual model deficiencies without biasing the posteriors. However, they entirely neglect potential spectral correlations introduced by the spectrograph or by syst-mematic deviations between data and model spread over large bandwidths. To address this limitation, we employed Gaussian processes (GPs) to model correlated noise between neighboring pixels. Building upon previous work (Czekala et al. 2015; Wang et al. 2020; de Regt et al. 2024, 2025; Rotman et al. 2025), we define the covariance matrix as

$$C_{i;m,n}={\sigma }_{i;n}^2{\delta }_{m,n}+a_i^2{\widetilde{{\sigma }_i}}^2\exp [-\frac{\mathrm{\Delta }{\lambda }_{m,n}^2}{2{\ell }_i^2}],$$(4)

where each C_{i; m,n} represents the total covariance between pixels m and n in the i-th spectrum. The matrix comprises two components: diagonal terms ${\sigma }_{i;n}^2$ δ_m,n that account for uncorrelated noise, and radial basis function (RBF) kernel terms that model spectral correlations. The RBF kernel is parameterized by an amplitude a_i and a length scale I_i, and is normalized by the squared median uncertainty ${\widetilde{{\sigma }_i}}^2$; ∆λ_m,n denotes the wavelength separation between data points m and n. Injecting this covariance matrix C_i into the standard log-likelihood, we obtain

$$\begin{array}{rl}& \ln {\mathcal{L}}^{\mathrm{G}\mathrm{P}}=\sum _i\ln {\mathcal{L}}_i^{\mathrm{G}\mathrm{P}}\\ & \ln {\mathcal{L}}_i^{\mathrm{G}\mathrm{P}}=-\frac{{\chi }_i^2}{2}-\frac{1}{\text{2}}\ln (|{\mathbf{C}}_{\mathbf{i}}|)-\frac{N_i}{2}\ln (2\pi )\\ \\ & \text{with}{\chi }_i^2=({\mathbf{d}}_{\mathbf{i}}-{\mathbf{m}}_{\mathbf{i}}{)}^T{{\mathbf{C}}_{\mathbf{i}}}^{-1}({\mathbf{d}}_{\mathbf{i}}-{\mathbf{m}}_{\mathbf{i}}).\end{array}$$(5)

This likelihood is very similar to the one presented in de Regt et al. (2025), except that we do not fit for a noise scaling factor. This choice is motivated by the fact that including a noise scaling term introduces strong biases in the model-data residuals for MIRI-LRS (see bottom panel of Fig. 1, bottom panels of Fig. 4, and discussion in Sect. 4.2). Equations (2) and (5) assume that each dataset is independent, a reasonable assumption given that each observation was taken with a different instrument at a different time. Each log(a_i) and log(l_i) are retrieved as free parameters during the inversion.

3 Results

We divided our analysis into three parts. First, we performed an exploration using different atmospheric models evaluated with a simple noise scaling (Eq. (2)). Then we carried out a more detailed investigation using GP-aided forward modeling (Eq. (5)), focusing on the best-fitting models identified in the first stage. Finally, we used the preferred solutions to rederive the bolometric luminosity using the full spectral energy distribution (SED) and, in combination with evolutionary models, infer the corresponding bulk properties. The aim of this strategy is to mitigate some of the limitations of current forward modeling approaches (such as model systematics and correlated noise) to derive more robust constraints on the atmospheric parameters of COCONUTS-2 b.

Similarly to Zhang et al. (2025), we masked five regions of the FLAMINGOS-2 observation and performed the fit on the following ranges: 1.05-1.12 μm, 1.18-1.33 μm, and 1.52-1.66 μm (see Fig. 1). Fluxes from wavelengths shorter than 1.05 μm were masked to avoid the area with low S/N near the edge of the detector and valleys between the Y, J, H, and K peaks due to the residuals from telluric correction and the relatively large flux uncertainties in the spectrum. We also ignored any variability and did not adjust for relative flux offsets between instruments as no variability has been reported for COCONUTS-2 b.

3.1 Model comparison

Figures 1 and 2 and Tables 3 and A.3 summarize the inversion results using the five different models. While most models manage to reproduce the overall shape of the combined observations, we observed significant data-model mismatches across all instruments. We observe an important dispersion between models and their retrieved posterior distributions (see Fig. 2), with almost no overlap. This huge dispersion suggests that the models are unable to robustly extract the physical parameters. Comparing the log-Bayes factors (see Table 3), the ATMO2020++ models (with and without PH₃) are statistically preferred among those tested. This result is in line with a previous atmospheric analysis of the target (Zhang et al. 2025), although the model with PH₃ seems to be preferred (ln B = 138). The PH₃ detection significance is discussed in Sect. 4.3.

The retrieved effective temperatures range from 450 to 540 K for models that converged within their parameter grids (excluding Sonora Diamondback, which converges at the grid edge with T_eff > 600 K). This is broadly consistent with previous expectations from evolutionary tracks, which predict T_eff = ${483}_{-53}^{+44}$ K (Zhang et al. 2025). The lower limit of T_eff > 600 K obtained with Sonora Diamondback probably explains the broader constraint observed in luminosity (see purple histogram in lower left panel of Fig. 2) due to the degeneracy with the radius. We note that this temperature is highly inconsistent with evolutionary predictions. The surface gravity is poorly constrained, with retrieved values between 4.5 and 5.0 dex, significantly higher than the evolutionary estimate of log(g) = ${4.19}_{-0.13}^{+0.18}$ dex. Notably, inversions using BT-Settl and Sonora Elf Owl reached the edges of their grids, with posteriors at >4.5 dex and <3.25 dex, respectively. However, comparisons between atmospheric and evolutionary log(g) estimates should be treated with caution as they trace different physical processes (e.g., atmospheric structure versus radius contraction; Petrus et al. 2025).

Except for Sonora Elf Owl, all grids that explored metallic-ity converged toward slightly supersolar values ([M/H] ∈ 0.20.4 dex). Since we used the same dataset, it is reasonable to assume that the observed differences in metallicity between the various inversions arise from systematic offsets between the models. Because metallicity is correlated with log(g) (Zhang et al. 2021b), the systematically higher surface gravities retrieved for most models in this setup may also bias the inferred metal-licity. Moreover, metallicity is highly sensitive to the wavelength range and datasets included in the inversion (Petrus et al. 2024). Sonora Elf Owl is the only grid that also probes the carbon-to-oxygen ratio (C/O), yielding a subsolar value of 0.317 ± 0.005.

Of the two grids that include cloud modeling, only Sonora Diamondback explores their effects explicitly through the sedimentation efficiency parameter of its silicate clouds (f_sed), which ranges from 1 (thick clouds) to 8 (thin clouds). Our best-fit value of f_sed = ${3.6}_{-0.2}^{+0.1}$ suggests a moderately cloudy atmosphere. However, since both cloudy models (BT-Settl and Sonora Diamondback) are among the least favored in our analysis, this constraint should be interpreted with caution as it may not provide accurate information about the actual cloud content and properties of COCONUTS-2 b.

The vertical (eddy) diffusion coefficient K_zz (in cm² s⁻¹) parameterizes the efficiency of atmospheric vertical mixing. In the Sonora Elf Owl grid, this parameter is explored logarithmically from 10² cm² s⁻¹ (inefficient mixing) to 10⁹ cm² s⁻¹ (efficient mixing). We retrieved a value of K_zz = (${9.8}_{-0.7}^{+0.5}$) × 10³ cm² s⁻¹, indicative of moderate vertical mixing. The vertical mixing can also be retrieved analytically from the inferred log(g) using Fig. 1 of Phillips et al. (2020) for the two ATMO2020++ grids. Both ATMO2020++ models suggest a more vigorous vertical mixing with K_zz = (2.5 ± 0.1) × 10⁵ cm².s⁻¹ and K_zz = (1.3 ± 0.1) × 10⁵ cm².s⁻¹ for the model with and without PH₃, respectively. This strong vertical mixing can also be seen in Fig. A.5, where all the key chemical abundances are quenched. However, as noted and explored in Mukherjee et al. (2022, 2024), K_zz is expected to vary significantly with pressure, depending on whether the local atmospheric layer is radiative or convective. This effect is not accounted for in this setup. In addition, the retrieved K_zz may partially encode viewing-geometry effects, further contributing to the degeneracy of this parameter (Tan & Showman 2021; Suárez et al. 2023; Petrus et al. 2025). Recent studies have investigated how mixing strength may vary in more complex 3D atmospheres (Visscher et al. 2010; Mukherjee et al. 2022; Lacy & Burrows 2023).

The top left panel of Fig. 2 compares the retrieved pressure-temperature (P-T) profiles for Sonora Diamondback (purple), Sonora Elf Owl (yellow), and ATMO2020++ (green). All the profiles agree well with each other.

Fig. 1 Top panel : forward modeling results of the combined COCONUTS-2 b observations. The black solid lines represent the spectroscopic data (from left to right: Gemini/FLAMINGOS-2, JWST/NIRSpec, and JWST/MIRI-LRS), while the colored lines represents each Rλ ∼ 100 model (classic). The pink and white stars indicate the photometric observations (WISE and Spitzer). The gray shaded regions indicate the masked wavelengths during the spectral inversion. Bottom panel : residuals for each fit. The dotted lines represents the ±σ (68%) confidence interval.

3.2 GP analysis

In Sect. 3.1, we identified ATMO2020++ and its variation ATMO2020++ (no PH₃) as the two preferred models; ATMO2020++ (no PH₃) is the only model that does not converge outside its grid range for any parameter. Using the GP framework presented in Sect. 2.2, we reanalyzed the combined dataset. Assuming that each dataset is independent, we fit for each covariance matrix separately, resulting in three length-scales and three amplitudes as additional extra-grid parameters.

This section focuses on the two ATMO2020++ models; however, we report the final parameters obtained with the other grids in Tables A.3 and A.4, and additional elements of discussion are provided in Appendix D. Notably, based on the various significance criteria (Table 4), the Sonora Elf Owl models become preferred when GPs are included, highlighting the importance of the GP framework in the modeling. However, these models still provide a comparatively poor fit, both visually (large residual offsets) and in terms of the inferred physical parameters (inconsistent with cooling model predictions). We therefore retain the focus of this section on the two ATMO2020++ grids. Final posteriors are compared with the classical approach and evolutionary models predictions in Fig. 3 and summary plots can be found in Figs. A.1, A.2, and 4. In particular, in Fig. 4, the GP approach significantly improves the fit to the MIRI-LRS observation, where the flux was previously underestimated. The effect is less noticeable for the other two spectra. In the first column of Table A.3, both GP-aided results exhibit a significant decrease in the relative log-Bayes factor, suggesting that accounting for correlated noise is relevant for this target. In this configuration, the ATMO2020++ model including PH₃ is again favored over the version without it, with a relative log-Bayes factor of ln B = 319 - 94 = 225 (see Table 4).

Focusing on ATMO2020++, Fig. 3 shows both a significant shift and broadening of the posterior distributions between the classical and GP approaches for all atmospheric parameters. We obtain T_eff = ${496}_{-3}^{+5}$ K and log(g) = ${4.30}_{-0.02}^{+0.04}$ dex; consistent with predictions from evolutionary models at 0.9 and 3.0σ, respectively, using the values reported in Table A.2 (${490}_{-4}^{+3}$ K and ${4.16}_{-0.04}^{+0.03}$ dex; see Sect. 3.3 for a more detailed description). The metallicity shifts markedly, from a supersolar [M/H] > 0.30 dex in the classical approach to solar [M/H] = $-{0.02}_{-0.02}^{+0.03}$ dex with the GP model, probably linked to the associated decrease in log(g). From our inversion, we infer a radius of R = ${1.03}_{-0.02}^{+0.01}$ R_Jup (consistent at 3.3σ), which translates into an estimated luminosity of log(L/L_⊙)= -6.163 ± 0.002 dex when combining the datasets. This luminosity estimate is also consistent with the value of -6.162 ± 0.006 dex obtained by Zhang et al. 2025 with the same model.

The retrieved hyperparameters for each observation listed in Table A.4 are consistent between the two ATMO2020++ models. These hyperparameters are also consistent across the other model grids (see Table A.4 and Appendix D). Among the datasets, MIRI-LRS exhibits the largest fitted correlation amplitude, with log(a)_MIRI = 1.24 ± 0.02.

The top row of Fig. A.3 presents a zoomed-in image of the three retrieved covariance matrices, normalized to highlight their noise structures. The correlated noise pattern is particularly visible in the FLAMINGOS-2 matrix with correlated noise extending across multiple spectral channels. Its block-like structure arises from the spectral windows used in the analysis.

The middle row displays the corresponding mean radial profiles, overplotted with the autocorrelation function (ACF) of the residuals. Comparing these profiles allows us to make an assessment of whether the retrieved correlation lengths adequately capture the actual noise structure in the data, assuming the model provides a good fit and the residuals are therefore dominated by noise. We find that the retrieved correlation lengths generally match this noise structure, although the correlation length may be slightly underestimated for the MIRI-LRS observation. Spectrographs with the finest wavelength sampling (FLAMINGOS-2 and NIRSpec) exhibit the largest correlation lengths, often spanning multiple pixels. In contrast, the retrieved correlation length for MIRI-LRS is nearly equal to the pixel size (~0.02 μm), suggesting more localized noise correlations.

The bottom row of Fig. A.3 compares the power spectral density (PSD) of the residuals, the injected observational uncertainties, and GP realizations. At low spectral frequencies, the PSD of the residuals (in green) is well captured by the GP (in red), and white noise (in blue) is not sufficient to explain the residual noise. At higher frequencies, for both the FLAMINGOS-2 and NIRSpec spectrographs, the PSD of the residuals is below that of the injected uncertainties and the GP noise, which suggests that the adopted noise model is too conservative at very small spectral scales.

Using the GP-aided spectral analysis on the combined dataset, we adopt the following parameters for COCONUTS-2 b: T_eff = ${496}_{-3}^{+5}$ K, log(g) = ${4.30}_{-0.02}^{+0.04}$ dex, [M/H] = $-{0.02}_{-0.02}^{+0.03}$ dex, R = ${1.03}_{-0.02}^{+0.01}$ R_Jup, and log(L/L_⊙) = -6.166 ± 0.002 dex (from Sect. 3.3).

Fig. 2 Top left panel : interpolated pressure-temperature (P-T) profiles for models providing a P-T grid (excluding BT-Settl) in Sect. 3.1. Remaining panels : posteriors of key parameters. The vertical dashed lines indicate the boundaries of the respective model grids when encountered during the inversion.

Fig. 3 Corner plot comparing forward modeling using the ATMO2020++ grid with (dark green) and without (light green) Gaussian processes (GPs). The black points indicate the mean predictions from evolutionary models with their associated error bars (see Sect. 3.3).

Fig. 4 Same as Fig. A.1, but for a comparison between the classical and GP-aided approaches on the JWST/MIRI-LRS spectrum using ATMO2020++ (in green) and ATMO2020++ (no PH3, in blue). Top panel : main molecular absorption features calculated using petitRADTRANS. In the other panels, the black lines show the observed spectrum; the colored lines show the corresponding model spectrum. The shaded regions around the model indicate the 1σ, 2σ, and 3σ dispersions from 200 random draws from the retrieved covariance matrix C. In the classical approach, the covariance matrix is assumed to be diagonal, C = ŝ2 diag(σ2), where ŝ2 is the noise scaling factor and σ the observational uncertainties. Zoom-in on the 5-11 μm region.

3.3 SED analysis

In Sects. 3.1 and 3.2, we constrain the bolometric luminosity of COCONUTS-2 b by propagating the posterior distributions of its effective temperature and radius using the Stefan-Boltzmann law (see the last column of Table A.3). This approach implicitly assumes that the bolometric flux can be accurately described by a single effective temperature. A more direct and less assumption-driven estimate can be obtained by directly integrating the full SED.

Because the wavelength coverage is not continuous, this SED must be completed. Previous studies (e.g., Filippazzo et al. 2015; Petrus et al. 2025; Kiman et al. 2026) have extrapolated the missing regions using Wien’s approximation at short wavelengths and the Rayleigh-Jeans tail at long wavelengths. This approach avoids introducing explicit atmospheric model assumptions (systematics) and is well suited for relatively warm L/T objects with smooth spectral energy distributions.

However, for very cold objects such as COCONUTS-2 b, whose emergent spectra are strongly shaped by molecular absorption, the assumption of smooth blackbody-like behavior outside the observed range may be less accurate. For this work we chose to use the posterior chain from our best-fit model (ATMO2020++ with GP; see Sect. 3.2) to generate a sample of model spectra over the widest possible wavelength range (0.230 μm). We then propagated this sample alongside the posterior distribution of the radius to derive a probability distribution for the bolometric luminosity. As previously, to account for potential model-data discrepancies, we also applied the corrective term given by Eq. (3) of Zhang et al. (2025), yielding

$${\text{log(L/L}}_{\odot })=-6.166\pm 0.002 \text{dex}.$$(6)

By repeating the procedure describe above on MIRI-LRS’s wavelength range, we can estimate the contribution of the MIRI-LRS observation to the total bolometric flux to be 41%. This value highlights the crucial role of the new dataset in constraining the bolometric luminosity of COCONUTS-2 b.

Finally, we used the ATMO2020++ (Phillips et al. 2020), Sonora Bobcat (Marley et al. 2021) and Sonora Diamondback hybrids (Morley et al. 2024) evolutionary models to rederive the physical properties of COCONUTS-2 b. We assumed a Gaussian distribution for the luminosity and the value we estimate (Eq. 6) and the Gaussian distribution for the age, based on the predictions from Kiman et al. (2026, 414 ± 23 Myr). The inferred parameters are summarized in Table A.2. While these parameters are fairly consistent across the models, there is an observable difference in radius and effective temperature between the ATMO2020 models and the Sonora Bobcat models (see Fig. A.6). After concatenating the chains for each parameter following the Monte Carlo sampling, evolution models predict that COCONUTS-2 b has T_eff = ${490}_{-4}^{+3}$ K, log(g) = ${4.16}_{-0.04}^{+0.03}$ dex, R = ${1.12}_{-0.01}^{+0.02}$ R_Jup, and M = 7.3 ± 0.3 M_Jup. These predictions are compared to those from the atmospheric models in Fig. 3.

4 Discussion

4.1 Atmosphere and formation of COCONUTS-2 b

In Sect. 3.1, we discuss the significant discrepancies in the atmospheric predictions from the different tested models. The cloud-free models (Sonora Elf Owl, ATMO2020++, and ATMO2020++ without PH₃) are significantly preferred based on Bayesian evidence, although the differences between them are also noticeable in the overall fit (see Fig. 1).

In particular, in the Y-band (FLAMINGOS-2) and N-band (MIRI), all models underestimate the observed flux. As suggested in Zhang et al. (2025), the discrepancies observed in the Y-band may be due to uncertainties in the modeling of the pressure-dependent Na and K lines. Although the preferred and most recent grids we use (ATMO2020++ and Sonora Elf Owl) include a self-consistent decrease in the abundance of these chemical species in the upper atmospheric layers through condensation and precipitation processes, we argue that increased precipitation or additional processes are needed to reduce the abundance of Na and/or K to match the observations. The flux underestimation in MIRI’s wavelength range is not seen in the GP inversions (Sect. 3.2, Fig. 4), suggesting that ATMO2020++ can reproduce this part of the spectrum reasonably well without invoking additional physical or chemical processes. This discrepancy between the classical and GP approaches arises from differences in the treatment of the noise as, in the classical approach, the noise marginalization effectively reduces the importance of the MIRI-LRS data during the inversion (see Sect. 4.2).

We pushed the atmospheric analysis using a GP-aided framework in Sect. 3.2 on the ATMO2020++ and ATMO2020++ (no PH₃) models, which resulted in a significant improvement in both the fit quality and consistency with evolutionary predictions. For both best-fit models, we retrieved metallicities consistent with the metallicity of the primary (0.00 ± 0.08 dex; Hojjatpanah et al. 2019) with [M/H] = $-{0.02}_{-0.02}^{+0.03}$ dex and [M/H] = -0.09 ± 0.02 dex for the models including and excluding PH₃, respectively. This general agreement suggests that COCONUTS-2 b has not experienced significant enrichment through core erosion or planetesimal accretion, and is instead compatible with a binary-like formation scenario. However, [M/H] is among the most degenerate parameters in our analysis analysis (Fig. 2 and Table A.3), with retrieved values spanning from supersolar to subsolar depending on the model employed. For instance, both Sonora models favor a more pronounced subsolar metallicity with [M/H] = -0.43 ± 0.01 dex and [M/H] = -0.30 ± 0.02 dex for Diamondback and Elf Owl, respectively, when using GP (see Table A.3). This could in turn suggest a different formation pathway (e.g., planetary, late accretion). Overall, this parameter is strongly coupled to the entire atmospheric composition, which means it is highly sensitive to the treatment of chemical processes, particularly non-equilibrium chemistry, which is nevertheless properly included in all the favored models.

4.2 Noise modeling

Accounting for all noise sources within a single inversion framework is a challenging task (e.g., Czekala et al. 2015; Gully-Santiago & Morley 2022). One of the simplest and least computationally expensive approaches is to include a noise scaling parameter in the Bayesian estimator, which can either be fitted directly (e.g., Zhang et al. 2025) or marginalized over (e.g., Ruffio et al. 2019). This scaling factor can be additive (e.g., Zhang et al. 2025) or, more commonly, multiplicative (e.g., Ruffio et al. 2019), and aims to compensate for missing and uncorrelated noise sources.

Marginalizing over this factor has been shown to significantly improve the robustness of retrieved posteriors in low (Landman et al. 2024; Denis et al. 2025) and in high S/N regimes (de Regt et al. 2024, 2025). However, our study demonstrates that such an approach is less suited to dealing with multi-modal datasets. Marginalizing over this scaling factor effectively equalizes the S/N contribution of each dataset. While MIRI-LRS intrinsically has the highest S/N and would therefore be expected to dominate the fit, the optimization process increases its noisescaling factor (see Table A.1), effectively reducing its weight and enhancing the contribution of the FLAMINGOS-2 and NIR-Spec spectra. This effect is clearly visible in Fig. 4: in the classical approach, the noise-scaling parameter allows the model to underestimate the MIRI-LRS flux, whereas this bias is not present in the GP approach, which does not marginalize over the noise-scaling factor.

For this study we used GPs to account for unmodeled spectral correlations, fitting for correlation lengths and amplitudes for each spectroscopic observations. We assessed the reliability of the retrieved correlation lengths from the ACF and PSD of the residuals (see Sect. 3.2 and Fig. A.3). In contrast, the correlation amplitudes depend on the overall scale of the residuals; their reliability is therefore best evaluated through injection tests as illustrated in Appendix C.

Accounting for correlated noise in the two ATMO2020++ models statistically improved both the fit quality and the consistency of most of the retrieved parameters (Table A.3). These correlations may arise from a combination of instrumental effects, and model systematics. In practice, the two contributions are difficult to disentangle from the residuals alone since the empirical ACF and PSD reflects both (residuals = observation - model). For the FLAMINGOS-2 and NIRSpec instruments, the residual autocorrelation and the GP reveal correlation lengths are both larger than (or on the order of) the instrumental line spread function (∼2-3 pixels, see middle row of Fig. A.3), indicating the presence of additional, more complex noise structures in the data that must be accounted for in the analysis.

More broadly, GP applications in exoplanet atmospheric characterization can extend beyond noise modeling. They have, for instance, been extensively used in transit spectroscopy to model systematic noise sources and correct for stellar contamination (e.g., Gibson et al. 2012). In the context of forward modeling, they can also be employed to remove local deviant spectra using the overall trends of the grids, therefore mitigating the effect of these spectra on the posterior distributions (Czekala et al. 2015). Moreover, GPs can propagate interpolation uncertainties into the final posteriors, as demonstrated in Czekala et al. (2015) with their Starfish⁷ framework. These last two aspects, while promising, were not explored for the present work. Future forward modeling efforts will benefit from incorporating these GP-based developments, to better characterize the noise structure and to more robustly disentangle observational systematics from model-driven discrepancies.

4.3 PH₃ analysis

Phosphine (PH₃) was expected to be the primary phosphorous molecule in the low-temperature atmospheres of brown dwarfs and giant exoplanets (Fegley & Lodders 1996; Visscher et al. 2006). However, most observational searches for PH₃ in these planetary-mass objects (e.g., Morley et al. 2018; Wallack et al. 2024) have provided only upper limits that are ∼100 times lower than the abundances predicted by atmosphere models (Beiler et al. 2024). In their paper Beiler et al. 2024 discussed the atmospheric mechanisms that could explain the observed underabundance of PH₃, including a vertical eddy diffusion coefficient that varies with altitude, incorrect chemical pathways, elements condensing out in forms such as NH₄H₂PO₄, or incorrect quenching approximations.

Throughout our analysis, the ATMO2020++ model including PH₃ remained statistically preferred compared to the same model without this molecule. The strongest features of PH₃ are located between 4 and 4.5 μm with a deep broad line at 4.3 μm. This deep 4.3 μm feature is clearly visible on the final fit (second and third panels from the top of Fig. A.2), but does not appear in the data at this location, clearly suggesting that PH₃ is indeed absent from the observation. On the other hand, the model without PH₃ manages to fit this region (fourth and fifth panels from the top of Fig. A.2), but consistently overestimate the flux between 4 μm and 4.2 μm (with or without GP) where CH₄ lines are present.

To evaluate whether this broader 4-4.2 μm region could be driving the preference for the model with PH₃ , we masked it and reran the inversions. In this case, ATMO2020++ without PH₃ was favored with ln B = 724 using the classical approach and ln B = 54 with the GP-aided approach. These results clearly indicate that the putative detection of PH₃ absorptions is highly sensitive to wavelength coverage and model systematics. Neither model successfully reproduces the NIRSpec data in this region implying that robust constraints on PH₃ detection are difficult to obtain with this approach.

To complement this analysis, we performed a crosscorrelation exploration of the NIRSpec spectrum using molecular templates. Templates for H₂O, CO, CH₄, CO₂, NH₃, and PH₃ were generated with petitRADTRANS⁸ (Mollière et al. 2019; Blain et al. 2024) using the median temperature and chemical profiles from the ATMO2020++ (GP) inversion (see Fig. A.5). We used the line-by-line opacity computation mode to generate spectra with only one molecule at a time. Both the observed spectrum and the templates were continuum-removed and normalized, and the template spectral resolution was matched to that of the observations. To avoid the 3.6-3.8 μm gap, we performed the cross-correlation on the two NIRSpec detectors separately. The resulting cross-correlation functions (CCF) are shown in Fig. A.7. H₂O, CO, CH₄, CO₂, and NH₃ are detected at significance levels above 4σ, while no significant signal is recovered for PH₃, further supporting its absence (or reduced abundance) in the atmosphere of COCONUTS-2 b.

5 Conclusion

In this work, we investigated the T9 planetary-mass companion COCONUTS-2 b. Using new MIRI-LRS observations in combination with existing spectro-photometric data within a self-consistent GP-aided framework, we refined the atmospheric properties of this object. From the GP-aided model analysis, we adopted the following parameters: T_eff = ${496}_{-3}^{+5}$ K, log(g) = ${4.30}_{-0.02}^{+0.04}$ dex, [M/H] = $-{0.02}_{-0.02}^{+0.03}$ dex, and R = ${1.03}_{-0.02}^{+0.01}$ R_Jup. These values are consistent with previous results (Zhang et al. 2021a, 2025; Kiman et al. 2026) and with the updated predictions from evolutionary models. The adopted metallicity, [M/H] = $-{0.02}_{-0.02}^{+0.03}$ dex, is consistent with that of the primary, supporting a binary-like formation scenario for COCONUTS-2 b. The wavelength extension provided by the MIRI-LRS observations contribute 41% of the bolometric flux. It enables a robust new constraint on the luminosity of COCONUTS-2 b: log(L/L_⊙) = -6.166 ± 0.002 dex.

However, systematic discrepancies between the data and atmospheric models in specific bands (Y and N), as well as inconsistencies in the retrieved metallicities and surface gravities across different model grids, remain the main limitations to providing a fully robust characterization of the object. These issues point to missing or incomplete chemistry in the current forward models, and highlight the need to incorporate improved alkali condensation and rainout, potentially more complex cloud structures, and more generally, diabatic processes. Such processes are increasingly being recognized as essential for interpreting the cold brown dwarf population revealed by JWST (e.g., Faherty et al. 2024; Alejandro Merchan et al. 2025; Kiman et al. 2026; Biller et al., in prep.). Accounting for them self-consistently will thus be equally important in preparation for the population of temperate near-water-ice-line planets expected from Gaia DR4 and the ELT, and, in the longer term, for extending these modeling efforts toward exo-Earth analogs.

We also find that classical noise-scaling schemes are ill-suited for heterogeneous multi-modal datasets with varying resolutions and S/N values. In contrast, parametric covariance models provide a robust alternative for small to moderate-size observation samples.

The aim of future retrieval analyses will be to constrain individual molecular abundances and explore more complex cloud prescriptions (blueCopeland et al., in prep.; Kühnle et al., subm.) using the data presented here and with the upcoming MIRI-MRS observations, providing further insight into the formation history of this object. Additionally, the high S/N MIRI-LRS spectrum will enable constraints on the isotope ratios (e.g., D/H, ¹²C/¹³C).

Acknowledgements

This work is heavily supported by a large collaboration of international people: Atmospheric modeling groups represented by P. Tremblin and M. W. Phillips (ATMO), F. Allard (deceased) (BT-Settl), J. J. Mang and C. V. Morley (Sonora). Collaborators who provided us with the data are represented by the M. Bonnefoy (MIRI-LRS), J. K. Faherty (NIRSpec), Z. Zhang (FLAMINGOS-2). This project has received funding from Agence Nationale de la Recherche (ANR) under grant ANR-23-CE31-0006-01 (MIRAGES). This project was provided with computing HPC and storage resources by GENCI at TGCC thanks to the grant 2024-15722 and 2025-15722 on the supercomputer Joliot Curie’s SKL and ROME partition. This work is based [in part] on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with GTO program 2124 and 3514. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. G.-D.M. acknowledges the partial support of the Deutsche Forschungsgemeinschaft (DFG) through grant “MA 9185/2-1”. This publication is based upon work from COST Action CA22133 “PLANETS” (https://costactionplanets.github.io), supported by COST (European Cooperation in Science and Technology).

References

Appendix A Additional figures and tables

Fig. A.1 Same as Fig. 4, but with Gemini/FLAMINGOS-2. Zoom-in on the Y- and J-bands.

Appendix B Models

We used five atmospheric grids with different physical and chemical makeup :

• BT-Settl (Allard et al. 2012) is a grid that explores complex cloud microphysics and includes non-equilibrium chemistry as well as vertical mixing. Clouds are simulated by dividing the atmosphere into multiple layers and computing the distribution and size of grains by comparing their characteristic times of condensation, coalescence, dispersion and sedimentation.

• Sonora Diamondback (Morley et al. 2024) model grid includes vertical mixing and parameterized clouds. Clouds are parameterized using Ackerman & Marley (2001) approach with the sedimentation efficiency parameter f_sed. This model assumes radiative-convective and chemical equilibrium. We used a custom version of this grid extended at lower temperature by rerunning the full forward model.

• Sonora Elf Owl (Mukherjee et al. 2024) is a cloudless grid that includes vertical mixing induced disequilibrium chemistry with subsolar to supersolar [M/H] and C/O. The atmospheric models have been computed using the open-source radiative-convective equilibrium model PICASO⁹ (Batalha et al. 2019; Marley et al. 2021; Mukherjee et al. 2023; JWST TECERST 2023). We used version v.2¹⁰ of the grid with a correction to the disequilibrium abundance of CO₂ and no PH₃ (Wogan et al. 2025).

• ATMO2020++ and ATMO2020++ (no PH₃) (Phillips et al. 2020; Leggett et al. 2021; Meisner et al. 2023) are extensions of the ATMO2020 model which additionally incorporate a non-adiabatic thermal structure of the atmospheres. These models are cloud-free but uses non-equilibrium chemical reactions of CO-CH₄ and N₂-NH₃, driving fingering convection in the atmosphere and changing the temperature gradient, to emulate their reddening effect.

Fig. A.2 Same as Fig. 4, but with JWST/NIRSpec. Zoom-in on the 3.84.7 μm region.

Fig. A.3 GP diagnostics. Top row: Zoom-in on the retrieved (normalized) covariance matrices from the GP-aided inversion using the ATMO2020++ model. The block-diagonal structure reflects the discontinuous spectral windows of the Gemini/FLAMINGOS-2 observations (see Sect. 3). Middle row: Autocorrelation of the residuals (data - model) in green, with the corresponding GP mean radial profile overplotted as a dotted green line. The vertical black lines indicate the GP correlation length-scales. Bottom row: Power spectral densities (PSD) of the residuals compared with random realizations from the injected instrumental noise (blue, ±1σ shaded) and from the fitted GP covariance (red, ±1σ shaded).

Appendix C GP injection test

This section describes the injection test we performed to probe the performance of the newly implemented GP.

To ensure that the correlation amplitude and length scale retrieved during our inversion are consistent with the correlated noise present in real spectra, we constructed a synthetic Gemini/FLAMINGOS-2 spectrum. This mock spectrum was interpolated from the ATMO2020++ grid using atmospheric parameters of T_eff = 483 K, log(g) = 4.19, and [M/H] = 0.0 (see Table 1).

The model spectrum was convolved and sampled at a spectral resolution of R_λ ~ 900 to match that of FLAMINGOS-2, resulting in a noise-free spectrum d. We then simulated observational noise using two components, as described by

$${\mathbf{d}}_{\mathbf{\text{noisy}}}=\mathbf{d}+{\mathbf{n}}_{\mathbf{\text{sky}}}+{\mathbf{n}}_{\mathbf{\text{inst}}},$$(C.1)

where n_sky and n_inst represent the sky and instrumental noise contributions, respectively, and d_noisy is the final synthetic noisy spectrum.

We model the sky noise as uncorrelated Gaussian noise, assuming an average S/N of 10

$${\mathbf{n}}_{\mathbf{\text{sky}}}\sim \mathcal{N}(0,{\sigma }^2),\text{with}{\sigma }^2={\left(\frac{\overline{d}}{\mathrm{S}/\mathrm{N}}\right)}^2,$$(C.2)

Fig. A.4 Injection test summary figure. Top panel: Comparison between the mock spectrum (black) and the fitted GP-aided model (red). Top subpanel: Corner plot of the grid parameters from the inversions with (red) and without (blue) GP. Green lines indicate the injected values. Bottom panel: Normalized residuals between the mock data and the fitted GP-aided model (black). Shaded regions represent the 1σ, 2σ, and 3σ dispersions from 200 random draws from the covariance matrix C.

where d̄ = mean(d). This component is the only one considered in the injected error bars.

To simulate the effect of correlated instrumental noise, we introduce a second component modeled with a GP characterized by a correlation length l = 10⁻² μm and amplitude a = σ (i.e., strong correlated noise). This is represented by

$${\mathbf{n}}_{\mathbf{\text{inst}}}\sim \mathcal{N}(0,\mathbf{C}),\text{with}{C}_{ij}=a^2\exp [-\frac{\mathrm{\Delta }{\lambda }_{i,j}^2}{2{\ell }^2}],$$(C.3)

where ∆λ_i,j denotes the wavelength separation between data points i and j. Finally, we inverted this mock data with the ATMO2020++ using 1000 living points and uniform priors, with (in red) and without GP (in blue). Results are summarized in Fig. A.4 and Table A.5.

The top panel of Fig. A.4 displays the fitted spectrum alongside the corner plots for both inversions. When spectral covariances are omitted, ForMoSA fails to recover the injected values, and the posterior uncertainties are noticeably underestimated across all parameters. In contrast, incorporating a GP model improves performance: although some bias remains, the true values lie within the broader posterior distributions. Additionally, the GP-aided inversion successfully retrieves the covariance hyperparameters, yielding final estimates of a = 1.0 ± 0.1 (true value: 1) and l = ${10.7}_{-0.1}^{+0.7}$ nm (true value: 10 nm). The bottom panel of Fig. A.4 presents the normalized residuals between the data and the model for the GP-aided inversion, clearly revealing correlated structures (in black). The red shaded contours represent the 1σ, 2σ, and 3σ dispersions of an ensemble of 200 random draws from the final (fitted) covariance matrix C. The inclusion of the GP kernel effectively captures both the structure and amplitude of the residuals.

Appendix D GP with additional atmospheric models

In Sect. 3.2, we focused our analysis on the two ATMO2020++ grids. Here, we provide an overview of the results obtained using the GP framework with the other grids (namely BT-Settl, Sonora Diamondback, and Sonora Elf Owl).

Fig. A.5 Key species mass fractions interpolated from ATMO2020++. The dashed lines correspond to the classical inversion, while the solid lines correspond to the GP-aided inversion.

According to all statistical criteria adopted in this work, Sonora Elf Owl provides the preferred fit among the tested models (see Table 4). Nevertheless, the overall quality of the fit remains poor, particularly in its ability to reproduce the MIRI observations. In addition, the retrieved atmospheric parameters, most notably T_eff (${549}_{-3}^{+2}$ K) and log(g) (< 3.25 dex), remain inconsistent with the predictions of all evolutionary models explored in this study (see Table A.2).

On the other hand, the GP significantly improves the consistency of the retrieved T_eff and log(g) for BT-Settl and Sonora Diamondback, with T_eff = 463 ± 3 K and log(g) = ${4.48}_{-0.02}^{+0.01}$ dex for BT-Settl, and T_eff = 515 ± 1 K and log(g) = ${4.500}_{-0.004}^{+0.003}$ dex for Sonora Diamondback.

Similar to ATMO2020++, the metallicity is found to be subsolar for all grids that explore this parameter, with [M/H] = -0.43 ± 0.01 dex and [M/H] = -0.30 ± 0.02 dex for Sonora Diamondback and Sonora Elf Owl, respectively. Contrary to the values retrieved with ATMO2020++ in this setup (see Sect. 4.1), these metallicities are not consistent with the stellar metallicity of 0.00 ± 0.08 dex (Hojjatpanah et al. 2019). Similarly, the C/O ratio remains subsolar, with C/O = 0.42 ± 0.01 retrieved using Sonora Elf Owl.

The vertical mixing inferred with Sonora Elf Owl (K_zz = $({3.2}_{-0.5}^{+0.6})$ × 10³ cm² s⁻¹) remains significantly lower than that predicted by the two ATMO2020++ grids (K_zz = (${2.5}_{-0.4}^{+0.2}$) × 10⁶ cm² s⁻¹ and K_zz = (${5.3}_{-0.5}^{+0.5}$) × 10⁶ cm² s⁻¹, with and without PH₃ , respectively).

With the exception of the BT-Settl grid, the inferred GP correlation lengths are remarkably consistent across the different atmospheric models (see Table A.4). In contrast, we observe larger variations in the GP correlation amplitudes. This trend suggests that the correlation lengths may be primarily driven by each instrument’s correlated noise pattern (see Sect. 4.2), while the correlation amplitudes may instead more closely reflect the GP compensating for each specific model mismatch.

Fig. A.6 Corner plot comparing the three families of evolutionary models (ATMO2020 in green, Sonora Bobcat in blue, and Sonora Diamondback in purple). The secondary bumps in the Sonora distributions usually corresponds to the different [M/H] nodes.

Fig. A.7 Cross-correlation functions of COCONUTS-2 b NIRSpec spectra (at full resolution) with H2O, CO, CH4, CO2, NH3, and PH3 molecular templates generated using petitRADTRANS. The S/N was computed using formula (1) of Houllé et al. (2021). Upper panel: CCFs obtained using the NRS1 detector where H2O, CH4, and NH3 are detected. Lower panel: CCFs obtained using the NRS2 detector where H2O, CO, CH4, and CO2 are detected.

All Tables

All Figures

Fig. 1 Top panel : forward modeling results of the combined COCONUTS-2 b observations. The black solid lines represent the spectroscopic data (from left to right: Gemini/FLAMINGOS-2, JWST/NIRSpec, and JWST/MIRI-LRS), while the colored lines represents each Rλ ∼ 100 model (classic). The pink and white stars indicate the photometric observations (WISE and Spitzer). The gray shaded regions indicate the masked wavelengths during the spectral inversion. Bottom panel : residuals for each fit. The dotted lines represents the ±σ (68%) confidence interval.

In the text

	Fig. 2 Top left panel : interpolated pressure-temperature (P-T) profiles for models providing a P-T grid (excluding BT-Settl) in Sect. 3.1. Remaining panels : posteriors of key parameters. The vertical dashed lines indicate the boundaries of the respective model grids when encountered during the inversion.
In the text

	Fig. 3 Corner plot comparing forward modeling using the ATMO2020++ grid with (dark green) and without (light green) Gaussian processes (GPs). The black points indicate the mean predictions from evolutionary models with their associated error bars (see Sect. 3.3).
In the text

Fig. 4 Same as Fig. A.1, but for a comparison between the classical and GP-aided approaches on the JWST/MIRI-LRS spectrum using ATMO2020++ (in green) and ATMO2020++ (no PH3, in blue). Top panel : main molecular absorption features calculated using petitRADTRANS. In the other panels, the black lines show the observed spectrum; the colored lines show the corresponding model spectrum. The shaded regions around the model indicate the 1σ, 2σ, and 3σ dispersions from 200 random draws from the retrieved covariance matrix C. In the classical approach, the covariance matrix is assumed to be diagonal, C = ŝ2 diag(σ2), where ŝ2 is the noise scaling factor and σ the observational uncertainties. Zoom-in on the 5-11 μm region.

In the text

	Fig. A.1 Same as Fig. 4, but with Gemini/FLAMINGOS-2. Zoom-in on the Y- and J-bands.
In the text

	Fig. A.2 Same as Fig. 4, but with JWST/NIRSpec. Zoom-in on the 3.84.7 μm region.
In the text

Fig. A.3 GP diagnostics. Top row: Zoom-in on the retrieved (normalized) covariance matrices from the GP-aided inversion using the ATMO2020++ model. The block-diagonal structure reflects the discontinuous spectral windows of the Gemini/FLAMINGOS-2 observations (see Sect. 3). Middle row: Autocorrelation of the residuals (data - model) in green, with the corresponding GP mean radial profile overplotted as a dotted green line. The vertical black lines indicate the GP correlation length-scales. Bottom row: Power spectral densities (PSD) of the residuals compared with random realizations from the injected instrumental noise (blue, ±1σ shaded) and from the fitted GP covariance (red, ±1σ shaded).

In the text

Fig. A.4 Injection test summary figure. Top panel: Comparison between the mock spectrum (black) and the fitted GP-aided model (red). Top subpanel: Corner plot of the grid parameters from the inversions with (red) and without (blue) GP. Green lines indicate the injected values. Bottom panel: Normalized residuals between the mock data and the fitted GP-aided model (black). Shaded regions represent the 1σ, 2σ, and 3σ dispersions from 200 random draws from the covariance matrix C.

In the text

	Fig. A.5 Key species mass fractions interpolated from ATMO2020++. The dashed lines correspond to the classical inversion, while the solid lines correspond to the GP-aided inversion.
In the text

	Fig. A.6 Corner plot comparing the three families of evolutionary models (ATMO2020 in green, Sonora Bobcat in blue, and Sonora Diamondback in purple). The secondary bumps in the Sonora distributions usually corresponds to the different [M/H] nodes.
In the text

Fig. A.7 Cross-correlation functions of COCONUTS-2 b NIRSpec spectra (at full resolution) with H2O, CO, CH4, CO2, NH3, and PH3 molecular templates generated using petitRADTRANS. The S/N was computed using formula (1) of Houllé et al. (2021). Upper panel: CCFs obtained using the NRS1 detector where H2O, CH4, and NH3 are detected. Lower panel: CCFs obtained using the NRS2 detector where H2O, CO, CH4, and CO2 are detected.

In the text