EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 689 (September 2024) A&A, 689 (2024) A212 Full HTML
Open Access
Issue
A&A
Volume 689, September 2024
Article Number A212
Number of page(s) 21
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/202450028
Published online 16 September 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Young brown dwarfs (BDs) and gas giant exoplanets share atmospheric and spectral characteristics, making BDs excellent analogues for studying exoplanet atmospheres (Chauvin et al. 2005; Currie et al. 2014; Faherty et al. 2016). These young, self-luminous objects, characterised by high effective temperatures and inflated radii, serve as benchmarks for understanding the complex processes shaping planetary atmospheres (Marley & Robinson 2015; Madhusudhan 2019). The advent of high-resolution spectroscopy (HRS) has transformed our ability to analyse these atmospheres, revealing details about their chemical compositions, thermal structures, and fundamental properties, such as rotational velocities and surface gravities (Snellen et al. 2014; Ruffio et al. 2021).

Recent advancements in atmospheric retrievals have enabled the quantitative characterisation of the atmospheres of young BDs and directly imaged exoplanets (Mollière et al. 2020; Zhang et al. 2021a, 2023; Wang et al. 2022; Landman et al. 2024; Inglis et al. 2024). Observations from state-of-the-art ground-based instruments like CRIRES+ (Dorn et al. 2023) and KPIC (Delorme et al. 2021), along with space-based facilities such as JWST (Rigby et al. 2023), have required more sophisticated models and retrieval frameworks to incorporate all physical and chemical processes in the atmosphere, in addition to the correct treatment of instrumental effects (e.g. Ruffio et al. 2024; de Regt et al. 2024).

Knowing elemental abundances, such as the carbon-to-oxygen ratio (C/O), allows us to trace the formation mechanisms and origins of exoplanets (Öberg et al. 2011). The local conditions within a protoplanetary disk determine the material available for planet formation and the chemical composition of the resulting objects. As the temperature decreases with distance from the star, volatile species such as water and CO freeze out, leading to chemical segregation within the disk (Lee et al. 2024). Consequently, the C/O in resulting planets is linked to their birth location within the disk. Previous studies have shown distinct C/Os in exoplanets, ranging from C/O ≈0.1 to C/Os close to unity for planets formed beyond the snow lines of water and CO (e.g. Öberg et al. 2011; Madhusudhan et al. 2017; Brewer et al. 2017; Line et al. 2021). Similarly, measurements of nearby low-mass stars and BDs indicate C/Os varying from ≈0.3 to ≈0.8 (Nissen 2013; Brewer & Fischer 2016). Precise measurements of C/O in stars are challenging due to the lack of molecular features in their spectra; the C/Os are typically inferred from the elemental abundances of C and O. Cooler objects, such as BDs and exoplanets, exhibit rich molecular spectra, allowing for a direct measurement of the C/O from molecular features in their spectra (Nissen 2013). However, challenges arise for ultracool dwarfs due to the limitations of existing atmospheric model grids in fitting their spectra and the degeneracies between the C/O and other parameters such as the temperature profile and the presence of clouds (Phillips et al. 2024).

Isotopic ratios, altered by chemical reactions and fraction-ation, may link the observed composition of objects to their formation environments. Objects that accrete material from the disk inherit its chemical composition, including isotopic ratios, potentially providing a link between the formation environment and the observed composition of the object (Mollière et al. 2022). Isolated BDs are thought to form through molecular cloud collapse and fragmentation, akin to star formation processes (Bate et al. 2002), or via disk fragmentation and ejection from extended disks around Sun-like stars (Stamatellos & Whitworth 2009). In contrast, massive gas giant planets, with masses ranging from 5 to 30 MJup, hereafter referred to as super-Jupiters (SJs), are believed to arise from combined processes of solid core accretion and subsequent gas capture (Pollack et al. 1996) or disk fragmentation (Mayer et al. 2002). Due to the uniqueness of these origins, chemical markers, such as the C/O and carbon isotopic ratios, may serve as indicators of their formation pathways. Precise measurements of these ratios in BDs and SJs are pivotal for understanding their formation mechanisms, offering insights into the broader phenomena of substellar and planetary formation (Nielsen et al. 2019; Vigan et al. 2021; Mollière et al. 2022; Zhang et al. 2021a).

In the Solar System, isotopic ratios serve as indispensable probes into the formation and evolutionary history of planets and smaller celestial bodies. Among these, the deuterium-to-hydrogen ratio (D/H) is particularly significant for unravelling the origins and historical dynamics of water (Aléon et al. 2022). This ratio is key to understanding various processes, including atmospheric loss and the potential roles of comets and asteroids in delivering water to planetary surfaces. The D/H isotope ratio varies across the Solar System, with a clear decreasing trend as a function of orbital distance: D/H ~1.6 × 10−2 on Venus (Donahue et al. 1982), ~1.5 × 10−4 on Earth (Hagemann et al. 1970), ~2.0 × 10−5 on the giant planets (Pierel et al. 2017), and ~4.0 × 10−5 on the ice giants (Feuchtgruber et al. 2013). Such measurements are crucial for piecing together the early environmental conditions and the subsequent development of habitable conditions on Earth-like planets (Goderis et al. 2016).

The carbon isotope ratio across the Solar System is mostly uniform, (12C/13C)⊙ ~ 89, with a terrestrial value of 89.3 ± 0.2 (Meija et al. 2016) and a solar value of 93.5 ± 3.1 (Lyons et al. 2018). In contrast, the local present-day interstellar medium (ISM) exhibits a lower value of (12C/13C)ISM = 68 ± 15 (Milam et al. 2005). The 12C/13C is a potential tracer of formation pathways, as it is altered by chemical reactions and fractionation processes during the formation stages (Mollière et al. 2022). Additional isotope ratios, such as 16O/18O,16O/17O, and 14N/15N, can provide complementary information and might be within the reach of current facilities (Gandhi et al. 2023; Mollière et al. 2022; Barrado et al. 2023).

The finding of a 13CO-rich atmosphere on the giant planet YSES 1b (Zhang et al. 2021a) sparked renewed interest in the 12C/13C as a formation tracer. A subsequent analysis of a young, isolated BD revealed a high 12C/13C (Zhang et al. 2021b), higher than that of the ISM and comparable to the solar value. Additional observations with CRIRES confirmed a high 12C/13C of 108 ± 10 in the young BD 2M0355 and hinted at the presence of C18O (Zhang et al. 2022). Recent HRS from Keck/KPIC has revealed solar-like values of 12C/13C in the HIP 55507 system and in the hot BD HD 984 B (Xuan et al. 2024; Costes et al. 2024). Space-based observations have also provided valuable measurements of the 12C/13C in the young, directly imaged planet VHS 1256 b (Gandhi et al. 2023) and a cool T dwarf (Hood et al. 2024). In addition to the 12C/13C, the 16O/18O has been tentatively constrained in the HIP 55507, pointing to homogeneous values in the system (Xuan et al. 2024). An analysis of JWST/NIRSPec M-band spectra (4.1-5.3 µm) of the young SJ VHS 1256 b resulted in a 12C/13C of 62 ± 2, in agreement with the value in the local ISM, within the uncertainties; the same analysis also found additional isotope ratios for oxygen, including 16O/18O and 16O/17O, that were slightly lower than the ISM value (Gandhi et al. 2023).

With the goal of testing the 12C/13C as a formation tracer, the ESO SupJup Survey (PI: Snellen) observed a large sample of isolated BDs and SJs with CRIRES+. An overview of the SupJup Survey is presented in de Regt et al. (2024), along with the initial results, which focused on the late-L dwarf DENIS J0255-4700. In this work we present the results of the atmospheric retrievals of three young BDs from the SupJup Survey, including precise chemical abundances of the main molecular and atomic spectrally active species, tight constraints on thermal profiles, and derived C/O and 12C/13C for each object.

Section 2 provides an overview of the sample of young BDs. Next, we describe the observations and data reduction (see Sect. 3). In Sects. 4 and 5, we present the forward modelling and retrieval framework. In Sect. 6, we report the results of the retrievals and discuss the implications of our findings in the context of young BDs and recent literature on isotope ratios. Finally, we summarise our conclusions in Sect. 8, highlighting the importance of HRS in characterising the atmospheres of young BDs and SJs and the potential of isotope ratios as tracers of their formation history.

2 Sample of young isolated brown dwarfs

Our sample includes three young (≲ 10 Myr) isolated BDs, with effective temperatures between 2300 and 2800 K and spectral types (SpTs) corresponding to late M-type objects (see Table 1). Due to ongoing gravitational contraction, these young BDs exhibit low surface gravity (log(𝑔) < 4.5) and low rotational velocities, typically υ sin i < 15 km s−1 (Stahler 1983; Baraffe et al. 2002; Cruz et al. 2009; Bouvier et al. 2014). The selection of these three objects is based on their shared properties with young SJs, including temperature, mass, and age (see Fig. 1). The detection of SJs with direct-imaging techniques is sensitive to self-luminous objects, making young and massive companions prime targets for detection and characterisation (e.g. McElwain et al. 2007; Vigan et al. 2008).

2.1 J1200

2MASS J12003792-7845082, hereafter referred to as J1200, is a late M-type BD with a circumstellar disk, as detailed by (Schutte et al. 2020). Located within the ε Cha association, J1200 is a member of this young (~3.7 Myr) stellar group. In their comprehensive analysis involving spectral energy distribution (SED) fitting and near-infrared spectroscopy, Schutte et al. (2020) establishes its SpT as M6γ with an effective temperature of approximately 2800 K and a mass ranging between 42 and 58 MJup. The assigned spectral subtype γ is indicative of low surface gravity, a commonly observed characteristic in young BDs.

Table 1

Fundamental parameters of the young BDs.

thumbnail Fig. 1

Colour-magnitude diagram of MJ vs. J − K. The objects observed in the SupJup Survey are indicated with three types of markers corresponding to isolated BDs, substellar companions, and hosts, which are either primary stars or BDs. The population of isolated cool dwarfs as catalogued in the UltracoolSheet1 is displayed as a reference. The colour of the points indicates the SpTs.

2.2 TWA 28

