EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 693 (January 2025) A&A, 693 (2025) A298 Full HTML
Open Access
Issue
A&A
Volume 693, January 2025
Article Number A298
Number of page(s) 20
Section Planets, planetary systems, and small bodies
DOI https://doi.org/10.1051/0004-6361/202451936
Published online 27 January 2025

© The Authors 2025

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

The discovery of widely separated substellar companions around young stars has challenged traditional theories of planetary formation and evolution (e.g. Chauvin et al. 2005; Marois et al. 2008; Kraus & Ireland 2012). Straddling the mass boundary between planets and brown dwarfs (BDs), these companions offer unique insights into the atmospheres of young, low-mass objects (e.g. Mayer et al. 2002; Fortney et al. 2005). Young BDs and directly imaged exoplanets share similar atmospheric properties, including high temperatures, low surface gravities, and comparable chemical compositions (e.g. Patience et al. 2012; Faherty et al. 2016; Palma-Bifani et al. 2023). To simplify terminology, we refer to companions with masses below 75 MJup as substellar companions, classifying a lower-mass subset (1−30 MJup) as super-Jupiters (SJs). The current sample of SJs includes young, directly imaged companions with orbital separations of ∼1–300 au, effective temperatures (Teff) of ∼800–2600 K, and low surface gravities (e.g. Lafrenière et al. 2007; Lagrange et al. 2009; Konopacky et al. 2013; Macintosh et al. 2015; Bohn et al. 2020; Mesa et al. 2023).

Spectroscopic studies provide critical insights into atmospheric composition, thermal structure, and dynamics, enabling the testing of formation theories (e.g. Spiegel & Burrows 2012).

Substellar companions are thought to form via two main pathways: core accretion (Pollack et al. 1996; Johansen & Bitsch 2019) and disc fragmentation or gravitational instability (Boss 1997; Mayer et al. 2002). In core accretion, solid cores form first, followed by the accretion of gas from the protoplanetary disc. While this process explains the formation of Solar System gas giants, it struggles to account for planets at separations greater than 30 au (Inaba et al. 2003; Alibert et al. 2005). Gravitational instability, on the other hand, involves the direct collapse of disc material and can produce more massive objects at distances of 20–60 au (Boss 2011; Vorobyov 2013). However, this mechanism fails to fully explain the population of substellar companions at 100–200 au (Dodson-Robinson et al. 2009). Subsequent dynamical interactions, such as disc-planet or planet-planet scattering, migration, or ejection, likely influence the final positions of these objects (Veras et al. 2009; Boley et al. 2012; Bitsch et al. 2019; Carter & Stamatellos 2023).

The detection and characterisation of widely separated sub- stellar companions rely on high-contrast imaging instruments capable of spatially resolving faint objects next to bright stars (e.g. Macintosh et al. 2014; Beuzit et al. 2019). Multi-wavelength photometric and spectroscopic observations reveal physical properties such as effective temperature, surface gravity, and chemical composition (e.g. Barman et al. 2011; Konopacky et al. 2013; Hoeijmakers et al. 2018). High-resolution spectroscopy (HRS) has revolutionised atmospheric studies, enabling precise retrievals of chemical composition, thermal structure, and cloud properties, as well as measurements of spin and surface gravity (e.g. Snellen et al. 2014; Ruffio et al. 2021; Whiteford et al. 2023; Landman et al. 2024).

The carbon isotopic ratio (12C/13C) has been proposed as a tracer of formation pathways for SJs and BDs. Like the C/O ratio, the 12C/13C is possibly sensitive to the formation environment, objects formed via accretion of solid material are expected to inherit the12 C/13 C of the solid material, while objects formed via gravitational instability are expected to have a 12C/13C similar to the local interstellar medium (Öberg & Bergin 2021; Zhang et al. 2021b). Isotope fractionation processes (e.g. Langer & Penzias 1993; Visser et al. 2009) in molecular clouds and protoplanetary discs can alter the CO isotopologue ratios, leading to variations in the 12C/13C across different regions of the disc (Yoshida et al. 2022; Lee et al. 2024). The first measurement of the 12C/13C in a substellar companion was reported for YSES 1 b 12C/13C=30−7+10$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 30_{ - 7}^{ + 10}$ (Zhang et al. 2021a), employing an infrared medium-resolution spectrograph. The12 C/13 C of YSES 1 b was found to be significantly lower than that of the local ISM ((12C/13C)ISM = 68 ± 15; Milam et al. 2005) and clearly lower than the Sun’s carbon isotope ratio ((12C/13C)⊙ = 93.5 ± 3.1; Lyons et al. 2018), suggesting that the object was enriched in 13C via the accretion of 13C-rich ices beyond the snow line. A subsequent study of the atmosphere of the young isolated BD 2M 0355 found a 12C/13C of 97−18+25$97_{ - 18}^{ + 25}$, a significantly higher value than that of YSES 1 b, suggesting that the two objects formed via different mechanisms (Zhang et al. 2021b). Additional data from the science verification of the upgraded CRyogenic high-resolution InfraRed Echelle Spectrograph (CRIRES+; Dorn et al. 2023) confirmed the detection of the 12C/13C with a value of 108 ± 1° and hinted at the detection of C18O in the atmosphere of 2M 0355 (Zhang et al. 2022). Recent work with the Keck Planet Imager and Characterizer (KPIC; Delorme et al. 2020) also measured the 12C/13C in the HIP 55507 system, revealing values of 98−22+28$98_{ - 22}^{ + 28}$ and 79−16+21$79_{ - 16}^{ + 21}$ for the companion HIP 55507 B and the host star HIP 55507 A, respectively (Xuan et al. 2024b), suggesting a homogeneous isotopic composition. Costes et al. 2024 also analysed the KPIC spectrum of HD 984 B and determined a C/O=0.50 ± 0.1 and a 12C/13C of 98−25+20$98_{ - 25}^{ + 20}$, suggesting that the object formed via gravitational collapse or disc instability. Measurements from space-based instruments, namely JWST, have also been used to measure the 12C/13C in the atmosphere of the young SJ VHS 1256 b, revealing a value of 62 ± 2 and a super-solar C/O ratio of 0.69 ± 0.01, in addition to placing constraints on 16O/18O and 16O/17O that indicate values lower than that of the local ISM for the oxygen isotope ratios yet a12 C/13 C consistent with the lower value of the ISM (Gandhi et al. 2023).

To expand the sample of SJs and BDs with measured12 C/13 C values and assess the role of isotope ratios as formation tracers, the ESO SupJup Survey (PI: Snellen; de Regt et al. 2024) aims to characterise the atmospheres of young substellar companions and isolated BDs using HRS with VLT/CRIRES+. The initial results of this survey were reported by de Regt et al. (2024), who analysed the atmosphere of a field L-dwarf was analysed and found a 12C/13C of 184−40+61$184_{ - 40}^{ + 61}$. Further characterisation of the12 C/13 C ratios in three young BDs with similar spectral types and ages as GQ Lup B was presented by González Picos et al. (2024) who found values of 81−19+28,79−14+20$81_{ - 19}^{ + 28},79_{ - 14}^{ + 20}$, and 114−33+69$114_{ - 33}^{ + 69}$ for TWA 28, 2M J0856, and 2M J1200, respectively. In the third paper of the series, Zhang et al. (2024) confirmed the presence of 13CO in the atmosphere of YSES 1 b, refining the 12C/13C to 81 ± 9 and reporting the detection of atmospheric CO and H2O in YSES 1 c.

In this work we analyse the atmosphere of GQ Lup B as part of the ongoing ESO SupJup survey, focusing on its temperature profile, chemical composition, spin, surface gravity, and the determination of the 12C/13C in its atmosphere and that of its host star, GQ Lup A. The paper is structured as follows: Sect. 2 provides an overview of the GQ Lup system; Sect. 3 details the observations and data reduction; Sect. 4 presents the atmospheric retrieval of GQ Lup B; Sect. 5 describes the modelling of the spectrum of the host star; Sect. 6 reports the results; Sect. 7 discusses the findings and compares them with previous studies; and Sect. 8 presents the conclusions.

Table 1

Properties of GQ Lup A and B.

2 The GQ Lup system

The GQ Lup system, situated in the Lupus I cloud (Tachihara et al. 1996), is a young stellar system comprising GQ Lup A, a K7Ve-type T Tauri star enveloped by a circumstellar disc (Dai et al. 2010; MacGregor et al. 2017), and a substellar companion, GQ Lup B, located at a wide angular separation of approximately 0.7″, equivalent to ~14° au (Neuhäuser et al. 2005). The properties of the system are listed in Table 1. Initial mass estimates for GQ Lup B ranged from 10 to 40 MJup (Marois et al. 2007; McElwain et al. 2007), covering a spectrum from planetary-mass objects to BDs; however, to date, its mass remains poorly constrained. Near-infrared spectroscopy has established GQ Lup B as a late M- to early L-type object, with Teff ≈ 2100–2700 K and log 𝑔 ≤ 4.5 (Neuhäuser et al. 2005; Marois et al. 2007; McElwain et al. 2007; Seifahrt et al. 2007; Lavigne et al. 2009; Stolker et al. 2021). The system’s estimated age, between 2 and 5 Myr (Donati et al. 2012), aligns with the observed low- surface gravity features (Seifahrt et al. 2007). Through HRS cross-correlation techniques, Schwarz et al. (2016) identified CO and H2O within GQ Lup B’s atmosphere and determined its rotational velocity to be υ sin i=5.3−0.9+1.0$i = 5.3_{ - 0.9}^{ + 1.0}$ km s−1, a slow spin characteristic of its young age. Studies focusing on the orbital dynamics of the GQ Lup system have estimated the inclination of GQ Lup B’s orbit to be i ≈ 60 deg and the semi-major axis to be a = 100–150 au (Janson et al. 2006; Ginski et al. 2014; Schwarz et al. 2016).

Further analysis combining optical, near-, and mid-infrared data confirmed GQ Lup B’s spectral type as M9 and highlighted low surface gravity features (Stolker et al. 2021). Additionally, the presence of mid-infrared excess in GQ Lup B’s spectral energy distribution and accretion signatures in the optical spectrum suggested the existence of a proto-lunar disc (Stolker et al. 2021), positioning this system as a prime target for investigations into satellite formation around substellar companions (Lazzoni et al. 2022). Recent JWST/MIRI (Gardner et al. 2023; Argyriou et al. 2023) of the GQ Lup system in the mid-infrared have confirmed the presence of a disc around GQ Lup B, offering insights into the disc’s extent and temperature (Cugno et al. 2024).

A newly released analysis of the K-band spectrum of GQ Lup B with KPIC has reported the detection of 13CO in its atmosphere (Xuan et al. 2024a). The authors of this study also provide estimates of the elemental composition and the carbon isotope ratio of a small sample of substellar companions with masses ranging from 10 to 30 MJup. Xuan et al. (2024a) report a 12C/13C of approximately 150 and a super-solar C/O of approximately 0.70 for GQ Lup B (see Sect. 7 for a comparison with our results).

A recent addition to the GQ Lup system is the discovery of 2MASS J15491331-3539118, a low-mass star bound to the GQ Lup system, designated GQ Lup C (Alcalá et al. 2020). Positioned at a projected separation of 16.0″, approximately 2400 au away, GQ Lup C adds to the architectural diversity of the GQ Lup system. Its near-infrared spectrum indicates M4 spectral type, with effective temperature around 3200 K, mass of approximately 157 MJup, and a notably low surface gravity of 3.74 ± 0.24 (Alcalá et al. 2020). The identification of GQ Lup C not only adds valuable context to the formation pathways within the GQ Lup system but also contributes to broader discussions on the formation of companions with extremely wide orbits around young stellar objects.

