11institutetext: Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France
11email<EMAIL_ADDRESS>22institutetext: Institute for
Advanced Study, Tsinghua University, Beijing 100084, People’s Republic of
China 33institutetext: Department of Astronomy, Tsinghua University, Beijing
100084, People’s Republic of China 44institutetext: Department of Earth and
Planetary Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo
113-0033, Japan 55institutetext: Department of Astronomy, University of
Michigan, 1085 S. University Ave., Ann Arbor, MI 48109, USA 66institutetext:
Institute for Astrophysical Research and Department of Astronomy, Boston
University, 725 Commonwealth Ave., Boston, MA 02215, USA 77institutetext:
Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, 47057
Duisburg, Germany 88institutetext: Institut für Astronomie und Astrophysik,
Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany
99institutetext: Physikalisches Institut, Universität Bern, Gesellschaftsstr.
6, CH-3012 Bern, Switzerland 1010institutetext: Max-Planck-Institut für
Astronomie, Königstuhl 17, 69117 Heidelberg, Germany 1111institutetext: Maison
de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay,
91191 Gif-sur-Yvette, France 1212institutetext: Université Côte d’Azur, OCA,
Lagrange CNRS, 06304 Nice, France 1313institutetext: Instituto de Física y
Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Av. Gran Bretaña
1111, Valparaíso, Chile 1414institutetext: Núcleo Milenio Formación Planetaria
- NPF, Universidad de Valparaíso, Av. Gran Bretaña 1111, Valparaíso, Chile
1515institutetext: Department of Astronomy, The University of Texas at Austin,
2515 Speedway, Stop C1400, Austin, TX 78712, USA 1616institutetext: LESIA,
Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités,
UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules
Janssen, 92195 Meudon, France

# Emission line variability of young 10–30 $\mathrm{M_{Jup}}$
companions††thanks: Based on observations collected at the European
Organisation for Astronomical Research in the Southern Hemisphere under ESO
program 099.C-0760.

I. The case of GQ Lup b and GSC 06214-00210 b

D. Demars 11 M. Bonnefoy 11 C. Dougados 11 Y. Aoyama 22 3 3 4 4 T.
Thanathibodee 55 6 6 G.-D. Marleau 77 8 8 9 9 10 10 P. Tremblin 1111 P.
Delorme 11 P. Palma-Bifani 1212 S. Petrus 1313 14 14 B. P. Bowler 1515 G.
Chauvin 1212 1 1 A.-M. Lagrange 1616 1 1

(Received …; accepted …)

###### Abstract

Context. Emission lines indicative of active accretion have been seen on a
handful of low-mass companions (M¡30 $\mathrm{M_{Jup}}$) to stars. Line
variability is ubiquitous on stellar accretors but has never been
characterized in detail on low-mass companions and can give insights on the
accretion mechanism at play.

Aims. We aim to characterize the short- to long-term $\mathrm{H\,\textsc{i}}$
Paschen $\upbeta$ emission line variability of two 10 to 30 $\mathrm{M_{Jup}}$
companions on wide orbits – GQ Lup b and GSC 06214-00210 b – to clarify the
accretion mechanisms at play.

Methods. We used J-band observations at R=1800–2360, obtained with VLT/SINFONI
in 2017, to record time series investigating hours-to-week variability of the
$\mathrm{H\,\textsc{i}}$ Paschen $\upbeta$ emission line (1.282 $\upmu$m).
Contrary to $\mathsf{H}_{\upalpha}$, it should be less affected by
chromospheric activity contamination. The photospheric emission is analyzed at
each epoch and removed with the FORMOSA forward modeling tool, using new grids
of ATMO models exploring different C/O and [M/H] values. The time series of
line profiles and intensities are compared to those of more massive accretors
and to predictions from the latest magnetospheric accretion and shock models.
To complement, we also re-investigate archival spectroscopic observations at
near-infrared wavelengths of each target to increase the time span up to a
decade and build a more comprehensive understanding of the variability
processes at play.

Results. For GQ Lup b, we find line variability on months and decade
timescales, while it is within noise level at shorter timescales. For GSC
06214-00210 b, we find line variability on timescales of tens of minutes up to
a decade. The line profiles of GSC 06214-00210 b are partly resolved in at
least one epoch. Both objects show $\mathrm{H\,\textsc{i}}$ Paschen $\upbeta$
flux variability that is moderate ($<50$ %) at timescales below their rotation
period, and more important at longer timescales (up to $\sim 1000$ % on decade
timescales), a behavior similar to that of classical T Tauri stars. The line
profiles of GQ Lup b are blue-shifted and can only be reproduced by
magnetospheric accretion models, while GSC 06214-00210 b’s are fairly well
reproduced by both magnetospheric accretion and shock models, except for one
epoch for which the shock model is highly favored. The companions have C/O
values broadly consistent with solar values.

Conclusions. While magnetospheric accretion is favored for GQ Lup b, higher
resolution ($\textrm{R}>10\,000$) observations are required to disentangle the
two (non-exclusive) emitting mechanisms. The similar variability behavior
observed in these low mass companions and in classical T Tauri stars may
support similar accretion mechanisms. The high amplitude of variability at
months and longer timescales found in both objects could be key to explain the
low yield of $\mathsf{H}_{\upalpha}$ imaging campaigns.

###### Key Words.:

Planets and satellites: formation, individual: GQ Lup b, GSC 06214-00210 b -
Accretion, accretion disks

##  1 Introduction

Gas accretion on forming giant planets and brown-dwarf companions is a
critical step determining their early physical evolution and internal
structure (e.g., Marley et al., 2007; Baraffe et al., 2009; Cumming et al.,
2018) and the dissipation of their angular momentum (e.g., Bryan et al., 2018;
Batygin, 2018). It may also drive disk–planet migration (e.g., Pierens &
Raymond, 2016) and the formation and water content of exo-moons (e.g., Heller
& Pudritz, 2015) which are possible habitats for life.

The theoretical understanding of the accretion processes is still in its
infancy. Companions should accrete part of their mass from circumplanetary
disks (hereafter CPDs, Miki, 1982) forming around them. Boundary-layer
accretion is proposed to develop if the CPD hits the planet surface (e.g.,
Dong et al., 2021; Fu et al., 2023, and ref. therein) producing a shock within
a thin boundary at the disk-planet interface. The most recent radiative-
hydrodynamic (RHD) simulations of planets embedded in circumplanetary-disks
predict that part of the mass should fall onto the planet embryo from the
circumstellar disk via high-altitude flows and producing an accretion shock at
the CPD surface or directly at the planet surface (e.g., Szulágyi et al.
(2014), or Fig. 1 in Szulágyi & Ercolano (2020) for a representation).

The picture of accretion may in fact not be fully captured by the RHD
approach. The strength of the dipole component of Jupiter’s current magnetic
field is only about 4 Gauss (Smith et al., 1974), but the interior of young
exoplanets and low-mass brown-dwarfs are predicted to be strongly convective
and ionized, and so be able to fuel $\sim$ kG magnetic fields (Reiners &
Christensen, 2010) sufficient to affect the gas flow (Cridland, 2018) and
truncate the CPD (Fendt, 2003). It thus remains to be known to what extent the
magnetospheric accretion framework applies in the planetary-mass regime.

The heated gas, either at the accretion shock or within the funnels, produces
a set of emission lines from the UV to the NIR whose intensity and spectral
profiles start to be predicted by accretion models (e.g., Thanathibodee et
al., 2019; Szulágyi & Ercolano, 2020; Aoyama et al., 2020; Marleau & Aoyama,
2022). These lines have now been found on a handful of young directly-imaged
young companions with mass below $\sim$30$\mathrm{M_{Jup}}$ (TWA 5 B, GQ Lup
B, CT Cha B, USco CTIO 108 B, DH Tau B, GSC 06214-00210 B, SR 12 C, PDS 70 b
and c, 2M0249c, Delorme 1(AB)b, HIP 77900B, HIP 78530B, USco
161031.9-16191305B, YSES 1b; Neuhäuser et al., 2000; Seifahrt et al., 2007;
Schmidt et al., 2008; Béjar et al., 2008; Bonnefoy et al., 2014; Bowler et
al., 2011; Santamaría-Miranda et al., 2018, 2019; Wagner et al., 2018; Zhou et
al., 2014; Haffert et al., 2019; Chinchilla et al., 2021; Petrus et al., 2020;
Zhang et al., 2021) and are estimated to be frequent ($\sim$50% frequency,
Bowler et al., 2017) below 15 Myr. These companions are found in systems
spanning a wide range of ages (1–40 Myr) and with diverse architectures
(within the cavity of circumstellar disks, in hierarchical systems, or
orbiting single or binary stars of various masses).

All the line detections have also been fortuitous and follow up-studies are
starting to understand their origin and constrain the accretion models. Zhou
et al. (2014) and Zhou et al. (2021) recorded measurements of the shock
emission of four companions in the UV (GQ Lup b, GSC 06214-00210 b, DH Tau b
and PDS 70 b, all within 10–30 $\mathrm{M_{Jup}}$ and $<$15 Myr). They could
measure accretion rates in the range $\mathrm{10^{-9}-10^{-11}}$
$\mathrm{M_{\odot}}\,\mathrm{yr^{-1}}$. They also argue that a significant
fraction of the energy released at the accretion shock is reprocessed at
$\mathsf{H}_{\upalpha}$. Line ratios measurements on Delorme 1 (AB)b ($\sim$14
$\mathrm{M_{Jup}}$, $\sim$40 Myr, Delorme et al., 2013) at near-infrared
wavelengths (Betti et al., 2022b, a) were found to be compatible with non-LTE
accretion-shock models for planets (Aoyama et al., 2018). High-resolution
observations of the Balmer lines of Delorme 1 (AB)b in the UV show line
profiles suggestive of magnetospheric accretion operating on that 40 Myr old
companion (Ringqvist et al., 2023). Constraints on the
$\mathsf{H}_{\upalpha}$/$\mathsf{H}_{\upbeta}$ line-ratio of PDS 70 b indicate
significant extinction in the line-of-sight due to surrounding dust or a CPD
(Hashimoto et al., 2020).

Variability of emission lines is observed across a wide stellar mass spectrum
(for a review, see Fischer et al., 2022), and traces episodic variability in
the accretion and ejection over the course of the stars evolution. Non-steady
accretion is predicted to occur on low-mass companions and has a significant
impact on their early physical evolution (Baraffe et al., 2009; Cumming et
al., 2018) and our ability to detect them (Brittain et al., 2020). Global
hydrodynamic and magneto-hydrodynamic simulations of protoplanets embedded in
a circumstellar disk indicate modulations on decade-long timescales (see
Gressel et al., 2013; Szulágyi & Ercolano, 2020). Instabilities or dead zones
developing through the CPD (Lubow & Martin, 2012; Zhu et al., 2016) should
modulate the accretion flow at Keplerian timescales (months to years). If
magnetospheric accretion operates, hot spots and funnel flows should produce
variability of the emission lines along the object rotation period alike on
stars (e.g., Sicilia-Aguilar et al., 2015). Variable extinction in the line of
sight should modulate the line emission further (Szulágyi & Ercolano, 2020;
Marleau et al., 2022).

The study of the emission line variability timescales and amplitude therefore
offers to probe both the accretion mechanism and the effect of the local
environment of young companions. Line variability has been detected at the
$\sim$year timescale on GQ Lup b and GSC 06214-00210 b (see below). Wolff et
al. (2017) found the $\mathsf{Pa}_{\upbeta}$ line of DH Tau b (a $\sim$1Myr
old $\sim 15$$\mathrm{M_{Jup}}$ companion, Itoh et al., 2005) to disappear
within $\sim 5$weeks. Eriksson et al. (2020) find the $\mathsf{H}_{\upalpha}$
line of Delorme 1(AB)b (a 30–40 Myr old, $\sim 14$$\mathrm{M_{Jup}}$ companion
Delorme et al., 2013) to vary over a few hours. Zhou et al. (2021) report
several measurements of the $\mathsf{H}_{\upalpha}$ line of PDS 70 b over days
to month timescales and find an upper limit of 30% on the variability.

We focus in this paper on the spectroscopic monitoring of the
$\mathsf{Pa}_{\upbeta}$ emission line of two low-mass companions – GQ Lup b
and GSC 06214-00210 b – with previously reported line variability over the
course of a semester (Lavigne et al., 2009). We combine these observations
with archival data to probe line variability from minutes to decade
timescales. The $\mathsf{Pa}_{\upbeta}$ line is well adapted to a monitoring
with standard adaptive-optics fed instruments from the ground and is more
intense above the continuum than the $\mathsf{Br}_{\upgamma}$ line in young
accreting late-type objects (Natta et al., 2004; Betti et al., 2022b).
Moreover, it is known to be less affected by chromospheric activity than
$\mathsf{H}_{\upalpha}$ in the stellar and sub-stellar mass regimes (Natta et
al., 2004; Antoniucci et al., 2011; Rigliaco et al., 2012; Manara et al.,
2013).

We present a compendium of the presently known properties of the two objects
in Section 2. The observations of GQ Lup b and GSC 06214-00210 b are detailed
in Section 3 followed by the data extraction strategy in Section 4. The
photospheric parameters, and emission line variability (timescale, intensity,
profiles) are presented in Section 5 and compared to models. We analyze how
the $\mathsf{Pa}_{\upbeta}$ line behavior relates to other accretors,
contextualize possible line profile variations, and discuss the implication of
line variability for the detectability of protoplanets in Section 6. Results
are summarized in Section 7.

##  2 Information on targets

###  2.1 GQ Lup system

GQ Lup b is a 10 to 40 $\mathrm{M_{Jup}}$ companion resolved by Neuhäuser et
al. (2005) at $\sim$0.7” ($\sim$108 au at $154.10\pm 0.69$ pc; Gaia
Collaboration, 2022) from a 2-5 Myr old K7Ve-type T Tauri star
($\textrm{M}_{\star}=1.02\pm 0.05$ $\mathrm{M_{\odot}}$; MacGregor et al.,
2017; Gaia Collaboration, 2022) member of the Lupus I cloud (Herbig, 1977).
The star is surrounded by a circumstellar disk inclined at 60∘ showing spiral
structures (e.g., van Holstein et al., 2021) in scattered light up to 0.84”
(129 au), while the disk shows smaller extension at millimeter and sub-
millimeter wavelengths both in continuum and 12CO (20–46 au; MacGregor et al.,
2017; Wu et al., 2017b; Long et al., 2020). Long et al. (2020) report a gap at
10 au in the dust component of the disk at millimeter wavelengths which might
be carved by a planet at this separation. It is unclear whether this possible
planet could explain the significant radial-velocity variation of GQ Lup A
noticed by Donati et al. (2012) on a 2 years time span. The star has a
reported spin axis inclination of $\sim$30∘ misaligned with the disk and is
experiencing dominantly pole-ward magnetospheric accretion that produces a set
of emission lines, including $\mathsf{Pa}_{\upbeta}$ (Broeg et al., 2007;
Donati et al., 2012; Frasca et al., 2017; Seifahrt et al., 2007).

GQ Lup b is also known to display a bright $\mathsf{H}_{\upalpha}$ emission
line found at four epochs spanning two decades: on April 10, 1999 HST/WFPC2
wide-field imaging data (F606W filter; $\uplambda_{cen}=605.4$,
$\textrm{FWHM}=220.5$ nm; Marois et al., 2007), in HST/WFC3-UVIS2 data (Zhou
et al., 2014) acquired on February 25, 2012 with a narrower filter (F656N,
$\uplambda_{\textrm{cen}}=656.2\leavevmode\nobreak\ nm,\textrm{FWHM}=187.1$
nm); on ground-based Magelan/MagAO imaging (Wu et al., 2017b) with a custom
$\mathsf{H}_{\upalpha}$ filter on April 16, 2015
($\uplambda_{\textrm{cen}}=656.0\leavevmode\nobreak\ nm,\textrm{FWHM}=6.3$
nm), and with VLT/MUSE medium-resolution (R=2516) optical spectroscopy on
April 19, 2019 (Stolker et al., 2021). Wu et al. (2017b) argues that the line
could be showing variability, but the different instrumental setup prevents
from conducting a robust homogeneous comparison of the reported
$\mathsf{H}_{\upalpha}$ fluxes. The MUSE spectrum shows the Ca II triplet in
emission ($\uplambda$8498, 8542, 8662 Å) but no $\mathsf{H}_{\upbeta}$ line.

The companion displays a $\mathsf{Pa}_{\upbeta}$ emission line first seen on
September 18, 2006 (EW=$-3.83\pm 0.12$ Å) from data acquired with VLT/SINFONI
(Seifahrt et al., 2007) at R=2360. The line is also detected in Gemini/NIFS
data (R$\sim$5000) collected on May 30, 2017 but at a significantly weaker
flux level (EW=$-0.46\pm 0.08$Å), supporting the observed variability at
$\mathsf{H}_{\upalpha}$. No $\mathsf{Br}_{\upgamma}$ emission line is reported
in K-band spectra obtained at four epochs (Neuhäuser et al., 2005; McElwain et
al., 2007; Seifahrt et al., 2007; Lavigne et al., 2009).