2MASS J11020983-3430355, also known as TWA 28, is an accreting substellar member of the dynamic TW Hydrae association (TWA), a young stellar group with an age ranging from 5 to 10 Myr (see Scholz et al. 2005 and Mamajek 2005). Analysis of medium-resolution near-infrared spectroscopy conducted by Venuti et al. (2019) point to an M9-type, corresponding to an effective temperature of Teff=2600 ± 70 K. The measured surface gravity log𝑔 = 4.1±0.3 is consistent with low-surface gravity, while the estimated rotational velocity is v sin i = 25 ± 10 km s−1. Additionally, the mass and radius are estimated via the accretion luminosity to 20.9 ± 5 MJup and 2.8 RJup, respectively. Low-resolution spectroscopy from Gaia Data Release 3 (Arenou et al. 2018) enabled the analysis of a large sample of ultracool dwarfs as reported in Cooper et al. (2024), where TWA 28 is identified as a low-gravity object with a SpT of M8.5γ and an effective temperature of 2382 ± 42 K, suggesting a lower temperature than previously measured.

Recent observations with the NIRSPec instrument aboard JWST, as presented in Manjavacas et al. (2024), compared the 0.97–5.3 µm spectra of TWA 28 to models of substellar atmospheres to constrain its atmospheric properties. The authors found that TWA 28 is best fit by models with effective temperatures of 2400–2600 K and a surface gravity of log(𝑔) = 4.0 ± 0.5. Additionally, the authors report the detection of excess flux in the >2.5 micron region, suggesting the presence of a disk around TWA 28, compatible with the accretion signatures found by Venuti et al. (2019).

2.3 J0856

2MASS J08561384-1342242, referred to as J0856, is a young (10 ± 3 Myr) M8γ BD with an estimated mass of 14.4 ± 1.4 MJup (Boucher et al. 2016). J0856 shows distinct low-gravity signatures and infrared excess flux, indicating the presence of a disk (Boucher et al. 2016). The properties of J0856 closely resemble those of TWA 28, hence providing a suitable comparison object.

Analysis of the SED by Morales-Gutiérrez et al. (2021) supports the presence of extended infrared excess emission. The authors modelled its SED with a central source of Teff ≈ 2500 K and a disk. The contribution from the disk becomes comparable to that of the BD photosphere at 4–5 µm and dominant at longer wavelengths.

In the study by Cooper et al. (2024), the parameters of J0856 parameters were refined to a SpT of M8.6β and an effective temperature of 2380 ± 32 K. The β spectral index indicates intermediate surface gravity, consistent with the object’s young age but not as low as the γ index.

3 Observations and data reduction

3.1 Observations

The observations presented in this work were conducted as part of the SupJup Survey (PI: Snellen, de Regt et al. 2024), employing the CRIRES+ instrument at the Very Large Telescope (VLT). CRIRES+, a state-of-the-art high-resolution slit-spectrograph equipped with adaptive optics, covers a broad wavelength range of the infrared spectrum (0.95–5.3 µm), at a resolving power R ~ 100 000 (Dorn et al. 2023). The data were collected March 3–4, 2023. Of the two observing nights, the first one experienced moderate seeing conditions (1.0–1.5″) and high humidity, while on the second night the seeing was stable (< 1.0″).

For this work, we opted for the K2166 wavelength setting, covering from 1.90 to 2.48 µm, specifically targeting regions rich in 12CO and, more specifically, 13CO absorption features. The observations employed the wide slit mode of 0.4″, achieving a resolving power of R ~ 50 000 as measured by the line shape of the telluric lines (see Appendix A). Each target’s total exposure time was calibrated to reach a S/N per pixel greater than 15 in the continuum (see Table 2). Observations were conducted in nodding mode, with individual exposures lasting 300s (total exposure time for each object outlined in Table 2), as this allows the correct subtraction of the sky thermal emission. The observations of telluric standard stars were interleaved with the science observations, with the goal of obtaining representative spectra of the Earth’s transmission. The telluric standards were observed at a similar airmass and with the same instrumental setup as the science targets.

Table 2

Observational properties of the young BDs.

3.2 Data reduction

The data reduction was performed with the open-source package excalibuhr2 (Zhang et al., in prep., also described in de Regt et al. 2024). The reduction started by tracing spectral orders on the flat-field images and determining the slit curvature via Fabry-Perot etalon frames, which provides an initial wavelength calibration. Subsequent steps included applying flat-fielding corrections, replacing bad pixels, fixing the slit curvature and tilt, and removing the sky emission through AB (or BA) subtraction. Frames from identical nodding positions were co-added into master frames for both A and B positions, which were then combined into a singular master frame. The conversion of calibrated images to one-dimensional spectra utilised the optimal extraction algorithm (Horne 1986). This procedure was similarly applied to standard star observations. The refinement of wavelength solutions for the extracted spectra was achieved through chi-square minimisation with telluric transmission templates provided by ESO’s Skycalc3 (Leschinski 2021).

The imprint of the Earth’s atmosphere on the observed spectra was corrected for with Molecfit (Smette et al. 2015) and observations of standard stars. Regions with saturated or deep telluric features were masked and ignored in the subsequent analysis (see Appendix A for more details). The telluric fit to the standard star also provides the instrumental throughput via a third-degree polynomial fit, which is used to correct the observed spectra for instrumental effects. The telluric-corrected spectra, containing 21 datasets (7 spectral orders with 3 detectors each), are divided by the mean flux of each order-detector pair, resulting in a normalised spectrum.

4 Forward modelling

To generate high-resolution spectra of BD atmospheres, we employed the petitRADTRANS code (Mollière et al. 2019). In this framework, the atmosphere is divided into a series of layers covering the pressure range from 102 to 10−5 bar, providing comprehensive coverage of the atmospheric vertical extent, including the region where spectral lines originate: the photosphere. Each layer is parameterised by a temperature and mass fractions of pertinent chemical species. At each layer, the radiative transfer equation is solved to obtain the outgoing flux as a function of wavelength. The resulting spectra are convolved with rotational and instrumental broadening (see Appendix A) kernels and resampled to the observed wavelength grid for direct comparison with the data.

Primary input for petitRADTRANS is the pressure-temperature profile, the opacities of spectrally active species, mass fractions of relevant chemical species and the surface gravity of the object. Additionally, the object’s radius can be utilised to scale the model to the observed flux.

4.1 Pressure-temperature profile

Self-luminous objects, devoid of significant external irradiation, exhibit a thermal profile expected to decrease with altitude. We parameterised the temperature as a function of pressure using values for the temperature gradient ∇Ti≡d ln Tid ln Pi$\nabla {T_i} \equiv {{{\rm{d}}\ln {T_i}} \over {{\rm{d}}\ln {P_i}}}$, where i = 1,2,3,4,5 corresponds to equally spaced points in pressure log-space. The temperature at each atmospheric layer was calculated as described in Zhang et al. (2023): T1=Tbottom at P1=102 barTj+1=exp(lnTj+(lnPj+1−lnPj)⋅∇Tj).$\eqalign{ & {T_1} = {T_{{\rm{bottom}}}}\quad {\rm{at}}\quad {P_1} = {10^2}\quad {\rm{bar}} \cr & {T_{j + 1}} = \exp \left( {\ln {T_j} + (\ln {P_{j + 1}} - \ln {P_j}) \cdot \nabla {T_j}} \right). \cr} $

Here, j = 1,…, 100 represents the atmospheric layers used in the radiative transfer calculation. The value of ∇Tj at each layer is determined with linear interpolation from the values at the knots ∇Ti. This approach facilitates the coupling of self-consistent radiative–convective-equilibrium (RCE) models with a ‘free’ parameterisation through custom prior probabilities. The uniform priors for ∇Ti· are set to replicate RCE thermal profiles (see Fig. 2). For this study, we considered the temperature profiles derived from state-of-the-art atmospheric models for M dwarfs (SPHINX; Iyer et al. 2023). The SPHINX profiles closely match those of other theoretical models in overlapping temperature ranges; however, the lack of high-temperature models of cool, substellar objects (Marley et al. 2021; Mukherjee et al. 2024) makes SPHINX the most suitable framework to derive physical constraints on the temperature profiles of the objects of interest in our study: late-M BDs. The prior probability distribution for each gradient ∇Ti and for the temperature at the bottom of the atmosphere is determined such that the resulting temperature profiles are consistent with the SPHINX models. In this way, the priors of the parameters are informed by RCE profiles, hence coupling the flexibility of the retrieval framework with the physical constraints of the models. Our parameterisation is expanded by including a free parameter that shifts the position of the intermediate pressure-temperature knots in the log-pressure space, denoted as log ∆P. This approach allows us to explore temperature profiles without biasing the photosphere’s location, which is expected to vary with the object’s surface gravity and effective temperature.

thumbnail Fig. 2

Thermal profiles (right) and temperature gradients (left) of the SPHINX model grid (https://zenodo.org/records/7416042) covering the range of temperatures between 2000 and 3000 K. The blue envelopes show the parameter space covered by the selected priors for the temperature gradients.

4.2 Chemical composition

The temperature range of 2000–3000 K straddles the boundary between planetary photospheres and those of low-mass stars. The spectra of objects within this temperature range are predominantly characterised by molecular lines (see Fig. 3), with some atomic lines as seen in stars from species such as Na, K, Ca, Fe, and Ti (Cushing et al. 2005).

In the K band (λ1.9–2.5 µm), the main spectral features arise from rotational–vibrational transitions of water and carbon monoxide. Dense water lines are pervasive across the observed wavelength range, while CO features are concentrated in the two reddest spectral orders of CRIRES+ (λ2.29–2.48 µm). Notable CO isotopologue overtones covered in our observations include 13CO(2,0) (λ2.345 µm) and 13CO(4,2) (λ2.403 µm) (Harrison & Marra 2017).

Opacities used in this work were calculated utilising the latest line lists from ExoMol (Tennyson et al. 2016) for H2O (Polyansky et al. 2017, 2018), HITEMP for the CO isotopo-logues (Rothman et al. 2010; Li et al. 2015), and for atomic opacities: Na, K (Allard et al. 2019), Mg, Ca, Ti, and Fe (Castelli & Kurucz 2003). Additionally, we included hydrogen fluoride (HF) as a line species (Wilzewski et al. 2016). Relevant continuum opacities encompass H2–H2, H2–He, H– (Dalgarno & Williams 1962; Chan & Dalgarno 1965; Gray 2022), and collision-induced absorption (Borysow et al. 1988).

Our retrieval framework treats the chemical abundance of each opacity source as an independent, free parameter. In contrast to equilibrium chemistry (EC), this approach avoids imposing physical constraints on the relative abundances among different species. This method enhances sensitivity to individual species and is less influenced by assumptions inherent in atmospheric models.

The number fractions, nX, of the main carbon- and oxygen-bearing species, CO and H2O respectively, are used to derive the atmospheric C/O for each object, defined as C/O=nCOnH2O+nCO=n12CO+n13COnH216O+nH218O+n12CO+n13CO.${\rm{C/O}} = {{{n_{{\rm{CO}}}}} \over {{n_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}} + {n_{{\rm{CO}}}}}} = {{{n_{12}}{\rm{CO}} + {n_{13}}{\rm{CO}}} \over {{n_{{{\rm{H}}_{\rm{2}}}^{16}{\rm{O}}}} + {n_{{{\rm{H}}_{\rm{2}}}^{18}{\rm{O}}}} + {n_{12}}{\rm{CO}} + {n_{13}}{\rm{CO}}}}.$(1)