3 Observations

We obtained high-resolution spectra of the GQ Lup system on February 28, 2023, as part of the ESO SupJup survey, using the CRIRES+ instrument at the VLT (Dorn et al. 2023). CRIRES+ is an upgraded version of CRIRES, featuring a long- slit, cross-dispersed echelle spectrograph equipped with the MACAO adaptive optics system (Paufique et al. 2004).

Observations were conducted in the 0.2″ narrow-slit mode, achieving a resolving power of ℛ ∼ 120 000, as determined from telluric line fits. This value exceeds the nominal ℛ = 100 000 under favourable observing conditions (see Appendix A), consistent with prior findings by Holmberg & Madhusudhan (2022). The slit was oriented at a position angle of 278.91 degrees (Ginski et al. 2014), enabling simultaneous observation of GQ Lup A and GQ Lup B (see Fig. 1). The K2166 wavelength setting was used, covering the K-band spectral range from 1.92 to 2.48 µm, with particular focus on the 13CO band heads near 2.345 µm. The total exposure time was 3 hours, split into 60 exposures of 180 seconds each, using an ABBA nodding sequence with a 5″ nod throw.

Sky conditions were photometric and stable, with an average seeing of 0.6″ and an airmass ranging from 1.5 to 1.0 (see Appendix A). To aid in removing Earth’s atmospheric transmission features, the standard star HD 1407 was observed immediately before GQ Lup using the same slit width and nodding sequence, though at a slightly higher airmass. To account for this discrepancy, we applied forward-modelling of the telluric transmission rather than directly removing it from the data (see Section 4.5).

Data reduction was performed using the Python package excalibuhr1 (Zhang et al. 2024). The reduction pipeline included the following steps:

  1. Flat-fielding: correcting for pixel-to-pixel variations in the detector using flat-field frames.

  2. Spectral order tracing: identifying the traces of the spectral orders by fitting a second-order polynomial to the spatial profile from the combined flat-field frame.

  3. Spectral order curvature tracing: determining the curvature of the spectral orders from Fabry-Pérot calibration frames and using it as an initial wavelength solution.

  4. Blaze extraction: extracting the blaze function from the flatfield frames and using it to correct the spectral orders.

  5. Nodding pair subtraction: subtracting the nodding pairs to remove the sky background.

  6. Nodding pair combination: generating a combined image for each nodding position by taking the median value of all frames in the same nodding position.

  7. Spectral extraction: extracting the spectra using optimal extraction (Horne 1986) with a Gaussian spatial profile. The extraction is performed at two different slit positions: one centred on GQ Lup A and the other on GQ Lup B. The size of the half-aperture for the extraction is set to 8 and 2 pixels for GQ Lup A and GQ Lup B, respectively, corresponding to 0.45″ and 0.11″ on the sky.

  8. Wavelength solution calibration: calibrating the wavelength solution using telluric lines present in the spectra and a telluric model matching the observing conditions obtained from Skycalc2. The telluric model is convolved with a Gaussian kernel to match the instrumental resolution. The best-fit wavelength solution is obtained by maximising the cross-correlation between the observed and telluric models shifted in wavelength by a second-order polynomial.

thumbnail Fig. 1

Calibrated 2D image of a single order-detector pair. The vertical axis is the spatial direction and the horizontal axis the spectral direction. The bright central source is GQ Lup A, and the fainter source at −0.71 arcseconds is GQ Lup B. The telluric lines are visible as vertical stripes in the image. The right panel shows the integrated line spread function and the spatial profile employed for the extraction of the data at the position of GQ Lup B.

4 Atmospheric retrieval of GQ Lup B

4.1 Linear model

The spectrum at the position of GQ Lup B was modelled using a linear framework, which assumes that the observed signal is a combination of the atmospheric model of the companion and the scattered starlight from GQ Lup A (Wilcomb et al. 2020; Ruffio et al. 2021; Xuan et al. 2024b). The model is expressed as: d=ϕM+n=ϕ0M0︸companion +ϕ1…iM1…i︸starlight +n,${\bf{d}} = \phi {\bf{M}} + {\bf{n}} = \underbrace {{\phi _0}{{\bf{M}}_0}}_{{\rm{companion }}} + \underbrace {{\phi _{1 \ldots i}}{{\bf{M}}_{1 \ldots i}}}_{{\rm{starlight }}} + {\bf{n}},$(1)

where d is the extracted spectrum at the position of GQ Lup B, M is the model matrix, ϕ is the vector of linear coefficients, and n is the noise. The coefficient ϕ0 represents the amplitude of the companion model, while ϕ1…i represent the amplitudes of the starlight model components.

The matrix M0 contains the atmospheric model of the companion, generated using petitRADTRANS based on a set of non-linear parameters ψ (see Section 4.4). The matrix M1…i represents the scattered starlight from GQ Lup A, derived from its observed spectrum (see Appendix B). The linear coefficients are determined at each likelihood evaluation using the procedure outlined in Eq. (3).

The likelihood function follows Ruffio et al. (2019) and is expressed as: ℒ(ψ,ϕ,s2)=1(2π)Nd| Σ0 |s2Nd                         ⋅exp(−12s2χ0,ϕ=ϕ˜,ψ2)                         ⋅exp(−12s2Δϕ⊤(M⊤Σ0−1M)Δϕ),$\eqalign{ & {\cal L}\left( {\psi ,\phi ,{s^2}} \right) = {1 \over {\sqrt {{{(2\pi )}^{{N_d}\left| {{\Sigma _0}} \right|{s^{2{N_d}}}}}} }} \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \exp \left( { - {1 \over {2{s^2}}}\chi _{0,\phi = \tilde \phi ,\psi }^2} \right) \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \exp \left( { - {1 \over {2{s^2}}}\Delta {\phi ^ \top }\left( {{M^ \top }\Sigma _0^{ - 1}M} \right)\Delta \phi } \right), \cr} $(2)

where Nd is the number of data points, Σ0 is the covariance matrix, s2 is the noise scaling factor, χ0,ϕ=ϕ˜,ψ2$\chi _{0,\phi = \tilde \phi ,\psi }^2$ is the optimal chi-squared value, and Δϕ represents the uncertainty in the linear parameters, ϕ=ϕ˜+Δϕ$\phi = \tilde \phi + \Delta \phi $. The optimal linear parameters ϕ and noise scaling factors s2 are obtained by minimising the chi-squared value: χ02=(d−ϕM)⊤Σ0−1(d−ϕM),∇ϕχ02=0⇒ϕ˜=(M⊤Σ0−1M)−1M⊤Σ0−1d,∇s2χ02=0⇒s2˜=1Nd(d−ϕ˜M)⊤Σ0−1(d−ϕ˜M)=χ0,ϕ=ϕ˜2Nd,$\eqalign{ & \chi _0^2 = {({\bf{d}} - \phi {\bf{M}})^ \top }\Sigma _0^{ - 1}({\bf{d}} - \phi {\bf{M}}), \cr & {\nabla _\phi }\chi _0^2 = 0 \Rightarrow \tilde \phi = {\left( {{{\bf{M}}^ \top }\Sigma _0^{ - 1}{\bf{M}}} \right)^{ - 1}}{{\bf{M}}^ \top }\Sigma _0^{ - 1}{\bf{d}}, \cr & {\nabla _{{s^2}}}\chi _0^2 = 0 \Rightarrow \widetilde {{s^2}} = {1 \over {{N_d}}}{({\bf{d}} - \tilde \phi {\bf{M}})^ \top }\Sigma _0^{ - 1}({\bf{d}} - \tilde \phi {\bf{M}}) = {{\chi _{0,\phi = \tilde \phi }^2} \over {{N_d}}}, \cr} $(3)

which in practise is solved using the non-negative least squares algorithm (Lawson & Hanson 1995)3. To marginalise over the linear parameters and noise scaling factor, we integrated the likelihood function across the parameter space (see Appendix D of Ruffio et al. 2019), which yields the log-likelihood function logℒ(ψ)=−Nd−Nϕ2log2π−12log| Σ0 |+logΓ(Nd−Nϕ+12)−12log| M⊤Σ0−1M |−Nd−Nϕ+12logχ0,ϕ=ϕ˜,ψ2,$\matrix{ {\log {\cal L}(\psi ) = } & { - {{{N_d} - {N_\phi }} \over 2}\log 2\pi - {1 \over 2}\log \left| {{{\rm{\Sigma }}_0}} \right|} \cr {} & { + \log {\rm{\Gamma }}\left( {{{{N_d} - {N_\phi } + 1} \over 2}} \right) - {1 \over 2}\log \left| {{M^ \top }{\rm{\Sigma }}_0^{ - 1}M} \right|} \cr {} & { - {{{N_d} - {N_\phi } + 1} \over 2}\log \chi _{0,\phi = \tilde \phi ,\psi }^2,} \cr } $(4)

where Nϕ is the number of linear parameters, and log Γ(n) denotes the logarithm of the gamma function4.

4.2 Correlated noise

Data uncertainties often exhibit correlations, potentially biasing parameter estimates and underestimating uncertainties (Kawahara et al. 2022). To account for correlated noise, we introduced a noise model that incorporates a covariance matrix Σ0 in the likelihood function. The covariance matrix is defined following de Regt et al. (2024) and González Picos et al. (2024): Σ0,ij=δijσi2+kij,${{\rm{\Sigma }}_{0,ij}} = {\delta _{ij}}\sigma _i^2 + {k_{ij}},$(5)

where i and j index data points, δij is the Kronecker delta, σi is the uncertainty of data point i, and kij represents the Gaussian process (GP) kernel. The GP kernel is given by: kij=a2(σi+σj2)2exp(−(xi−xj)22l2),${k_{ij}} = {a^2}{\left( {{{{\sigma _i} + {\sigma _j}} \over 2}} \right)^2}\exp \left( { - {{{{\left( {{x_i} - {x_j}} \right)}^2}} \over {2{l^2}}}} \right),$(6)

where a and l are free parameters representing the kernel amplitude and length scale, respectively.

4.3 Starlight model

The signal at the position of the companion is dominated by starlight, as illustrated in Fig. 2. The on-axis starlight is scattered across the slit due to atmospheric effects and telescope optics (Vigan et al. 2008). To model the scattered starlight at the position of the companion, we employed a technique involving a linear combination of different versions of the observed on-axis spectrum of GQ Lup A, adjusted with different velocity shifts and scaling factors (Landman et al. 2024).

Our approach involves constructing a model for the starlight using shifted versions of the observed on-axis spectrum of GQ Lup A. We considered a range of pixel shifts denoted by Ns = 2s + 1, where s ranges from −3 to 3 in steps of 1 pixel, equivalent to a maximum shift of 9 km/s, which sufficiently covers the observed broadening of the spectrum due to scattering (see Appendix B for more details).

Moreover, we accounted for low-frequency variations in the starlight model, which may stem from various sources such as atmospheric conditions, telescope optics, and detector characteristics. To model these variations, we adopted the spline decomposition method introduced in Ruffio et al. (2023). With each shifted component of the model, we generated a cubic spline decomposition with Nk = 10 knots, corresponding to approximately 200 pixels in width. The complete starlight model is represented as follows: Mstarlight=[ ms=−3,k=0…ms=3,k=9 ],${{\bf{M}}_{{\rm{starlight}}}} = \left[ {{{\bf{m}}_{s = - 3,k = 0}} \ldots {{\bf{m}}_{s = 3,k = 9}}} \right]{\rm{,}}$(7)