The set of available near-infrared spectra revealed an inconsistency in the
spectral slope of the SINFONI spectrum of Seifahrt et al. (2007) that we re-
investigate in Section 4. The compendium of 0.6–5 $\upmu$m spectro-photometry
on the companion allowed Stolker et al. (2021) to derive
$\textrm{T}_{\textrm{eff}}$ of 2700 K consistent with the optical spectral
type of M9 (e.g, Luhman et al., 2003) and the luminosity-based “hot-start”
mass of 30 $\mathrm{M_{Jup}}$. The analysis also confirms a significant in-
line extinction ($\textrm{Av}=2.3$ mag) and an excess emission long-wards of
2.5 $\upmu$m that could be caused by infalling material from the
circumplanetary disk or a CPD. ALMA observations have failed to detect this
CPD thus far (MacGregor et al., 2017; Wu et al., 2017b) and no intrinsic
linear polarization is detected at the location of the companion (van Holstein
et al., 2021). High-resolution spectroscopy (VLT/CRIRES, R$\sim$100 000) of GQ
Lup b centered on the 12CO overtone at K-band allowed to measure a projected
rotational velocity ($5.3^{+0.9}_{-1.0}$ km/s; Schwarz et al., 2016), in the
lower range of values measured on young 10-50 $\mathrm{M_{Jup}}$ objects
(Kurosawa et al., 2006; Bryan et al., 2018, 2020) which sets an upper limit on
the rotation period of 82h.

GQ Lup b’s orbital motion has been monitored since its original discovery
(Ginski et al., 2014; Schwarz et al., 2016) and now covers $\sim 1$ % of its
estimated orbit. Tentative constraints of the mutual inclination with the
circumstellar disk indicate the companion might be on a polar orbit (Stolker
et al., 2021) with a semi-major axis of $\textrm{a}=117^{+24}_{-23}$ au and
low eccentricity possibly crossing the circumstellar disk.

A second accreting low-mass stellar companion (2MASS 15491331-3539118;
projected separation $\sim 2400$ au; $M\sim 0.15$ $\mathrm{M_{\odot}}$, GQ Lup
C) was recently found to be bound to the star. GQ Lup c bears its own detached
disk roughly aligned with the one of the primary star and with the large-scale
dust filament of Lupus I containing GQ Lup (Alcalá et al., 2020; Lazzoni et
al., 2020).

In summary, this system shows companions of different masses and provide a
unique mean to better understand the early formation and dynamical evolution
of planetary systems.

###  2.2 GSC 06214-00210 system

GSC 06214-00210 b (hereafter GSC06214b) is a $14.5\pm 2$ $\mathrm{M_{Jup}}$
co-moving companion (Kraus et al., 2008; Ireland et al., 2011; Pearce et al.,
2019) resolved on a wide orbit (2.2” separation, corresponding to a projected
physical separation of 240 au at $108.8\pm 0.3$ pc; Gaia Collaboration, 2022)
from the weak-line T Tauri star GSC 06214-00210 (K$5\pm 1$;
$\textrm{M}_{\star}=0.80\pm 0.11$ $\mathrm{M_{\odot}}$; Bowler et al., 2014;
Pearce et al., 2019).

The star is reported to be a member of the Upper Scorpius sub-group of the
Scorpius-Centaurus (Sco-Cen) OB association (Preibisch et al., 1998), with a
median age of $\sim 10$ Myr (Pecaut et al., 2012; Feiden, 2016). We notice
that the kinematic analysis tool BANYAN $\Sigma$ (Gagné et al., 2018b) now
suggests a possible alternative membership (58.8% probability) to the older
Upper-Scorpius Centaurus (UCL) sub-group of Sco-Cen (16 Myr; Pecaut et al.,
2012) based on the GAIA-DR3 astrometry. This possible membership to UCL is
better consistent with the revised isochronal age of $24^{+7}_{-5}$ Myr
reported in Pearce et al. (2019) and based on the Baraffe et al. (2015)
evolutionary tracks and their $16.9^{+2.9}_{-1.9}$ Myr age inferred using the
magnetized evolutionary models of Feiden (2016) and a prior on the age of Sco-
Cen. We adopt this later value in our study.

The star is known to have no reported excess emission, a weak
$\mathsf{H}_{\upalpha}$ emission line, and no $\mathsf{Pa}_{\upbeta}$ emission
(Preibisch et al., 1998; Bowler et al., 2011). Conversely, medium-resolution
1.17–1.83 $\upmu$m spectroscopy (R=3800) of the system with OSIRIS at Keck in
2010 (Bowler et al., 2011) revealed a bright $\mathsf{Pa}_{\upbeta}$ emission
line (EW$=-11.4\pm 0.3$ Å) on GSC06214b. HST/WFC3 photometry also reveals the
companion to emit at $\mathsf{H}_{\upalpha}$ (Zhou et al., 2014). The
companion has a continuum excess emission noticeable at K- and L’- bands
(1.8-4.2 $\upmu$m; Bowler et al., 2011; van Holstein et al., 2021) which
indicates it likely bears a CPD. The CPD has remained elusive at sub-
millimeter wavelengths (Bowler et al., 2015; Wu et al., 2017a) indicative of a
compact detached disk. A likely intrinsic polarization signal (van Holstein et
al., 2021) was however detected with VLT/SPHERE. van Holstein et al. (2021)
argue that the CPD must have a inclination of $\sim 48^{\circ}$ to explain the
measured degree of polarization. They used the possible constraints on the CPD
inclination and the measured rotational velocity of $6.1^{+4.9}_{-3.8}$ km/s
of GSC06214b (Bryan et al., 2018) to narrow down the companion rotation period
to 22–77 h.

An additional 0.9–2.5 $\upmu$m spectrum of GSC06214b was obtained with GNIRS
on June 23, 2011 by Lachapelle et al. (2015) with a resolution lower than
OSIRIS (R$\sim$1000). It shows a weaker $\mathsf{Pa}_{\upbeta}$ emission
(EW$=-4.4\pm 0.5$ Å) than in Bowler et al. (2011) and a
$\mathsf{Br}_{\upgamma}$ emission line (EW$=-0.24\pm 0.05$ Å), making
GSC06214b the third companion with mass below 20 $\mathrm{M_{Jup}}$ to display
this line (Zhang et al., 2021; Betti et al., 2022b). K-band (1.965–2.381
$\upmu$m) OSIRIS spectra obtained on June 26, 2012 at R=3800 do not show the
reported $\mathsf{Br}_{\upgamma}$ line. Bryan et al. (2018) do not mention the
detection of this line in their R=25000 Keck/NIRSPEC spectrum (2.03-2.38
$\upmu$m) of the target obtained on June 3, 2015 either. But this line is re-
detected in archival VLT/SINFONI spectra (R$\sim$5000) obtained on May 28,
2014 and July 07, 2014 (Palma-Bifani et al., in preparation) further hinting
for significant year-to-year variability of accretion lines on that companion.
All sets of near-infrared spectra allow to classify this companion as a
M9-M9.5 low-gravity object (Lavigne et al., 2009; Bowler et al., 2014).

GSC06214b offers to probe the accretion processes at the Deuterium-burning
boundary in a more evolved system, thus complementing GQ Lup b’s analysis.

##  3 Observations and data reduction

###  3.1 Monitoring sequences with VLT/SINFONI

We monitored GQ Lup b and GSC06214b with the Integral Field Spectrograph (IFS)
SINFONI. SINFONI is a now decommissioned instrument previously installed at
VLT/UT4 Cassegrain focus and composed of a custom adaptive optics module
(MACAO) feeding an IFS (SPIFFI). SPIFFI splits the field-of-view into 32
horizontal slices (slitlets) and re-arranges them to form a pseudo-slit which
is dispersed by a grating on a 2k$\times$2k HAWAII-2RG detector (Eisenhauer et
al., 2003; Bonnet et al., 2004). The instrument was operated with pre-optics
and gratings in the J-band (1.1 – 1.35 $\upmu$m). For GQ Lup b, the setup
provided a 0.8” square field of view (FoV) composed of $12.5\times 25$ mas
spaxels, with a resolving power of
$R_{\uplambda}=\frac{\uplambda}{\Delta\uplambda}=2360$. For GSC06214b, the
setup provided a 3.2” square field of view (FoV) composed of $50\times 100$
mas spaxels, with a resolving power of
$R_{\uplambda}=\frac{\uplambda}{\Delta\uplambda}=1800$. They translate into
Line Spread Function (LSF) widths of 127 and 167 km/s for GQ Lup b and
GSC06214b respectively.

MACAO used the host stars for wavefront sensing. The de-rotator at the
Cassegrain focus was used to align the star and companions horizontally in the
FoV at acquisition. A telescope offset was then applied to place the host
stars PSF cores outside of the FoV and allow for deeper exposures. Small
additional dithering were applied along the vertical direction of the FoV in-
between each exposure to ensure that the companion signal did not fall all the
time on bad pixels. A blank-sky exposure (with identical DIT) was performed at
the end of the dithering sequence to ensure a proper removal of the OH-
emission lines and thermal background. We then moved the star back in the FoV
and recorded an additional short-exposure to record a PSF, to be used for
calibration purposes.

GSC06214b was observed on April 25, May 4, June 1, and August 18, 2017 as part
of our program (099.C-0760). GQ Lup b was observed on April 25, 30, 2017 and
June 1, 2017 as part of the same program (UT dates). We also decided to re-
process SINFONI data of GQ Lup b obtained on 2006-09-18 (program 077.C-0264)
and presented in Seifahrt et al. (2007) for consistency and given that the
published SINFONI spectrum has been shown to be discrepant with respect to
other studies (McElwain et al., 2007; Lavigne et al., 2009). These
observations used the same setup and followed an identical data-extraction
strategy as for the 2017 observations. The log of observations is detailed in
Table LABEL:tab:obs.

###  3.2 SINFONI cube reconstruction and calibration

SINFONI records bi-dimensional raw frames of the 32 dispersed slitlets. We
used the Toolkit for Exoplanet deTection and chaRacterization with IfS
(hereafter TExTRIS) (Petrus et al., 2021, Bonnefoy et al., in prep) to correct
them from static noise caused by the detector electronics (Bonnefoy et al.,
2014).

The SINFONI data handling pipeline v3.0.0 (Abuter et al., 2006) was then used
to reconstruct calibrated data cubes encoding the spatial and spectral
dimensions (X, Y, $\uplambda$) from these raw data. The pipeline uses
calibration frames obtained at day time to apply basic cosmetic steps on the
raw frames and correct them from the distortion introduced by optical elements
and the detector. It identifies the position of slitlets on the raw frames and
builds a reference 2D map associating the pixels to given wavelengths. The sky
emission was evaluated and removed through the field subtracting cubes of
reconstructed exposures on sky. Sky OH- emission lines residuals were then
evaluated and subtracted using the method of Davies (2007). Because residuals
hot/cosmics pixels were found in the cubes, they were corrected for using the
lacosmic Python package on each individual wavelength frame (van Dokkum,
2001).

The wavelength solution relies on day-time calibration frames obtained with an
Argon lamp and has proven to be inaccurate on other datasets (Petrus et al.,
2021; Zhang et al., 2021). We therefore improved the pipeline wavelength
calibration with routines from TExTRIS, comparing the many telluric absorption
lines contained in each spaxel to a model generated for the observing
conditions using the ESO Skycalc tool (Noll, S. et al., 2012; Jones, A. et
al., 2013). The method is applicable whenever the spaxels contain enough
stellar light to infer the telluric absorptions. We evaluated the shifts in
the wavelength interval 1.1–1.3 $\upmu$m using cross-correlations of spaxels
with a signal-to-noise ratio (S/N) above 10 at $\mathsf{Pa}_{\upbeta}$. A
master median-combined map of wavelength shift was then created and each
individual cube was corrected from the median of the shifts of individual
spaxels meeting the S/N criterion. The corrected wavelength shifts at each
epochs are reported in Table LABEL:tab:obs, for an average wavelength
calibration uncertainty of $\sim 0.2$ Å and $\sim 0.3$ Å for GQ Lup b and
GSC06214b respectively. Because of the relatively low S/N in the science
target cubes, we used the wavelength solution found with the telluric standard
stars (relying on the same original calibration as the science target) as the
solution for the whole sequence.

###  3.3 New investigation of Keck/OSIRIS data

We also include in the analysis the individual spectra of GSC06214b obtained
on July 9, 2010 with the OH-Suppressing Infrared Imaging Spectrograph (OSIRIS;
Larkin et al., 2006) mounted on Keck II. OSIRIS is a lenslet-based IFS
operated with the Keck adaptive optics system. The instrument produced data
cubes sampling a rectangular FoV ($0.8\times 3.2$”) with 0.05x0.05′′ spaxels.
The J${}_{\textrm{bb}}$ band grating provided 1.18-1.416 $\upmu$m spectra at
R$\sim$3800\. Eight 5-minutes long exposures were recorded. The stacked data
cubes are presented in Bowler et al. (2011) and the companion is sufficiently
separated from the star to directly extract its spectrum within a circular
aperture. The data analysis is described in more details in Bowler et al.
(2011). We re-analyzed here the individual data cubes in an attempt to detect
$\mathsf{Pa}_{\upbeta}$ variability from one exposure to the next.

##  4 Spectral extraction

###  4.1 Halo removal

Both GQ Lup b and GSC06214b spectra are contaminated by the PSF wings of the
primary stars located outside of the FoV. The data cubes obtained at the end
of each sequence of observations with the host stars in the FoV do not probe
the extended wings of the PSF and therefore do not allow for direct
subtraction of the contamination. We therefore had to develop a specific
procedure described below and applied to each individual SINFONI data cube on
both targets.

Before halo removal, the cube rectangular spaxels were linearly interpolated
to produce a regular spatial sampling of 12.5 mas in the horizontal and
vertical directions. The contamination was then evaluated and corrected in
each individual frame of the data cubes, with a CLEAN-like algorithm (Högbom,
1974). The algorithm fits a two component flux model made of: (i) a third-
order bi-dimensional gradient, approximating the primary star PSF wings, (ii)
a PSF model placed at the companion location and scaled to the companion’s
flux. The latter step requires the definition of a PSF model. We used the
primary star data cubes (normalized to 1) for that purpose in the 2017 data.
Data cubes of the primary star were not available in the 2006 data and we
chose to use those corresponding to the telluric standard star instead. The
two flux components were successively fitted and subtracted from the original
frames, for a total of 4 iterations, which is a trade-off between
$\chi^{2}_{r}$ and S/N optimization. See Fig. 1 for a visual representation,
and Fig. 11 for residuals.

The host star gradient subtraction allows us to subtract the large-scale
structures but can not remove more complex features such as the diffraction
pattern caused by the telescope spiders. This component is most important on
GQ Lup b and was then approximated as a 1-d Gaussian profile in the direction
perpendicular to the spider axis and subtracted in the gradient-subtracted
cubes.

SINFONI shows a static dark horizontal area in the lower half part of all data
cubes (see Fig 1) likely caused by persistence appearing in the lamp flat-
fields taken at day time. This was discovered during data reduction and
affected the spectra profiles when the companion was located in this region.
We decided to exclude the corresponding cubes and spectra from the analysis.

Figure 1: Illustration of the halo subtraction, for the first data cube of the
GQ Lup b 2017-04-30 epoch (median along the wavelength axis). We also report
the location of the diffraction pattern of the telescope spiders, along with
the identified area with improperly calibrated instrument transmission that
affects the lower-half part of the images. Figure 2: Companion spectra. Top
row: GQ Lup b spectra. Bottom row: GSC 06214-00210 b spectra. Left column:
full spectra. Right column: zoom on the $\mathsf{Pa}_{\upbeta}$ line. In both
cases, the spectra continuum are normalized to 1 as the mean value on the
1.27–1.28 $\upmu$m and 1.29–1.30 $\upmu$m ranges. The 2010-07-09 (OSIRIS)
epoch of GSC06214b was degraded to SINFONI resolution (R=1800) in the left
column, but is shown at full resolution (R=3800) in the right column. The
spectral slopes differences are the results of different Strehl Ratios between
epochs, despite the AO correction. Our re-extraction of GQ Lup b spectrum
initially presented in Seifahrt et al. (2007) produces a pseudo-continuum
short-wards of 1.24 $\upmu$m in better agreement with the more recent SINFONI
epochs and observations produced by NIFS.

Figure 3: $\mathsf{Pa}_{\upbeta}$ profiles, corrected from the atmospheric
contribution and baseline subtracted (see Sects. 5.2 & 5.3). Each spectrum is
labeled with the corresponding UT time at the start of the integration. First
row: GQ Lup b lines for each epoch. Second and third rows: same for GSC06214b,
with the addition of the individual OSIRIS lines from Bowler et al. (2011)
(2010-07-09 epoch).

###  4.2 Extraction and flux calibration