We note that the C/O might change as a function of altitude if there are vertical gradients in the abundances of the species. In this work, we assumed constant-with-altitude abundances and, hence, a constant C/O. The atmospheric C/O might differ from the bulk C/O of the object, as cloud formation can result in the depletion of oxygen in the deep atmosphere (Line et al. 2015). At the hot temperatures of our objects, most condensates are expected to be in the gas phase; however, there are a few species, such as Al2O3 and Mg2SiO4, that can form condensates at temperatures above 2000 K (Helling & Casewell 2014; Woitke et al. 2020). The presence of clouds can affect the slope of the continuum and the shape of the temperature profile. Our K-band observations cover a small wavelength region and are not sensitive to the overall shape of the continuum (see Sect. 3.2); hence, we do not expect the presence of clouds to significantly affect our results.

As a proxy for metallicity, we employed the carbon abundance relative to hydrogen scaled to solar composition as [C/H]=log10(nCnH)−log10(nCnH)⊙,$[{\rm{C/H}}] = {\log _{10}}({{{n_{\rm{C}}}} \over {{n_{\rm{H}}}}}) - {\log _{10}}{({{{n_{\rm{C}}}} \over {{n_{\rm{H}}}}})_ \odot },$(2)

where the solar value is log10(nCnH)⊙=−3.57${\log _{10}}{\left( {{{{n_{\rm{C}}}} \over {{n_{\rm{H}}}}}} \right)_ \odot } = - 3.57$ (Asplund et al. 2021). In a similar manner, the fluorine abundance is calculated from HF as [F/H]=log10(nHFnH)−log10(nHFnH)⊙,$[{\rm{F/H}}] = {\log _{10}}\left( {{{{n_{{\rm{HF}}}}} \over {{n_{\rm{H}}}}}} \right) - {\log _{10}}{\left( {{{{n_{{\rm{HF}}}}} \over {{n_{\rm{H}}}}}} \right)_ \odot },$(3)

where the solar value is log10(nFnH)⊙=−7.6±0.25${\log _{10}}{\left( {{{{n_{\rm{F}}}} \over {{n_{\rm{H}}}}}} \right)_ \odot } = - 7.6 \pm 0.25$ (Maiorca et al 2014; Asplund et al. 2021).

To ascertain carbon isotope ratios, we retrieved the number fractions of 12CO and 13CO independently. The 12C/13C is then calculated as 12C/13C=n12COn13CO,$^{12}{\rm{C}}{/^{13}}{\rm{C}} = {{{n_{12}}{\rm{CO}}} \over {{n_{13}}{\rm{CO}}}},$(4)

similarly, the 16O/18O ratios are calculated as 16O/18O=nH2160nH2180.$^{16}{\rm{O}}{/^{18}}{\rm{O}} = {{{n_{{{\rm{H}}_2}{{16}_0}}}} \over {{n_{{{\rm{H}}_2}{{18}_0}}}}}.$(5)

To robustly confirm the presence of molecular species, we conducted a cross-correlation check following the methodology outlined in (Zhang et al. 2021a). The cross-correlation function (CCF) for each species X was calculated as CCFX(v)=∑i=1Nλ(di−mi,X˜)(mi−mi,X˜)σi2,${\rm{CC}}{{\rm{F}}_X}(v) = \sum\limits_{i = 1}^{{N_\lambda }} {{{({d_i} - {m_{i,\tilde X}})({m_i} - {m_{i,\tilde X}})} \over {\sigma _i^2}}} ,$(6)

where mi,X˜${m_{i,\tilde X}}$ is the model without the species X, and mi is the model containing all species included in the retrieval. The CCF is calculated for shifted versions of the model spectrum in velocity space covering υ = (−1000, 1000) km s−1 with steps of 1 km s−1. The uncertainties σi are calculated from the diagonal values of the best-fit covariance matrix, which accounts for the uncertainty of each data point and correlated noise among nearby pixels (see Sect. 5.2).

thumbnail Fig. 3

Opacity sources included in the forward modelling. The opacities are shown for a temperature of 2400 K and at the retrieved mixing ratios of TWA 28 (see Table 3).
thumbnail Fig. 4

Best-fitting models retrieved for the three objects. The spectra correspond to the three detectors of a single order of CRIRES+. This wavelength range contains several 12CO and 13CO lines. The observed data are shown in black, and the corresponding best-fit models are in blue (J1200), green (TWA 28), and red (J0856). Residuals are shown in the bottom panel. The average scaled uncertainty is indicated on the right side of each panel. The data and the models are displayed in the rest frame of each object. The full wavelength coverage used in the retrievals is shown in Appendix D.

4.3 Line profile

Several physical mechanisms contribute to the observed line spread function. Intrinsic rotation of the object induces Doppler-shifted lines across different regions of the observed disk, leading to an overall spectrum with broadened spectral lines.

To account for this, we utilised the fastRotBroad broadening kernel from PyAstronomy4 (Czesla et al. 2019). This kernel is employed to broaden our model lines based on the projected rotational velocity υ sin i, where i denotes the inclination of the system, and εlimb represents the linear limb-darkening coefficient. It is important to note that fastRotBroad assumes a constant-with-wavelength broadening kernel for computational efficiency, a valid approximation within our small wavelength range (Δλ ≈ 50 nm), showing a negligible difference (<0.5%) with direct disk integration (Carvalho & Johns-Krull 2023).

To correct for the shift in line positions induced by the radial velocity (υrad) of the object, we shift the forward-modelled spectra via linear interpolation to match the observed Doppler shift. Furthermore, to incorporate the instrumental broadening introduced by the slit spectrograph, we convolve the lines with a Gaussian kernel. The width of this kernel is determined by the full width at half maximum (FWHM), which corresponds to the spectral resolution of the instrument (R ~ 50 000; see Appendix A).

4.4 Veiling

Classical T Tauri stars often exhibit veiling of photospheric lines, characterised by shallower absorption features than expected for their SpT (Hartigan et al. 1989; Stempels & Piskunov 2003; Sullivan & Kraus 2022; Sousa et al. 2023). Veiling is attributed to the presence of a circumstellar disk, which can introduce an additional continuum source, frequently observed as an infrared excess (McClure et al. 2013).

In our analysis, veiling is modelled as an additional continuum component to the observed flux, with wavelength slope and amplitude as free parameters. This component is included in the forward model and retrieved simultaneously with other parameters. The veiling continuum is expected to be stronger for objects with accretion signatures. We do not observe any emission lines in our sample that would be indicative of active accretion, such as Brγ emission (Beck et al. 2010).

The veiling continuum is modelled as a power-law function of the form d(λ)=ϕ[mspec(λ)+mveil(λ)]${\bf{d}}(\lambda ) = \phi [{{\bf{m}}_{{\rm{spec}}}}(\lambda ) + {{\bf{m}}_{{\rm{veil}}}}(\lambda )]$(7) =ϕ[ mspec(λ)+r0(λλ0)α ],$ = \phi \left[ {{{\bf{m}}_{{\rm{spec}}}}(\lambda ) + {r_0}{{\left( {{\lambda \over {{\lambda _0}}}} \right)}^\alpha }} \right],$(8)

where ϕ is a normalising constant fitted for each order-detector pair (see Sect. 5.1), r0 is the veiling factor at the shortest wavelength of our dataset (λ = 1.90 µm), and α is the power-law index.

In this parameterisation, the wavelength-dependent K-band veiling factor, rk, is rk(λ)=r0(λλ0)α,${{\bf{r}}_{\bf{k}}}(\lambda ) = {r_0}{\left( {{\lambda \over {{\lambda _0}}}} \right)^\alpha },$(9)

which is a quantity often reported in the literature, with values ranging from 0 (no veiling) to of the order of 5 for very active T Tauri stars (Rei et al. 2018; Alcalá et al. 2021; Sousa et al. 2023).

5 Atmospheric retrieval framework

Atmospheric retrievals aim to solve the inverse problem of finding the best-fitting parameters and their uncertainties. We employed a Bayesian inference framework for our retrieval process, allowing us to quantify the uncertainties in the derived atmospheric parameters. We used the nested sampling algorithm (Skilling 2004), implemented through MultiNest (Feroz et al. 2009) with the Python wrapper PyMultiNest (Buchner 2016). MultiNest is adept at efficiently exploring and sampling from complex, high-dimensional parameter spaces, effectively managing parameter degeneracies. The sampling employs a set of live points drawn from the hyper-dimensional prior parameter space, defining an iso-likelihood contour. Through repeated likelihood evaluations (refer to Sect. 5.1), a fraction of live points is iteratively replaced to shrink the parameter space to regions of higher likelihood. We made use of the importance nested sampling technique in constant efficiency mode, which provides reasonably accurate evidence (𝒵) estimates at significantly higher efficiency than traditional nested sampling (Cameron & Pettitt 2014; Feroz et al. 2019). Our retrievals operate with 400 live points at a constant efficiency of 5%, with a tolerance of ∆ log𝒵 = 0.5, following the recommendations in Feroz et al. (2019). Additional retrievals with different numbers of live points and efficiencies were performed to ensure convergence and consistency of the results. We note that the results presented here are consistent with the number of live points varying from 200 to 1000 and an efficiency of 1–5% (see Dittmann 2024 for a detailed discussion on MultiNest’s hyperparameters).

5.1 Likelihood function

The likelihood function, as defined in Ruffio et al. (2019) and employed in Wang et al. (2022); Landman et al. (2024), is expressed for a data vector dik of shape (Ni, Nk) = (21, 2048), corresponding to the combination of 7 spectral orders and 3 detectors (totalling 21 order-detector pairs) with 2048 pixels each. A model Mψ,ik with non-linear parameters ψ and amplitudes ϕi is used. The log-likelihood is computed as follows: lnℒ=∑iln ℒi$\ln {\cal L} = \sum\limits_i {\ln } \,{{\cal L}_i}$(10) lnℒi=−12(Nkln(2π)+ln(|Σ0,i|)+Nkln(si2)+1si2riTΣ0,i−1ri),$\ln {{\cal L}_i} = - {1 \over 2}\left( {{N_k}\ln (2\pi ) + \ln (|{\Sigma _{0,i}}|) + {N_k}\ln (s_i^2) + {1 \over {s_i^2}}{\bf{r}}_i^T\Sigma _{0,i}^{ - 1}{{\bf{r}}_i}} \right),$(11)