where the linear coefficients of the starlight model components are independent for each order-detector (Nod = 21) and nodding position (Nab = 2) and are optimised at every evaluation of the log-likelihood function (see Eq. (4)). The total number of linear parameters for a single order-detector pair of a single nodding position is Nϕ = Ns × Nk = 90, and the total number of linear parameters for all order-detector pairs and nodding positions is Nϕ = Nod × Nab × Ns × Nk = 3780.

thumbnail Fig. 2

Best-fit model of the spectrum at the position of GQ Lup B. Top: observed spectrum (black), joint best-fit model (green) and components of the starlight model and companion (blue and brown, respectively). Centre: data with the best-fit starlight model subtracted compared to the companion model. Bottom: residuals of the best-fit model. The shaded region represents the root-squared diagonal of the covariance matrix with the optimal noise scaling factor.

4.4 Atmospheric model

To accurately characterise the high-resolution spectra of GQ Lup B, we employed the radiative transfer code petitRADTRANS5 (v2.7; Mollière et al. 2019). The outgoing spectrum is calculated by solving the radiative transfer equation within a 1D plane-parallel atmosphere, where each layer is characterised by its temperature, pressure, and mass fractions of relevant chemical species. The atmospheric model is defined by several key parameters, including the temperature profile, chemical composition, and surface gravity of the object.

4.4.1 Chemical composition

The atmospheres of substellar objects, such as GQ Lup B, exhibit rich molecular compositions, that prominently feature water (H2O) and carbon monoxide (CO). Our atmospheric model incorporates a variety of line species, continuum, and collision- induced absorption (CIA) to accurately reproduce the spectral characteristics observed in GQ Lup B (see Table 2 for a comprehensive list of species and their corresponding references).

We adopted a free composition approach, treating the volume mixing ratio (VMR) of each species as a free parameter instead of relying on chemical equilibrium calculations. This method provides greater flexibility, allowing us to explore a broader range of chemical compositions, which can result in better model fits (Gandhi et al. 2023; de Regt et al. 2024; González Picos et al. 2024).

From the retrieved abundances, expressed as volume mixing ratios nX , we can infer relevant atmospheric properties and derive elemental and isotopic ratios of interest. For instance, the C/O ratio is derived from the abundances of key carbon- and oxygenbearing species: C/O=n12CO+n13COnH2O+n12CO+n13CO.${\rm{C}}/{\rm{O}} = {{{n_{12}}_{_{{\rm{CO}}}} + {n_{13}}_{_{{\rm{CO}}}}} \over {{n_{{{\rm{H}}_2}{\rm{O}}}} + {n_{12}}{\rm{CO}} + {n_{13}}_{_{{\rm{CO}}}}}}.$(8)

Similarly, the ratio of12 CO/13 CO provides insights into the isotopic composition of carbon: 12C/13C=n12COn13CO.$^{12}{\rm{C}}{/^{13}}{\rm{C}} = {{{n_{12}}_{_{{\rm{CO}}}}} \over {{n_{13}}_{_{{\rm{CO}}}}}}.$(9)

Furthermore, we can estimate the atmospheric metallicity using the carbon-to-hydrogen ratio normalised to solar metallicity: [C/H]=log10(nCnH)−log10(nCnH)⊙,$[{\rm{C}}/{\rm{H}}] = {\log _{10}}\left( {{{{n_{\rm{C}}}} \over {{n_{\rm{H}}}}}} \right) - {\log _{10}}{\left( {{{{n_{\rm{C}}}} \over {{n_{\rm{H}}}}}} \right)_ \odot },$(10)

where the solar value is log10(nCnH)⊙=−3.54${\log _{10}}{\left( {{{{n_{\rm{C}}}} \over {{n_{\rm{H}}}}}} \right)_ \odot } = - 3.54$ (Asplund et al. 2021).

Table 2

Opacity species and references.

4.4.2 Temperature profile

For widely separated companions like GQ Lup B, which receive minimal irradiation from their host stars, the internal heat flux primarily dictates the temperature profile. As one moves from high to low pressures within the atmosphere, the temperature typically decreases with altitude, following a gradient, d log T d log P=∇${{{\rm{d}}\log T} \over {{\rm{d}}\log P}} = \nabla $, that reflects the dominant energy transfer mechanism.

In the denser, deeper regions of the atmosphere, energy transport is primarily governed by convective processes, which ideally follow the adiabatic temperature gradient (∇ = ∇ad; Baraffe et al. 2002). As the altitude increases and the atmosphere becomes less dense, the mean free path of photons grows, causing some wavelength regions to become optically thin. This shift allows radiation to emerge as the dominant mechanism of energy transfer. The transition between these energy transport regimes defines the radiative-convective equilibrium (RCE; Marley & Robinson 2015) boundary.

Various parameterisations for temperature profiles have been proposed, ranging from physically motivated models to more flexible, free-form models. Physically motivated models attempt to describe the temperature structure using parameters with direct physical interpretations (Madhusudhan & Seager 2009; Mollière et al. 2020; Line et al. 2012). On the other hand, traditional planetary science models often retrieve temperature at each discrete atmospheric layer (Rodgers 2000; Irwin et al. 2008). However, in exoplanetary atmospheres, where data may be sparse and noisy, parameterised models are often preferred for their computational speed and flexibility (Line et al. 2015). Such models typically specify temperature at discrete pressure levels (Pi) with corresponding temperatures Ti, and interpolate between these levels using spline functions (Line et al. 2015).

To balance physical accuracy with computational efficiency, recent studies have introduced hybrid approaches that combine the physical basis of self-consistent models with the flexibility of free-form models. For instance, Zhang et al. (2023) proposed a parameterisation that incorporates constraints on temperature gradients from self-consistent models.

In our study, we explored two variations of the temperature parameterisation proposed by Zhang et al. (2023). The first model, termed the static gradients (SG) model, utilises a set of temperature gradients (∇i) measured at equally spaced pressure levels on logarithmic space (Pi). The second model, referred to as the dynamic gradients (DG) model, extends the SG model by allowing the positions of the measured gradients to vary. In the DG model we added three free parameters that determine the pressure level of the RCE boundary (PRCE), the spacing between the upper (Δ log Ptop) and lower pressure levels (Δ log Pbot). The resulting pressure levels, from the bottom to the top of the atmosphere, are defined as: P0=102bar,P1=PRCE−2Δ log Pbot,P2=PRCE−Δ log Pbot,P3=PRCE,P4=PRCE+Δ log Ptop,P5=PRCE+2Δ log Ptop,P6=10−5bar.$\matrix{ {} \hfill & {{P_0} = {{10}^2}{\rm{bar,}}} \hfill \cr {} \hfill & {{P_1} = {P_{{\rm{RCE}}}} - 2{\rm{\Delta }}\log {P_{bot}},} \hfill \cr {} \hfill & {{P_2} = {P_{{\rm{RCE}}}} - {\rm{\Delta }}\log {P_{{\rm{bot}}}},} \hfill \cr {} \hfill & {{P_3} = {P_{{\rm{RCE}}}},} \hfill \cr {} \hfill & {{P_4} = {P_{{\rm{RCE}}}} + {\rm{\Delta }}\log {P_{{\rm{top}}}},} \hfill \cr {} \hfill & {{P_5} = {P_{RCE}} + 2{\rm{\Delta }}\log {P_{{\rm{top}}}},} \hfill \cr {} \hfill & {{P_6} = {{10}^{ - 5}}{\rm{bar}}{\rm{.}}} \hfill \cr } $(11)

An additional constraint is imposed for the DG model to ensure that the largest temperature gradient occurs at the RCE boundary. The temperature gradients at every layer ∇j , with j = 0,…, 50, are obtained via linear interpolation of the retrieved parameters ∇j, with i = 0,…, 6, and the temperature at every layer is calculated according to Tj=Tj−1(PjPj−1)∇j,${T_j} = {T_{j - 1}}{\left( {{{{P_j}} \over {{P_{j - 1}}}}} \right)^{{\nabla _j}}},$(12)

where Tj=0 is the temperature at the bottom of the atmosphere and is retrieved as a free parameter. This approach enables us to explore a wide range of temperature profiles while maintaining physical consistency and computational efficiency.

4.4.3 Broadening

Accurately modelling line broadening mechanisms is crucial for HRS. In addition to intrinsic line broadening, the rotation of the object, the turbulence in the atmosphere, and the resolving power of the instrument contribute to the observed line profile.

Rotational broadening is influenced by υ sin i, where v is the rotational velocity and i is the inclination angle of the rotation axis relative to the line of sight. To account for this, we applied a convolution of the model spectrum with a rotational broadening kernel (Gray 2022). We used the fastRotBroad function from PyAstronomy6 to perform this convolution for each spectral order. This method achieves accuracy within 0.5% of the more computationally intensive approach of direct disc integration (Carvalho & Johns-Krull 2023).

Instrumental broadening accounts for changes in the line profile caused by the Earth’s atmosphere and the optical system of the instrument, leading up to the final detector measurement. We modelled this effect using a Gaussian profile, with its standard deviation determined by the full width at half maximum of the telluric lines observed in the spectrum of the standard star (see Sect. 3). For spectra of objects with υ sin i values greater than the instrumental full width at half maximum, the exact shape of the instrumental line profile has less of an impact, making the use of a Gaussian kernel sufficient for modelling instrumental broadening in such cases.

The total broadening of the model spectrum is obtained by convolution with the rotational broadening kernel and the instrumental line shape. The high resolving power of CRIRES+ enables precise measurements of the rotational velocity for slow rotators like GQ Lup B. Previous studies (Schwarz et al. 2016) reported a υ sin i=5.3−0.9+1.0$i = 5.3_{ - 0.9}^{ + 1.0}$ km/s for GQ Lup B using the older CRIRES instrument. Our updated retrievals provide a more precise υ sin i measurement (see Eq. (20)).

4.5 Telluric transmission

The observed spectrum contains telluric absorption features originating from the Earth’s atmosphere. The contribution on the telluric lines from the stellar signal is already included in the starlight model, as we utilised the observed on-axis spectrum as the basis for the starlight model. However, we must account for the telluric absorption features in the companion spectrum to accurately reproduce the observed spectrum.

To obtain a high-fidelity telluric model, we employed Molecfit to generate synthetic telluric spectra. We utilised the spectrum of a standard star observed immediately before the observations of GQ Lup B to fit a telluric model from Molecfit. The best-fit telluric model corresponds to the atmospheric conditions, including airmass and precipitable water vapour content, at the time of the standard star observations.

To account for changes in the telluric lines between the science and standard star observations, we scaled the telluric model to match the average precipitable water vapour content during the observations of GQ Lup B, which remained approximately constant (refer to Appendix A for details). Additionally, we modelled the scaled best-fit telluric model with one free parameter to fit the effective airmass of the telluric lines in the observed spectrum of GQ Lup B. Saturated telluric lines were masked out during the fitting process, namely we ignored telluric lines with transmission values below 1% and a region of 20 pixels around them. This procedure aims to accurately incorporate telluric features into our forward model while simultaneously propagating uncertainties in the retrieved parameters.

4.6 Model comparison

To find the model that best-fits the spectrum of GQ Lup B, we performed a Bayesian model comparison using the loglikelihood function defined in Eq. (4) and we report the goodness of fit of different models in terms of the chi-square statistic (Cochran 1952). We compared the SG and DG temperature profiles, as well as a model without 13CO to evaluate the necessity of including 13CO in the model. The Bayesian evidence is calculated as part of the retrieval process, allowing us to calculate the Bayes factor B12 between two competing models 𝓜1and 𝓜2(Kass & Raftery 1995), B12=p(d∣ℳ1)p(d∣ℳ2)=p(d∣ℳ1)p(d∣ℳ2)p(ℳ1)p(ℳ2)=Z1Z2p(ℳ1)p(ℳ2),${B_{12}} = {{p\left( {{\bf{d}}\mid {{\cal M}_1}} \right)} \over {p\left( {{\bf{d}}\mid {{\cal M}_2}} \right)}} = {{p\left( {{\bf{d}}\mid {{\cal M}_1}} \right)} \over {p\left( {{\bf{d}}\mid {{\cal M}_2}} \right)}}{{p\left( {{{\cal M}_1}} \right)} \over {p\left( {{{\cal M}_2}} \right)}} = {{{Z_1}} \over {{Z_2}}}{{p\left( {{{\cal M}_1}} \right)} \over {p\left( {{{\cal M}_2}} \right)}},$(13)