The PSF model estimated from the primary (or telluric standard) star was not a
good enough match to the companion PSF to extract spectra from profile
fitting. Therefore, the companions spectra were extracted within a circular
aperture on the contamination-corrected data cubes of constant radius with
wavelengths. The size of the aperture was chosen to conserve the spectral
slope and to avoid integrating noise due to the imperfect removal of the
stellar halo in the background-limited area. This procedure is similar as the
one adopted in Bonnefoy et al. (2014) and proven to provide a robust
extraction of the spectral features and slope on J-band SINFONI data. The
adopted radii were ranging from 50 to 75mas. We retrieve spectra of various
S/N depending on the observing conditions. For GQ Lup b, the S/N ranges from
20 to 30 in the 2006-09-18 and 2017-06-01 epochs, while it ranges from 35 to
50 in the 2017-04-25 and 2017-04-30 epochs. For GSC06214b, the S/N is
relatively stable, within the 15-30 range for all epochs.

The atmospheric transmission was evaluated computing the ratio between the
extracted spectra of telluric standard stars and their synthetic spectra:
black-bodies scaled to the 2MASS JHK, and WISE W1 and W2 photometry of the
objects (Cutri et al., 2003; Cutri & et al., 2012) collected with the
VOSA111http://svo2.cab.inta-csic.es/theory/vosa/ tool (Bayo et al., 2008). The
$\mathsf{Pa}_{\upbeta}$ absorption line of the standards were evaluated and
removed beforehand as follows. The telluric standard star spectrum was
corrected from telluric absorptions by a model created with the SkyCalc tool
(Noll, S. et al., 2012; Jones, A. et al., 2013). The line profile was then
fitted in the by a Voigt profile function, then subtracted from the spectrum.
The spectrum was then applied the same telluric absorption model. This ensures
that the line shape evaluation is not impacted by weak telluric features in
the line wings and that it does not bias the companion emissions at these
wavelengths.

We performed an absolute flux calibration of the individual spectra using the
telluric standard star as a spectro-photometric reference (2MASS J-band
photometry). The use of the telluric standard stars was favored over the use
of the host stars because of their potential photometric variability. This is
the case for GQ Lup A which is a known high-amplitude variable star ($\Delta
J=0.87$, Batalha et al., 2001; Broeg et al., 2007). The remaining J-band
magnitudes inferred from the individual epochs were averaged to remove part of
the uncertainty on the photometric calibration due to variation in observing
conditions and AO performance between the telluric standard stars and the
companion observations. We decided to exclude the 2006 epoch on GQ Lup b from
this process because the companion and the telluric standard were observed
under poor and variable adaptive-optics conditions.

We find a magnitude J = $16.21^{+0.26}_{-0.21}$ on GSC06214b. The value is
well consistent with the one measured by Lachapelle et al. (2015) and Bowler
et al. (2011) (J=$16.19\pm 0.04$ and J$=16.25\pm 0.04$, respectively) derived
from the contrast between the primary and the companion, and the tabulated
2MASS J-band magnitude of the system. We decided to use the latter value.

We find J = $14.38^{+0.13}_{-0.12}$ for GQ Lup b. Lavigne et al. (2009) found
J = 15.13 $\pm$ 0.44 but the measurement accuracy was limited by the partial
saturation of the host star in their data and the photometry relies on the
2MASS magnitude of the host star which is known to be variable, as previously
mentioned. McElwain et al. (2007) found J = 14.90 $\pm$ 0.11 using a similar
procedure as ours and this photometry was used in all contemporaneous studies
on the system (e.g., Stolker et al., 2021). But this measurement was performed
on OSIRIS data at high airmass with a narrower wavelength range and at a lower
signal-to-noise ratio. Given the amplitude of variability of the primary star
($\Delta$J = 0.87, Broeg et al., 2007), the discrepancy between published
values and ours is not surprising. The SINFONI photometry makes the J-K color
of GQ Lup b ($1.0\pm 0.2$ mag, assuming the K-band magnitude of GQ Lup bre-
evaluated by Neuhäuser et al., 2008) in better agreement with that of other
companions of the same luminosity, spectral type, and age range (e.g., Fig. 4
of Gagné et al., 2018a). Furthermore, our value is also consistent with the
magnitude correction (0.30 to 0.66 mag) found by Stolker et al. (2021) to
adjust with a single atmospheric model the J-band scaled to the McElwain et
al. (2007) photometry with the rest of the companion spectral energy
distribution. We therefore choose to adopt the J-band magnitude inferred from
the SINFONI data later on (see Sect. 6.3 for a discussion on possible
photometric variations).

Spectra are also corrected for the Barycentric Earth Radial Velocity (BERV),
the primary radial velocities, and the companion keplerian motions. The BERV
is computed following Piskunov & Valenti (2002) for each individual data cube.
GQ Lup b radial velocity is $-0.8\pm 0.4$ km/s (Schwarz et al., 2016).
GSC06214b radial velocity is $-3.98\pm 0.88$ km/s, it was computed from a
combination of its orbital velocity ($0.61\pm 0.09$ km/s, Pearce et al.
(2019)) and the radial velocity of the primary star ($-4.59\pm 0.88$ km/s,
GAIA DR2). Radial velocities are given in the heliocentric system.

For all the following analysis, error bars were computed from the dispersion
around $\mathsf{Pa}_{\upbeta}$, on the residual photosphere-subtracted and
baseline-corrected (polynomial approximation of residual baselines, see Sect.
5.3) spectra. For velocity conversion (e.g., right panels of Fig. 2), we used
the $\mathsf{Pa}_{\upbeta}$ rest-frame wavelength gathered on NIST (Kramida et
al., 2022), $\uplambda_{0}=1.2821578$ $\upmu$m (vacuum, Chang & Deming, 1996).

###  4.3 Photospheric contribution

Both GQ Lup b and GSC06214b exhibit a rich spectrum of molecular absorption
lines and shallow KI doublets at 1.169/1.177 $\upmu$m and 1.243/1.253 $\upmu$m
characteristic of young M-L transition objects (e.g., Bonnefoy et al., 2014).
Molecular absorption lines overlap with the $\mathsf{Pa}_{\upbeta}$ line,
affecting the emission line profiles. The spectra were therefore corrected for
the photospheric features using our forward modeling tool ForMoSA (Petrus et
al., 2020, 2023).

ForMoSA compares pre-computed grids of atmospheric models with observations
using the Nested Sampling Bayesian inference algorithm (Skilling, 2004). We
used as input the ATMO models (Tremblin et al., 2015, 2016) which explore
effective temperature and surface gravity in the expected range of values for
the late-M companions. The grids are described in more details in Table C.1 of
Petrus et al. (2023) and allow to explore non-solar C/O and [M/H] values, as
well as the adiabatic index $\upgamma$. ForMoSA also allows for fitting of the
radial velocity (RV), rotational broadening (vsin(i)), limb-darkening
coefficient ($\upepsilon$), and interstellar extinction (Av). We did not fit
the photospheric contribution on each individual spectrum, but rather on the
averaged spectrum at each epoch and used the RV shift to compensate for
residual wavelength shifts from one epoch to another that may arise from
imperfect wavelengths re-calibration.

The resulting parameter corner plots are presented in Fig. 4. The
1D-histograms were used to find the resulting best values for each parameter,
bi-modal Gaussian fitting (different $\upsigma$ for the upper and lower
halves). The free parameters are the peak location $\upmu_{0}$ and the two
sigmas $\upsigma_{1}$ and $\upsigma_{2}$. The retained value is $\upmu_{0}$
with the corresponding error bars. The dispersion of values between epochs
makes it clear that the Bayesian process on individual epochs provides
posterior distributions whose extents are not representative of the true error
bars on the values. We therefore estimate the adopted value as the mean of all
epochs, and error bars from the dispersion of the ensemble of posteriors of
the different epochs. Values at the edge of the grid (upper/lower limits) were
ignored when computing the adopted value. The results are presented in Table 1
and discussed in the following section.

Table 1: Atmospheric parameters retrieved with ForMoSA

| Date | Teff | log(g) | [M/H] | $\upgamma$ | C/O | Radius | $A_{v}$ | RV
---|---|---|---|---|---|---|---|---|---
| YYYY-MM-DD | (K) |  |  |  |  |  ($\mathrm{R_{Jup}}$) |  | (km/s)
GQ Lup b | 2006-09-18 | $2587$ | $4.16$ |  $>0.59$111Unconstrained parameter, fit converged at the edge of the parameter range. |  $>1.05$111Unconstrained parameter, fit converged at the edge of the parameter range. | $0.57$ | $4.45$ | $2.38$ | $3.03$
2017-04-25 | $2671$ | $4.16$ | $0.30$ |  $>1.05$111Unconstrained parameter, fit converged at the edge of the parameter range. | $0.42$ | $4.08$ | $2.13$ | $5.59$
2017-04-30 | $2655$ | $4.22$ | $0.17$ |  $>1.05$111Unconstrained parameter, fit converged at the edge of the parameter range. | $0.33$ | $4.07$ | $1.91$ | $13.59$
222Epoch ignored when computing the adopted value2017-06-01 | $2711$ | $3.92$ | $0.30$ |  $>1.05$111Unconstrained parameter, fit converged at the edge of the parameter range. | $0.55$ |  $>5.00$111Unconstrained parameter, fit converged at the edge of the parameter range. | $4.32$ | $13.11$
adopted | $2638^{+33}_{-51}$ | $4.18^{+0.04}_{-0.02}$ | $0.23^{+0.06}_{-0.06}$ | ¿1.05 | $0.44^{+0.13}_{-0.11}$ | $4.20^{+0.25}_{-0.13}$ | $2.14^{+0.24}_{-0.23}$ | -
GSC06214b | 2017-04-25_set1 | $2274$ | $4.30$ | $-0.13$ | $1.03$ | $0.50$ | $1.45$ | $0.45$ | $-18.06$
2017-04-25_set2 | $2230$ | $4.10$ | $-0.33$ | $1.03$ |  $<0.30$111Unconstrained parameter, fit converged at the edge of the parameter range. | $1.47$ |  $<0.00$111Unconstrained parameter, fit converged at the edge of the parameter range. | $-8.78$
2017-05-04 | $2317$ | $4.30$ | $-0.15$ | $1.03$ | $0.42$ | $1.56$ | $1.50$ | $-12.73$
222Epoch ignored when computing the adopted value2017-06-01 | $2269$ | $4.33$ | $-0.18$ | $1.02$ | $0.39$ | $2.45$ |  $>4.99$111Unconstrained parameter, fit converged at the edge of the parameter range. | $-19.58$
2017-08-19 | $2189$ | $4.36$ | $0.00$ | $1.01$ | $0.64$ | $1.58$ |  $<0.00$111Unconstrained parameter, fit converged at the edge of the parameter range. | $-6.51$
2010-07-09 | $2232$ | $4.17$ | $-0.21$ | $1.04$ | $0.36$ | $1.37$ |  $<0.00$111Unconstrained parameter, fit converged at the edge of the parameter range. | $-11.90$
adopted | $2248^{+69}_{-59}$ | $4.24^{+0.11}_{-0.15}$ | $-0.16^{+0.17}_{-0.17}$ | $1.03^{+0.01}_{-0.02}$ | $0.48^{+0.16}_{-0.12}$ | $1.49^{+0.10}_{-0.12}$ | $0.97^{+0.53}_{-0.53}$ | -

222Adopted values are obtained as described in Sect. 4.3.

Figure 4: ForMoSA resulting photosphere parameters. Left: GQ Lup b results.
Right: GSC06214b results. The Radius, Av, radial velocity (RV) and limb-
darkening coefficients are omitted from these plots for the sake of clarity.
Figure 5: Visual comparison of the photosphere fitting with ForMoSA, for each
epoch of both objects, with a zoom on $\mathsf{Pa}_{\upbeta}$ in the right
column. Grayed out wavelength ranges were excluded from fitting as the model
either did not include corresponding species (FeH absorption band,
$\mathsf{Pa}_{\upbeta}$) or data were too noisy ($>$1.34 $\upmu$m).

##  5 Results

In this section, we present the resulting spectra and atmospheric parameters,
along with their lines properties (profiles and variability) and modeling.

###  5.1 Companion J-band spectra

The mean spectra for each epoch are presented in Fig. 2, along with a zoom on
the $\mathsf{Pa}_{\upbeta}$ lines. $\mathsf{Pa}_{\upbeta}$ lines of individual
data cubes are presented in Fig. 3. The spectra all have a relatively similar
spectral slope, with a discrepancy for the 2017-06-01 epoch for both objects.
We argue that this is due to the lower Strehl ratios during that night – the
AO performance was significantly worse than for other epochs – which is
affecting mostly the blue part of the spectrum, leading to larger flux loss at
smaller wavelengths. Overall, our spectral slope is in agreement with the
already published spectra (Lavigne et al., 2009; Bowler et al., 2011), but the
GQ Lup b spectra differ from the Seifahrt et al. (2007), which was already
pointed out by Lavigne et al. (2009). This may be explained by residual
primary flux in their spectrum, which is in line with the wider
$\mathsf{Pa}_{\upbeta}$ line profile (See Fig. 2, top-right panel). We observe
spectroscopic features typical of late-M/early-L dwarfs for both objects at
all epochs: the two KI doublets (1.17 and 1.25 $\upmu$m), the FeH absorption
band (1.17–1.22 $\upmu$m) and the numerous H2O absorption lines (1.13–1.16
$\upmu$m and $>$1.32 $\upmu$m, see Fig. 2). The spectra are also composed of
numerous less prominent molecular features all over the band, which are best
described by photosphere models (see Section 5.2).

###  5.2 Atmospheric and physical parameters

We present here the results from atmospheric modeling of the J-Band spectra
with ForMoSA (Sect. 4.3). Best fits for both sources can be found in Table 1.
A visual comparison of data and fitted spectra can be found in Fig. 5. A
comparison of the various atmospheric parameters can be found in Fig. 4, in
the form of corner plots. Some parameters converged at the edge of the
parameter range, such as the adiabatic index $\upgamma$, either consistently
or occasionally. The spectral slopes of both companions on June 1, 2017 are
discrepant with the other epochs and the retrieved atmospheric parameters (at
the very least $\textrm{T}_{\textrm{eff}}$ and Av) are likely non
representative.

####  5.2.1 Inferred values and correlations

Correlations are observed between $\textrm{T}_{\textrm{eff}}$ and log g. The
behavior is similar to that observed on the spectral analysis of the young
late-M companion HIP 78530B with ForMoSA in Petrus et al. (2020) (albeit they
used BT-SETTL atmospheric models). The correlation might be caused by the
increase in strength of the K I and FeH absorptions with both decreasing
$\textrm{T}_{\textrm{eff}}$ and increasing surface gravities at the M-L
transition (Rice et al., 2010; Bonnefoy et al., 2014).

On the other hand, the log g and C/O are anti-correlated. This is due to the
sensitivity to these two parameters of the numerous FeH and H2O absorptions at
J-band. Finally, both sources show either no correlations between [M/H] and
C/O, or a positive one. We do not notice strong degeneracies between
$\upgamma$ and the others parameters, which suggests that the limited
exploration of that parameter does not bias the estimate of the other
parameters.

For GQ Lup b, atmospheric parameters are consistent for all epochs:
$\textrm{log g}\sim 4.2$, $\textrm{T}_{\textrm{eff}}=2638^{+33}_{-51}$ K
(beware that error bars do not include systematics of the model) and
$\upgamma>1.05$. The log g and $\textrm{T}_{\textrm{eff}}$ are both well
consistent with the values found in Stolker et al. (2021) adjusting the
VLT/MUSE optical spectrum of the companion with BT-SETTL models. We however
find with ForMoSA a higher radius than the estimates of Stolker et al. (2021)
(found by adjusting a synthetic spectra on the companion’s spectro-
photometry). This might be due to the increased J-band apparent flux that we
measure in Section 4.2. Part of the difference could also be caused by the use
of different atmospheric models (see below).

For GSC06214b, we find $\textrm{log g}\sim 4.24$,
$\textrm{T}_{\textrm{eff}}=2248^{+69}_{-59}$ K and
$\upgamma=1.03^{+0.01}_{-0.02}$. Lachapelle et al. (2015) found
$\textrm{T}_{\textrm{eff}}=2300\pm 200$ K and $\textrm{log g}=\textrm{3.0--
4.0}$ using different models (BT-SETTL and DRIFT-PHOENIX, Allard et al., 2011;
Helling et al., 2008; Witte et al., 2009, 2011) and spectra, but which remain
broadly consistent with our values.

We present the first estimates of C/O and [M/H] for both companions. Both
objects have C/O values that are compatible with solar value
($\textrm{C/O}_{\odot}=0.55$) within $1\upsigma$. However, they have different
behavior in terms of metallicity: GSC06214b has a solar metallicity (within
1$\upsigma$) while the one of GQ Lup b appears slightly super-solar
($>$3$\upsigma$). We interpret these abundances in section 6.1.

The extinction of GQ Lup b is consistent with the value found by
($\textrm{Av}\sim 2.3$ mag, Stolker et al., 2021). It is above the values
reported on GQ Lup A and C ($\textrm{Av}\sim 1.0$ mag, McElwain et al., 2007;
Alcalá et al., 2020). The extinction of GSC06214b is within the range of
values measured on the primary star ($\textrm{Av}=0.7\pm 0.5$ mag, Barenfeld
et al., 2016; Anders et al., 2022).

####  5.2.2 Use of evolutionary tracks