where Σ0 is the covariance matrix, ri = dik − ϕimik represents the residuals for each order-detector pair, and si is the uncertainty scaling factor for the total covariance matrix Σi=si2Σ0,i${\Sigma _i} = s_i^2{\Sigma _{0,i}}$.

The optimal parameters ϕ˜i${\tilde \phi _i}$ and s˜i2$\tilde s_i^2$ are are calculated at every likelihood evaluation using a least square solver5 and the uncertainty scaling formula, respectively: ϕ˜i=(MΣ0−1M)−1MΣ0−1d${\tilde \phi _i} = {(M\Sigma _0^{ - 1}M)^{ - 1}}M\Sigma _0^{ - 1}{\bf{d}}$(12) s˜i2=1NkriTΣ0,i−1ri|ϕ=ϕ˜.$\tilde s_i^2 = {1 \over {{N_k}}}{\bf{r}}_i^T\Sigma _{0,i}^{ - 1}{{\bf{r}}_i}{|_{\phi = \tilde \phi }}.$(13)

5.2 Correlated noise

Uncertainties in the data play a crucial role in the retrieval process. Recent findings by de Regt et al. (2024) underscore the significance of correlated noise in CRIRES+ spectra on the retrieved parameters. When uncertainties are correlated among nearby pixels, the covariance matrix Σ0 is not diagonal. Gaussian processes (GPs) provide a powerful means to model such correlated noise (Kawahara et al. 2022). In this study, we employed a radial basis function kernel to model correlated noise in the data, defined as kij=a2σeff,ij2exp(−12(xi−xj)2l2),${k_{ij}} = {a^2}\sigma _{{\rm{eff}},ij}^2\exp \left( { - {1 \over 2}{{{{({x_i} - {x_j})}^2}} \over {{l^2}}}} \right),$(14)

where a represents the amplitude, σeff,ij2$\sigma _{{\rm{eff}},ij}^2$ the effective variance, calculated as the average variance per order-detector pair, l the length-scale in units of wavelength, and xi the wavelength of the i-th pixel. The resulting covariance matrix for each order-detector pair is the sum of the diagonal covariance matrix and the GP kernel: Σ0,ij=δijσi2+kij.${\Sigma _{0,ij}} = {\delta _{ij}}\sigma _i^2 + {k_{ij}}.$(15)

During retrieval, the amplitude a and length scale l of the GP kernel are treated as free parameters. The effective variance is derived from the data, and the diagonal covariance matrix corresponds to the uncertainty in the data. The GP kernel is computed at each likelihood evaluation and added to the covariance matrix.

6 Results

Our best-fitting models successfully reproduce the observed spectral features within the specified uncertainties as shown in Fig. 4 and quantified via the reduced chi-square (without the noise model, i.e. no uncertainty scaling) in Appendix C. The chemical abundances and temperature profiles for the three objects are broadly comparable, with variations attributed to differences in their effective temperatures, surface gravities, spin, and the S/N of our observations. A summary of the retrieved parameters is provided in Table 3 and a detailed discussion of the results is presented in the following sections.

thumbnail Fig. 5

(a) Best-fit pressure-temperature profiles. The shaded regions indicate the 1-, 2-, and 3-sigma regions. On the right axis we overplot the integrated contribution function. (b) Posterior distributions for C/O (top) and 12C/13C (bottom). The values for the solar C/O = 0.59 ± 0.08 (Asplund et al. 2021), and the isotope ratios for the Solar System, 89.3 ± 0.2 (Meija et al. 2016), and the ISM, (12C/13C)ISM= 68 ± 16 (Milam et al. 2005), are marked.
Table 3

Retrieved parameters and their uncertainties.

thumbnail Fig. 6

Posterior distributions of the surface gravity and carbon-to-oxygen ratio relative to the solar value for each object. The contours of the 2D histogram represent the 1, 2, and 3σ confidence intervals, with dashed lines indicating the median values. The Pearson correlation coefficient is shown in the top-right corner, where a value of r = 0 indicates no correlation and r = 1 indicates a perfect positive correlation.

6.1 Chemical composition

Our retrievals provide the best-fit volume-mixing ratios (VMRs) of spectrally active species in the atmospheres of the three objects (see Table 3). The derived C/Os according to Eq. (1) are C/OJ1200=0.62−0.02+0.02,C/OTWA28=0.61−0.02+0.02,C/OJ0856=0.58−0.01+0.01,$\eqalign{ & {\rm{C/}}{{\rm{O}}_{{\rm{J1200}}}} = 0.62_{ - 0.02}^{ + 0.02}, \cr & {\rm{C/}}{{\rm{O}}_{{\rm{TWA28}}}} = 0.61_{ - 0.02}^{ + 0.02}, \cr & {\rm{C/}}{{\rm{O}}_{{\rm{J0856}}}} = 0.58_{ - 0.01}^{ + 0.01}, \cr} $

where the uncertainties indicate the 1σ intervals. The C/Os are comparable among the three targets and close to the solar value of 0.59 ± 0.08 (Asplund et al. 2021; see Fig. 5b) and similar to the C/O of the young BD 2M0355 (0.56 ± 0.02; Zhang et al. 2021b).

The derived posterior distributions of [C/H] and surface gravity for each object are displayed in Fig. 6. Using the carbon-to-hydrogen ratio as a proxy for metallicity, we find that the retrieved metallicities are broadly consistent with the solar value. However, their interpretation is complicated by the high correlation between surface gravities and metallicities (see Sect. 7.4).

In addition to the main molecular species, we detect 13CO, HF, Na, and Ca in the atmospheres of the three objects. In TWA 28 and J0856, we retrieve the VMR of H218O and tentatively detect it in J1200 (see Sect. 7.3). We also report the presence of Ti in the two hottest targets, J1200 and TWA 28, and provide an upper limit for J0856 (see Fig. C.3).

To assess the contribution to the best-fit model and quantify the detection significance of the minor isotopologues, we ran a series of retrievals with the species removed from the model. The significance of the detection is determined via a Bayesian model comparison, where the Bayes factor, Bm (Kass & Raftery 1995) is calculated as the ratio of the marginal likelihoods of the models with and without the species lnBm=ln Zm−ln ZmX˜.$\ln {B_} = \ln {cal Z_} - \ln {cal Z_{{_{\tilde X}}}}.$(16)

The Bayes factor is then converted to a detection significance, with values greater than 3 indicating a significant detection (Benneke & Seager 2013).

We find that the detection of 13CO is strong for TWA 28 and J0856, with 3.9σ and 4.7σ significance, respectively, while the detection of H218O is moderate with 2.1σ and 3.3σ significance for TWA 28 and J0856, respectively. The best-fit model of J1200 exhibits weak evidence for the presence of 13CO with a 1.8σ significance. The best-fit model for J1200 with H218O is disfavoured at 1.6σ; however, the lowest reduced chi-squared (χr) corresponds to the model with both minor isotopologues (see Appendix C). Given the lack of strong evidence of H218O in J1200, we report the retrieved VMR as an upper limit.

The CCFs support the presence of the minor isotopologues in the atmospheres of the three objects (see Fig. 7). The S/Ns of the CCFs are consistent with the Bayesian significances for the strong detections, with a slightly higher detection significance in the cross-correlation analysis. The CCFs of J1200 show a moderate detection of 13CO and H218O at 3.4σ and 2.6σ level, respectively, higher than the Bayesian significance. We note that in the retrievals without the minor isotopologues, the atmospheric and noise parameters can be adjusted to compensate for the absence of the species; however, the retrieved parameters are consistent to within 1σ confidence across all retrievals of the same target. The best fit to the data for all three targets is achieved when the minor isotopologues are included in the model, as reflected in the lowest reduced chi-squared values (see Appendix C).

thumbnail Fig. 7

CCFs for 13CO and H218O. The CCFs are calculated for each order-detector pair and summed over all orders and detectors. The CCFs are converted to S/N by dividing them by the standard deviation of the CCFs away from the peak (|υrad| > 100 km s−1). The S/N of the peaks of the CCFs is shown in the legend. The residuals between the observed and modelled CCFs – the auto-correlation function (ACF) – are shown in the bottom panel.

6.2 Thermal profile

The best-fitting temperature profiles (see Fig. 5a) exhibit remarkable similarities among all three objects, consistent with their similar SpTs and effective temperatures. However, the photospheric region of each object shows slight variations reflecting their individual temperature profiles and surface gravities. The observed spectral lines originate from a region spanning 0.01 to 3 bar, as indicated by the integrated emission contribution function and the narrower uncertainties within this region (see Fig. 5a).

Moreover, the retrieved thermal profiles align well with literature values for the effective temperature of each object (see Table 1). The inclusion of a free parameter for the spacing of the temperature knots avoids potential biases on the location of the photosphere due to fixed pressure knots and provides a wider range of temperature profiles. The retrieved log ∆P values for TWA 28 and J0856 highlight the ability of this parameterisa-tion to capture the location of the photosphere accurately. In contrast, the retrieved log ∆P for J1200 is constrained around zero, indicating that even when the fixed pressure knots can already capture the photosphere’s location, the additional parameter offers a more accurate representation of the uncertainties in the temperature profile (see Table 3 and Fig. 5a).

6.3 Rotational velocity

For TWA 28, we report a lower value than the previous estimate of 25 ± 10 km s−1 (Venuti et al. 2019), υsini TWA28=11.6−0.1+0.1 km s−1.$\upsilon \sin {i_{{\rm{TWA}}28}} = 11.6_{ - 0.1}^{ + 0.1}{\rm{km}}\,{{\rm{s}}^{ - 1}}.$

TWA 28 is the fastest rotator of our sample, in agreement with the observed broadening of the spectral lines (see Fig. 4 and Appendix D). The retrieved υ sin i values for J0856 and J1200 are consistent with the expected slow rotation for young BDs (Vos et al. 2020): υsini J1200=5.3−0.2+0.2 km s−1,υsini J0856=4.5−0.1+0.2 km s−1,$\eqalign{ & \upsilon \sin {i_{{\rm{J1200}}}} = 5.3_{ - 0.2}^{ + 0.2}{\rm{km}}\,{{\rm{s}}^{ - 1}}, \cr & \upsilon \sin {i_{{\rm{J}}0856}} = 4.5_{ - 0.1}^{ + 0.2}{\rm{km}}\,{{\rm{s}}^{ - 1}}, \cr} $