where p(d|𝓜1) and p(d|𝓜2) are the marginal likelihoods of models 𝓜1 and 𝓜2, respectively, and p(𝓜1) and p(𝓜2) are the prior probabilities of the models. The Bayes factor quantifies the relative probability of the models given the data, with a value of B12 > 1 indicating that model 𝓜1 is more probable than model 𝓜2. Assuming that the prior probabilities of the models are equal, the Bayes factor is equivalent to the ratio of the marginal likelihoods, or in logarithmic form, the difference in the log-evidence values In B12 = log Z1− log Z2. The Bayes factor can be converted to a σ value using the methodology outlined in Benneke & Seager (2013), providing a quantitative measure of the significance of the model comparison.

5 Modelling the spectrum of GQ Lup A

The spectrum of the host star, GQ Lup A, is crucial for understanding the properties of the entire system. Here, we focus on modelling the spectrum of GQ Lup A to measure its carbon isotope ratio. Instead of modelling the spectrum with a parameterised model as in the case of GQ Lup B, we use a grid of precomputed PHOENIX models extended to different carbon isotope ratios to fit the observed spectrum of GQ Lup A. This approach allows us to determine the carbon isotope ratio of GQ Lup A while accounting for the veiling effect.

GQ Lup A’s K-band spectrum is predominantly characterised by CO features and atomic lines. Previous studies have identified it as a young, magnetically active star (Donati et al. 2012). The presence of strong Zeeman broadening in atomic lines indicates significant magnetic activity, providing opportunities to measure the star’s magnetic field strength. For our analysis, we focus on the last two spectral orders, ranging from 2.20 to 2.48 µm, where CO lines dominate the spectrum and the influence of Zeeman-broadened atomic lines is small. Additionally, previous high-resolution spectroscopic studies have revealed substantial veiling of photospheric lines in GQ Lup A (Sousa et al. 2023), a common feature among young stars (McClure et al. 2013; Rei et al. 2018).

5.1 Isotopic PHOENIX models

Models of stellar atmospheres are a valuable tool for interpreting the observed spectra of stars. The PHOENIX model grid (Husser et al. 2013) provides a wide range of stellar models with varying parameters, including effective temperature, surface gravity, and metallicity. However, the standard PHOENIX grid does not include models with different carbon isotope ratios. To address this limitation, we generated a custom grid of PHOENIX models with varying 12C/13C from 16 to 301 in steps of 15. For our purposes, we used models with fundamental parameters similar to GQ Lup A (see Table 1).

5.2 Linear veiling model

Veiling reduces the depth of absorption lines, making them appear shallower than those of a main sequence star of the same spectral type (Hartigan et al. 1989), with the veiling factor (rk) typically increasing with wavelength (Fischer et al. 2011). To model the spectrum of GQ Lup A, we used a linear combination of a PHOENIX model for the star (MA) and a veiling model parameterised as a slowly varying function of wavelength (Mveiling). The veiling component was fitted to the data using a spline model with ten knots, following a similar approach used for the starlight model. The linear model is expressed as: d=ϕM+n=ϕAMA+ϕweiling Mveiling +n.$d = \phi M + n = {\phi _A}{M_A} + {\phi _{{\rm{weiling }}}}{M_{{\rm{veiling }}}} + n.$(14)

During the likelihood calculation, the amplitude of the linear parameters is fitted, providing a robust estimation of each source’s contribution. The veiling factor rk, defined as the ratio of the veiling amplitude to the amplitude of the model spectrum of GQ Lup A, is calculated as: rk=ϕveilingϕA.${r_k} = {{{\phi _{{\rm{veiling}}}}} \over {{\phi _A}}}.$(15)

To determine the best-fit model of the spectrum of GQ Lup A, we fitted the observed spectrum with the model spectrum of GQ Lup A and the veiling model. We set up a retrieval framework with six free parameters, covering the stellar properties (radial velocity, spin, and carbon isotope ratio), the effective airmass of the observations, and the parameters of the GP kernel. The likelihood function is defined similarly to Eq. (4), with the amplitude of the PHOENIX model and the veiling model as the linear components.

5.3 Retrieval of GQ Lup A

We retrieved the best-fit parameters of GQ Lup A using the linear model described in Section 5.2. To avoid degeneracies between veiling and the stellar parameters, we fixed the effective temperature to Teff=4300 K as reported in Table 1. In addition, the surface gravity and metallicity were also fixed to log g= 4.0 and [Fe/H] = 0.0, respectively. We modelled the noise of the observed spectrum using a GP kernel following the approach outlined in Section 4.2. The model spectra were broadened according to the instrumental resolution of CRIRES+ and the rotational velocity v sin i of GQ Lup A, which we fitted for as part of the retrieval process. Additionally, we forward-modelled the telluric transmission spectrum by allowing the airmass of the observations to vary as a free parameter.

6 Results

6.1 The atmosphere of GQ Lup B

The atmospheric retrieval process provides us with a best-fit model characterising the atmosphere of GQ Lup B, along with its rotational and radial velocities. This model effectively reproduces observed spectral features across the entire wavelength range up to the reported uncertainty level. However, the derived surface gravity of GQ Lup B faces challenges due to degeneracies with the temperature profile, resulting in less constrained estimates. We present the main results from the retrieval with the DG model, which is favoured over the SG with a Bayesian evidence of BDG/SG= 8.53, equivalent to a 4.5σ preference (Benneke & Seager 2013).

6.1.1 Chemical composition

Table C.1 summarises the abundances of the primary opacity sources in GQ Lup B’s atmosphere. From the retrieved abundances of H2O, 12CO and 13CO we derived the C/O ratio, C/O=0.50−0.01+0.01,${\rm{C}}/{\rm{O}} = 0.50_{ - 0.01}^{ + 0.01},$(16)

which is comparable to the solar value of C/O = 0.59 ± 0.08 (Asplund et al. 2021). The presence of 13CO is confirmed at 5.8σ and clearly detected in the cross-correlation with the best- fit model (S/N=6.8; see the top panel of Fig. 4). The derived12 C/13 C from the retrieved abundances of 12CO and 13CO is 12C/13C=53−6+7,$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 53_{ - 6}^{ + 7},$(17)

which is broadly consistent with the lower value of the local ISM of 12C/13C = 69 ± 15 (Milam et al. 2005) and lower than the Sun’s value of 12C/13C = 93.5 ±3.1 (Lyons et al. 2018).

We report the metallicity of GQ Lup B as the carbon- to-hydrogen ratio normalised to the solar value as defined in Eq. (10), [C/H]=0.50−0.17+0.16,$[{\rm{C}}/{\rm{H}}] = 0.50_{ - 0.17}^{ + 0.16},$(18)

we note that [C/H] is strongly degenerate with log g, exhibiting a Pearson correlation coefficient of r = 0.977, which indicates that the retrieved metallicity is highly dependent on the surface gravity of the object. In the best-fit model with the SG temperature profile, the retrieved metallicity is [C/H]=−0.12−0.10+0.10,$[{\rm{C}}/{\rm{H}}] = - 0.12_{ - 0.10}^{ + 0.10}$, with r = 0.95, corresponding to the lower log g retrieved in this model.

We constrained the abundances of HF, Ca, Na and Ti in the atmosphere of GQ Lup B, with values comparable to those reported for three young BDs with similar spectral types, ages, and temperatures to GQ Lup B (González Picos et al. 2024). We report the non-detection of the following species: H218O, C18O, Mg, K and Fe. Considering the temperature and wavelength range of the observations, no additional molecules are expected to be relevant sources of opacity in the atmosphere (Iyer et al. 2023).

thumbnail Fig. 3

Temperature profile of GQ Lup B retrieved using the SG and DG models. Left: temperature profile. Right: temperature gradient. Profiles from the ATMO grid (Tremblin et al. 2015) with Teff=2300 K and log g=4.0 are shown for two distinct models with adiabatic (Phillips et al. 2020) and diabatic convection (Petrus et al. 2021).

6.1.2 Temperature profile

The retrieved temperature profile and its gradient are shown in the bottom panels of Fig. 3. Both pressure-temperature (PT) parameterisations are consistent in the upper atmosphere and near the radiative-convective boundary. However, slight differences arise in the deeper atmospheric layers.

The more complex DG model is favoured over the simpler SG model, as supported by both frequentist and Bayesian analyses. From the frequentist perspective, the chi-square difference is Δχ2 = 11.24, which exceeds the critical value of 9.49 at the 5% significance level (α = 0.05) and corresponds to a p-value of 0.024. Bayesian evidence further supports this preference, with a Bayes factor of BDG/SG = 8.53, equivalent to a 4.5σ preference (Benneke & Seager 2013). We note that the retrieved parameters of the covariance matrix are identical (<1σ) for the two models, and hence the Bayesian evidence comparison should not be affected by the covariance matrix. These results underscore the DG model’s improved fit and greater flexibility in capturing the temperature profile.

The retrieval constrains the position and temperature gradient of the RCE boundary in GQ Lup B’s atmosphere, enabling the derivation of TRCE, PRCE=0.19−0.02+0.02 bar, ∇RCE=0.32−0.02+0.01,TRCE=2319.91−43.80+36.28 K,$\eqalign{ & {P_{{\rm{RCE}}}} = 0.19_{ - 0.02}^{ + 0.02}{\rm{ bar, }} \cr & {\nabla _{{\rm{RCE}}}} = 0.32_{ - 0.02}^{ + 0.01}{\rm{,}} \cr & {T_{{\rm{RCE}}}} = 2319.91_{ - 43.80}^{ + 36.28}{\rm{K}}, \cr} $(19)

this highlights the pressure and temperature range of the atmosphere that our observations probe. The temperature gradient at the RCE boundary aligns with the expected value for the adiabatic cooling of a hydrogen and helium gas mixture (∇ad ≈ 0.28–0.38; Cox & Giuli 1968; Hansen et al. 2004).

thumbnail Fig. 4

Cross-correlation functions (CCFs) of the data and the best- fit models of 13CO for GQ Lup B (top) and GQ Lup A (bottom). The CCFs are shown in black, with the auto-correlation functions of the templates shown in brown and orange. The CCFs are calculated between the data residuals and the model residuals following the methodology from Zhang et al. (2021a) (also in de Regt et al. 2024).

6.1.3 Rotational and radial velocity

The best-fit projected rotational velocity of GQ Lup B is υsini=5.56−0.02+0.02 km/s,$\upsilon \sin i = 5.56_{ - 0.02}^{ + 0.02}{\rm{km}}/{\rm{s}},$(20)

which is consistent with the previous measurement of vsini=5.3−0.9+1.0$v\sin i = 5.3_{ - 0.9}^{ + 1.0}$ km/s using the old CRIRES (Schwarz et al. 2016), confirming a slow spin characteristic of its young age. The retrieved radial velocity of GQ Lup B, corrected for barycentric motion is υradB=−0.49−0.02+0.02 km/s−υsys=2.03−0.04+0.04 km/s,$\upsilon _{{\rm{rad}}}^B = - 0.49_{ - 0.02}^{ + 0.02}{\rm{km}}/{\rm{s}} - {\upsilon _{{\rm{sys}}}} = 2.03_{ - 0.04}^{ + 0.04}{\rm{km}}/{\rm{s}},$(21)