Using our estimated $\textrm{T}_{\textrm{eff}}$ from ForMoSA, recent ”hot-
start” evolution models (Phillips et al., 2020) predict a mass of $15.5\pm
0.5$ $\mathrm{M_{Jup}}$ at the age of the system (Sect. 2.2) for GSC06214b,
consistent with previous estimates based on the bollometric luminosity. Our
estimated surface gravity is broadly consistent with expectations at the age
of UCL ($\textrm{log g}=4.07\pm 0.01$) and Upper Scorpius ($\textrm{log
g}=4.04\pm 0.02$) and does not allow for discrimination between the two
possible system memberships. Predicted radii range from 1.8 to 2.0
$\mathrm{R_{Jup}}$, e.g. $\sim 20$% larger than the value found with
atmospheric models alone. The mismatch amplitude is consistent with the
predictions of the same atmospheric models for VHS1256b, a cooler
($\textrm{T}_{\textrm{eff}}\sim 1400$ K) companion of the same mass range as
GSC06214b and might point to a systematic overestimation of the absolute near-
infrared fluxes by the ATMO models across a broad $\textrm{T}_{\textrm{eff}}$
range (Petrus et al., 2023). This is, however, not an issue for removing the
photospheric contribution since the radii are only used to scale the synthetic
spectrum to the companion continuum emission level.

The use of GQ Lup b’s $\textrm{T}_{\textrm{eff}}$ from ForMoSA as input for
the ATMO evolutionary tracks predict a mass of $22^{+2}_{-3}$
$\mathrm{M_{Jup}}$, a radius in the range $2.65-3.3$ $\mathrm{R_{Jup}}$, and
$\textrm{log g}\sim 3.8-3.9$ dex at the system age (Section 2.1). Both the
predicted radius and log g are inconsistent with those inferred from ForMoSA
and the deviation is at odd with that observed for GSC06214b. The use of the
older AMES-DUSTY tracks (Chabrier et al., 2000) do not change the conclusions.
This could arise from uncertainties in the evolutionary tracks not capturing
the physics of the interior of this more massive companion (e.g., deuterium
burning). However, the radius is degenerate with the extinction in our fit and
higher extinction than estimated would naturally lead to a better agreement
between ForMoSA and evolutionary model predictions. Stolker et al. (2021) find
radii of $3.55-4.13$ $\mathrm{R_{Jup}}$ and a higher extinction (2.3-2.7 mag),
which remain slightly inconsistent with the Phillips et al. (2020)
evolutionary models. The large radii and inconsistent log g derived from
atmospheric models could also be produced because of additional flux
contribution such as veiling increasing the flux and requiring a larger
dilution factor (Sect. 6.3). A similar mismatch has also been observed on
other companions of the same $\textrm{T}_{\textrm{eff}}$ range in Upper
Scorpius (Petrus et al., 2020) and unresolved binarity has been proposed as an
explanation. Should GQ Lup b be an equal-luminosity equal
$\textrm{T}_{\textrm{eff}}$ binary, the apparent J-band flux would be deviant
overestimated by a factor two with respect to a single object. In such a case,
the corresponding radius of individual components found with ForMoSA
($3^{+0.2}_{-0.1}$ $\mathrm{R_{Jup}}$) would match well the evolutionary model
predictions for the fitted $\textrm{T}_{\textrm{eff}}$.

Figure 6: Gaussian-fitting resulting parameters for the
$\mathsf{Pa}_{\upbeta}$ line of GSC06214b. Left: correlations between
intrinsic FWHM and the line peak location $\upmu_{0}$, arrows represent an
upper limit on the intrinsic FWHM. Right: correlations between the
$\mathsf{Pa}_{\upbeta}$ line flux and the line peak location $\upmu_{0}$. The
intrinsic FWHM is computed as the square root of the quadratic difference
between the measured line FWHM and the LSF FWHM. The $\upmu_{0}$ parameter is
relative to the line expected wavelength in vacuum (1.2821578 $\upmu$m).

###  5.3 Analysis of $\mathsf{Pa}_{\upbeta}$ time series

$\mathsf{Pa}_{\upbeta}$ emission line profiles for individual data cubes are
presented in Fig. 3. We performed Gaussian fitting of the lines to infer their
width, position and flux (see Table 3).

We find the $\mathsf{Pa}_{\upbeta}$ line profiles to be partly resolved in the
2006 epoch of GQ Lup b (1–1.5 times the LSF width, 1–3 $\upsigma$), and mostly
unresolved in the 2017 epochs of GSC06214b (1–1.2 times the LSF width, 1–2
$\upsigma$). Although, we do find the line to be partly resolved in the
2010-07-09 epoch (OSIRIS) (1.2–1.9 the LSF width, 4–8 $\upsigma$). See Table 3
for details. We also find the peak of the $\mathsf{Pa}_{\upbeta}$ emission
lines in GQ Lup b to be slightly blue-shifted. For GSC06214b, they are mostly
centered at zero velocity and show no significant asymmetries (with respect to
the line center), except for the 2017-06-01 epoch of GSC06214b for which the
line has a wider red-wing in comparison to the blue wing. This is reminiscent
of emission lines formed in magnetospheric accretion columns in T Tauri stars,
which show prominent blue/red asymmetry (Edwards et al., 1994). Shock induced
accretion lines also produce lines with a red-shifted peak, hence a wider red
wing. These will be explored in more details in Sect. 6.2.

We studied both the line profiles and fluxes variability. Note that studying
integrated fluxes is equivalent to studying EWs, as they were all computed on
a similar continuum baseline. Integrated fluxes were computed by Gaussian
fitting of the line profiles after local baseline correction (polynomial,
n=2), for which main results are reported in Table 3 and shown in Fig. 6.

As illustrated in Fig. 6, we find, for GSC06214b, possible correlations
between the line intrinsic width and its peak velocity
($\textrm{FWHM}_{\textrm{in}}$–$\upmu_{0}$, left panel), and between the line
flux and its peak velocity (F($\mathsf{Pa}_{\upbeta}$)–$\upmu_{0}$, right
panel). We computed Pearson’s correlations coefficients and p-values using an
MCMC approach almost equal to the pymccorrelation tool (Curran, 2015; Privon
et al., 2020), with the added constraint that the intrinsic FWHM distributions
cannot be negative. For the first correlation, ignoring upper limits we find
$r=0.39^{+0.19}_{-0.22}$ and $\textrm{log p}=-0.88^{+0.59}_{-0.87}$. Note that
while they were excluded, upper limits still appear in the correct region for
the trend. For the second correlation (excluding the 2010-07-09 epoch), we
find $r=0.52^{+0.10}_{-0.12}$ and $\textrm{log p}=-1.43^{+0.52}_{-0.61}$. The
$\textrm{FWHM}_{\textrm{in}}$–$\upmu_{0}$ correlation, if true, could be the
sign of either (i) red-shifted absorption within the line (as absorption
increases, the line centroid will be moved to the blue) or (ii) increasing
red-wing contribution as the line flux increases. The
F($\mathsf{Pa}_{\upbeta}$)–$\upmu_{0}$ correlation, if true could be explained
in the same way as the (ii) explanation for the peak-velocity/FWHM
correlation. Possible implications are discussed in Sect. 6.2.

From Fig. 3, we see unambiguous EW variability from an epoch to another for
both sources. There is no clear sign of line profile variability within the
SINFONI sequences – at the exception of a possible red-wing variability on the
GSC06214b 2017-06-01 epoch. However, we do see line flux variability within
multiple sequences of GSC06214b, in particular the 2010-07-09 one.

Figure 7 shows the relative $\mathsf{Pa}_{\upbeta}$ flux variation, computed
as (Fmax-Fmin)/Fmin for each individual pair of lines, as a function of their
time difference. The marker’s color reports the significance of the
measurement, i.e., the $\upsigma$-distance between the two lines used to
compute this value. While a small and noisy denominator may artificially
increase the relative variation, the measurement will appear gray as the two
values would not be statistically far apart. In other words, gray measurements
are compatible with no relative flux variation (given our uncertainties),
while colored points highlight an actual variation (given the $\upsigma$
level). GQ Lup b shows overall close to no flux variability ($3\upsigma$
confidence level, $<50$%) at the hour timescales. However, we do see
significant variability ($>100$%) at the $\sim$month and longer timescales.
Note that the $\sim$year timescale corresponds to the variation between the
2006-09-18 and the Lavigne et al. (2009) (2007-05-30) epochs. GSC06214b, on
the other hand, shows definitive flux variability around the 30-min and
$\sim$hour timescales ($>5\upsigma$ significance). We also do see clear
($>5\upsigma$) line variability at the weeks, months and years timescales. The
line flux variability at the $\sim$hour timescale appears within multiple
datasets: the OSIRIS epoch (2010-07-09), the 2017-05-04 epoch and between the
two sets of the 2017-04-25 epochs. In both objects, the amplitude of
variability stays below 50–60% for timescales smaller than their estimated
rotation periods, while stronger variability amplitude (up to a factor 10) is
observed at longer timescales. We discuss these findings in the more general
context of variable young stars in Sect. 6.2.

Figure 7: $\mathsf{Pa}_{\upbeta}$ flux variability for both sources, as a
function of the timescale, with respect to the lowest-flux point of the pair.
For GQ Lup b, are included data points from this study as well as the spectrum
from Lavigne et al. (2009). For GSC06214b, are included data points from this
study as well as those from the individual Bowler data cubes (Bowler et al.,
2011) and the mean spectrum from Lachapelle et al. (2015). The color of the
points represents the significance of the measurement, i.e., the
$\upsigma$-distance between the two lines used to compute this value. The
horizontal and vertical bars correspond respectively to the bin size, and the
spread (weighted by significance) within that bin. The red-shaded regions mark
the average observed continuum variability among the SINFONI sequences (see
Table LABEL:tab:obs). The orange shaded regions represent the expected
rotational period of the objects (see Sect. 2).

Figure 8: Comparison of the best fitting models for both GQ Lup b (first row)
and GSC06214b (2nd and 3rd rows). Colorbars represent the reduced $\chi^{2}$
($\upchi^{2}_{r}$) of the Shock (blue) and MagAcc (red) models. For each
epoch, all best models (within the $\upchi^{2}_{r}$ range given in the
colorbar) are over-plotted with a color scaling with the reduced $\chi^{2}$
color. The fitting limits and number of parameters of each model yield a
number of degree of freedom equal to 14 for the Shock model and 9 for the
MagAcc model.

Figure 9: $\upchi^{2}_{r}$ maps illustrating the accretion models degeneracies
when fitting the $\mathsf{Pa}_{\upbeta}$ lines. Left panel:
$\textrm{T}_{\textrm{max}}$ vs $\dot{\textrm{M}}$ degeneracy for the MagAcc
model (GQ Lup b, epoch: 2006-09-18). Right panel: $\textrm{n}_{0}$ vs
$\textrm{v}_{0}$ degeneracy for the Shock model (GSC06214b, epoch: 2010-07-09
– OSIRIS)

###  5.4 Comparison to protoplanet accretion models

$\mathsf{Pa}_{\upbeta}$ emission in the spectra of low-mass free-floating
brown-dwarfs of the same mass, age, and spectral-type range as the two
companions studied here have been investigated by Natta et al. (2004) and
Manara et al. (2013). The $\mathsf{Pa}_{\upbeta}$ line intensity is shown to
be strongly correlated with the accretion luminosity, while the chromospheric
activity does not contribute significantly to the line intensity for accretion
rates above $10^{-11}-10^{-12}$ M⊙ yr-1 at ages 3–10 Myr. Furthermore, none of
the non-accreting sources display $\mathsf{Pa}_{\upbeta}$ emission lines. We
therefore assume here that most, if not all of the $\mathsf{Pa}_{\upbeta}$
emission, is produced by accretion. In this section, we explore two accretion
models in the planetary mass regime predicting $\mathsf{Pa}_{\upbeta}$ line
flux and profiles, and compare their predictions to our time-series
observations, in order to infer the corresponding physical parameters.

####  5.4.1 Description of accretion models

Thanathibodee et al. (2019) proposed a version of the magnetospheric accretion
models developed for Classical T Tauri Stars (CTTS) (see the review of
Hartmann et al., 2016, and ref. therein) extended to the planetary mass range.
It assumes that the central object’s dipolar magnetic field is intense enough
at young ages to truncate the inner CPD, lift material from the CPD and
channel it towards the central object. Material then falls down to the
planetary surface at free fall velocity along accretion funnels. In these
funnels, the material is heated to temperatures high enough ($\sim 10^{4}$ K)
to trigger atomic Hydrogen line emission ($\mathsf{H}_{\upalpha}$,
$\mathsf{Pa}_{\upbeta}$, $\mathsf{Br}_{\upgamma}$, …). The free-falling
material then creates a shock at the base of the funnel flows when it hits the
planetary surface, which triggers hot continuum emission that can be absorbed
by the infalling gas, leading to red-shifted absorption features within the
emission line profiles. Red-shifted absorption also occurs when the accretion
flow absorbs the photosphere emission itself. Blue-shifted absorption may also
be present due to jets or outflows, but these are not accounted for by these
models. The magnetospheric accretion model from Thanathibodee et al. (2019) is
described by: the accretion column launching radius
$\mathrm{R}_{\mathrm{in}}$, its 0pt, its maximum temperature
$\textrm{T}_{\textrm{max}}$, the accretion rate $\dot{\textrm{M}}$, the
inclination, the planet mass M, the planet radius $\textrm{R}_{\textrm{p}}$,
and the planet effective temperature $\textrm{T}_{\textrm{eff}}$. The model is
axis-symmetric and assumes that the large scale dipolar magnetic field is
aligned with the central object rotation axis. The model neglects line
emission contribution from the shock at the base of the accretion funnels.
This model will be referred to as the MagAcc model hereafter.

Conversely, Aoyama et al. (2018) proposed that accretion lines, rather than
being triggered within accretion funnels, would be triggered within an
accretion shock either at the planet or the CPD surface. In such a scenario,
material flows from the circumstellar disk down to the planetary or CPD
surface. This may happen for planet embedded within the circumstellar disk of
their host star, through polar-inflow (Marleau et al., 2022) or, as was
proposed for GQ Lup b, if material from the circumstellar disk channels up to
the planet immediate environment (Stolker et al., 2021). Alternatively, this
model could also describe line emission formed in the shock at the base of
magnetospheric accretion columns. This model considers two parameters: the
infall velocity $\textrm{v}_{0}$ and the pre-shock density $\textrm{n}_{0}$.
We add an additional flux scaling parameter, $\textrm{S}_{\textrm{eff}}$, the
effective shock surface (i.e., $4\pi R_{p}^{2}\textrm{f}_{\textrm{fill}}$ if
it takes place at the planetary surface). See details in next section. This
model will be referred to as the Shock model hereafter.

####  5.4.2 Model fitting

We used both aforementioned models (MagAcc and Shock, Thanathibodee et al.,
2019; Aoyama et al., 2018) to try and retrieve physical parameters of the
accretion mechanisms by fitting the $\mathsf{Pa}_{\upbeta}$ line of both GQ
Lup b and GSC06214b.

Both models were scaled to the distance of the objects (Sect. 2). For the
MagAcc model, this was performed by applying a factor
($\textrm{R}_{\textrm{p}}$/$\textrm{D})^{2}$, with D the distance and
$\textrm{R}_{\textrm{p}}$ the planet radius. For the Shock model, we applied a
factor $\textrm{S}_{\textrm{eff}}$/(4$\uppi\textrm{D}^{2}$). The fits were
performed on photosphere-subtracted and baseline-corrected
$\mathsf{Pa}_{\upbeta}$ emission line profiles, with continuum-subtracted
models. This allows accretion models fitting not to be biased by continuum
discrepancies, and to focus solely on the line profile. The models are then
convolved to the instruments resolution (R=2360 and R=1800 for GQ Lup b and
GSC06214b SINFONI data respectively, R=3800 for OSIRIS).

The fits were performed on the mean emission line of each epoch, computing the
$\chi^{2}_{r}$ on the whole grid of models. Fiducial results are presented in
Fig. 9. The GQ Lup b lines are only reproduced by the MagAcc models: this is
because the Shock models cannot reproduce a blue shifted emission line peak,
as is clearly apparent from the first row of Fig. 9. Note that the red-shifted
absorption in the 2017 epochs of GQ Lup b and fitted by the models is mostly
likely caused by noise in the data. However, GSC06214b lines are well
reproduced by both models. Maps of the best $\chi^{2}_{r}$ are shown in Fig.
9. We only represent the $\textrm{T}_{\textrm{eff}}$–$\dot{\textrm{M}}$ axis
in the figure for the MagAcc model, because others parameters are mostly
unconstrained.