The rotational velocity of low-mass objects is a key parameter in understanding their evolution (Bryan et al. 2020). Young objects (< 1000 Myr) are expected to exhibit slow rotation due to the conservation of angular momentum during the contraction phase (Bouvier et al. 2014). The observed values for υ sin i are are consistent with slow rotation; however, the inclination of the system is unknown, and the true rotational velocity may be higher for a face-on system.

7 Discussion

7.1 Hydrogen fluoride

Hydrogen fluoride, marked by its 1–0 vibrational transitions with prominent absorption lines in the 2.3–2.5 µm region (see Appendix B), is clearly detected in our observations. With precedence in solar observations (Hall & Noyes 1969) and various giant and dwarf stars (Wallace & Hinkle 1996), the inclusion of HF as a line species stems from its detection in the analysis of high S/N data from a nearby BD (de Regt et al., in prep.). The absolute abundances of fluorine with respect to hydrogen are calculated from the VMRs of HF (see Fig. 8b), [F/H] J1200=−0.18−0.22+0.20,${[{\rm{F/H}}]_{{\rm{J1200}}}} = - 0.18_{ - 0.22}^{ + 0.20},$(17) [F/H] TWA28=0.01−0.21+0.20,${[F/H]_{TWA28}} = 0.01_{ - 0.21}^{ + 0.20},$(18) [F/H] J0856=0.05−0.15+0.14,${[{\rm{F/H}}]_{{\rm{J0856}}}} = 0.05_{ - 0.15}^{ + 0.14},$(19)

which are in excellent agreement with the solar fluorine abundance (see Fig. 8b) as measured from a sunspot spectrum (Maiorca et al. 2014; Asplund et al. 2021).

thumbnail Fig. 8

(a) Posterior distributions of the oxygen isotope ratio 16O/18O for the three objects. The solar value of 525 ± 21 (Lyons et al. 2018) and the ISM value of 557 ± 30 (Wilson 1999) are marked. (b) Posterior distributions of the abundance of fluorine with respect to hydrogen normalised to the solar value. The solar fluorine abundance is overplotted as a Gaussian distribution with a 1σ scatter (Maiorca et al. 2014).

7.2 Chemical equilibrium

In this work we adopted an approach where the mixing ratios of line species were treated as independent parameters, a method we refer to as ‘free composition’. This approach does not enforce constraints on the relative abundances, assuming constant-with-altitude mixing ratios. Given the limited vertical extent probed by the observed spectral lines, this assumption is generally valid for many scenarios. Nonetheless, it is insightful to contrast our retrieved abundances against those predicted by EC. The EC profiles for most species display minimal variation across the temperature and pressure ranges of our subjects, except for titanium (Ti), which shows a notable decrease in concentration with rising temperatures (Fig. 9). Using FastChem6 (Kitzmann et al. 2024), we compared our retrieved mixing ratios against EC predictions, considering elemental solar abundances adjusted of the retrieved metallicity and scaling the abundance of12 CO to match the retrieved C/O. To adjust for any potential systematic offsets due to the metallicity-surface gravity degeneracy, we normalised the VMRs relative to VMR(12CO).

Our findings indicate that the vertical distributions of the primary molecular species (12CO and H2O) are relatively stable over the pressure-temperature range of the three objects. In general, the retrieved values align with the EC predictions to within a 1σ confidence level, with the significant exception of calcium (Ca). While EC predicts comparable mixing ratios for sodium (Na) and Ca, our analysis find the retrieved VMRs of Na to be consistent with EC, while the Ca abundances are retrieved at significantly higher values than EC. However, the contribution of Ca mostly originates from a small region at the bluest order of CRIRES+ (see Fig. 3), where telluric lines are prevalent (see Fig. A.1); hence, the interpretation of the Ca abundance is not straightforward. For Ti, the retrieved mixing ratios broadly align to EC in the atmosphere’s deepest regions while for HF the retrieved mixing ratios are in the range of 10−8 to 10−7, resulting in close agreement with EC.

thumbnail Fig. 9

Retrieved VMRs of line species compared to EC values of the retrieved temperature profiles of each target. The dashed lines indicate the values from EC at the retrieved C/O and metallicity. The shaded regions represent the retrieved 1σ confidence region around the best fit. The integrated contribution function is plotted in the left column and indicated in the right column as a dark shaded region.

7.3 Isotopic composition

7.3.1 Carbon isotope

The obtained 12C/13C for J1200, TWA 28, and J0856 (see Fig. 5b), 12C/13C J1200=114−33+69,12C/13C TWA28=81−19+28,12C/13C J0856=79−14+20,$\eqalign{ & ^{12}{\bf{C}}{/^{13}}{{\bf{C}}_{{\bf{J}}1200}} = 114_{ - 33}^{ + 69}, \cr & ^{12}{\rm{C}}{/^{13}}{{\rm{C}}_{{\rm{TWA28}}}} = 81_{ - 19}^{ + 28}, \cr & ^{12}{\bf{C}}{/^{13}}{{\bf{C}}_{{\rm{J0856}}}} = 79_{ - 14}^{ + 20}, \cr} $

are consistent among the three objects within 1 σ uncertainties. The retrieved values are in agreement with the ISM given the scatter around the present-day value (12C/13C)ISM = 68 ± 16 (Milam et al. 2005) and show no significant enrichment in 13CO. In the context of young substellar objects, our results are in good agreement with other studies that have reported 12C/13C values consistent with the ISM (Xuan et al. 2024; Gandhi et al. 2023). The retrieved 1σ uncertainties cover the solar value (12C/13C)⊙~ 89 (Meija et al. 2016) and partially overlap with objects of higher 12C/13C, ~100 (Zhang et al. 2021b; Costes et al. 2024; Hood et al. 2024). In contrast to the first object of the SupJup Survey: DENIS J0255 (12C/13C=184−40+61$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 184_{ - 40}^{ + 61}$; de Regt et al. 2024), the retrieved 12C/13C values are consistent with the ISM, indicating no significant enrichment or depletion of13 CO in the atmospheres of these objects. As noted in de Regt et al. (2024), the age of the objects may play a significant role in the isotopic composition, with younger objects exhibiting a lower 12C/13C due to the 13CO enrichment in the ISM over time (Milam et al. 2005; Romano et al. 2017).

7.3.2 Oxygen isotope

The detection of H218O in TWA 28 and J0856, and the marginal detection in J1200, provides additional insights into the iso-topic composition of these objects. From the retrieved H216O and H218O VMRs, we calculated the 16O/18O ratio for each object, 16O/18OJ1200>159,16O/18OTWA28=205−62+140,16O/18OJ0856=141−28+42,$\eqalign{ & ^{16}{\rm{O}}{/^{18}}{{\rm{O}}_{{\rm{J}}1200}} > 159, \cr & ^{16}{\rm{O}}{/^{18}}{{\rm{O}}_{{\rm{TWA}}28}} = 205_{ - 62}^{ + 140}, \cr & ^{16}{\rm{O}}{/^{18}}{{\rm{O}}_{{\rm{J}}0856}} = 141_{ - 28}^{ + 42}, \cr} $

where the lower limit for J1200 is due to its marginal detection. The retrieved 16O/18O ratios are lower than the solar value of 525 ± 21 (Lyons et al. 2018) and the ISM value of 557 ± 30 (Wilson 1999), indicating a potential enrichment of 18O in the atmospheres of these objects. The 16O/18O ratios reported for HIP 55507 B (16O/18O=288−70+125$^{16}{\rm{O}}{{\rm{/}}^{18}}{\rm{O = }}288_{ - 70}^{ + 125}$; Xuan et al. 2024) also exhibited a lower value than the ISM (i.e. a higher 18O abundance). Gandhi et al. (2023) showcased the excellent capabilities of JWST to measure the 16O/18O ratio in the atmospheres of young substellar objects, with a reported value of 16O/18O=425−28+33$^{16}{\rm{O}}{{\rm{/}}^{18}}{\rm{O = }}425_{ - 28}^{ + 33}$ for the young SJ VHS 1256 b. The 16O/18O retrieved in this work for J0856, with the caveat that H218O is only detected at 3.3σ, is one of the lowest values reported for young substellar objects, indicating a potential deviation of 18O in the atmosphere of this object from the ISM value.

7.4 Surface gravity

Surface gravity significantly influences the observed spectra by affecting line shapes, the temperature profile, and the overall continuum (Marley & Robinson 2015; Mukherjee et al. 2024). Its inverse relationship with atmospheric scale height impacts the pressure range of the photosphere, a key property determining spectral features. We report low surface gravities for the three objects (see Table 3), consistent with expectations for young objects (Baraffe et al. 2002; Allers et al. 2007; Bonnefoy et al. 2014).

The posterior distributions are highly correlated7 (see Fig. 6), making it difficult to precisely constrain surface gravities and metallicities. Determining surface gravity from K-band observations is challenging due to the scarcity of reliable gravity-sensitive features (Zhang et al. 2021b). This limitation could be mitigated with additional data at shorter (e.g. J- or H-band) wavelengths. Furthermore, incorporating prior knowledge of precise dynamical mass and radius measurements could help resolve existing degeneracies (Stolker et al. 2020).

7.5 Veiling

The best-fit models for TWA 28 and J0856 show no evidence of veiling in their spectra, that is, r0 ≈ 0.0. However, for J1200, the retrieved veiling factor is r0=0.14−0.06+0.06${r_0} = 0.14_{ - 0.06}^{ + 0.06}$, suggesting the presence of a weak veiling continuum. The wavelength dependence of the veiling factor is loosely constrained by the data, favouring a power-law index of α=0.74−0.44+0.58$\alpha = 0.74_{ - 0.44}^{ + 0.58}$, which increases towards redder wavelengths.

Veiling is a widely observed phenomenon in the near-infrared spectra of many young stellar objects, particularly accreting sources (Folha & Emerson 1999; McClure et al. 2013; Sullivan & Kraus 2022; Sousa et al. 2023). However, it is less common in substellar objects, with only a few detections in young BDs (e.g. White & Basri 2003). The physical origin of veiling is still debated, with potential sources including accretion, chromospheric activity, and magnetic fields (Hartigan et al. 1989; Folha & Emerson 1999; Stempels & Piskunov 2003; Rei et al. 2018).

There is a potential link between mass accretion rate and near-infrared veiling; more massive objects with higher accretion rates are expected to exhibit stronger veiling (White & Basri 2003; Sousa et al. 2023). The presence of veiling in J1200 may indicate ongoing accretion, which is plausible given its young age, the detection of a disk around this object (Schutte et al. 2020), and its higher mass compared to TWA 28 and J0856 (see Table 1).