where we have used the radial velocity of the host star as the systemic velocity, vsys =−2.52−0.04+0.04${v_{{\rm{sys }}}} = - 2.52_{ - 0.04}^{ + 0.04} $ km/s (see Eq. (24)). Our radial velocity measurement is in good agreement with the earlier measurement of 2.0 ± 0.4 km/s by Schwarz et al. (2016). We note that our reported uncertainty does not account for potential variability due to magnetospheric accretion as reported in the literature, which could introduce an additional uncertainty source of approximately 0.4 km/s (Donati et al. 2012).

6.1.4 Surface gravity

The surface gravity of the best-fit model of GQ Lup B is logg=3.83−0.18+0.17.$\log g = 3.83_{ - 0.18}^{ + 0.17}.$(22)

This value should be interpreted with caution as it is subject to degeneracies with the temperature profile and the absolute abundances or metallicity. We performed a test retrieval applying a high-pass filter to the data and the model to remove the continuum, which led to similar results to the nominal retrieval.

6.2 The spectrum of GQ Lup A

The veiling model effectively captures the continuum excess and the decreased line depths present in the observed spectrum, providing a robust estimation of the veiling factor at each wavelength and enabling the retrieval of the isotope ratio with PHOENIX models. The best-fit model of GQ Lup A is shown in Fig. 5, with the residuals indicating a good fit to the data to within <5%. The contribution from the veiling continuum and the PHOENIX model are shown in Appendix D.

6.2.1 Veiling factor

The excess continuum observed in the spectrum of GQ Lup A dominates over the PHOENIX continuum. The strength of veiling increases with wavelength, with a veiling factor of 2.78 ± 0.13 and 3.85 ± 0.19 at 2346 nm and 2448 nm, respectively. In Fig. 6 we show that the veiling factor increases linearly with wavelength over the observed range, consistent with a previous measurement at 2260 nm (Sousa et al. 2023).

6.2.2 Rotational and radial velocity

The retrieved projected rotational velocity of GQ Lup A υsini=5.37±0.05 km/s,$\upsilon \sin i = 5.37 \pm 0.05{\rm{km}}/{\rm{s}},$(23)

indicates a slow rotator, in line with expected values for young stars (Bouvier et al. 2014). Previous measurements of the rotational velocity of GQ Lup A have reported values ranging from 5 to 7 km/s (Guenther et al. 2005; Donati et al. 2012; Schwarz et al. 2016). The presence of stellar spots and magnetic activity in GQ Lup A can lead to variations in the observed v sin i , which may explain the range of values reported in the literature. Donati et al. (2012) conducted a detailed study of the magnetic field of GQ Lup A on spectropolarimetric data and reported a rotational velocity of v sin i = 5.0 ± 1.0 km/s. Our measurement is consistent with this value, thus providing a refined measurement of the projected rotational velocity of GQ Lup A. We report a radial velocity of υrad=−2.52±0.04 km/s${\upsilon _{{\rm{rad}}}} = - 2.52 \pm 0.04{\rm{km}}/{\rm{s}}$(24)

which is broadly consistent with previous measurements of the systemic velocity (Donati et al. 2012; Schwarz et al. 2016).

6.2.3 Carbon isotope ratio

The best-fit carbon isotope ratio for GQ Lup A is 12C/13C=51−8+10,$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 51_{ - 8}^{ + 10},$(25)

with uncertainties that include a systematic error of Δ12C/13C=5, arising from the level of telluric masking (0.10–0.50) and the selection of Teff for the PHOENIX model within the 4200– 4400 K range. The significance of this measurement is underscored by the Bayes factor, comparing models with and without 13CO in the PHOENIX grid, resulting in In B with/without = 15.1, equivalent to a 5.8σ preference for the inclusion of 13CO.

Additionally, the cross-correlation function of the residuals between the data and the GQ Lup A model further corroborates a strong detection of 13CO (S/N=7.8; see the bottom panel of Fig. 4).

thumbnail Fig. 5

Best-fit model of GQ Lup A. a. Spectrum with the veiling continuum subtracted and divided by the telluric model. The observed spectrum is shown in black, the PHOENIX model with optimal linear amplitudes in orange and the 13CO lines from the best-fit model in purple. The residuals of the are shown in the bottom panel, the shaded grey region indicates the root-squared diagonal of the covariance matrix with the GP kernel. b. Posterior distributions of the carbon isotope ratio, projected rotational velocity and radial velocity of GQ Lup A.
thumbnail Fig. 6

Veiling factor of GQ Lup A as a function of wavelength. The veiling factors for each spectral order derived from the fitted linear amplitudes are shown in orange, with error bars representing the standard deviation between the three detectors of each order. The measurement from Sousa et al. (2023) is included to fit the three points with a linear function of slope 0.011 nm−1.

7 Discussion

7.1 Comparison with previous studies

The analysis of the KPIC spectrum of GQ Lup B, as presented in Xuan et al. (2024a) (hereafter X24), reports a 12C/13C value for GQ Lup B that is inconsistent with our findings at a significance level of ≥3σ: 12C/13CX24=153−31+43,12C/13Cthis work =53−6+7,$\eqalign{ & ^{12}{\rm{C}}{/^{13}}{{\rm{C}}_{{\rm{X}}24}} = 153_{ - 31}^{ + 43}, \cr & ^{12}{\rm{C}}{/^{13}}{{\rm{C}}_{{\rm{this work }}}} = 53_{ - 6}^{ + 7}, \cr} $

both studies confirm the presence of 13CO in the atmosphere of GQ Lup B; however, the derived 12C/13C values differ significantly. Differences in wavelength coverage and spectral resolution between KPIC and CRIRES+ likely contribute to this discrepancy. X24 employs the three reddest spectral orders of KPIC, spanning wavelength ranges of (2290–2340 nm), (2360–2410 nm), and (2440–2490 nm) at a resolution of 𝓡 ∼ 35 000. In contrast, our analysis utilises the full wavelength range of CRIRES+ K2166, covering seven spectral orders between 1920 nm and 2472 nm (see detailed ranges in Fig. C.1) at 𝓡 ∼ 120 000. Furthermore, our observations achieve a higher signal-to-noise ratio (S/N), with an estimated S/N ∼ 20 for the companion spectrum compared to S/N ∼ 12 for KPIC. Our retrieval independently fits the abundances of individual species without enforcing equilibrium chemistry, whereas X24 derives species abundances from free parameters of C/O and metallicity under the assumption of chemical equilibrium.

X24 reports a super-adiabatic temperature gradient around the photosphere of GQ Lup B, exceeding predictions from RCE (see Fig. E.1). In contrast, our retrieved temperature profile shows a high photospheric gradient of approximately 0.32, though significantly lower than 0.50 as retrieved by X24, and we do not find evidence of a super-adiabatic region. These discrepancies in temperature profiles may contribute to differences in the derived abundances and 12C/13C values, as discussed in Section 7.2.

The C/O ratio retrieved by X24 (0.70 ± 0.01) is higher than the solar value and our result (0.50 ± 0.01). A higher abundance of 12CO could inflate both the C/O and 12C/13C values while maintaining a constant 13CO abundance. Our independent fits to the main opacity sources allow us to avoid assumptions of chemical equilibrium, providing an alternative perspective. Further studies are required to reconcile these differences, including a detailed comparison of individual abundances and the interaction between the pressure-temperature profile and surface gravity.

7.2 Atmospheric retrieval with fixed thermal profile from Xuan et al. (2024)

The absorption line depths of atmospheric species are sensitive to the temperature profile at different altitudes. Spectral lines of different species form at distinct atmospheric layers depending on the temperature structure and opacity. As a result, species probing different line depths sample varying pressure regions. The 12C/13C is particularly sensitive to the temperature profile because 13CO lines originate deeper in the atmosphere than 12CO lines. These regions can be identified using the contribution function (Mollière et al. 2019), which calculates the relative line contrast of a species across pressures, indicating relevant altitudes where absorption lines originate (see Fig. C.3).

To evaluate whether the temperature profiles retrieved in this work and X24 explain the 12C/13C discrepancy, we performed a retrieval using the fixed temperature profile from X24 (FPTX24). The posterior distributions for our nominal retrieval (DG model) and FPTX24 are shown in Fig. E.1. A preference for the DG model is observed, with an 8.6σ improvement in fit quality over FPTX24. The FPTX24 retrieval yields systematically lower abundances and a reduced surface gravity: loggFPTX24=2.28−0.12+0.14<3.83−0.18+0.17=loggDG.$\log {g_{{\rm{FPTX}}24}} = 2.28_{ - 0.12}^{ + 0.14} < 3.83_{ - 0.18}^{ + 0.17} = \log {g_{{\rm{DG}}}}.$

The lower surface gravity in FPTX24 stems from differences in the temperature profile. The 12C/13C retrieved with FPTX24 (68−10+12$68_{ - 10}^{ + 12}$) is higher than our nominal value (53−6+7$53_{ - 6}^{ + 7}$) but consistent within uncertainties. Similarly, the FPTX24 C/O ratio (0.57 ± 0.02) exceeds the DG model’s value (0.50 ± 0.01). These results suggest that distinct temperature profiles impact species abundances and spectral interpretations. However, the extent of the discrepancy does not fully account for the difference in12 C/13C values reported by X24 and this work.

7.3 The carbon isotope ratio of the GQ Lup system

The 12C/13C has been proposed as a potential tracer of formation pathways of BDs and SJs. In the presence of ice accretion, giant planets might be enriched in 12C/13C compared to its hoststar (Zhang et al. 2021a). The measured 12C/13C ratios in the GQ Lup system suggest that GQ Lup B is not significantly enriched in 13C compared to GQ Lup A. This result is consistent with the formation of GQ Lup B through gravitational collapse or disc instability (Stamatellos & Whitworth 2009; Stolker et al. 2021). The 12C/13C of GQ Lup A/B is slightly lower than the local ISM ∼ 69 (Wilson 1999) but broadly consistent within the 1σ scatter (12C/13C)ISM = 68 ± 15 (Milam et al. 2005). The measured12 C/13 C of GQ Lup B lies between the reported 13CO atmosphere of YSES 1 b (30−7+10$30_{ - 7}^{ + 10}$; Zhang et al. 2021a) and VHS 1256 b (62 ± 2; Gandhi et al. 2023). Additional measurements of the12 C/13 C of young SJs and BDs are required to understand the diversity of formation pathways and the role of 12C/13C in the formation of substellar objects.

To assess the role of 12C/13C as a formation tracer it is critical to measure the 12C/13C of the host star and the companion. A high 13CO abundance can be interpreted as a signature of ice accretion during the formation of the object (Zhang et al. 2021a); however, it should be compared to the 12C/13C of the host star to understand the formation history of the system and whether there is any actual enrichment due to accretion processes. Recent work by Xuan et al. (2024b) presented a homogeneous isotopic composition for the binary system HIP 555077, where the 12C/13C of the host star and the companion are consistent with each other. We report a similar consistency in the 12C/13C of GQ Lup A and GQ Lup B, suggesting that the system formed from the same parent cloud that was somewhat enriched in 13CO compared to the present-day local ISM. This result is cohesive with the view that the 12C/13C of young objects should have a higher abundance of 13CO compared to older objects due to the enrichment of 13CO in the ISM over time (e.g. Prantzos et al. 1996; Milam et al. 2005). In the first results of the SupJup survey, de Regt et al. (2024) reported the carbon isotope ratio of an old T dwarf (DENIS 0255, 12C/13C=184−61+40$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 184_{ - 61}^{ + 40}$), which is significantly higher than the ISM. On the other hand, the 12C/13C of three young BDs (age ≤ 10 Myr) as presented in González Picos et al. (2024) are consistent with the ISM, with values ranging from ≈80–120.