In the case of the MagAcc scenario, there is a clear degeneracy between the
accretion column temperature and the accretion rate. This behavior was already
seen in the context of magnetospheric accretion on T Tauri stars (Muzerolle et
al., 2001) and the degeneracy could be lifted with very high resolution
observations and/or fitting multiple accretion lines simultaneously. This
shows that in the magnetospheric accretion framework, we cannot derive an
accurate measurement of the accretion rate, but only a lower limit. For GQ Lup
b, this leads to $\dot{\textrm{M}}\gtrsim 10^{-11}$ $\mathrm{M_{\odot}}$/yr
for the (bright) 2006-09-18 epoch and $\dot{\textrm{M}}\gtrsim
10^{-11.5}$$\mathrm{M_{\odot}}$/yr for the other epochs. For GSC06214b, this
leads to a lower limit of $\dot{\textrm{M}}\gtrsim
10^{-12}$$\mathrm{M_{\odot}}$/yr for all epochs. Zhou et al. (2014) derived
values of $10^{-10.8}$ $\mathrm{M_{\odot}}$/yr and $10^{-9.3}$
$\mathrm{M_{\odot}}$/yr from fitting the UV-optical excess emission for
GSC06214b and GQ Lup b respectively, which are consistent with our lower
limits.

GSC06214b lines are also well reproduced with the Shock model. In this
scenario, fits are mostly degenerated over the parameter space, except for the
2010-07-09 (OSIRIS) epoch for which high velocity models are strongly favored,
as shown in Fig. 9. The high shock velocity models correspond to infall
velocities $v_{ff}>150$ km/s, which imply a shock location at $R<2.4$
$\mathrm{R_{Jup}}$, consistent with the estimated radius of the object (Sect.
5.2). These solutions would indicate that the shock is located at the
planetary surface. However, lower infall velocities in the range 30–170 km/s
cannot be fully excluded. They are consistent with a shock located anywhere
within 60 $\mathrm{R_{Jup}}$, so either at the protoplanet surface or on the
CPD at larger scales. In the 2010-07-09 (OSIRIS) epoch, our best fit values of
$\textrm{S}_{\textrm{eff}}$, $\textrm{v}_{0}$ and $\textrm{n}_{0}$ ($1.5\times
10^{-2}$ S${}_{\textrm{jup}}$, 180 km/s and $10^{20}$ $\textrm{m}^{-3}$) yield
a mass accretion rate of $\dot{\textrm{M}}$ $\sim$ $10^{-9}$
$\mathrm{M_{\odot}}$/yr. Assuming a shock at the planetary surface, the
accretion rate can be given as $\dot{\textrm{M}}=S\mu\nu_{0}v_{0}$, with
$\upmu$ the mean molecular weight ($3.84\times 10^{-27}$ kg as per Aoyama et
al., 2018). Inverting the equation, we find the shock surface given as a
fraction of the planetary surface (with an estimated radius of 1.8–2.0
$\mathrm{R_{Jup}}$), this corresponds to a filling factor of 0.38–0.47 %. This
is of the same order (0.3%) as reported for the older analogue companion
Delorme 1 (AB)b and using similar models on Balmer lines (Ringqvist et al.,
2023).

##  6 Discussion

In this section, we first interpret the derived atmospheric abundances of both
objects. We then discuss the implications brought by our results for the
origin of the $\mathrm{H\,\textsc{i}}$ line emission in planetary mass objects
and the possible impact of variable veiling. We finally extrapolate how
emission line variability as seen here could impact protoplanet detection
surveys.

###  6.1 Atmospheric abundances as formation tracers

The star forming regions Sco-Cen and Lupus I are expected to have nearly solar
metallicity (e.g., Bubar et al., 2011; Biazzo et al., 2017; Murphy et al.,
2021, and ref. therein) and it has been shown that nearby solar-type stars
with $\textrm{[M/H]}=0$ have C/O distributed around the solar value (Suárez-
Andrés et al., 2018; Tautvaišienė et al., 2022). We note however a
spectroscopic determination of $\textrm{[M/H]}=-0.06\pm 0.01$ and $-0.35\pm
0.01$ dex for GSC06214A and GQ Lup A, respectively (Swastik et al., 2021). The
value found for GQ Lup A is at odd with expectations and implies that GQ Lup b
is significantly enriched in metals
($\Delta\textrm{[M/H]}=0.61^{+0.05}_{-0.10}$) with respect to its host star
while GSC06214b metallicity remains compatible within 1$\sigma$ with that of
GSC06214A. Swastik et al. (2021) accounted for veiling affecting the depth of
stellar absorption lines of GQ Lup A in their inversion of the high-resolution
optical spectra of the source. They chose in particular the less contaminated
and emission-free region 5900–5965 Å for the analysis. However, the
$\mathsf{H}_{\upalpha}$ line of GQ Lup A appears as the most intense of their
sample and remaining effects of veiling at considered wavelengths might lead
to an underestimation of the [M/H] value for that particular star.

The metallicity of GSC06214b and GQ Lup b are both compatible with the
envelope of values predicted for planet population synthesis models based on
the core-accretion framework assuming a pristine solar-composition (see Fig. 9
of Petrus et al., 2021). Disk-instability models predict a wide range of
metallicity values, including solar (Boley et al., 2011). If the companions
formed inside the circumstellar disks surrounding their host stars, the
nearly-solar C/O implies they should have accreted a substantial fraction of
solids. In summary, the derived abundances based on our J-band data do not
permit to firmly conclude on the formation scenario of either objects, but an
in-situ formation via turbulent fragmentation remains plausible (i.e.,
fragmentation of the solar nebula into the primary and its companion, Padoan &
Nordlund, 1999, 2004).

###  6.2 Constraints on accretion processes in planetary mass objects

It is shown in Sect. 5.4.2 that GQ Lup b lines are only reproduced by the
MagAcc model, while GSC06214b lines are well reproduced by both MagAcc and
Shock models.

####  6.2.1 MagAcc and Shock models relative contributions

Although we studied the MagAcc and the Shock scenarios independently, they are
not mutually exclusive. The magnetospheric accretion scenario also involves a
shock at the base of the accretion columns on the planetary surface, while the
Shock model makes no assumption as to where the accretion shock takes place.
Therefore, the Shock model could very well explain the line emission at the
base of a magnetospherically induced funnel.

The relative contributions of the accretion funnels and the accretion shock to
the $\mathsf{Pa}_{\upbeta}$ emission line may vary with time. The MagAcc
induced lines are usually asymmetric with a larger blue wing whereas Shock-
induced lines tend to be more symmetric with a slightly stronger red wing. A
time-modulation of the MagAcc-Shock contributions could then be the key to
explain the correlations observed between the line centroid, its width and
flux (Fig. 6).

Figure 2 of Aoyama et al. (2021) shows the relative contributions of the
funnel flows and the planetary shock $\mathsf{H}_{\upalpha}$ emission,
assuming a filling factor of 1–10 % and planet mass of 6 $\mathrm{M_{Jup}}$
(corresponding to PDS 70 b) depending on the accretion rate. Our higher mass
range and different line may quantitatively change the accretion rate
threshold at which both emissions become comparable, but this will hardly
change the qualitative picture. For the mass accretion rate of GSC06214b
($10^{-10.8}$ $\mathrm{M_{\odot}}$/yr, Zhou et al., 2014), both MagAcc and
Shock contributions are about the same order of magnitude. For the mass
accretion rate of GQ Lup b ($10^{-9.3}$ $\mathrm{M_{\odot}}$/yr, Zhou et al.,
2014), the contribution from the accretion funnels largely dominates, which
could explain why only the MagAcc model reproduces the blue-shifted lines of
GQ Lup b.

####  6.2.2 Others explanations for line variability

While we focused the MagAcc and Shock scenarii, other physical mechanisms have
been proposed to modulate the accretion lines properties. Marleau et al.
(2022) studied the effects of gas and dust extinction on the
$\mathsf{H}_{\upalpha}$ line intensity profile. According to their Fig. 5, for
the gas extinction to alter significantly the $\mathsf{H}_{\upalpha}$ flux
($A_{R}>1$) for an object of $\sim 15$ $\mathrm{M_{Jup}}$, the accretion rate
must be at least $10^{-5}$ $\mathrm{M_{Jup}}$/yr, which is far above the
accretion rates of GQ Lup b or GSC06214b. Figures 6, 7 and 8 in Marleau et al.
(2022) show intrinsic $\mathsf{H}_{\upalpha}$ line profiles for the same Shock
model as in this work, with and without gas extinction, for various accretion
flow geometries, planet masses and accretion rates. Some of these
$\mathsf{H}_{\upalpha}$ lines show significant asymmetries (in their intrinsic
profile), including an absorption of the red-wing of the line. If the behavior
in the $\mathsf{Pa}_{\upbeta}$ line is similar, the correlations observed in
Fig. 6 could be explained by such gas self-absorption effects. As for dust
extinction, they find that it should be negligible for masses above a few
$\mathrm{M_{Jup}}$.

Szulágyi & Ercolano (2020) performed 3D simulations of an accreting
protoplanet embedded in a circumstellar disk and studied the lines variation
due to the variable shock-front (as the planet orbits around the star) and
density variations in the protoplanet vicinity. They find the
$\mathsf{Pa}_{\upbeta}$ line amplitude variability to be about $\sim 8-41$%
(depending on the modeling hypothesis) for a 10 $\mathrm{M_{Jup}}$ planet, at
keplerian timescales. The different environments of GSC06214b and GQ Lup b
(non-embedded) make this scenario implausible for these two objects and the
proposed timescales do not match our observations.

####  6.2.3 Line variability on similar objects

Variability of PMCs has already been seen at various timescales. In the
$\mathsf{Pa}_{\upbeta}$ study of DH Tau b by Wolff et al. (2017), the authors
report 4 epochs separated by days and weeks. According to their Table 1, they
find the $\mathsf{Pa}_{\upbeta}$ line to disappear within $\sim 5$ weeks
($3.8\upsigma$) and a tentative variability within a week ($80\pm 50$%,
$1.7\upsigma$). Eriksson et al. (2020) conducted a multi-epoch
$\mathsf{H}_{\upalpha}$ study of Delorme 1 (AB)b. Applying the same methods as
in Fig. 7 to their Table B.1 data, we find $\mathsf{H}_{\upalpha}$ EW
variability spanning $\sim 30\pm 20$% along the $\sim 2$ hours sequence, with
no difference between the $<1$ h and $>1$ h timescales. In comparison, we find
the variability amplitude to be $\sim 7\pm 7$% and $\sim 20\pm 17$% for GQ Lup
b and GSC06214b respectively for timescales $<2$h (note that error bars
represent the spread of the variability, not to an uncertainty), compatible
with Delorme 1 (AB)b.

Figure 10: Correlations between equivalent widths (EW) of both potassium
absorption doublets and the $\mathsf{Pa}_{\upbeta}$ line. The equivalent width
is positive for absorption lines and negative for emission lines. Gray lines
correspond to a slope equal to -1, where the EW variability can be explained
by sole continuum variability.

####  6.2.4 Similarities with CTTS

Both GQ Lup b and GSC06214b show maximum variability for timescales above
their estimated rotation periods (see Fig 7). This behavior is also seen in
CTTS (Costigan et al., 2014; Zsidi et al., 2022). The $\mathsf{Pa}_{\upbeta}$
line profiles of both companions also show strong similarities with the
centrally peaked line profiles observed in solar mass accreting young stars
(CTTS, see Folha & Emerson, 2001). The line profile of GSC06214b on May 04,
2017 and August 19, 2017 matches surprisingly well the one of DG Tau observed
at October 10, 1994 and December 15, 1995 by Folha & Emerson (2001). The blue-
shifted $\mathsf{Pa}_{\upbeta}$ line of GQ Lup b is also well reproduced by
their observations of DD Tau. Both objects are sub-solar-mass CTTS.
Magnetospheric accretion models have been very successful at reproducing these
line profiles (Muzerolle et al., 2001). These similarities suggest that the
mechanisms proposed for the formation of the $\mathrm{H\,\textsc{i}}$ emission
lines in CTTS (i.e., magnetospheric accretion) are at play down to the
planetary mass regime. The estimates of mass accretion rates for GSC06214b and
GQ Lup b are also in line with the global decrease of mass accretion rates
with object mass which extends down to isolated $\sim$15$\mathrm{M_{Jup}}$
objects (e.g., Muzerolle et al., 2005; Joergens et al., 2013).

Ringqvist et al. (2023) report recent high-resolution UV spectra (2 epochs
taken within 1h) of the $13\pm 5$ $\mathrm{M_{Jup}}$ companion Delorme 1 (AB)b
where they find line blue-shifted profiles indicative of magnetospheric
accretion. While the profiles of these lines cannot be directly compared to
those of $\mathsf{Pa}_{\upbeta}$, this further suggests this mechanism could
be universal for young 10-30 $\mathrm{M_{Jup}}$ companions with ages from 3 to
40 Myr. Both Delorme 1 (AB)b and GSC06214b have estimated filling factors of
the accretion shock of the order $\sim$1%. These values are consistent with
those found in CTTS (e.g., Alencar & Batalha, 2002; Batalha et al., 2002;
Pittman et al., 2022) although the methodologies used to derive these factors
on stars can be different.

Magnetospheric accretion holds only if the objects interior can generate
powerful magnetic fields. Following Eq. 1 from Hartmann et al. (2016) with the
values found from evolutionary tracks (Sect. 5.2, GQ Lup b: $\textrm{M}=22$
$\mathrm{M_{Jup}}$, $\textrm{R}=3$ $\mathrm{R_{Jup}}$,
$\dot{\textrm{M}}=10^{-9.3}$ $\mathrm{M_{\odot}}$/yr; GSC06214b:
$\textrm{M}=15$ $\mathrm{M_{Jup}}$, $\textrm{R}=2$ $\mathrm{R_{Jup}}$,
$\dot{\textrm{M}}=10^{-10.8}$ $\mathrm{M_{\odot}}$/yr), the lower limits on
their rotational period (5h for GQ Lup b, 20h for GSC06214b) and $\upxi=0.7$,
they would require a minimum dipolar magnetic field of $\sim 70$ G and $\sim
150$ G respectively. These increase to $\sim 750$ G and $1.7$ kG for the
rotational period upper limits (80h). Magnetic field strengths on young
companions are just starting to be investigated (Katarzyński et al., 2016;
Cendes et al., 2022) but no measurement exists yet. Recent radio and
$\mathsf{H}_{\upalpha}$ observations of low-mass brown-dwarfs and planetary-
mass objects in the same mass range as our companions suggest that the
magnetic field strength could reach several k-Gauss at early ages (Kao et al.,
2018). Time-series of line profiles at higher resolutions (R$>$10000) should
further help to place constraints on the MagAcc and Shock models (Marleau et
al., 2022). More direct measurements of the magnetic fields strength on
accreting free-floating planets or companions on wider-orbits with spectro-
polarimeters (e.g., SPIRou) might be possible and would also help to confirm
the MagAcc scenario.

###  6.3 Evidence for continuum variability?

In our analysis of the $\mathsf{Pa}_{\upbeta}$ line variability, we have
assumed that the continuum emission remained constant. Photospheric emission
of young dwarfs at the M-L transition is not expected to show high-amplitude
modulation due to non-homogeneous cloud coverage (e.g., Artigau, 2018; Vos et
al., 2022). It remains possible that variable veiling or extinction due to
accretion and the circum-”planetary” environment or magnetic spots affect the
J-band and could change the contrast between the lines and their surrounding
continua, thus producing at least a portion the observed variability.
Carpenter et al. (2002) find a 0.15mag variability at J-band of the faintest
members of Chameleon whose colors (J-K$\sim$0.8-1.0) are close to the one of
GQ Lup b and correspond to spectral type in the mid- to late-M range (Cao et
al., 2022).

Zhou et al. (2014) measured the shock contribution from 1D slab models
(continuum excess produced by the accretion shock) and the
$\mathsf{H}_{\upalpha}$ luminosity of both companions using HST data. Their
model indicate a negligible contribution of the slab for GSC06214b beyond 1
$\upmu$m and a contamination at the percent level for GQ Lup b at J-band.
However, the shock contribution is expected to be linked to the line
intensity, this could work as a variable continuum component at J band that
could explain observed $\mathsf{Pa}_{\upbeta}$ equivalent width variations.
Photometric variability of GQ Lup b is known to be negligible at K-band
($\upsigma_{\textrm{K}}=0.073$ mag, Neuhäuser et al., 2008), but no monitoring
has been conducted at J-band where the effect of veiling could be more
important.

Veiling would both decrease the prominent $\mathrm{K\,\textsc{i}}$ doublets
and the $\mathsf{Pa}_{\upbeta}$ emission lines. Figure 10 shows the EW of
$\mathsf{Pa}_{\upbeta}$ and the KI absorption lines (corresponding values are
given in Table 2). A variation along a line with slope -1 is expected if the
observed EW variability is only due to a variation in the continuum level. GQ
Lup b lines show no clear correlations in 2017, although there might be a
slight trend in 2006. Because the 2006 epoch and the 2017 epochs are not
following the same slope, their EW difference cannot be explained only by
continuum variability. To explain the EW variation observed within the 2006
epoch, a variation of the continuum level of $\sim 30-50$ % within 2 hours
would be required, corresponding to changes of J-band magnitudes of $0.3-0.4$
mag, compatible with the dispersion of values reported in Table LABEL:tab:obs.