7.6 Correlated noise

Our analysis revealed significant correlation in the data as reflected by the retrieved parameters of the GP kernel. We observe a tentative relation between the GP amplitude, a, and length scale, l, and the spin of the objects. Specifically, slow rotators, characterised by narrower spectral features, exhibit a smaller GP length scale and a larger amplitude compared to fast rotators. The retrieved values, ordered from the slowest to the fastest rotator, are presented below: J0856 (υsini≈4.5 km s−1)     a=2.15−0.02+0.02     l=0.0278−0.0004+0.0001≈3.6 pixelsJ1200(υsini≈5.3 km s−1)    a=1.91−0.02+0.02   l=0.0317−0.0010+0.0002≈4.1 pixelsTWA28 (υsini≈11.6 km s−1)   a=1.06−0.01+0.01  l=0.0410−0.0003+0.0002≈5.3 pixels.$\eqalign{ & J0856\quad (\upsilon \sin i \approx 4.5{\rm{km}}\,{{\rm{s}}^{ - 1}}) \cr & \,\,\,\,\,a = 2.15_{ - 0.02}^{ + 0.02} \cr & \,\,\,\,\,l = 0.0278_{ - 0.0004}^{ + 0.0001} \approx 3.6{\rm{pixels}} \cr & J1200(\upsilon \sin i \approx 5.3{\rm{km }}{{\rm{s}}^{ - 1}}) \cr & \,\,\,\,a = 1.91_{ - 0.02}^{ + 0.02} \cr & \,\,\,l = 0.0317_{ - 0.0010}^{ + 0.0002} \approx 4.1{\rm{pixels}} \cr & {\rm{TWA28}}\quad (\upsilon \sin i \approx 11.6{\rm{km}}\,{{\rm{s}}^{ - 1}}) \cr & \,\,\,a = 1.06_{ - 0.01}^{ + 0.01} \cr & \,\,l = 0.0410_{ - 0.0003}^{ + 0.0002} \approx 5.3{\rm{pixels}}. \cr} $

In this context, the GP amplitude serves as a dimensionless parameter adjusting the covariance matrix’s off-diagonal elements, while the length scale, expressed in nanometers (nm), reflects the correlation distance over which the correlation between two points decreases to 1/e of its initial value. This scale generally matches the width of the spectral lines within the observed spectra. Employing the CRIRES+sampling rate, the length scale is converted to a number of pixels. Accurately accounting for correlated noise is essential for obtaining reliable uncertainties in the retrieved parameters. We performed a test retrieval without the GP kernel to assess the impact of correlated noise on the retrieved parameters. The results show that the uncertainties in the retrieved parameters are underestimated when correlated noise is not accounted for (see Fig. E.1). This consideration is particularly pertinent to HRS, where correlated noise can substantially alter the interpretation of the retrieved parameters (de Regt et al. 2024).

8 Conclusions

We analysed high-resolution CRIRES+ data for three young BDs, gaining insights into their atmospheric composition, including carbon and oxygen isotopes, thermal structure, rotational velocities, and the presence of veiling. The retrieved thermal profiles of the three objects have similar structures, reflecting their comparable SpTs. The surface gravities and rotational velocities are consistent with those of young, low-gravity objects, indicating the slow rotation typical of young BDs still undergoing contraction.

The retrieved metallicities are broadly consistent with solar values, but their interpretation is complicated by the strong correlation between surface gravity and metallicity. The uncertainties in surface gravity are likely underestimated due to this correlation (see Fig. 6). We improved the existing methodology for atmospheric retrievals by incorporating correlated noise, which is crucial for HRS (see Appendix E). Our results highlight the importance of understanding and accounting for these correlations when interpreting atmospheric parameters, especially in the context of degeneracies between surface gravity and metallicity.

Precise measurements of the 12C/13C isotope ratios are challenging, yet we demonstrate that HRS can constrain the isotopic composition of BDs. The measured 12C/13C in our sample are consistent with that of the ISM within 1σ, considering the scatter of the ISM value (12C/13C)ISM= 68 ± 15. This supports the idea that isolated BDs share a common formation mechanism with stars, likely through fragmentation processes. These findings add valuable data points for assessing the role of the 12C/13C as a formation tracer and suggest that oxygen isotopes might be accessible with current and future K-band HRS. The measured C/Os are consistent with the solar value and show no significant enrichment in carbon or oxygen. Additional measurements of the C/O in young BDs are needed to assess the role of this parameter in the formation and evolution of substellar objects. Looking ahead, expanding our sample size to include more isolated BDs and substellar companions, as planned in the ESO SupJup Survey, will enhance our understanding of the 12C/13C as a potential formation tracer and deepen our insights into the diverse atmospheres of BDs and SJs, including new measurements of the C/O across objects of different ages and SpTs. Additional isotope ratios might be within the reach of state-of-the-art high-resolution spectrographs, such as CRIRES+, and current space-based facilities, such as JWST (e.g. Gandhi et al. 2023; Barrado et al. 2023). In the future, high-resolution spec-trographs on next-generation ground-based telescopes, such as the METIS instrument at the E-ELT (Brandl et al. 2021), will be powerful tools for characterising the atmospheres and measuring isotope ratios of cool BDs and directly imaged exoplanets.

Acknowledgements