thumbnail Fig. 7

Carbon-to-oxygen ratio of GQ Lup B compared to other young substellar companions (<100 Myr). The solar value and its uncertainty are shown as a grey band. The size of the markers is proportional to the mass. The data used to create this figure and their uncertainties can be found in Petrus et al. (2021), Zhang et al. (2021a), Lacour et al. (2021), Dupuy et al. (2023), Gandhi et al. (2023), Hoch et al. (2023), Brown-Sevilla et al. (2023), Mesa et al. (2023), Palma-Bifani et al. (2023), Palma-Bifani et al. (2024), Landman et al. (2024) and Costes et al. (2024).

7.4 The carbon-to-oxygen ratio

The carbon-to-oxygen ratio (C/O) has been proposed as a chemical indicator of formation pathways for substellar objects (e.g. Öberg et al. 2011; Fortney 2012; Hoch et al. 2023). Objects formed within protoplanetary discs may exhibit a C/O ratio exceeding the solar value (around 0.59) due to water freeze-out onto dust grains, enriching the disc’s gas-phase C/O (Öberg et al. 2011; Mordasini et al. 2016; Brewer et al. 2017; Cridland et al. 2019). Conversely, a C/O ratio closer to the solar value might suggest formation through mechanisms like gravitational collapse or disc instability, where the C/O reflects the bulk composition of the parent cloud (Boss 2011).

However, interpreting atmospheric C/O in substellar objects remains challenging due to potential condensation of oxygenbearing molecules (e.g. Line et al. 2015; Calamari et al. 2022). This effect is particularly pronounced for cooler objects than GQ Lup B, such as VHS 1256 b (Gandhi et al. 2023), where a fraction of the oxygen is expected to be sequestered in clouds, leading to an enhanced C/O ratio in the gas phase.

Our analysis reveals a C/O ratio for GQ Lup B of 0.50 ± 0.01, consistent with previous measurements from medium-resolution spectroscopy (0.44−0.11+0.13$0.44_{ - 0.11}^{ + 0.13}$; Demars et al. 2023). This finding aligns with a formation scenario via gravitational collapse or disc instability for GQ Lup B. Furthermore, our measurement for GQ Lup B falls within the range observed for other young sub- stellar objects (see Fig. 7), such as YSES 1 b (0.52−0.03+0.04$0.52_{ - 0.03}^{ + 0.04}$; Zhang et al. 2021a) and HD 984 B (0.50−0.01+0.01$0.50_{ - 0.01}^{ + 0.01}$; Costes et al. 2024). We note that the reported value from Xuan et al. (2024a) (0.70 ± 0.01) is higher than the one reported here, as discussed in the previous section, this may be attributed to a different model interpretation. While C/O alone may not be a definitive formation tracer, it provides valuable insights into the atmospheric chemistry and complements isotope ratio studies in piecing together the formation history of substellar objects.

7.5 Non-adiabatic convection in GQ Lup B

The retrieved temperature profile of GQ Lup B shows a slight deviation from the adiabatic temperature gradient in the lower atmosphere, suggesting the presence of non-adiabatic convection (Tremblin et al. 2015, 2019; Petrus et al. 2021). This result underscores the importance of considering non-adiabatic convection in the interpretation of spectra of SJs and BDs, as it can significantly affect the temperature structure and the derived abundances. The presence of non-adiabatic temperature profiles has been reported in other objects, such as the cool Y dwarf WISE 0359 (Leggett & Tremblin 2023).

7.6 Sensitivity of the surface gravity to the temperature profile

Atmospheric retrievals using different parameterisations of the PT profile reveal that small differences in temperature gradients can significantly influence the retrieved surface gravity. The DG model offers a more flexible representation of the atmospheric temperature structure by allowing the photosphere’s position to be fitted, avoiding biases associated with predefined pressure levels. In contrast, the SG model constrains the surface gravity to log g = 3.14 ± 0.12, whereas the DG model retrieves a higher surface gravity of log g = 3.83 ± 0.18. The DG model is preferred due to its more adaptable temperature structure and superior fit to the data, as evidenced by a 4.5σ preference and a lower reduced chi-squared value. The surface gravity retrieved with the DG model aligns well with previous literature measurements, suggesting it better captures the true atmospheric temperature structure of GQ Lup B (Seifahrt et al. 2007; Lavigne et al. 2009), also consistent with the expected surface gravity for young substellar objects (Baraffe et al. 2002).

Using evolutionary models, Stolker et al. (2021) and Xuan et al. (2024a) estimated a surface gravity of log g ≈ 3.8. While this value agrees with our retrieval, uncertainties in the system’s age and limitations in evolutionary models may introduce discrepancies in the predicted surface gravity.

7.7 Veiling of the GQ LupA spectrum

The K-band veiling factor of GQ Lup A is consistent with a rapid increase with wavelength including the measurement from Sousa et al. (2023) (see Fig. 6). The veiling contribution is an essential parameter for the accurate modelling of the spectrum of GQ Lup A, as it affects the depth of the absorption lines (see Fig. D.1).

Previous attempts to link the observed veiling to disc emission have found black body temperatures that are too high for dust to survive in the inner disc (Fischer et al. 2011; Alcalá et al. 2021). Alternative descriptions of veiling invoke the presence of hot spots on the stellar surface or accretion emission (Hartigan et al. 1989; McClure et al. 2013; Kidder et al. 2021). Observations at longer wavelengths (e.g JWST/MIRI; Cugno et al. 2024) might be able provide insights into the origin of veiling in GQ Lup A.

8 Conclusions

In this study we analysed the carbon isotope ratio (12C/13C) in the GQ Lup system and the carbon-to-oxygen ratio (C/O) of its companion, GQ Lup B, to investigate their formation pathways. Our findings suggest that GQ Lup B formed through gravitational collapse or disc instability, as indicated by the homogeneous12 C/13 C of the system and the lack of enhanced carbon elemental abundance. The key conclusions of this study are as follows:

  • Carbon isotope ratio: the 12C/13C of GQ Lup B (12C/13C=53−6+7$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 53_{ - 6}^{ + 7}$) is consistent with that of GQ Lup A (12C/13C=51−8+10$^{12}{\rm{C}}{/^{13}}{\rm{C}} = 51_{ - 8}^{ + 10}$), indicating a shared origin from the same parent cloud with no evidence of 13CO enrichment. This supports the idea of formation through gravitational collapse or disc instability, rather than significant solid accretion.

  • Comparison with other systems: the 12C/13C of GQ Lup B falls within the range observed for other young sub- stellar objects but is slightly lower than the ISM value ((12C/13C)ISM = 68 ± 15). This suggests potential diversity in 12C/13C among substellar objects of different ages and spectral types.

  • Carbon-to-oxygen ratio: the C/O ratio of GQ Lup B (0.50 ± 0.01) is consistent with the lower bound of the solar value. The absence of the higher C/O ratios predicted for formation within a protoplanetary disc supports an origin via gravitational collapse or disc instability. However, interpreting the C/O ratio in the context of planet formation is challenging due to the complex interplay between cloud condensation and atmospheric abundances. This underscores the importance of further studies using isotope ratios as complementary tracers of formation pathways.

  • Temperature profiles: we introduce a novel retrieval method that uses dynamic pressure-temperature knots to provide flexible temperature profiles and minimise biases in retrieved parameters such as surface gravity. The retrieved temperature profile of GQ Lup B suggests the possible presence of non-adiabatic convection in the lower atmosphere, as evidenced by a reduced temperature gradient.

  • Impact on isotope ratios: the accurate determination of 12C/13C depends on precise temperature profiles, as the 13CO and 12CO lines arise from different atmospheric altitudes. This study demonstrates how variations in temperature profiles can lead to differing interpretations of spectral features, which affects retrieved surface gravities, abundances, and isotope ratios.

  • Veiling in GQ Lup A: veiling observed in the K -band spectrum of GQ Lup A indicates an additional emitting source or photospheric spots on the stellar surface. Veiling was modelled as a wavelength-dependent excess continuum, with the veiling factor fitted independently for each order-detector pair. Observations at longer wavelengths (>2.5 µm) could help determine the physical origin of veiling. Results show that the veiling factor increases with wavelength, consistent with previous measurements by Sousa et al. (2023).

Our results highlight the importance of combining isotopic and elemental ratios to decipher the complex formation histories of substellar objects. Further measurements spanning different ages and spectral types are crucial for advancing our understanding of substellar formation pathways. This study contributes to the broader effort of characterising young planetary and substellar systems by leveraging isotope ratios. We stress the necessity of including isotope measurements of host stars to fully understand the formation history of widely separated companions. Specifically, we present a method for measuring the carbon isotope ratio of a K-type host star, incorporating veiling into the analysis.

The next generation of high-resolution spectrographs will enable the study of substellar companions at smaller separations from their host stars and at fainter magnitudes. The methods outlined in this paper for accounting for starlight contamination can be applied to extract atmospheric properties from such systems. Upcoming mid-infrared spectrographs such as METIS (Brandl et al. 2021) at the European extremely large telescope (Gilmozzi & Spyromilio 2007) will enable carbon isotope measurements at longer wavelengths. In this regime, understanding veiling and the effects of disc emission on the observed spectrum will be critical for studying young systems.

Acknowledgements

We thank J. Xuan for providing the temperature profile of GQ Lup B retrieved with KPIC. We thank P. Hauschildt for generating the isotope-dependent PHOENIX models. We are grateful to the anonymous referee for their constructive comments. 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 pro- gramme(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), peti- tRADTRANS (Mollière et al. 2019), PyAstronomy (Czesla et al. 2019), Astropy (Astropy Collaboration 2022), corner (Foreman-Mackey 2016).

Appendix A Observing conditions

thumbnail Fig. A.1

Observing conditions over a 3-hour science exposure sequence. Airmass, integrated water vapour, and seeing are plotted in the top, middle, and bottom panels, respectively. Average values are indicated by dashed black lines. Pink data points represent standard star observations, taken immediately before the science sequence. The retrieved effective airmass, with 1σ uncertainty indicated by the line width, is overlaid in the airmass panel.

Appendix B Linear starlight model

The starlight model is constructed using a linear combination of the observed spectrum of GQ Lup A, where individual components are shifted versions of the on-axis spectrum. In addition, we accounted for low-frequency variations in the scattered starlight using a spline decomposition method, as presented in Ruffio et al. 2023. The linear model is constructed and fitted for each order-detector pair and nodding position separately, hence providing an efficient method for modelling the starlight spectrum for each chunk of the data without the need to combine data with potential different systematics and noise properties.

thumbnail Fig. B.1

Components of the starlight model. Left panel: Broadening kernel. The black spectrum represents the combined spectrum with equally weighted components. Right panel: Spline decomposition of the starlight model. The upper part displays the resulting spectra corresponding to different amplitudes, and the lower part the individual spline components. The linear amplitudes used to generate the spectra in the upper part are indicated by markers in the lower part.

Appendix C Extended results of the atmospheric retrieval of GQ Lup B

Table C.1

Description of the free parameters of the retrieval, prior ranges, and best-fit values.

thumbnail Fig. C.1

Same as the top panel of Fig. 2, but for all spectral orders.
thumbnail Fig. C.2

Posterior distributions of the best-fit parameters of GQ Lup B. The uncertainties indicate the 16th and 84th percentiles of the distributions.
thumbnail Fig. C.3