For GSC06214b, the trend appears between all epochs except 2017-08-19 and
2010-07-09. The $\mathsf{Pa}_{\upbeta}$ EW variability between those might
therefore be explained by continuum variability. Although, in a similar
fashion as for GQ Lup b, the $\mathsf{Pa}_{\upbeta}$ EW variation between the
2017-08-19/2010-07-09 and other 2017 epochs cannot be explained only be
continuum variability.

Given the current constraints, we assume that both objects have EW variability
that is driven by intrinsic line flux variability. Accurate photometric
monitoring of both companions at J-band simultaneous to
$\mathsf{Pa}_{\upbeta}$ emission line variability campaign would be required
to fully confirm this hypothesis.

###  6.4 Implications for protoplanet detection

Several surveys have been conducted to detect accreting planets nested in
circumstellar disks using narrow-band imaging at $\mathsf{H}_{\upalpha}$
(Mendigutía et al., 2018; Wagner et al., 2018; Cugno et al., 2019; Uyama et
al., 2020; Zurlo et al., 2020; Huélamo et al., 2022; Follette et al., 2022),
or integral field spectroscopy targeting $\mathsf{H}_{\upalpha}$ (Haffert et
al., 2019; Xie et al., 2020) and $\mathsf{Pa}_{\upbeta}$ emissions (Uyama et
al., 2017, 2021). Thus far, $\mathsf{H}_{\upalpha}$ point source detections
have been claimed in four systems (LkCa15, AB Aur, PDS70, CS Cha; Sallum et
al., 2015; Huélamo et al., 2018; Currie et al., 2022; Zhou et al., 2021;
Follette et al., 2022). The nature of the detection remains unclear in two of
these systems (LkCa15, AB Aur; Whelan et al., 2015; Follette et al., 2022;
Zhou et al., 2022). With the exception of Delorme 1 (AB)b, which does not seem
to be in a protoplanetary-disk, no detection has ever been achieved at
$\mathsf{Pa}_{\upbeta}$.

PDS 70 b and c are known to be affected by significant in-line extinction
(Hashimoto et al., 2020; Wang et al., 2021) caused by surrounding material
while they are otherwise found within a large cavity. Extinction is therefore
a good candidate to explain the poor yield of $\mathsf{H}_{\upalpha}$ and
$\mathsf{Pa}_{\upbeta}$ deep imaging searches, especially in disks with
significant amount of gas and dust. However, most are single-epoch
observations and emission line variability might also provide an additional
explanation for the low yield of these observations.

Xie et al. (2020) report in their Fig. 11 the $\mathsf{H}_{\upalpha}$
detection limit achieved for 30min of effective VLT/MUSE observations of PDS
70. Their observations of PDS 70 are the deepest presented in the paper and
probe deeper apparent line fluxes than all the narrow-band observations
presented in Cugno et al. (2019) and Huélamo et al. (2022). They also report
the measured apparent $\mathsf{H}_{\upalpha}$ emission line flux of the two
accreting planets.

Both PDS 70 b & c, unlike GQ Lup b and GSC06214b, are embedded and interacting
with the circumstellar disk. It is therefore reasonable to assume that the
majority of variability happens at keplerian timescales. Given the variety of
ages and environments of known accreting objects, our $\mathsf{Pa}_{\upbeta}$
variability analysis might only serve as a first order approximation of the
general $\mathsf{H}_{\upalpha}$ variability. Should our observed
$\mathsf{Pa}_{\upbeta}$ variations still be representative of
$\mathsf{H}_{\upalpha}$ variability of others companions, the 100 to 400%
variability amplitude at year to decade timescales measured for GSC06214b and
GQ Lup b would imply that at other epochs, both PDS 70 planets might have been
emitting below the VLT/MUSE detection threshold. The various
$\mathsf{H}_{\upalpha}$ flux measurements of PDS70 b shown in Zhou et al.
(2021) tend to indicate that such amplitude of variability is plausible for
the companion at the year timescale. Our findings clearly advocate new visits
of these young systems to maximize the detection probability. Similarly,
variability might have also helped detection surveys.

##  7 Conclusions

We presented the results of a J-band spectroscopic monitoring campaign of GQ
Lup b and GSC06214b at medium resolving power ($\textrm{R}=1800-2360$) with
VLT/SINFONI. We re-detect the $\mathsf{Pa}_{\upbeta}$ emission line at each
epoch. We report below our main conclusions:

  * •

We confirm that GQ Lup b is a $\sim 2640$ K, $22^{+2}_{-3}$ $\mathrm{M_{Jup}}$
companion, and GSC06214b a $\sim 2250$ K, $15.5\pm 0.5$ $\mathrm{M_{Jup}}$
companion, given photosphere fitting and evolutionary tracks. Note that error
bars may not be representative of systematic uncertainties of the models.

  * •

We measure C/O values on both objects compatible with solar value (GQ Lup b:
C/O = $0.44^{+0.13}_{-0.11}$, GSC06214b: C/O = $0.48^{+0.16}_{-0.12}$,
C/O${}_{\odot}=0.55$). We find GSC06214b’s [M/H] to be compatible with solar
value ($-0.16\pm 0.17$), while GQ Lup b’s is slightly super-solar ($0.23\pm
0.06$). The uncertainties on the atmospheric parameters are driven by small
systematic errors in the data revealed by our multi-epoch observations, while
Bayesian inference at single epoch provide too optimistic estimates of the
error bars.

  * •

For both objects, the $\mathsf{Pa}_{\upbeta}$ intensity variability is found
to be moderate ($<50$ %) at timescales below their rotation period, and more
important at longer timescales (up to $\sim 1000$ % for a decade), a behavior
similar to that of CTTS. We studied whether variable veiling could change the
contrast of the line above the continuum. We find that such flux component
does not explain the long-term variability of either object but we cannot
exclude it might explain the shorter term evolution of the lines.

  * •

Both companions had their maximum recorded line brightness at first
(historical) epochs. Their variability amplitude can be up to $\sim 1000$ % at
the decade timescale. Should that amplitude extend to $\mathsf{H}_{\upalpha}$,
this might explain the low yield of $\mathsf{H}_{\upalpha}$ and
$\mathsf{Pa}_{\upbeta}$ (spectro)-imaging campaigns for protoplanet detection.

  * •

The blue-shifted $\mathsf{Pa}_{\upbeta}$ lines of GQ Lup b are reproduced only
by the protoplanet magnetospheric accretion models. Those of GSC06214b are
well reproduced by both magnetospheric accretion and shock models at
quiescence, while the shock model is favored when the line is more intense
(2010 epoch). We find filling factors of the accretion shock below 1% at that
epoch and compatible with the value found on the older companion Delorme 1
(AB)b.

  * •

We identify a tentative correlation between GSC06214b’s emission line center
and its width, which could be explained either by asymmetries in
magnetospheric accretion-induced lines, or by gas extinction in shock-induced
lines.

Our study presents the first monitoring of emission lines of young accreting
companions at the edge of the deuterium-burning boundary. It illustrates how
variability timescales could be used to investigate whether the phenomena at
play are a down-sizing of those happening on protostars or more specific ones.
Our results advocate multiple visits of young systems on a year-span at least
to maximize the detection probability of young protoplanets. Observations at
higher resolution ($\textrm{R}>10000$, e.g. ERIS) would now be required to
better constrain the line profiles and emitting mechanism.

###### Acknowledgements.