D.G.P. and I.S. acknowledge NWO grant OCENW.M.21.010. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme(s) 1110.C-4264(F). This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grant no. EINF-4556. This research has made use of NASA’s Astrophysics Data System. This research has made use of adstex (https://github.com/yymao/adstex). Software: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), petitRADTRANS (Mollière et al. 2019), PyAstronomy (Czesla et al. 2019), Astropy (Astropy Collaboration 2022), corner (Foreman-Mackey 2016).

Appendix A Telluric correction with Molecfit

The open-source software Molecfit (Smette et al. 2015) is a tool for correcting astronomical observations for the transmission of the Earth’s atmosphere. Molecfit fits a telluric model to the observed data based on airmass, precipitable water vapour, the thermal profile of the Earth’s atmosphere at the time of the observations and the abundances of the main opacity sources (H2O, CO2, CO, CH4, and N2O). During the fitting process, the wavelength solution and the continuum of the model are adjusted to match the observed spectrum. Additionally, Molecfit fits for the instrumental profile of the spectrograph, which is crucial for HRS. The best-fit telluric model from the standard star is then scaled to the airmass and precipitable water vapour of the science observations. Lastly, the telluric model is divided from the observed data to obtain the telluric-corrected spectrum (see Fig. A.1). The refined wavelength solution and the continuum of the telluric model are used to calibrate the data. The fitted continuum on the standard star is divided by the black body function corresponding to the effective temperature of the star to obtain the wavelength-dependent throughput. From the fitting of the line spread function with a Lorentzian profile we obtain an empirical estimate of the spectral resolution of the instrument, R = λΔλ= cΔv=cFWHM=50,000±1,5000,$R\, = \,{\lambda \over {\Delta \lambda }} = \,{c \over {\Delta v}} = {c \over {{\rm{FWHM}}}} = 50,000 \pm 1,5000,$(A.1)

which is in perfect agreement with the nominal resolution of CRIRES+ in the wide slit mode (R ~ 50,000) and is equivalent to a resolution element of Δv ~ 3 km s−1.

thumbnail Fig. A.1

Telluric correction for TWA 28. Each panel shows a different spectral order. The corrected spectrum (black) is the result of dividing the data by the telluric model (red). Telluric lines deeper than 0.65 (continuum normalised to 1.0) with respect to the continuum are masked for the atmospheric retrieval.

Appendix B Hydrogen fluoride

thumbnail Fig. B.1

Top: Retrieved HF line depths for J1200. Bottom: Best-fitting models with HF (magenta) and without HF (blue), plotted around three regions with strong HF lines.
thumbnail Fig. B.2

Top: Retrieved HF line depths for TWA 28. Bottom: Best-fitting models with HF (magenta) and without HF (green), plotted around three regions with strong HF lines.
thumbnail Fig. B.3

Top: Retrieved HF line depths for J0856. Bottom: Best-fitting models with HF (magenta) and without HF (red), plotted around three regions with strong HF lines.

Appendix C Extended corner plots

thumbnail Fig. C.1

Posterior distributions of the retrieved parameters of J1200. The reduced chi-squared values for the models with and without 13CO and H218O are shown in the top-right corner along with the Bayes factors, Bm, and the corresponding significance levels.
thumbnail Fig. C.2

Posterior distributions of the retrieved parameters of TWA 28. The reduced chi-squared values for the models with and without 13CO and H218O are shown in the top-right corner along with the Bayes factors, Bm, and the corresponding significance levels.
thumbnail Fig. C.3

Posterior distributions of the retrieved parameters of J0856. The reduced chi-squared values for the models with and without 13CO and H218O are shown in the top-right corner along with the Bayes factors, Bm, and the corresponding significance levels.

Appendix D Best-fitting spectra

thumbnail Fig. D.1

The best-fitting models of the three objects for the first (bluest) spectral order. The central detector is masked due to the presence of strong telluric lines.
thumbnail Fig. D.2

The best-fitting models of the three objects for the second spectral order.
thumbnail Fig. D.3

The best-fitting models of the three objects for the third spectral order. The Brackett-γ line (λ2166 nm) of the standard star affects the central detector of this spectral order for J1200. It is masked at the preprocessing stage.
thumbnail Fig. D.4

The best-fitting models of the three objects for the fourth spectral order.
thumbnail Fig. D.5

The best-fitting models of the three objects for the sixth spectral order. This region contains several strong telluric lines that are masked at the preprocessing stage.

Appendix E Correlated noise comparison

thumbnail Fig. E.1

Posterior distributions of the retrieved parameters of J0856: including correlated noise with a GP (red) and without accounting for correlated noise (cyan).

References

  1. Alcalá, J. M., Gangi, M., Biazzo, K., et al. 2021, A&A, 652, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Aléon, J., Lévy, D., Aléon-Toppani, A., et al. 2022, Nat. Astron., 6, 458 [CrossRef] [Google Scholar]
  3. Allard, N. F., Spiegelman, F., Leininger, T., & Molliere, P. 2019, A&A, 628, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Allers, K. N., Jaffe, D. T., Luhman, K. L., et al. 2007, ApJ, 657, 511 [Google Scholar]
  5. Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  8. Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563 [CrossRef] [EDP Sciences] [Google Scholar]
  9. Barrado, D., Mollière, P., Patapis, P., et al. 2023, Nature, 624, 263 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 332, L65 [NASA ADS] [CrossRef] [Google Scholar]
  11. Beck, T. L., Bary, J. S., & McGregor, P. J. 2010, ApJ, 722, 1360 [NASA ADS] [CrossRef] [Google Scholar]
  12. Benneke, B., & Seager, S. 2013, ApJ, 778, 153 [Google Scholar]
  13. Bonnefoy, M., Chauvin, G., Lagrange, A. M., et al. 2014, A&A, 562, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Borysow, J., Frommhold, L., & Birnbaum, G. 1988, ApJ, 326, 509 [NASA ADS] [CrossRef] [Google Scholar]
  15. Boucher, A., Lafrenière, D., Gagné, J., et al. 2016, ApJ, 832, 50 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bouvier, J., Matt, S. P., Mohanty, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 433 [Google Scholar]
  17. Brandl, B., Bettonvil, F., van Boekel, R., et al. 2021, The Messenger, 182, 22 [NASA ADS] [Google Scholar]
  18. Brewer, J. M., & Fischer, D. A. 2016, ApJ, 831, 20 [CrossRef] [Google Scholar]
  19. Brewer, J. M., Fischer, D. A., & Madhusudhan, N. 2017, AJ, 153, 83 [NASA ADS] [CrossRef] [Google Scholar]
  20. Bryan, M. L., Ginzburg, S., Chiang, E., et al. 2020, ApJ, 905, 37 [NASA ADS] [CrossRef] [Google Scholar]
  21. Buchner, J. 2016, PyMultiNest: Python interface for MultiNest, Astrophysics Source Code Library, [record ascl:1606.005] [Google Scholar]
  22. Cameron, E., & Pettitt, A. 2014, Stat. Sci., 29, 397 [CrossRef] [Google Scholar]
  23. Carvalho, A., & Johns-Krull, C. M. 2023, RBAAS, 7, 91 [Google Scholar]
  24. Castelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, 210, eds. N. Piskunov, W. W. Weiss, & D. F. Gray, A20 [NASA ADS] [Google Scholar]
  25. Chan, Y. M., & Dalgarno, A. 1965, Proc. Phys. Soc., 85, 227 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chauvin, G., Lagrange, A. M., Dumas, C., et al. 2005, A&A, 438, L25 [CrossRef] [EDP Sciences] [Google Scholar]
  27. Cooper, W. J., Smart, R. L., Jones, H. R. A., & Sarro, L. M. 2024, MNRAS, 527, 1521 [Google Scholar]
  28. Costes, J. C., Xuan, J. W., Vigan, A., et al. 2024, A&A, 686, A294 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Cruz, K. L., Kirkpatrick, J. D., & Burgasser, A. J. 2009, AJ, 137, 3345 [NASA ADS] [CrossRef] [Google Scholar]
  30. Currie, T., Daemgen, S., Debes, J., et al. 2014, ApJ, 780, L30 [Google Scholar]
  31. Cushing, M. C., Rayner, J. T., & Vacca, W. D. 2005, ApJ, 623, 1115 [Google Scholar]
  32. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
  33. Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages, Astrophysics Source Code Library, [record ascl:1906.010] [Google Scholar]
  34. Dalgarno, A., & Williams, D. A. 1962, ApJ, 136, 690 [NASA ADS] [CrossRef] [Google Scholar]
  35. de Regt, S., Gandhi, S., Snellen, I. A. G., et al. 2024, A&A, 688, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Delorme, J.-R., Jovanovic, N., Echeverri, D., et al. 2021, J. Astron. Telescopes Instrum. Syst., 7, 035006 [NASA ADS] [Google Scholar]
  37. Dittmann, A. J. 2024, arXiv e-prints, [arXiv:2404.16928] [Google Scholar]
  38. Donahue, T. M., Hoffman, J. H., Hodges, R. R., & Watson, A. J. 1982, Science, 216, 630 [NASA ADS] [CrossRef] [Google Scholar]
  39. Dorn, R. J., Bristow, P., Smoker, J. V., et al. 2023, A&A, 671, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016, ApJS, 225, 10 [Google Scholar]
  41. Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  42. Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10 [Google Scholar]
  43. Feuchtgruber, H., Lellouch, E., Orton, G., et al. 2013, A&A, 551, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Folha, D. F. M., & Emerson, J. P. 1999, A&A, 352, 517 [NASA ADS] [Google Scholar]
  45. Foreman-Mackey, D. 2016, J. Open Source Softw., 1, 24 [Google Scholar]
  46. Gandhi, S., de Regt, S., Snellen, I., et al. 2023, ApJ, 957, L36 [NASA ADS] [CrossRef] [Google Scholar]
  47. Goderis, S., Brandon, A. D., Mayer, B., & Humayun, M. 2016, Geochim. Cos-mochim. Acta, 191, 203 [NASA ADS] [CrossRef] [Google Scholar]
  48. Gray, D. F. 2022, The observation and analysis of stellar photospheres (Cambridge University Press) [Google Scholar]
  49. Hagemann, R., Nief, G., & Roth, E. 1970, Tellus, 22, 712 [Google Scholar]
  50. Hall, D. N., & Noyes, R. W. 1969, ApJ, 4, 143 [Google Scholar]
  51. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  52. Harrison, T. E., & Marra, R. E. 2017, ApJ, 843, 152 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hartigan, P., Hartmann, L., Kenyon, S., Hewett, R., & Stauffer, J. 1989, ApJS, 70, 899 [NASA ADS] [CrossRef] [Google Scholar]
  54. Helling, C., & Casewell, S. 2014, A&A Rev., 22, 80 [NASA ADS] [CrossRef] [Google Scholar]
  55. Hood, C. E., Mukherjee, S., Fortney, J. J., et al. 2024, arXiv e-prints [arXiv:2402.05345] [Google Scholar]
  56. Horne, K. 1986, PASP, 98, 609 [Google Scholar]
  57. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  58. Inglis, J., Wallack, N. L., Xuan, J. W., et al. 2024, AJ, 167, 218 [NASA ADS] [CrossRef] [Google Scholar]
  59. Iyer, A. R., Line, M. R., Muirhead, P. S., Fortney, J. J., & Gharib-Nezhad, E. 2023, ApJ, 944, 41 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [Google Scholar]
  61. Kawahara, H., Kawashima, Y., Masuda, K., et al. 2022, ApJS, 258, 31 [NASA ADS] [CrossRef] [Google Scholar]
  62. Kitzmann, D., Stock, J. W., & Patzer, A. B. C. 2024, MNRAS, 527, 7263 [Google Scholar]
  63. Landman, R., Stolker, T., Snellen, I. A. G., et al. 2024, A&A, 682, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Lee, S., Nomura, H., & Furuya, K. 2024, ApJ, 969, 41 [NASA ADS] [CrossRef] [Google Scholar]
  65. Leschinski, K. 2021, SkyCalc_ipy: SkyCalc wrapper for interactive Python, Astrophysics Source Code Library, [record ascl:2109.007] [Google Scholar]
  66. Li, G., Gordon, I. E., Rothman, L. S., et al. 2015, ApJS, 216, 15 [NASA ADS] [CrossRef] [Google Scholar]
  67. Line, M. R., Teske, J., Burningham, B., Fortney, J. J., & Marley, M. S. 2015, ApJ, 807, 183 [NASA ADS] [CrossRef] [Google Scholar]
  68. Line, M. R., Brogi, M., Bean, J. L., et al. 2021, Nature, 598, 580 [NASA ADS] [CrossRef] [Google Scholar]
  69. Lyons, J. R., Gharib-Nezhad, E., & Ayres, T. R. 2018, Nat. Commun., 9, 908 [Google Scholar]
  70. Madhusudhan, N. 2019, ARA&A, 57, 617 [NASA ADS] [CrossRef] [Google Scholar]
  71. Madhusudhan, N., Bitsch, B., Johansen, A., & Eriksson, L. 2017, MNRAS, 469, 4102 [NASA ADS] [CrossRef] [Google Scholar]
  72. Maiorca, E., Uitenbroek, H., Uttenthaler, S., et al. 2014, ApJ, 788, 149 [NASA ADS] [CrossRef] [Google Scholar]
  73. Mamajek, E. E. 2005, ApJ, 634, 1385 [NASA ADS] [CrossRef] [Google Scholar]
  74. Manjavacas, E., Tremblin, P., Birkmann, S., et al. 2024, AJ, 167, 168 [NASA ADS] [CrossRef] [Google Scholar]
  75. Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279 [Google Scholar]
  76. Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85 [NASA ADS] [CrossRef] [Google Scholar]
  77. Mayer, L., Quinn, T., Wadsley, J., & Stadel, J. 2002, Science, 298, 1756 [NASA ADS] [CrossRef] [Google Scholar]
  78. McClure, M. K., Calvet, N., Espaillat, C., et al. 2013, ApJ, 769, 73 [NASA ADS] [CrossRef] [Google Scholar]
  79. McElwain, M. W., Metchev, S. A., Larkin, J. E., et al. 2007, ApJ, 656, 505 [Google Scholar]
  80. Meija, J., Coplen, T. B., Berglund, M., et al. 2016, Pure Appl. Chem., 88, 293 [Google Scholar]
  81. Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126 [Google Scholar]
  82. Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67 [Google Scholar]
  83. Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131 [Google Scholar]
  84. Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74 [CrossRef] [Google Scholar]
  85. Morales-Gutiérrez, S., Nagel, E., & Barragan, O. 2021, MNRAS, 506, 5361 [CrossRef] [Google Scholar]
  86. Mukherjee, S., Fortney, J. J., Morley, C. V., et al. 2024, ApJ, 963, 73 [NASA ADS] [CrossRef] [Google Scholar]
  87. Nielsen, E., De Rosa, R., Macintosh, B., et al. 2019, in AAS/Division for Extreme Solar Systems Abstracts, 51, 100.02 [Google Scholar]
  88. Nissen, P. E. 2013, A&A, 552, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [Google Scholar]
  90. Phillips, M. W., Liu, M. C., & Zhang, Z. 2024, ApJ, 961, 210 [NASA ADS] [CrossRef] [Google Scholar]
  91. Pierel, J. D. R., Nixon, C. A., Lellouch, E., et al. 2017, AJ, 154, 178 [NASA ADS] [CrossRef] [Google Scholar]
  92. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  93. Polyansky, O. L., Kyuberis, A. A., Lodi, L., et al. 2017, MNRAS, 466, 1363 [NASA ADS] [CrossRef] [Google Scholar]
  94. Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597 [NASA ADS] [CrossRef] [Google Scholar]
  95. Rei, A. C. S., Petrov, P. P., & Gameiro, J. F. 2018, A&A, 610, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001 [NASA ADS] [CrossRef] [Google Scholar]
  97. Romano, D., Matteucci, F., Zhang, Z. Y., Papadopoulos, P. P., & Ivison, R. J. 2017, MNRAS, 470, 401 [NASA ADS] [CrossRef] [Google Scholar]
  98. Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spec. Radiat. Transf., 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
  99. Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200 [Google Scholar]
  100. Ruffio, J.-B., Konopacky, Q. M., Barman, T., et al. 2021, AJ, 162, 290 [NASA ADS] [CrossRef] [Google Scholar]
  101. Ruffio, J.-B., Perrin, M. D., Hoch, K. K. W., et al. 2024, AJ, 168, 73 [NASA ADS] [CrossRef] [Google Scholar]
  102. Scholz, A., Jayawardhana, R., & Brandeker, A. 2005, ApJ, 629, L41 [NASA ADS] [CrossRef] [Google Scholar]
  103. Schutte, M. C., Lawson, K. D., Wisniewski, J. P., et al. 2020, AJ, 160, 156 [Google Scholar]
  104. Skilling, J. 2004, AIP Conf. Ser., 735, 395 [Google Scholar]
  105. Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Snellen, I. A. G., Brandl, B. R., De Kok, R. J., et al. 2014, Nature, 509, 63 [NASA ADS] [CrossRef] [Google Scholar]
  107. Sousa, A. P., Bouvier, J., Alencar, S. H. P., et al. 2023, A&A, 670, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Stahler, S. W. 1983, ApJ, 274, 822 [Google Scholar]
  109. Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392, 413 [Google Scholar]
  110. Stempels, H. C., & Piskunov, N. 2003, A&A, 408, 693 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182 [EDP Sciences] [Google Scholar]
  112. Sullivan, K., & Kraus, A. L. 2022, ApJ, 928, 134 [NASA ADS] [CrossRef] [Google Scholar]
  113. Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al. 2016, J. Mol. Spectrosc., 327, 73 [Google Scholar]
  114. Venuti, L., Stelzer, B., Alcalá, J. M., et al. 2019, A&A, 632, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Vigan, A., Langlois, M., Moutou, C., & Dohlen, K. 2008, A&A, 489, 1345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72 [EDP Sciences] [Google Scholar]
  117. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  118. Vos, J. M., Biller, B. A., Allers, K. N., et al. 2020, AJ, 160, 38 [NASA ADS] [CrossRef] [Google Scholar]
  119. Wallace, L., & Hinkle, K. 1996, ApJS, 107, 312 [NASA ADS] [CrossRef] [Google Scholar]
  120. Wang, J., Kolecki, J. R., Ruffio, J.-B., et al. 2022, AJ, 163, 189 [CrossRef] [Google Scholar]
  121. White, R. J., & Basri, G. 2003, ApJ, 582, 1109 [Google Scholar]
  122. Wilson, T. L. 1999, Rep. Progr. Phys., 62, 143 [CrossRef] [Google Scholar]
  123. Wilzewski, J. S., Gordon, I. E., Kochanov, R. V., Hill, C., & Rothman, L. S. 2016, J. Quant. Spec. Radiat. Transf., 168, 193 [NASA ADS] [CrossRef] [Google Scholar]
  124. Woitke, P., Helling, C., & Gunn, O. 2020, A&A, 634, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Xuan, J. W., Wang, J., Finnerty, L., et al. 2024, ApJ, 962, 10 [NASA ADS] [CrossRef] [Google Scholar]
  126. Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021a, Nature, 595, 370 [NASA ADS] [CrossRef] [Google Scholar]
  127. Zhang, Y., Snellen, I. A. G., & Mollière, P. 2021b, A&A, 656, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Zhang, Y., Snellen, I. A. G., Brogi, M., & Birkby, J. L. 2022, RNAAS, 6, 194 [NASA ADS] [Google Scholar]
  129. Zhang, Z., Mollière, P., Hawkins, K., et al. 2023, AJ, 166, 198 [NASA ADS] [CrossRef] [Google Scholar]
7

Correlation is quantified by the Pearson coefficient, calculated using scipy.stats.pearsonr, where a value of r = 1 indicates perfect positive correlation and r = 0 indicates no correlation.

All Tables

Table 1

Fundamental parameters of the young BDs.

In the text
Table 2

Observational properties of the young BDs.

In the text
Table 3

Retrieved parameters and their uncertainties.

In the text

All Figures

thumbnail Fig. 1

Colour-magnitude diagram of MJ vs. J − K. The objects observed in the SupJup Survey are indicated with three types of markers corresponding to isolated BDs, substellar companions, and hosts, which are either primary stars or BDs. The population of isolated cool dwarfs as catalogued in the UltracoolSheet1 is displayed as a reference. The colour of the points indicates the SpTs.
In the text
thumbnail Fig. 2

Thermal profiles (right) and temperature gradients (left) of the SPHINX model grid (https://zenodo.org/records/7416042) covering the range of temperatures between 2000 and 3000 K. The blue envelopes show the parameter space covered by the selected priors for the temperature gradients.
In the text
thumbnail Fig. 3

Opacity sources included in the forward modelling. The opacities are shown for a temperature of 2400 K and at the retrieved mixing ratios of TWA 28 (see Table 3).
In the text
thumbnail Fig. 4

Best-fitting models retrieved for the three objects. The spectra correspond to the three detectors of a single order of CRIRES+. This wavelength range contains several 12CO and 13CO lines. The observed data are shown in black, and the corresponding best-fit models are in blue (J1200), green (TWA 28), and red (J0856). Residuals are shown in the bottom panel. The average scaled uncertainty is indicated on the right side of each panel. The data and the models are displayed in the rest frame of each object. The full wavelength coverage used in the retrievals is shown in Appendix D.
In the text
thumbnail Fig. 5

(a) Best-fit pressure-temperature profiles. The shaded regions indicate the 1-, 2-, and 3-sigma regions. On the right axis we overplot the integrated contribution function. (b) Posterior distributions for C/O (top) and 12C/13C (bottom). The values for the solar C/O = 0.59 ± 0.08 (Asplund et al. 2021), and the isotope ratios for the Solar System, 89.3 ± 0.2 (Meija et al. 2016), and the ISM, (12C/13C)ISM= 68 ± 16 (Milam et al. 2005), are marked.
In the text
thumbnail Fig. 6

Posterior distributions of the surface gravity and carbon-to-oxygen ratio relative to the solar value for each object. The contours of the 2D histogram represent the 1, 2, and 3σ confidence intervals, with dashed lines indicating the median values. The Pearson correlation coefficient is shown in the top-right corner, where a value of r = 0 indicates no correlation and r = 1 indicates a perfect positive correlation.
In the text
thumbnail Fig. 7

CCFs for 13CO and H218O. The CCFs are calculated for each order-detector pair and summed over all orders and detectors. The CCFs are converted to S/N by dividing them by the standard deviation of the CCFs away from the peak (|υrad| > 100 km s−1). The S/N of the peaks of the CCFs is shown in the legend. The residuals between the observed and modelled CCFs – the auto-correlation function (ACF) – are shown in the bottom panel.
In the text
thumbnail Fig. 8

(a) Posterior distributions of the oxygen isotope ratio 16O/18O for the three objects. The solar value of 525 ± 21 (Lyons et al. 2018) and the ISM value of 557 ± 30 (Wilson 1999) are marked. (b) Posterior distributions of the abundance of fluorine with respect to hydrogen normalised to the solar value. The solar fluorine abundance is overplotted as a Gaussian distribution with a 1σ scatter (Maiorca et al. 2014).
In the text
thumbnail Fig. 9

Retrieved VMRs of line species compared to EC values of the retrieved temperature profiles of each target. The dashed lines indicate the values from EC at the retrieved C/O and metallicity. The shaded regions represent the retrieved 1σ confidence region around the best fit. The integrated contribution function is plotted in the left column and indicated in the right column as a dark shaded region.
In the text
thumbnail Fig. A.1

Telluric correction for TWA 28. Each panel shows a different spectral order. The corrected spectrum (black) is the result of dividing the data by the telluric model (red). Telluric lines deeper than 0.65 (continuum normalised to 1.0) with respect to the continuum are masked for the atmospheric retrieval.
In the text
thumbnail Fig. B.1

Top: Retrieved HF line depths for J1200. Bottom: Best-fitting models with HF (magenta) and without HF (blue), plotted around three regions with strong HF lines.
In the text
thumbnail Fig. B.2

Top: Retrieved HF line depths for TWA 28. Bottom: Best-fitting models with HF (magenta) and without HF (green), plotted around three regions with strong HF lines.
In the text
thumbnail Fig. B.3

Top: Retrieved HF line depths for J0856. Bottom: Best-fitting models with HF (magenta) and without HF (red), plotted around three regions with strong HF lines.
In the text
thumbnail Fig. C.1

Posterior distributions of the retrieved parameters of J1200. The reduced chi-squared values for the models with and without 13CO and H218O are shown in the top-right corner along with the Bayes factors, Bm, and the corresponding significance levels.
In the text
thumbnail Fig. C.2

Posterior distributions of the retrieved parameters of TWA 28. The reduced chi-squared values for the models with and without 13CO and H218O are shown in the top-right corner along with the Bayes factors, Bm, and the corresponding significance levels.
In the text
thumbnail Fig. C.3

Posterior distributions of the retrieved parameters of J0856. The reduced chi-squared values for the models with and without 13CO and H218O are shown in the top-right corner along with the Bayes factors, Bm, and the corresponding significance levels.
In the text
thumbnail Fig. D.1

The best-fitting models of the three objects for the first (bluest) spectral order. The central detector is masked due to the presence of strong telluric lines.
In the text
thumbnail Fig. D.2

The best-fitting models of the three objects for the second spectral order.
In the text
thumbnail Fig. D.3

The best-fitting models of the three objects for the third spectral order. The Brackett-γ line (λ2166 nm) of the standard star affects the central detector of this spectral order for J1200. It is masked at the preprocessing stage.
In the text
thumbnail Fig. D.4

The best-fitting models of the three objects for the fourth spectral order.
In the text
thumbnail Fig. D.5

The best-fitting models of the three objects for the sixth spectral order. This region contains several strong telluric lines that are masked at the preprocessing stage.
In the text
thumbnail Fig. E.1

Posterior distributions of the retrieved parameters of J0856: including correlated noise with a GP (red) and without accounting for correlated noise (cyan).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.