Contribution function of the full best-fit model of GQ Lup B (top panel). The individual contribution functions of 12CO and 13CO are displayed for the last two spectral orders on the bottom panels.

Appendix D Veiling in the spectrum of GQ Lup A

thumbnail Fig. D.1

Observed spectrum of GQ Lup A (black) for the spectral order centred at 2345 nm, compared with a PHOENIX template (Teff = 4300 K, log g = 4.0) that includes instrumental and rotational broadening. The mismatch between the observed line depths and the model without veiling is evident. The second panel shows the best-fit linear model (green), which incorporates a veiling continuum (pink) alongside the PHOENIX model (orange). Residuals of the best-fit model are presented in the third panel.

Appendix E Posterior comparison with fixed temperature profile from X24

thumbnail Fig. E.1

Comparison of posterior distributions for selected best-fit parameters of GQ Lup B, highlighting differences between our fiducial model (with DG) and the fixed PT profile model from Xuan et al. 2024a (X24). Dashed lines in the corner plot represent the 16th and 84th percentiles of the distributions. Derived C/O and 12C/13C values are displayed in the top-right corner, with the dashed grey probability density indicating the values reported by X24. The middle-right panel shows the retrieved temperature profiles, while the rightmost panel presents the corresponding temperature gradients.

References

  1. Alcalá, J. M., Majidi, F. Z., Desidera, S., et al. 2020, A&A, 635, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Alcalá, J. M., Gangi, M., Biazzo, K., et al. 2021, A&A, 652, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Allard, N. F., Spiegelman, F., Leininger, T., & Molliere, P. 2019, A&A, 628, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111 [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. Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65 [Google Scholar]
  10. Batalha, C., Lopes, D. F., & Batalha, N. M. 2001, ApJ, 548, 377 [CrossRef] [Google Scholar]
  11. Benneke, B., & Seager, S. 2013, ApJ, 778, 153 [Google Scholar]
  12. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bitsch, B., Izidoro, A., Johansen, A., et al. 2019, A&A, 623, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, MNRAS, 492, 431 [Google Scholar]
  15. Boley, A. C., Payne, M. J., & Ford, E. B. 2012, ApJ, 754, 57 [Google Scholar]
  16. Borysow, J., Frommhold, L., & Birnbaum, G. 1988, ApJ, 326, 509 [NASA ADS] [CrossRef] [Google Scholar]
  17. Boss, A. P. 1997, Science, 276, 1836 [Google Scholar]
  18. Boss, A. P. 2011, ApJ, 731, 74 [Google Scholar]
  19. 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 (Tucson: University odf Arizona Press), 433 [Google Scholar]
  20. Brandl, B., Bettonvil, F., van Boekel, R., et al. 2021, The Messenger, 182, 22 [NASA ADS] [Google Scholar]
  21. Brewer, J. M., Fischer, D. A., & Madhusudhan, N. 2017, AJ, 153, 83 [NASA ADS] [CrossRef] [Google Scholar]
  22. Broeg, C., Schmidt, T. O. B., Guenther, E., et al. 2007, A&A, 468, 1039 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Brown-Sevilla, S. B., Maire, A. L., Mollière, P., et al. 2023, A&A, 673, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Calamari, E., Faherty, J. K., Burningham, B., et al. 2022, ApJ, 940, 164 [NASA ADS] [CrossRef] [Google Scholar]
  25. Carter, E. J., & Stamatellos, D. 2023, MNRAS, 525, 1912 [NASA ADS] [CrossRef] [Google Scholar]
  26. Carvalho, A., & Johns-Krull, C. M. 2023, Res. Notes Am. Astron. Soc., 7, 91 [Google Scholar]
  27. Castelli, F., & Kurucz, R. L. 2003, IAU Symp., 210, A20 [Google Scholar]
  28. Chan, Y. M., & Dalgarno, A. 1965, Proc. Phys. Soc., 85, 227 [NASA ADS] [CrossRef] [Google Scholar]
  29. Chauvin, G., Lagrange, A. M., Dumas, C., et al. 2005, A&A, 438, L25 [CrossRef] [EDP Sciences] [Google Scholar]
  30. Cochran, W. G. 1952, Annal. Math. Stat., 23, 315 [CrossRef] [Google Scholar]
  31. Costes, J. C., Xuan, J. W., Vigan, A., et al. 2024, A&A, 686, A294 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (Cambridge: Cambridge Scientific Publishers) [Google Scholar]
  33. Cridland, A. J., Eistrup, C., & van Dishoeck, E. F. 2019, A&A, 627, A127 [Google Scholar]
  34. Cugno, G., Patapis, P., Banzatti, A., et al. 2024, ApJ, 966, L21 [NASA ADS] [CrossRef] [Google Scholar]
  35. Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, Astrophysics Source Code Library [record ascl:1906.010] [Google Scholar]
  36. Dai, Y., Wilner, D. J., Andrews, S. M., & Ohashi, N. 2010, AJ, 139, 626 [NASA ADS] [CrossRef] [Google Scholar]
  37. Dalgarno, A., & Williams, D. A. 1962, ApJ, 136, 690 [NASA ADS] [CrossRef] [Google Scholar]
  38. Delorme, J.-R., Jovanovic, N., Echeverri, D., et al. 2020, SPIE Conf. Ser., 11447, 114471P [NASA ADS] [Google Scholar]
  39. Demars, D., Bonnefoy, M., Dougados, C., et al. 2023, A&A, 676, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. de Regt, S., Gandhi, S., Snellen, I. A. G., et al. 2024, A&A, 688, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Dodson-Robinson, S. E., Veras, D., Ford, E. B., & Beichman, C. A. 2009, ApJ, 707, 79 [Google Scholar]
  42. Donati, J. F., Gregory, S. G., Alencar, S. H. P., et al. 2012, MNRAS, 425, 2948 [Google Scholar]
  43. Dorn, R. J., Bristow, P., Smoker, J. V., et al. 2023, A&A, 671, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Dupuy, T. J., Liu, M. C., Evans, E. L., et al. 2023, MNRAS, 519, 1688 [Google Scholar]
  45. Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016, ApJS, 225, 10 [Google Scholar]
  46. Fischer, W., Edwards, S., Hillenbrand, L., & Kwan, J. 2011, ApJ, 730, 73 [Google Scholar]
  47. Foreman-Mackey, D. 2016, J. Open Source Softw., 1, 24 [Google Scholar]
  48. Fortney, J. J. 2012, ApJ, 747, L27 [NASA ADS] [CrossRef] [Google Scholar]
  49. Fortney, J. J., Marley, M. S., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Astron. Nachr., 326, 925 [NASA ADS] [CrossRef] [Google Scholar]
  50. Gandhi, S., de Regt, S., Snellen, I., et al. 2023, ApJ, 957, L36 [NASA ADS] [CrossRef] [Google Scholar]
  51. Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001 [NASA ADS] [CrossRef] [Google Scholar]
  52. Gilmozzi, R., & Spyromilio, J. 2007, The Messenger, 127, 3 [NASA ADS] [Google Scholar]
  53. Ginski, C., Schmidt, T. O. B., Mugrauer, M., et al. 2014, MNRAS, 444, 2280 [Google Scholar]
  54. González Picos, D., Snellen, I. A. G., de Regt, S., et al. 2024, A&A, 689, A212 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Gray, D. F. 2022, The Observation and Analysis of Stellar Photospheres (Cambridge: Cambridge University Press) [Google Scholar]
  56. Guenther, E. W., Neuhäuser, R., Wuchterl, G., et al. 2005, Astron. Nachr., 326, 958 [CrossRef] [Google Scholar]
  57. Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar Interiors: Physical Principles, Structure, and Evolution (Berlin: Springer) [Google Scholar]
  58. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  59. Hartigan, P., Hartmann, L., Kenyon, S., Hewett, R., & Stauffer, J. 1989, ApJS, 70, 899 [NASA ADS] [CrossRef] [Google Scholar]
  60. Herbig, G. H. 1962, Adv. Astron. Astrophys., 1, 47 [Google Scholar]
  61. Hoch, K. K. W., Konopacky, Q. M., Theissen, C. A., et al. 2023, AJ, 166, 85 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hoeijmakers, H. J., Schwarz, H., Snellen, I. A. G., et al. 2018, A&A, 617, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Holmberg, M., & Madhusudhan, N. 2022, AJ, 164, 79 [NASA ADS] [CrossRef] [Google Scholar]
  64. Horne, K. 1986, PASP, 98, 609 [Google Scholar]
  65. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  66. Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Inaba, S., Wetherill, G. W., & Ikoma, M. 2003, Icarus, 166, 46 [Google Scholar]
  68. Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, J. Quant. Spec. Radiat. Transf., 109, 1136 [NASA ADS] [CrossRef] [Google Scholar]
  69. Iyer, A. R., Line, M. R., Muirhead, P. S., Fortney, J. J., & Gharib-Nezhad, E. 2023, ApJ, 944, 41 [NASA ADS] [CrossRef] [Google Scholar]
  70. Janson, M., Brandner, W., Henning, T., & Zinnecker, H. 2006, A&A, 453, 609 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Johansen, A., & Bitsch, B. 2019, A&A, 631, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [Google Scholar]
  73. Kawahara, H., Kawashima, Y., Masuda, K., et al. 2022, ApJS, 258, 31 [NASA ADS] [CrossRef] [Google Scholar]
  74. Kidder, B., Mace, G., López-Valdivia, R., et al. 2021, ApJ, 922, 27 [NASA ADS] [CrossRef] [Google Scholar]
  75. Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398 [Google Scholar]
  76. Kraus, A. L., & Ireland, M. J. 2012, ApJ, 745, 5 [Google Scholar]
  77. Lacour, S., Wang, J. J., Rodet, L., et al. 2021, A&A, 654, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Lafrenière, D., Doyon, R., Marois, C., et al. 2007, ApJ, 670, 1367 [Google Scholar]
  79. Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Landman, R., Stolker, T., Snellen, I. A. G., et al. 2024, A&A, 682, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Langer, W. D., & Penzias, A. A. 1993, ApJ, 408, 539 [NASA ADS] [CrossRef] [Google Scholar]
  82. Lavigne, J.-F., Doyon, R., Lafrenière, D., Marois, C., & Barman, T. 2009, ApJ, 704, 1098 [Google Scholar]
  83. Lawson, C. L., & Hanson, R. J. 1995, Solving Least Squares Problems (USA: SIAM) [CrossRef] [Google Scholar]
  84. Lazzoni, C., Desidera, S., Gratton, R., et al. 2022, MNRAS, 516, 391 [NASA ADS] [CrossRef] [Google Scholar]
  85. Lee, S., Nomura, H., & Furuya, K. 2024, ApJ, 969, 41 [NASA ADS] [CrossRef] [Google Scholar]
  86. Leggett, S. K., & Tremblin, P. 2023, ApJ, 959, 86 [NASA ADS] [CrossRef] [Google Scholar]
  87. Li, G., Gordon, I. E., Rothman, L. S., et al. 2015, ApJS, 216, 15 [NASA ADS] [CrossRef] [Google Scholar]
  88. Line, M. R., Zhang, X., Vasisht, G., et al. 2012, ApJ, 749, 93 [Google Scholar]
  89. Line, M. R., Teske, J., Burningham, B., Fortney, J. J., & Marley, M. S. 2015, ApJ, 807, 183 [NASA ADS] [CrossRef] [Google Scholar]
  90. Lyons, J. R., Gharib-Nezhad, E., & Ayres, T. R. 2018, Nat. Commun., 9, 908 [Google Scholar]
  91. MacGregor, M. A., Wilner, D. J., Czekala, I., et al. 2017, ApJ, 835, 17 [Google Scholar]
  92. Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci., 111, 12661 [NASA ADS] [CrossRef] [Google Scholar]
  93. Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64 [Google Scholar]
  94. Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24 [NASA ADS] [CrossRef] [Google Scholar]
  95. Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279 [Google Scholar]
  96. Marois, C., Macintosh, B., & Barman, T. 2007, ApJ, 654, L151 [NASA ADS] [CrossRef] [Google Scholar]
  97. Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348 [Google Scholar]
  98. Mayer, L., Quinn, T., Wadsley, J., & Stadel, J. 2002, Science, 298, 1756 [NASA ADS] [CrossRef] [Google Scholar]
  99. McClure, M. K., Calvet, N., Espaillat, C., et al. 2013, ApJ, 769, 73 [NASA ADS] [CrossRef] [Google Scholar]
  100. McElwain, M. W., Metchev, S. A., Larkin, J. E., et al. 2007, ApJ, 656, 505 [Google Scholar]
  101. Mesa, D., Gratton, R., Kervella, P., et al. 2023, A&A, 672, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126 [Google Scholar]
  103. Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67 [Google Scholar]
  104. Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131 [Google Scholar]
  105. Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41 [Google Scholar]
  106. Neuhäuser, R., Guenther, E. W., Wuchterl, G., et al. 2005, A&A, 435, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Neuhäuser, R., Mugrauer, M., Seifahrt, A., Schmidt, T. O. B., & Vogt, N. 2008, A&A, 484, 281 [CrossRef] [EDP Sciences] [Google Scholar]
  108. Öberg, K. I., & Bergin, E. A. 2021, Phys. Rep., 893, 1 [Google Scholar]
  109. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [Google Scholar]
  110. Palma-Bifani, P., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Palma-Bifani, P., Chauvin, G., Borja, D., et al. 2024, A&A, 683, A214 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Patience, J., King, R. R., De Rosa, R. J., et al. 2012, A&A, 540, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Paufique, J., Biereichel, P., Donaldson, R., et al. 2004, SPIE Conf. Ser., 5490, 216 [NASA ADS] [Google Scholar]
  114. Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59 [EDP Sciences] [Google Scholar]
  115. Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  117. Polyansky, O. L., Kyuberis, A. A., Lodi, L., et al. 2017, MNRAS, 466, 1363 [NASA ADS] [CrossRef] [Google Scholar]
  118. Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597 [NASA ADS] [CrossRef] [Google Scholar]
  119. Prantzos, N., Aubert, O., & Audouze, J. 1996, A&A, 309, 760 [NASA ADS] [Google Scholar]
  120. Rei, A. C. S., Petrov, P. P., & Gameiro, J. F. 2018, A&A, 610, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Rodgers, C. D. 2000, Inverse Methods For Atmospheric Sounding: Theory And Practice (Singapore: World Scientific) [CrossRef] [Google Scholar]
  122. Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spec. Radiat. Transf., 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
  123. Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200 [Google Scholar]
  124. Ruffio, J.-B., Konopacky, Q. M., Barman, T., et al. 2021, AJ, 162, 290 [NASA ADS] [CrossRef] [Google Scholar]
  125. Ruffio, J.-B., Horstman, K., Mawet, D., et al. 2023, AJ, 165, 113 [NASA ADS] [CrossRef] [Google Scholar]
  126. Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Seifahrt, A., Neuhäuser, R., & Hauschildt, P. H. 2007, A&A, 463, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Nature, 509, 63 [Google Scholar]
  129. Sousa, A. P., Bouvier, J., Alencar, S. H. P., et al. 2023, A&A, 670, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174 [Google Scholar]
  131. Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392, 413 [Google Scholar]
  132. Stolker, T., Haffert, S. Y., Kesseli, A. Y., et al. 2021, AJ, 162, 286 [NASA ADS] [CrossRef] [Google Scholar]
  133. Tachihara, K., Dobashi, K., Mizuno, A., Ogawa, H., & Fukui, Y. 1996, PASJ, 48, 489 [NASA ADS] [CrossRef] [Google Scholar]
  134. Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17 [CrossRef] [Google Scholar]
  135. Tremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144 [Google Scholar]
  136. Veras, D., Crepp, J. R., & Ford, E. B. 2009, ApJ, 696, 1600 [Google Scholar]
  137. Vigan, A., Langlois, M., Moutou, C., & Dohlen, K. 2008, A&A, 489, 1345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  139. Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Vorobyov, E. I. 2013, A&A, 552, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Whiteford, N., Glasse, A., Chubb, K. L., et al. 2023, MNRAS, 525, 1375 [NASA ADS] [CrossRef] [Google Scholar]
  142. Wilcomb, K. K., Konopacky, Q. M., Barman, T. S., et al. 2020, AJ, 160, 207 [Google Scholar]
  143. Wilson, T. L. 1999, Rep. Prog. Phys., 62, 143 [Google Scholar]
  144. 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]
  145. Wu, Y.-L., Sheehan, P. D., Males, J. R., et al. 2017, ApJ, 836, 223 [CrossRef] [Google Scholar]
  146. Xuan, J. W., Hsu, C.-C., Finnerty, L., et al. 2024a, ApJ, 970, 71 [CrossRef] [Google Scholar]
  147. Xuan, J. W., Wang, J., Finnerty, L., et al. 2024b, ApJ, 962, 10 [NASA ADS] [CrossRef] [Google Scholar]
  148. Yoshida, T. C., Nomura, H., Furuya, K., Tsukagoshi, T., & Lee, S. 2022, ApJ, 932, 126 [NASA ADS] [CrossRef] [Google Scholar]
  149. Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021a, Nature, 595, 370 [NASA ADS] [CrossRef] [Google Scholar]
  150. Zhang, Y., Snellen, I. A. G., & Mollière, P. 2021b, A&A, 656, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  151. Zhang, Y., Snellen, I. A. G., Brogi, M., & Birkby, J. L. 2022, Res. Notes Am. Astron. Soc., 6, 194 [Google Scholar]
  152. Zhang, Z., Mollière, P., Hawkins, K., et al. 2023, AJ, 166, 198 [NASA ADS] [CrossRef] [Google Scholar]
  153. Zhang, Y., González Picos, D., de Regt, S., et al. 2024, AJ, 168, 246 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Properties of GQ Lup A and B.