We thank the ESO staff for the help during the preparation of the observations
and their execution. We are grateful to B. Tessore and J. Bouvier for useful
discussions on their magnetospheric accretion models. We acknowledge support
in France from the French National Research Agency (ANR) through project grant
ANR-20-CE31-0012 and the Programmes Nationaux de Planetologie et de Physique
Stellaire (PNP and PNPS). This project has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (COBREX; grant agreement 885 593). G.-D. M. acknowledges
the support of the DFG priority program SPP 1992 “Exploring the Diversity of
Extrasolar Planets” (MA 9185/1), and from the Swiss National Science
Foundation under grant 200021_204847 “PlanetsInTime”. Parts of this work have
been carried out within the framework of the NCCR PlanetS supported by the
Swiss National Science Foundation. This publication makes use of VOSA,
developed under the Spanish Virtual Observatory (https://svo.cab.inta-csic.es)
project funded by MCIN/AEI/10.13039/501100011033/ through grant
PID2020-112949GB-I00. VOSA has been partially updated by using funding from
the European Union’s Horizon 2020 Research and Innovation Programme, under
Grant Agreement nr. 776403 (EXOPLANETS-A).

## References

  * Abuter et al. (2006) Abuter, R., Schreiber, J., Eisenhauer, F., et al. 2006, New A Rev., 50, 398
  * Alcalá et al. (2020) Alcalá, J. M., Majidi, F. Z., Desidera, S., et al. 2020, A&A, 635, L1
  * Alencar & Batalha (2002) Alencar, S. H. P. & Batalha, C. 2002, ApJ, 571, 378
  * Allard et al. (2011) Allard, F., Homeier, D., & Freytag, B. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West, 91
  * Anders et al. (2022) Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91
  * Antoniucci et al. (2011) Antoniucci, S., García López, R., Nisini, B., et al. 2011, A&A, 534, A32
  * Aoyama et al. (2018) Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84
  * Aoyama et al. (2021) Aoyama, Y., Marleau, G.-D., Ikoma, M., & Mordasini, C. 2021, ApJ, 917, L30
  * Aoyama et al. (2020) Aoyama, Y., Marleau, G.-D., Mordasini, C., & Ikoma, M. 2020, arXiv e-prints, arXiv:2011.06608
  * Artigau (2018) Artigau, É. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 94
  * Baraffe et al. (2009) Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Barenfeld et al. (2016) Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142
  * Batalha et al. (2002) Batalha, C., Batalha, N. M., Alencar, S. H. P., Lopes, D. F., & Duarte, E. S. 2002, ApJ, 580, 343
  * Batalha et al. (2001) Batalha, C., Lopes, D. F., & Batalha, N. M. 2001, ApJ, 548, 377
  * Batygin (2018) Batygin, K. 2018, AJ, 155, 178
  * Bayo et al. (2008) Bayo, A., Rodrigo, C., Barrado Y Navascués, D., et al. 2008, A&A, 492, 277
  * Béjar et al. (2008) Béjar, V. J. S., Zapatero Osorio, M. R., Pérez-Garrido, A., et al. 2008, ApJ, 673, L185
  * Betti et al. (2022a) Betti, S. K., Follette, K. B., Ward-Duong, K., et al. 2022a, ApJ, 941, L20
  * Betti et al. (2022b) Betti, S. K., Follette, K. B., Ward-Duong, K., et al. 2022b, ApJ, 935, L18
  * Biazzo et al. (2017) Biazzo, K., Frasca, A., Alcalá, J. M., et al. 2017, A&A, 605, A66
  * Boley et al. (2011) Boley, A. C., Helled, R., & Payne, M. J. 2011, ApJ, 735, 30
  * Bonnefoy et al. (2014) Bonnefoy, M., Chauvin, G., Lagrange, A. M., et al. 2014, A&A, 562, A127
  * Bonnet et al. (2004) Bonnet, H., Conzelmann, R., Delabre, B., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5490, Advancements in Adaptive Optics, ed. D. Bonaccini Calia, B. L. Ellerbroek, & R. Ragazzoni, 130–138
  * Bowler et al. (2015) Bowler, B. P., Andrews, S. M., Kraus, A. L., et al. 2015, ApJ, 805, L17
  * Bowler et al. (2017) Bowler, B. P., Kraus, A. L., Bryan, M. L., et al. 2017, AJ, 154, 165
  * Bowler et al. (2014) Bowler, B. P., Liu, M. C., Kraus, A. L., & Mann, A. W. 2014, ApJ, 784, 65
  * Bowler et al. (2011) Bowler, B. P., Liu, M. C., Kraus, A. L., Mann, A. W., & Ireland, M. J. 2011, ApJ, 743, 148
  * Brittain et al. (2020) Brittain, S. D., Najita, J. R., Dong, R., & Zhu, Z. 2020, ApJ, 895, 48
  * Broeg et al. (2007) Broeg, C., Schmidt, T. O. B., Guenther, E., et al. 2007, A&A, 468, 1039
  * Bryan et al. (2018) Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, Nature Astronomy, 2, 138
  * Bryan et al. (2020) Bryan, M. L., Ginzburg, S., Chiang, E., et al. 2020, ApJ, 905, 37
  * Bubar et al. (2011) Bubar, E. J., Schaeuble, M., King, J. R., Mamajek, E. E., & Stauffer, J. R. 2011, AJ, 142, 180
  * Cao et al. (2022) Cao, L., Pinsonneault, M. H., Hillenbrand, L. A., & Kuhn, M. A. 2022, ApJ, 924, 84
  * Carpenter et al. (2002) Carpenter, J. M., Hillenbrand, L. A., Skrutskie, M. F., & Meyer, M. R. 2002, AJ, 124, 1001
  * Cendes et al. (2022) Cendes, Y., Williams, P. K. G., & Berger, E. 2022, AJ, 163, 15
  * Chabrier et al. (2000) Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464
  * Chang & Deming (1996) Chang, E. S. & Deming, D. 1996, Sol. Phys., 165, 257
  * Chinchilla et al. (2021) Chinchilla, P., Béjar, V. J. S., Lodieu, N., Zapatero Osorio, M. R., & Gauza, B. 2021, A&A, 645, A17
  * Costigan et al. (2014) Costigan, G., Vink, J. S., Scholz, A., Ray, T., & Testi, L. 2014, MNRAS, 440, 3444
  * Cridland (2018) Cridland, A. J. 2018, A&A, 619, A165
  * Cugno et al. (2019) Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A&A, 622, A156
  * Cumming et al. (2018) Cumming, A., Helled, R., & Venturini, J. 2018, MNRAS, 477, 4817
  * Curran (2015) Curran, P. A. 2015, MCSpearman: Monte Carlo error analyses of Spearman’s rank test, Astrophysics Source Code Library, record ascl:1504.008
  * Currie et al. (2022) Currie, T., Lawson, K., Schneider, G., et al. 2022, Nature Astronomy, 6, 751
  * Cushing et al. (2005) Cushing, M. C., Rayner, J. T., & Vacca, W. D. 2005, ApJ, 623, 1115
  * Cutri & et al. (2012) Cutri, R. M. & et al. 2012, VizieR Online Data Catalog, II/311
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Davies (2007) Davies, R. I. 2007, MNRAS, 375, 1099
  * Delorme et al. (2013) Delorme, P., Gagné, J., Girard, J. H., et al. 2013, A&A, 553, L5
  * Donati et al. (2012) Donati, J. F., Gregory, S. G., Alencar, S. H. P., et al. 2012, MNRAS, 425, 2948
  * Dong et al. (2021) Dong, J., Jiang, Y.-F., & Armitage, P. J. 2021, ApJ, 921, 54
  * Edwards et al. (1994) Edwards, S., Hartigan, P., Ghandour, L., & Andrulis, C. 1994, AJ, 108, 1056
  * Eisenhauer et al. (2003) Eisenhauer, F., Abuter, R., Bickert, K., et al. 2003, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, 4841, 1548
  * Eriksson et al. (2020) Eriksson, S. C., Asensio Torres, R., Janson, M., et al. 2020, A&A, 638, L6
  * Feiden (2016) Feiden, G. A. 2016, A&A, 593, A99
  * Fendt (2003) Fendt, C. 2003, A&A, 411, 623
  * Fischer et al. (2022) Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2022, arXiv e-prints, arXiv:2203.11257
  * Folha & Emerson (2001) Folha, D. F. M. & Emerson, J. P. 2001, A&A, 365, 90
  * Follette et al. (2022) Follette, K. B., Close, L. M., Males, J. R., et al. 2022, arXiv e-prints, arXiv:2211.02109
  * Frasca et al. (2017) Frasca, A., Biazzo, K., Alcalá, J. M., et al. 2017, A&A, 602, A33
  * Fu et al. (2023) Fu, Z., Huang, S., & Yu, C. 2023, ApJ, 945, 165
  * Gagné et al. (2018a) Gagné, J., Allers, K. N., Theissen, C. A., et al. 2018a, ApJ, 854, L27
  * Gagné et al. (2018b) Gagné, J., Mamajek, E. E., Malo, L., et al. 2018b, ApJ, 856, 23
  * Gaia Collaboration (2022) Gaia Collaboration. 2022, VizieR Online Data Catalog, I/355
  * Ginski et al. (2014) Ginski, C., Schmidt, T. O. B., Mugrauer, M., et al. 2014, MNRAS, 444, 2280
  * Gressel et al. (2013) Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135
  * Hashimoto et al. (2020) Hashimoto, J., Aoyama, Y., Konishi, M., et al. 2020, AJ, 159, 222
  * Heller & Pudritz (2015) Heller, R. & Pudritz, R. 2015, ApJ, 806, 181
  * Helling et al. (2008) Helling, C., Ackerman, A., Allard, F., et al. 2008, MNRAS, 391, 1854
  * Herbig (1977) Herbig, G. H. 1977, ApJ, 214, 747
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Huélamo et al. (2022) Huélamo, N., Chauvin, G., Mendigutía, I., et al. 2022, A&A, 668, A138
  * Huélamo et al. (2018) Huélamo, N., Chauvin, G., Schmid, H. M., et al. 2018, A&A, 613, L5
  * Ireland et al. (2011) Ireland, M. J., Kraus, A., Martinache, F., Law, N., & Hillenbrand, L. A. 2011, ApJ, 726, 113
  * Itoh et al. (2005) Itoh, Y., Hayashi, M., Tamura, M., et al. 2005, ApJ, 620, 984
  * Joergens et al. (2013) Joergens, V., Bonnefoy, M., Liu, Y., et al. 2013, A&A, 558, L7
  * Jones, A. et al. (2013) Jones, A., Noll, S., Kausch, W., Szyszka, C., & Kimeswenger, S. 2013, A&A, 560, A91
  * Kao et al. (2018) Kao, M. M., Hallinan, G., Pineda, J. S., Stevenson, D., & Burgasser, A. 2018, ApJS, 237, 25
  * Katarzyński et al. (2016) Katarzyński, K., Gawroński, M., & Goździewski, K. 2016, MNRAS, 461, 929
  * Kramida et al. (2022) Kramida, A., Yu. Ralchenko, Reader, J., & and NIST ASD Team. 2022, NIST Atomic Spectra Database (ver. 5.10), [Online]. Available: https://physics.nist.gov/asd [2022, November 30]. National Institute of Standards and Technology, Gaithersburg, MD.
  * Kraus et al. (2008) Kraus, A. L., Ireland, M. J., Martinache, F., & Lloyd, J. P. 2008, ApJ, 679, 762
  * Kurosawa et al. (2006) Kurosawa, R., Harries, T. J., & Littlefair, S. P. 2006, MNRAS, 372, 1879
  * Lachapelle et al. (2015) Lachapelle, F.-R., Lafrenière, D., Gagné, J., et al. 2015, ApJ, 802, 61
  * Larkin et al. (2006) Larkin, J., Barczys, M., Krabbe, A., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6269, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. I. S. McLean & M. Iye, 62691A
  * Lavigne et al. (2009) Lavigne, J.-F., Doyon, R., Lafrenière, D., Marois, C., & Barman, T. 2009, ApJ, 704, 1098
  * Lazzoni et al. (2020) Lazzoni, C., Gratton, R., Alcalá, J. M., et al. 2020, A&A, 635, L11
  * Long et al. (2020) Long, D. E., Zhang, K., Teague, R., & Bergin, E. A. 2020, ApJ, 895, L46
  * Lubow & Martin (2012) Lubow, S. H. & Martin, R. G. 2012, ApJ, 749, L37
  * Luhman et al. (2003) Luhman, K. L., Stauffer, J. R., Muench, A. A., et al. 2003, ApJ, 593, 1093
  * MacGregor et al. (2017) MacGregor, M. A., Wilner, D. J., Czekala, I., et al. 2017, ApJ, 835, 17
  * Manara et al. (2013) Manara, C. F., Testi, L., Rigliaco, E., et al. 2013, A&A, 551, A107
  * Marleau & Aoyama (2022) Marleau, G.-D. & Aoyama, Y. 2022, Research Notes of the American Astronomical Society, 6, 262
  * Marleau et al. (2022) Marleau, G. D., Aoyama, Y., Kuiper, R., et al. 2022, A&A, 657, A38
  * Marley et al. (2007) Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541
  * Marois et al. (2007) Marois, C., Macintosh, B., & Barman, T. 2007, ApJ, 654, L151
  * McElwain et al. (2007) McElwain, M. W., Metchev, S. A., Larkin, J. E., et al. 2007, ApJ, 656, 505
  * Mendigutía et al. (2018) Mendigutía, I., Oudmaijer, R. D., Schneider, P. C., et al. 2018, A&A, 618, L9
  * Miki (1982) Miki, S. 1982, Progress of Theoretical Physics, 67, 1053
  * Murphy et al. (2021) Murphy, S. J., Joyce, M., Bedding, T. R., White, T. R., & Kama, M. 2021, MNRAS, 502, 1633
  * Muzerolle et al. (2001) Muzerolle, J., Calvet, N., & Hartmann, L. 2001, ApJ, 550, 944
  * Muzerolle et al. (2005) Muzerolle, J., Luhman, K. L., Briceño, C., Hartmann, L., & Calvet, N. 2005, ApJ, 625, 906
  * Natta et al. (2004) Natta, A., Testi, L., Muzerolle, J., et al. 2004, A&A, 424, 603
  * Neuhäuser et al. (2000) Neuhäuser, R., Guenther, E. W., Petr, M. G., et al. 2000, A&A, 360, L39
  * Neuhäuser et al. (2005) Neuhäuser, R., Guenther, E. W., Wuchterl, G., et al. 2005, A&A, 435, L13
  * Neuhäuser et al. (2008) Neuhäuser, R., Mugrauer, M., Seifahrt, A., Schmidt, T. O. B., & Vogt, N. 2008, A&A, 484, 281
  * Noll, S. et al. (2012) Noll, S., Kausch, W., Barden, M., et al. 2012, A&A, 543, A92
  * Padoan & Nordlund (1999) Padoan, P. & Nordlund, Å. 1999, ApJ, 526, 279
  * Padoan & Nordlund (2004) Padoan, P. & Nordlund, Å. 2004, ApJ, 617, 559
  * Pearce et al. (2019) Pearce, L. A., Kraus, A. L., Dupuy, T. J., et al. 2019, AJ, 157, 71
  * Pecaut et al. (2012) Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154
  * Petrus et al. (2020) Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2020, A&A, 633, A124
  * Petrus et al. (2021) Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59
  * Petrus et al. (2023) Petrus, S., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, L9
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38
  * Pierens & Raymond (2016) Pierens, A. & Raymond, S. N. 2016, MNRAS, 462, 4130
  * Piskunov & Valenti (2002) Piskunov, N. E. & Valenti, J. A. 2002, A&A, 385, 1095
  * Pittman et al. (2022) Pittman, C. V., Espaillat, C. C., Robinson, C. E., et al. 2022, AJ, 164, 201
  * Preibisch et al. (1998) Preibisch, T., Guenther, E., Zinnecker, H., et al. 1998, A&A, 333, 619
  * Privon et al. (2020) Privon, G. C., Ricci, C., Aalto, S., et al. 2020, ApJ, 893, 149
  * Reiners & Christensen (2010) Reiners, A. & Christensen, U. R. 2010, A&A, 522, A13
  * Rice et al. (2010) Rice, E. L., Barman, T., Mclean, I. S., Prato, L., & Kirkpatrick, J. D. 2010, ApJS, 186, 63
  * Rigliaco et al. (2012) Rigliaco, E., Natta, A., Testi, L., et al. 2012, A&A, 548, A56
  * Ringqvist et al. (2023) Ringqvist, S. C., Viswanath, G., Aoyama, Y., et al. 2023, A&A, 669, L12
  * Sallum et al. (2015) Sallum, S., Follette, K. B., Eisner, J. A., et al. 2015, Nature, 527, 342
  * Santamaría-Miranda et al. (2018) Santamaría-Miranda, A., Cáceres, C., Schreiber, M. R., et al. 2018, MNRAS, 475, 2994
  * Santamaría-Miranda et al. (2019) Santamaría-Miranda, A., Cáceres, C., Schreiber, M. R., et al. 2019, MNRAS, 488, 5852
  * Schmidt et al. (2008) Schmidt, T. O. B., Neuhäuser, R., Seifahrt, A., et al. 2008, A&A, 491, 311
  * Schwarz et al. (2016) Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74
  * Seifahrt et al. (2007) Seifahrt, A., Neuhäuser, R., & Hauschildt, P. H. 2007, A&A, 463, 309
  * Sembach & Savage (1992) Sembach, K. R. & Savage, B. D. 1992, ApJS, 83, 147
  * Sicilia-Aguilar et al. (2015) Sicilia-Aguilar, A., Fang, M., Roccatagliata, V., et al. 2015, A&A, 580, A82
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405
  * Smith et al. (1974) Smith, E. J., Davis, L., J., Jones, D. E., et al. 1974, Science, 183, 305
  * Stolker et al. (2021) Stolker, T., Haffert, S. Y., Kesseli, A. Y., et al. 2021, AJ, 162, 286
  * Suárez-Andrés et al. (2018) Suárez-Andrés, L., Israelian, G., González Hernández, J. I., et al. 2018, A&A, 614, A84
  * Swastik et al. (2021) Swastik, C., Banyal, R. K., Narang, M., et al. 2021, AJ, 161, 114
  * Szulágyi & Ercolano (2020) Szulágyi, J. & Ercolano, B. 2020, ApJ, 902, 126
  * Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65
  * Tautvaišienė et al. (2022) Tautvaišienė, G., Mikolaitis, Š., Drazdauskas, A., et al. 2022, ApJS, 259, 45
  * Thanathibodee et al. (2019) Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, 885, 94
  * Tremblin et al. (2016) Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19
  * Tremblin et al. (2015) Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17
  * Uyama et al. (2020) Uyama, T., Norris, B., Jovanovic, N., et al. 2020, Journal of Astronomical Telescopes, Instruments, and Systems, 6, 045004
  * Uyama et al. (2017) Uyama, T., Tanigawa, T., Hashimoto, J., et al. 2017, AJ, 154, 90
  * Uyama et al. (2021) Uyama, T., Xie, C., Aoyama, Y., et al. 2021, AJ, 162, 214
  * van Dokkum (2001) van Dokkum, P. G. 2001, PASP, 113, 1420
  * van Holstein et al. (2021) van Holstein, R. G., Stolker, T., Jensen-Clem, R., et al. 2021, A&A, 647, A21
  * Vos et al. (2022) Vos, J. M., Faherty, J. K., Gagné, J., et al. 2022, ApJ, 924, 68
  * Wagner et al. (2018) Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8
  * Wang et al. (2021) Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148
  * Whelan et al. (2015) Whelan, E. T., Huélamo, N., Alcalá, J. M., et al. 2015, A&A, 579, A48
  * Witte et al. (2011) Witte, S., Helling, C., Barman, T., Heidrich, N., & Hauschildt, P. H. 2011, A&A, 529, A44
  * Witte et al. (2009) Witte, S., Helling, C., & Hauschildt, P. H. 2009, A&A, 506, 1367
  * Wolff et al. (2017) Wolff, S. G., Ménard, F., Caceres, C., et al. 2017, AJ, 154, 26
  * Wu et al. (2017a) Wu, Y.-L., Close, L. M., Eisner, J. A., & Sheehan, P. D. 2017a, AJ, 154, 234
  * Wu et al. (2017b) Wu, Y.-L., Sheehan, P. D., Males, J. R., et al. 2017b, ApJ, 836, 223
  * Xie et al. (2020) Xie, C., Haffert, S. Y., de Boer, J., et al. 2020, A&A, 644, A149
  * Zhang et al. (2021) Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021, Nature, 595, 370
  * Zhou et al. (2021) Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244
  * Zhou et al. (2014) Zhou, Y., Herczeg, G. J., Kraus, A. L., Metchev, S., & Cruz, K. L. 2014, ApJ, 783, L17
  * Zhou et al. (2022) Zhou, Y., Sanghi, A., Bowler, B. P., et al. 2022, ApJ, 934, L13
  * Zhu et al. (2016) Zhu, Z., Ju, W., & Stone, J. M. 2016, ApJ, 832, 193
  * Zsidi et al. (2022) Zsidi, G., Manara, C. F., Kóspál, Á., et al. 2022, A&A, 660, A108
  * Zurlo et al. (2020) Zurlo, A., Cugno, G., Montesinos, M., et al. 2020, A&A, 633, A119

##  Appendix A Supplementary material

Figure 11 contains a visual representation of the extraction process for both
objects, along with residuals. Table 2 contains the measured equivalent widths
of the two KI doublets and the $\mathsf{Pa}_{\upbeta}$ line of each individual
cube. Table 3 contains the resulting fitted Gaussian parameters of the
$\mathsf{Pa}_{\upbeta}$ line. The $\mathsf{Pa}_{\upbeta}$ integrated fluxes
used for Fig. 7 were computed from these parameters as
$I_{\textrm{tot}}=\sigma I_{0}\sqrt{2\pi}$. The error bar on the integrated
flux accounts for correlations between parameters:

$\left(\frac{\textrm{err}(I_{\textrm{tot}})}{I_{\textrm{tot}}}\right)^{2}=\left(\frac{\textrm{err}(\sigma)}{\sigma}\right)^{2}+\left(\frac{\textrm{err}(I_{0})}{I_{0}}\right)^{2}+\frac{2\,\textrm{cov}(\sigma,I_{0})}{\sigma
I_{0}}.$ (1)

Figure 11: Visual representation of the last iteration of the extraction for
both GQ Lup b and GSC06214b, on epochs 2017-04-30 and 2017-08-19 respectively.
As mentioned in Sect. 4, the horizontal spider is subtracted for GQ Lup b, but
not for GSC06214b as it does not extend below the companion. Are also marked
the spiders locations for both targets and some persistence effects in the
GSC06214b data. Table 2: Equivalent Widths

Source | Date |  KI 1.169 $\upmu$m |  KI 1.177 $\upmu$m |  KI 1.243 $\upmu$m |  KI 1.253 $\upmu$m | $\mathsf{Pa}_{\upbeta}$
---|---|---|---|---|---|---
|  | Å | Å | Å | Å | Å
GQ Lup b | 2006-09-17 – 23:22:53 | $2.07\pm 0.31$ | $2.81\pm 0.32$ | $1.81\pm 0.30$ | $1.62\pm 0.25$ | $-3.65\pm 0.21$
GQ Lup b | 2006-09-17 – 23:39:49 | $2.42\pm 0.37$ | $2.31\pm 0.36$ | $1.95\pm 0.25$ | $2.04\pm 0.21$ | $-3.34\pm 0.19$
GQ Lup b | 2006-09-17 – 23:45:21 | $2.40\pm 0.40$ | $2.79\pm 0.35$ | $1.22\pm 0.27$ | $1.31\pm 0.24$ | $-3.15\pm 0.20$
GQ Lup b | 2006-09-18 – 00:02:15 | $2.48\pm 0.35$ | $2.46\pm 0.33$ | $1.99\pm 0.28$ | $2.10\pm 0.21$ | $-3.34\pm 0.20$
GQ Lup b | 2006-09-18 – 00:07:41 | $1.94\pm 0.34$ | $2.55\pm 0.34$ | $2.31\pm 0.30$ | $1.65\pm 0.25$ | $-3.00\pm 0.20$
GQ Lup b | 2006-09-18 – 00:24:36 | $2.87\pm 0.51$ | $3.26\pm 0.47$ | $0.97\pm 0.44$ | $1.30\pm 0.39$ | $-2.96\pm 0.28$
GQ Lup b | 2006-09-18 – 00:30:08 | $0.74\pm 0.68$ | $1.80\pm 0.67$ | $1.47\pm 0.55$ | $1.82\pm 0.36$ | $-3.09\pm 0.31$
GQ Lup b | 2006-09-18 – 00:47:01 | $0.78\pm 1.25$ | $3.63\pm 1.22$ | $3.39\pm 0.82$ | $1.97\pm 0.66$ | $-4.40\pm 0.58$
GQ Lup b | 2006-09-18 – 00:52:33 | $5.36\pm 1.78$ | $4.16\pm 1.53$ | $3.12\pm 1.04$ | $2.67\pm 0.99$ | $-3.94\pm 0.73$
GQ Lup b | 2017-04-25 – 05:21:16 | $1.74\pm 0.24$ | $2.18\pm 0.21$ | cosmic | $1.67\pm 0.18$ | $-0.58\pm 0.16$
GQ Lup b | 2017-04-25 – 05:26:40 | cosmic | $3.19\pm 0.18$ | $1.21\pm 0.22$ | $1.57\pm 0.15$ | $-0.77\pm 0.13$
GQ Lup b | 2017-04-25 – 05:37:30 | $1.34\pm 0.20$ | $2.66\pm 0.16$ | $1.70\pm 0.19$ | $1.32\pm 0.17$ | $-0.44\pm 0.13$
GQ Lup b | 2017-04-25 – 05:42:55 | $1.50\pm 0.19$ | $2.72\pm 0.16$ | $1.42\pm 0.16$ | $1.38\pm 0.13$ | $-0.71\pm 0.12$
GQ Lup b | 2017-04-30 – 05:34:35 | $1.94\pm 0.17$ | $2.78\pm 0.15$ | $1.35\pm 0.17$ | $1.29\pm 0.15$ | $-0.44\pm 0.13$
GQ Lup b | 2017-04-30 – 05:56:16 | $1.94\pm 0.20$ | $2.47\pm 0.20$ | $1.77\pm 0.19$ | $1.81\pm 0.17$ | $-0.58\pm 0.15$
GQ Lup b | 2017-04-30 – 05:50:51 | $1.73\pm 0.23$ | $2.39\pm 0.23$ | $2.05\pm 0.23$ | $1.35\pm 0.17$ | $-0.56\pm 0.15$
GQ Lup b | 2017-04-30 – 05:40:00 | $1.59\pm 0.19$ | $2.58\pm 0.15$ | $1.68\pm 0.20$ | $1.47\pm 0.14$ | $-0.43\pm 0.14$
GQ Lup b | 2017-06-01 – 06:29:29 | $2.16\pm 0.36$ | $2.33\pm 0.38$ | $1.42\pm 0.23$ | $0.96\pm 0.23$ | $-0.76\pm 0.22$
GQ Lup b | 2017-06-01 – 06:34:55 | $2.94\pm 0.42$ | $2.25\pm 0.35$ | $1.82\pm 0.28$ | $1.14\pm 0.25$ | $-0.96\pm 0.23$
GQ Lup b | 2017-06-01 – 06:45:44 | $0.61\pm 0.42$ | $0.71\pm 0.43$ | $1.06\pm 0.32$ | $1.35\pm 0.33$ | $-0.82\pm 0.29$
GQ Lup b | 2017-06-01 – 06:51:09 | $1.87\pm 0.54$ | $1.11\pm 0.47$ | $1.44\pm 0.35$ | $1.31\pm 0.29$ | $-1.44\pm 0.24$
GSC06214b | 2017-04-25 – 06:21:07 | $3.16\pm 0.38$ | $5.09\pm 0.35$ | $3.55\pm 0.35$ | $3.58\pm 0.29$ | $-2.79\pm 0.29$
GSC06214b | 2017-04-25 – 06:32:12 | $3.48\pm 0.33$ | $4.50\pm 0.33$ | $3.91\pm 0.31$ | $3.59\pm 0.30$ | $-3.18\pm 0.26$
GSC06214b | 2017-04-25 – 06:48:51 | cosmic | $4.44\pm 0.40$ | $3.10\pm 0.36$ | $2.75\pm 0.33$ | $-3.29\pm 0.27$
GSC06214b | 2017-04-25 – 07:35:50 | $2.37\pm 0.44$ | $3.38\pm 0.39$ | $1.98\pm 0.39$ | $2.92\pm 0.31$ | $-2.15\pm 0.27$
GSC06214b | 2017-04-25 – 07:46:53 | $2.37\pm 0.43$ | $3.47\pm 0.43$ | $1.84\pm 0.43$ | $2.35\pm 0.39$ | $-2.19\pm 0.32$
GSC06214b | 2017-04-25 – 08:03:31 | $2.91\pm 0.48$ | $4.19\pm 0.43$ | $2.61\pm 0.37$ | $2.81\pm 0.37$ | $-1.72\pm 0.32$
GSC06214b | 2017-05-04 – 04:29:49 | $3.39\pm 0.35$ | $5.05\pm 0.34$ | $2.55\pm 0.36$ | $2.64\pm 0.29$ | $-3.38\pm 0.26$
GSC06214b | 2017-05-04 – 04:40:54 | $3.56\pm 0.30$ | $4.37\pm 0.30$ | $2.99\pm 0.36$ | $2.78\pm 0.30$ | $-3.16\pm 0.26$
GSC06214b | 2017-05-04 – 04:57:38 | $4.22\pm 0.37$ | $5.42\pm 0.26$ | $2.80\pm 0.33$ | $3.07\pm 0.26$ | $-3.88\pm 0.21$
GSC06214b | 2017-06-01 – 01:54:27 | $3.21\pm 0.44$ | $4.74\pm 0.47$ | $2.98\pm 0.31$ | $2.79\pm 0.31$ | $-2.73\pm 0.28$
GSC06214b | 2017-06-01 – 02:16:29 | $3.07\pm 0.36$ | $4.62\pm 0.35$ | $2.95\pm 0.32$ | $3.32\pm 0.26$ | $-3.21\pm 0.28$
GSC06214b | 2017-06-01 – 02:27:32 | $3.74\pm 0.36$ | $3.99\pm 0.38$ | $3.09\pm 0.30$ | $3.39\pm 0.31$ | $-2.85\pm 0.28$
GSC06214b | 2017-06-01 – 02:44:12 | $3.81\pm 0.40$ | $4.40\pm 0.32$ | $3.05\pm 0.30$ | $3.69\pm 0.30$ | $-3.60\pm 0.27$
GSC06214b | 2017-08-19 – 01:28:50 | $3.16\pm 0.41$ | $4.76\pm 0.36$ | $2.92\pm 0.37$ | $3.08\pm 0.29$ | $-4.29\pm 0.29$
GSC06214b | 2017-08-19 – 01:39:28 | $3.16\pm 0.43$ | $4.44\pm 0.41$ | $2.92\pm 0.40$ | $3.17\pm 0.32$ | $-4.12\pm 0.27$
GSC06214b | 2017-08-19 – 01:55:28 | $2.62\pm 0.61$ | $5.07\pm 0.53$ | $3.27\pm 0.50$ | $2.70\pm 0.37$ | $-4.14\pm 0.43$
GSC06214b | 2010-07-09 – 08:28:56 | - | - | $3.08\pm 1.05$ | $4.76\pm 0.92$ | $-13.28\pm 1.08$
GSC06214b | 2010-07-09 – 08:34:42 | - | - | $1.18\pm 0.98$ | $2.30\pm 0.69$ | $-11.05\pm 0.72$
GSC06214b | 2010-07-09 – 08:45:15 | - | - | $3.90\pm 1.10$ | $3.92\pm 0.74$ | $-10.36\pm 0.79$
GSC06214b | 2010-07-09 – 08:51:02 | - | - | $2.23\pm 1.20$ | $2.71\pm 0.85$ | $-9.13\pm 0.78$
GSC06214b | 2010-07-09 – 08:59:40 | - | - | $2.14\pm 0.96$ | $4.04\pm 0.91$ | $-10.25\pm 0.79$
GSC06214b | 2010-07-09 – 09:05:26 | - | - | $3.09\pm 1.17$ | $3.69\pm 0.84$ | $-13.89\pm 0.90$
GSC06214b | 2010-07-09 – 09:12:50 | - | - | $3.77\pm 0.80$ | $1.12\pm 0.82$ | $-9.37\pm 0.76$
GSC06214b | 2010-07-09 – 09:18:36 | - | - | $6.39\pm 1.21$ | $4.93\pm 1.01$ | $-14.88\pm 1.17$

333Measured equivalent widths of the various $\mathrm{K\,\textsc{i}}$ and
$\mathsf{Pa}_{\upbeta}$ lines, following the method described in Sembach &
Savage (1992) and Cushing et al. (2005). Date is given in YYYY-MM-DD –
hh:mm:ss format. Possible correlations between the lines are shown in Fig. 10.

Table 3: Gaussian fitting results

Source | Date |  I0 | $\upsigma$ | $\upmu_{0}$ | Int. Flux | Intrinsic FWHM
---|---|---|---|---|---|---
| YYYY-MM-DD – hh:mm:ss |  $10^{-16}$ W/m2/$\upmu$m | km/s | km/s |  $10^{-19}$ W/m2 | km/s
GQ Lup b | 2006-09-17 – 23:22:53 | $25.0\pm 1.3$ | $63.5\pm 3.7$ | $-27.8\pm 3.7$ | $17.0\pm 0.9$ | $79.1\pm 16.5$
GQ Lup b | 2006-09-17 – 23:39:49 | $19.4\pm 1.2$ | $79.5\pm 6.0$ | $-19.0\pm 5.9$ | $16.5\pm 1.1$ | $137.5\pm 19.4$
GQ Lup b | 2006-09-17 – 23:45:21 | $20.5\pm 1.3$ | $72.1\pm 5.3$ | $-20.4\pm 5.3$ | $15.8\pm 1.0$ | $112.6\pm 18.9$
GQ Lup b | 2006-09-18 – 00:02:15 | $22.4\pm 0.5$ | $67.9\pm 1.8$ | $-22.5\pm 1.8$ | $16.3\pm 0.4$ | $97.0\pm 7.1$
GQ Lup b | 2006-09-18 – 00:07:41 | $19.0\pm 0.9$ | $75.9\pm 4.1$ | $-16.8\pm 4.0$ | $15.4\pm 0.7$ | $125.6\pm 13.7$
GQ Lup b | 2006-09-18 – 00:24:36 | $25.3\pm 2.2$ | $54.5\pm 5.4$ | $-17.1\pm 5.4$ | $14.8\pm 1.3$ | $18.8\pm 87.2$
GQ Lup b | 2006-09-18 – 00:30:08 | $22.5\pm 1.5$ | $65.4\pm 4.9$ | $-21.2\pm 4.9$ | $15.8\pm 1.0$ | $87.2\pm 20.5$
GQ Lup b | 2006-09-18 – 00:47:01 | $29.5\pm 3.3$ | $61.5\pm 7.9$ | $-1.3\pm 7.9$ | $19.4\pm 2.2$ | $69.3\pm 39.0$
GQ Lup b | 2006-09-18 – 00:52:33 | $45.5\pm 6.9$ | $30.9\pm 5.5$ | $-25.6\pm 5.4$ | $15.1\pm 2.3$ | $0.0\pm 58.6$
GQ Lup b | 2017-04-25 – 05:21:16 | $5.7\pm 1.3$ | $38.3\pm 9.9$ | $-19.2\pm 9.9$ | $2.3\pm 0.5$ | $0.0\pm 80.3$
GQ Lup b | 2017-04-25 – 05:26:40 | $4.4\pm 0.9$ | $63.8\pm 15.0$ | $-53.8\pm 14.8$ | $3.0\pm 0.6$ | $80.4\pm 66.3$
GQ Lup b | 2017-04-25 – 05:37:30 | $4.7\pm 1.3$ | $46.9\pm 15.1$ | $-41.4\pm 15.1$ | $2.4\pm 0.7$ | $0.0\pm 101.4$
GQ Lup b | 2017-04-25 – 05:42:55 | $4.9\pm 0.8$ | $43.5\pm 7.9$ | $-2.6\pm 7.9$ | $2.3\pm 0.4$ | $0.0\pm 71.1$
GQ Lup b | 2017-04-30 – 05:34:35 | $4.8\pm 1.3$ | $47.0\pm 14.2$ | $-49.0\pm 14.2$ | $2.4\pm 0.6$ | $0.0\pm 98.2$
GQ Lup b | 2017-04-30 – 05:56:16 | $3.9\pm 1.1$ | $45.0\pm 14.3$ | $-34.3\pm 14.3$ | $1.9\pm 0.5$ | $0.0\pm 98.6$
GQ Lup b | 2017-04-30 – 05:50:51 | $5.5\pm 1.5$ | $40.4\pm 12.6$ | $-37.1\pm 12.6$ | $2.4\pm 0.6$ | $0.0\pm 91.7$
GQ Lup b | 2017-04-30 – 05:40:00 | $4.1\pm 1.4$ | $35.7\pm 13.9$ | $-22.4\pm 13.9$ | $1.6\pm 0.5$ | $0.0\pm 96.8$
GQ Lup b | 2017-06-01 – 06:29:29 | $9.5\pm 2.4$ | $49.3\pm 14.5$ | $-55.9\pm 14.5$ | $5.0\pm 1.3$ | $0.0\pm 99.1$
GQ Lup b | 2017-06-01 – 06:34:55 | $7.4\pm 1.9$ | $52.2\pm 15.7$ | $-14.2\pm 15.7$ | $4.1\pm 1.1$ | $0.0\pm 103.6$
GQ Lup b | 2017-06-01 – 06:45:44 | $11.2\pm 2.6$ | $56.6\pm 15.1$ | $-54.0\pm 15.0$ | $6.8\pm 1.6$ | $40.1\pm 118.4$
GQ Lup b | 2017-06-01 – 06:51:09 | $11.6\pm 2.3$ | $48.2\pm 11.0$ | $-69.2\pm 10.9$ | $6.0\pm 1.2$ | $0.0\pm 85.1$
GSC06214b | 2017-04-25 – 06:21:07 | $3.2\pm 0.3$ | $54.7\pm 6.7$ | $-14.2\pm 6.7$ | $1.9\pm 0.2$ | $0.0\pm 74.3$
GSC06214b | 2017-04-25 – 06:32:12 | $4.0\pm 0.2$ | $62.4\pm 2.8$ | $-7.3\pm 2.7$ | $2.7\pm 0.1$ | $0.0\pm 47.2$
GSC06214b | 2017-04-25 – 06:48:51 | $3.8\pm 0.3$ | $65.6\pm 5.5$ | $-10.1\pm 5.3$ | $2.7\pm 0.2$ | $0.0\pm 66.9$
GSC06214b | 2017-04-25 – 07:35:50 | $3.2\pm 0.3$ | $62.5\pm 7.9$ | $-9.7\pm 7.7$ | $2.1\pm 0.2$ | $0.0\pm 80.8$
GSC06214b | 2017-04-25 – 07:46:53 | $3.2\pm 0.5$ | $63.7\pm 11.9$ | $-19.8\pm 11.6$ | $2.2\pm 0.3$ | $0.0\pm 100.8$
GSC06214b | 2017-04-25 – 08:03:31 | $2.5\pm 0.2$ | $57.5\pm 4.5$ | $-11.5\pm 4.4$ | $1.5\pm 0.1$ | $0.0\pm 60.2$
GSC06214b | 2017-05-04 – 04:29:49 | $3.6\pm 0.2$ | $74.2\pm 3.8$ | $-2.4\pm 3.6$ | $2.9\pm 0.1$ | $52.9\pm 29.9$
GSC06214b | 2017-05-04 – 04:40:54 | $3.7\pm 0.4$ | $78.6\pm 10.1$ | $7.7\pm 9.1$ | $3.1\pm 0.3$ | $81.0\pm 54.2$
GSC06214b | 2017-05-04 – 04:57:38 | $4.5\pm 0.2$ | $75.7\pm 3.9$ | $8.2\pm 3.6$ | $3.7\pm 0.2$ | $63.5\pm 25.7$
GSC06214b | 2017-06-01 – 01:54:27 | $4.4\pm 0.5$ | $66.4\pm 8.3$ | $-5.7\pm 8.1$ | $3.1\pm 0.3$ | $0.0\pm 83.3$
GSC06214b | 2017-06-01 – 02:16:29 | $4.3\pm 0.3$ | $64.9\pm 5.1$ | $-0.3\pm 5.0$ | $3.0\pm 0.2$ | $0.0\pm 64.7$
GSC06214b | 2017-06-01 – 02:27:32 | $4.1\pm 0.2$ | $72.7\pm 4.4$ | $5.8\pm 4.1$ | $3.2\pm 0.2$ | $40.0\pm 44.1$
GSC06214b | 2017-06-01 – 02:44:12 | $4.8\pm 0.3$ | $78.5\pm 7.1$ | $13.4\pm 6.5$ | $4.0\pm 0.3$ | $80.0\pm 38.9$
GSC06214b | 2017-08-19 – 01:28:50 | $6.3\pm 0.2$ | $73.5\pm 3.2$ | $-5.9\pm 3.0$ | $4.9\pm 0.2$ | $46.9\pm 28.0$
GSC06214b | 2017-08-19 – 01:39:28 | $5.8\pm 0.2$ | $75.3\pm 4.0$ | $0.3\pm 3.7$ | $4.7\pm 0.2$ | $60.8\pm 27.6$
GSC06214b | 2017-08-19 – 01:55:28 | $5.5\pm 0.4$ | $83.5\pm 8.0$ | $5.7\pm 7.1$ | $4.9\pm 0.4$ | $104.4\pm 35.6$
GSC06214b | 2010-07-09 – 08:28:56 | $17.4\pm 0.9$ | $61.1\pm 3.7$ | $7.5\pm 3.7$ | $11.4\pm 0.6$ | $120.2\pm 10.4$
GSC06214b | 2010-07-09 – 08:34:42 | $19.6\pm 1.0$ | $44.5\pm 2.7$ | $7.2\pm 2.7$ | $9.4\pm 0.5$ | $69.1\pm 9.7$
GSC06214b | 2010-07-09 – 08:45:15 | $18.1\pm 1.3$ | $52.6\pm 4.5$ | $8.1\pm 4.5$ | $10.2\pm 0.8$ | $95.4\pm 13.8$
GSC06214b | 2010-07-09 – 08:51:02 | $17.2\pm 0.6$ | $42.0\pm 1.8$ | $8.5\pm 1.8$ | $7.8\pm 0.3$ | $59.8\pm 6.9$
GSC06214b | 2010-07-09 – 08:59:40 | $16.9\pm 0.8$ | $56.6\pm 3.1$ | $5.3\pm 3.1$ | $10.2\pm 0.5$ | $107.4\pm 9.1$
GSC06214b | 2010-07-09 – 09:05:26 | $18.1\pm 1.4$ | $61.1\pm 5.6$ | $12.0\pm 5.6$ | $11.8\pm 0.9$ | $120.2\pm 15.7$
GSC06214b | 2010-07-09 – 09:12:50 | $16.0\pm 0.6$ | $55.2\pm 2.5$ | $9.8\pm 2.5$ | $9.4\pm 0.4$ | $103.2\pm 7.5$
GSC06214b | 2010-07-09 – 09:18:36 | $19.6\pm 1.1$ | $63.7\pm 4.2$ | $17.9\pm 4.1$ | $13.4\pm 0.8$ | $127.7\pm 11.6$

444Gaussian fitting results of the $\mathsf{Pa}_{\upbeta}$ lines. Date is
given in YYYY-MM-DD – hh:mm:ss format. I0 is the scaling parameter of the
gaussian.