In the text
Table 2

Opacity species and references.

In the text
Table C.1

Description of the free parameters of the retrieval, prior ranges, and best-fit values.

In the text

All Figures

thumbnail Fig. 1

Calibrated 2D image of a single order-detector pair. The vertical axis is the spatial direction and the horizontal axis the spectral direction. The bright central source is GQ Lup A, and the fainter source at −0.71 arcseconds is GQ Lup B. The telluric lines are visible as vertical stripes in the image. The right panel shows the integrated line spread function and the spatial profile employed for the extraction of the data at the position of GQ Lup B.
In the text
thumbnail Fig. 2

Best-fit model of the spectrum at the position of GQ Lup B. Top: observed spectrum (black), joint best-fit model (green) and components of the starlight model and companion (blue and brown, respectively). Centre: data with the best-fit starlight model subtracted compared to the companion model. Bottom: residuals of the best-fit model. The shaded region represents the root-squared diagonal of the covariance matrix with the optimal noise scaling factor.
In the text
thumbnail Fig. 3

Temperature profile of GQ Lup B retrieved using the SG and DG models. Left: temperature profile. Right: temperature gradient. Profiles from the ATMO grid (Tremblin et al. 2015) with Teff=2300 K and log g=4.0 are shown for two distinct models with adiabatic (Phillips et al. 2020) and diabatic convection (Petrus et al. 2021).
In the text
thumbnail Fig. 4

Cross-correlation functions (CCFs) of the data and the best- fit models of 13CO for GQ Lup B (top) and GQ Lup A (bottom). The CCFs are shown in black, with the auto-correlation functions of the templates shown in brown and orange. The CCFs are calculated between the data residuals and the model residuals following the methodology from Zhang et al. (2021a) (also in de Regt et al. 2024).
In the text
thumbnail Fig. 5

Best-fit model of GQ Lup A. a. Spectrum with the veiling continuum subtracted and divided by the telluric model. The observed spectrum is shown in black, the PHOENIX model with optimal linear amplitudes in orange and the 13CO lines from the best-fit model in purple. The residuals of the are shown in the bottom panel, the shaded grey region indicates the root-squared diagonal of the covariance matrix with the GP kernel. b. Posterior distributions of the carbon isotope ratio, projected rotational velocity and radial velocity of GQ Lup A.
In the text
thumbnail Fig. 6

Veiling factor of GQ Lup A as a function of wavelength. The veiling factors for each spectral order derived from the fitted linear amplitudes are shown in orange, with error bars representing the standard deviation between the three detectors of each order. The measurement from Sousa et al. (2023) is included to fit the three points with a linear function of slope 0.011 nm−1.
In the text
thumbnail Fig. 7

Carbon-to-oxygen ratio of GQ Lup B compared to other young substellar companions (<100 Myr). The solar value and its uncertainty are shown as a grey band. The size of the markers is proportional to the mass. The data used to create this figure and their uncertainties can be found in Petrus et al. (2021), Zhang et al. (2021a), Lacour et al. (2021), Dupuy et al. (2023), Gandhi et al. (2023), Hoch et al. (2023), Brown-Sevilla et al. (2023), Mesa et al. (2023), Palma-Bifani et al. (2023), Palma-Bifani et al. (2024), Landman et al. (2024) and Costes et al. (2024).
In the text
thumbnail Fig. A.1

Observing conditions over a 3-hour science exposure sequence. Airmass, integrated water vapour, and seeing are plotted in the top, middle, and bottom panels, respectively. Average values are indicated by dashed black lines. Pink data points represent standard star observations, taken immediately before the science sequence. The retrieved effective airmass, with 1σ uncertainty indicated by the line width, is overlaid in the airmass panel.
In the text
thumbnail Fig. B.1

Components of the starlight model. Left panel: Broadening kernel. The black spectrum represents the combined spectrum with equally weighted components. Right panel: Spline decomposition of the starlight model. The upper part displays the resulting spectra corresponding to different amplitudes, and the lower part the individual spline components. The linear amplitudes used to generate the spectra in the upper part are indicated by markers in the lower part.
In the text
thumbnail Fig. C.1

Same as the top panel of Fig. 2, but for all spectral orders.
In the text
thumbnail Fig. C.2

Posterior distributions of the best-fit parameters of GQ Lup B. The uncertainties indicate the 16th and 84th percentiles of the distributions.
In the text
thumbnail Fig. C.3

Contribution function of the full best-fit model of GQ Lup B (top panel). The individual contribution functions of 12CO and 13CO are displayed for the last two spectral orders on the bottom panels.
In the text
thumbnail Fig. D.1

Observed spectrum of GQ Lup A (black) for the spectral order centred at 2345 nm, compared with a PHOENIX template (Teff = 4300 K, log g = 4.0) that includes instrumental and rotational broadening. The mismatch between the observed line depths and the model without veiling is evident. The second panel shows the best-fit linear model (green), which incorporates a veiling continuum (pink) alongside the PHOENIX model (orange). Residuals of the best-fit model are presented in the third panel.
In the text
thumbnail Fig. E.1

Comparison of posterior distributions for selected best-fit parameters of GQ Lup B, highlighting differences between our fiducial model (with DG) and the fixed PT profile model from Xuan et al. 2024a (X24). Dashed lines in the corner plot represent the 16th and 84th percentiles of the distributions. Derived C/O and 12C/13C values are displayed in the top-right corner, with the dashed grey probability density indicating the values reported by X24. The middle-right panel shows the retrieved temperature profiles, while the rightmost panel presents the corresponding temperature gradients.
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.