# Characterizing the Protolunar Disk of the Accreting Companion GQ Lupi B

Tomas Stolker, Sebastiaan Y. Haffert, Aurora Y. Kesseli, Rob G. van Holstein, Yuhiko Aoyama, Jarle Brinchmann, Gabriele Cugno, Julien H. Girard, Gabriel-Dominique Marleau, Michael R. Meyer, et al.

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

► To cite this version:

Tomas Stolker, Sebastiaan Y. Haffert, Aurora Y. Kesseli, Rob G. van Holstein, Yuhiko Aoyama, et al.. Characterizing the Protolunar Disk of the Accreting Companion GQ Lupi B. The Astronomical Journal, 2021, 162, 10.3847/1538-3881/ac2c7f. insu-03705258

HAL Id: insu-03705258

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

Submitted on 2 Aug 2025

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

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

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

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

# Characterizing the Protolunar Disk of the Accreting Companion GQ Lupi B*

Tomas Stolker¹, Sebastiaan Y. Haffert²,<fnref n="16" />, Aurora Y. Kesseli¹, Rob G. van Holstein¹,³, Yuhiko Aoyama⁴,⁵,⁶, Jarle Brinchmann¹,⁷, Gabriele Cugno⁸, Julien H. Girard⁹, Gabriel-Dominique Marleau¹⁰,¹¹,¹², Michael R. Meyer¹³, Julien Milli¹⁴, Sascha P. Quanz⁸, Ignas A. G. Snellen¹, and Kamen O. Todorov¹⁵

¹ Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands; stolker@strw.leidenuniv.nl

² Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

³ European Southern Observatory, Alonso de Córdova 3107, Casilla 19001, Vitacura, Santiago, Chile

⁴ Department of Earth and Planetary Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

⁵ Institute for Advanced Study, Tsinghua University, Beijing 100084, People's Republic of China

⁶ Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China

⁷ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal

⁸ Institute for Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland

⁹ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁰ Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany

¹¹ Physikalisches Institut, Universität Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

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

¹³ Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109-1107, USA

¹⁴ Université Grenoble Alpes, IPAG, F-38000 Grenoble, France

¹⁵ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090 GE Amsterdam, The Netherlands
Received 2021 August 6; revised 2021 September 20; accepted 2021 October 1; published 2021 December 6

## Abstract

GQ Lup B is a young and accreting, substellar companion that appears to drive a spiral arm in the circumstellar disk of its host star. We report high-contrast imaging observations of GQ Lup B with VLT/NACO at 4–5 µm and medium-resolution integral field spectroscopy with VLT/MUSE. The optical spectrum is consistent with an M9 spectral type, shows characteristics of a low-gravity atmosphere, and exhibits strong Hα emission. The H–M' color is ≳1 mag redder than field dwarfs with similar spectral types, and a detailed analysis of the spectral energy distribution (SED) from optical to mid-infrared wavelengths reveals excess emission in the L', NB4.05, and M' bands. The excess flux is well described by a blackbody component with T<sub>disk</sub> ≈ 460 K and R<sub>disk</sub> ≈ 65 R<sub>J</sub> and is expected to trace continuum emission from small grains in a protolunar disk. We derive an extinction of A<sub>V</sub> ≈ 2.3 mag from the broadband SED with a suspected origin in the vicinity of the companion. We also combine 15 yr of astrometric measurements and constrain the mutual inclination with the circumstellar disk to 84 ± 9 deg, indicating a tumultuous dynamical evolution or a stellar-like formation pathway. From the measured Hα flux and the estimated companion mass, M<sub>p</sub> ≈ 30 M<sub>J</sub>, we derive an accretion rate of M̄ ≈ 10<sup>-6.5</sup> M<sub>J</sub> yr<sup>-1</sup>. We speculate that the disk is in a transitional stage in which the assembly of satellites from a pebble reservoir has opened a central cavity while GQ Lup B is in the final stages of its formation.

*Unified Astronomy Thesaurus concepts: Accretion (14); Brown dwarfs (185); Direct imaging (387); Planet formation (1241); Natural satellite formation (1425); High angular resolution (2167)*

## 1. Introduction

Directly imaged planets and brown dwarfs are an intriguing population of low-mass objects that have been discovered at large distances from their star (e.g., Oppenheimer & Hinkley 2009; Bowler 2016). Their super-Jupiter masses and wide orbits challenge planet formation theories due to the long timescales associated with a bottom-up formation of rocky cores and subsequent gas accretion in circumstellar disks (CSDs). Some of these objects may instead have formed through a stellar-like formation pathway—in particular those with a high companion-to-star mass ratio and/or on wide orbits.

Large-scale surveys have revealed that planets are typically found at smaller separations around intermediate-mass stars, whereas brown dwarfs are detected on larger orbits around lower-mass stars (e.g., Nielsen et al. 2019; Vigan et al. 2021). This may suggest a stellar formation scenario for substellar-mass objects on wide orbits, which is in line with eccentricity constraints inferred

from their orbital dynamics (Bowler et al. 2020). The chemical composition of their atmospheres provides another glimpse into their formation history (e.g., Öberg et al. 2011; Madhusudhan et al. 2014). Atmospheric retrievals are starting to reveal signatures of the physical and chemical processes that are at play during formation, for example, through the spectral inference of the C/O ratio and metallicity, which could point to nonsolar chemical compositions (e.g., Lee et al. 2013; Mollière et al. 2020).

The youngest (i.e., ≲10 Myr) of the directly imaged, planetary-mass companions (PMCs) provide a direct window to the formation of giant planets and brown dwarfs and their circumplanetary/substellar disk characteristics. The PDS 70 planets have become the signpost for empirical studies on embedded protoplanets (e.g., Keppler et al. 2018; Müller et al. 2018). These Jupiter-mass planets have opened a gap at 20–30 au in the CSD while accreting gas and dust from their environment (e.g., Haffert et al. 2019; Stolker et al. 2020a; Wang et al. 2020). Other accreting PMCs (see Table 2 in Wu et al. 2017a, for an overview) typically orbit further away from their star (≈100 au), and the connection with the CSD is less clear (e.g., Ireland et al. 2011; Kraus et al. 2014). The spectral

<footnote>* Based on observations collected under ESO programmes 0101.C-0502(B), 0102.C-0649(A), and 0103.C-0524(A).</footnote>

<footnote n="16">¹⁶ NASA Hubble Fellow.</footnote>
energy distributions (SEDs) of some of these objects are unusually red, which can be evidence of a dusty disk (e.g., Bailey et al. 2013; Wu et al. 2015)—in line with expectations from isolated, planetary, and substellar objects—that is expected to serve as the formation site of satellites. The study of the accretion and disk characteristics of both populations (i.e., embedded in versus detached from the CSD) may reveal clues about possible differences in their formation pathways and the processes by which giant planets and brown dwarfs accumulate a gaseous envelope.

One such directly imaged substellar object is GQ Lup B, which was discovered by Neuhäuser et al. (2005) orbiting a K7Ve-type T Tauri star (Herbig 1977) with a well-studied CSD (e.g., Dai et al. 2010; McClure et al. 2012). The mass of GQ Lup B has been debated ever since, with inferred masses of ~10–40 $M_J$ (e.g., Marois et al. 2007; McElwain et al. 2007; Seifahrt et al. 2007). Its near-infrared (NIR) spectrum has been classified as mid M to early L (McElwain et al. 2007) with $T_{\text{eff}} \approx 2500$ K and $\log g \approx 4$ (see overview in Table 1 by Lavigne et al. 2009). GQ Lup B was also detected with optical photometry. This revealed $H\alpha$ emission that has been linked with accretion (Marois et al. 2007; Zhou et al. 2014; Wu et al. 2017b), in line with the detection of $Pa\beta$ emission by Seifahrt et al. (2007). In addition to the atmosphere, the orbit has been analyzed by several authors (Janson et al. 2006; Neuhäuser et al. 2008; Ginski et al. 2014) with the most recent constraints pointing to an inclination of ~60 deg and semimajor axis of 100–150 au (Schwarz et al. 2016). The companion has a low projected spin velocity, which could suggest that it is still gaining angular momentum (Schwarz et al. 2016) at the age of 2–5 Myr (Donati et al. 2012). Recently, a second companion, GQ Lup C, was detected at a projected separation of ~2400 au (Alcalá et al. 2020). This low-mass star is also surrounded by a dusty accretion disk, as inferred from UV and IR excess, with a somewhat similar orientation to the disk of GQ Lup (Lazzoni et al. 2020). Apart from GQ Lup B and C, who are orbiting at a large separation, there is also a gap detected in the CSD at ~10 au, which could be evidence of a hidden planet on a solar-system scale (Long et al. 2020).

In this work, we report the first detection of GQ Lup B with 4–5 µm imaging and with optical spectroscopy. Mid-infrared (MIR) wavelengths are sensitive to effects of clouds and carbon chemistry but are also a powerful probe for thermal emission from circumplanetary/substellar material. With the optical spectrum, we search for emission lines that carry accretion signatures. We analyze the 4–5 µm photometry and optical spectroscopy in combination with archival, medium-resolution JHK-band spectra to detect and characterize excess emission from the disk around GQ Lup B and constrain the atmospheric and extinction properties.

## 2. Observations

### 2.1. Mid-infrared Imaging with VLT/NACO

We observed GQ Lup with the Very Large Telescope (VLT) on Cerro Paranal in Chile as part of the Mid-InfraRed Atmospheric Characterization of Long-period Exoplanets and Substellar companions (MIRACLES) survey (ESO program ID: 0102.C-0649(A)). The program aims at the systematic characterization of directly imaged planets and brown dwarfs at 4–5 µm (Stolker et al. 2020b) by using the adaptive optics-assisted imaging capabilities of NACO (Lenzen et al. 2003;

Rousset et al. 2003). The observations were carried out with the NB4.05 (Brα; $\lambda_0 = 4.05$ µm, $\Delta\lambda = 0.06$ µm$^{17}$) and M' ($\lambda_0 = 4.78$ µm, $\Delta\lambda = 0.59$ µm) filters on the nights of 2019 March 14 and 8, respectively. A description of the observing strategy can be found in Stolker et al. (2020b), but we will provide a few specific details here.

With the NB4.05 filter, we used a detector integration time (DIT) of 1 s and obtained a total of 1400 frames. In addition, we obtained 600 frames with a shorter DIT of 0.2 s for calibration purposes. The parallactic rotation was only 11.6 deg but sufficient for a detection of GQ Lup B. The atmospheric conditions during the observations were stable with a seeing of 0.′′78 ± 0.′′04 and a coherence time of 5.0 ± 0.8 ms. The angular resolution, as measured from the point-spread function (PSF), was 116 mas (1 FWHM), and the stellar flux in the unsaturated exposures varied by 2.8%.

Similarly, we observed GQ Lup with the M' filter but with a DIT of 50 ms to prevent saturation by the high thermal background emission. We obtained a total of 34,400 frames with a comparable seeing, 0.′′85 ± 0.′′04, and coherence time, 4.3 ± 0.3 ms as the NB4.05 observations. The total parallactic rotation was 15.6 deg, and the angular resolution was 137 mas. The stellar flux varied by 3.9% across the full stack of frames.

### 2.2. Integral Field Spectroscopy (IFS) with VLT/MUSE

GQ Lup was also observed with the MUSE integral field spectrograph on the night of 2019 April 19 (ESO program ID: 0103.C-0524(A)). MUSE operates in the visible (4800–9300 Å) with a resolving power of $R = 1740–3450$ (Bacon et al. 2010). The instrument is mounted on Unit Telescope 4 (UT4) of the VLT and therefore benefits from the adaptive optics facility (AOF; Arsenault et al. 2008; Ströbele et al. 2012). We used the narrow-field mode, which leverages the laser tomography adaptive optics system (GALACSI) of the AOF and samples the field with $(25 \text{ mas})^2 \text{ spaxel}^{-1}$ to enable high-angular resolution observations.

The observations were obtained with field tracking and comprised 12 exposures of GQ Lup with a DIT of 139 s. For each subsequent exposure, the derotator offset was increased by 90 deg to reduce detector artifacts and improve the sampling of the line-spread function (LSF). The atmospheric conditions were excellent with a seeing of 0.′′32 ± 0.′′04 and a coherence time of 16 ± 2 ms, but the target was observed at a somewhat high airmass of 1.18–1.40. Earlier in the night, a single exposure with a DIT of 240 s was obtained of the standard star GD 108, which is a B-type subdwarf (Greenstein 1969). The conditions were a bit poorer but still good with a seeing of 0.′′64 and a coherence time of 9 ms. The calibration target was however observed at a smaller airmass of 1.07.

## 3. Data Reduction

### 3.1. NACO

#### 3.1.1. Image Processing

The processing and calibration of the NACO data sets were done with version 0.8.3 of PynPoint$^{18}$ (Stolker et al. 2019), which applies an implementation of full-frame principal component analysis (PCA; Amara & Quanz 2012;

<footnote>$^{17}$ The FWHM of the NB4.05 filter has been incorrectly listed on the ESO website as 0.02 µm, so we revised the value to 0.06 µm.</footnote>

<footnote>$^{18}$ https://pynpoint.readthedocs.io</footnote>
<table><caption>**Table 1**
Photometry of GQ Lup B</caption><thead><tr><th>Filter</th><th>Contrast (mag)</th><th>GQ Lup (mag)</th><th>Apparent Magnitude (mag)</th><th>Absolute Magnitude (mag)</th><th>Flux (W m<sup>-2</sup> µm<sup>-1</sup>)</th></tr></thead><tbody><tr><td>NACO NB4.05</td><td>6.54 ± 0.05</td><td>5.75 ± 0.04</td><td>12.29 ± 0.06</td><td>6.36 ± 0.06</td><td>(4.80 ± 0.28) × 10<sup>-16</sup></td></tr><tr><td>NACO M'</td><td>6.50 ± 0.05</td><td>5.47 ± 0.06</td><td>11.97 ± 0.08</td><td>6.03 ± 0.08</td><td>(3.52 ± 0.25) × 10<sup>-16</sup></td></tr></tbody></table>

<table><caption>**Table 2**
Astrometry of GQ Lup B</caption><thead><tr><th>MJD</th><th>Instrument</th><th>Filter</th><th>Separation (mas)</th><th>Position Angle (deg)</th></tr></thead><tbody><tr><td>58345.1</td><td>SPHERE</td><td>B_H</td><td>711.6 ± 2.4</td><td>278.27 ± 0.24</td></tr><tr><td>58551.3</td><td>NACO</td><td>M'</td><td>720.3 ± 3.4</td><td>277.90 ± 0.18</td></tr><tr><td>58557.2</td><td>NACO</td><td>NB4.05</td><td>715.0 ± 2.8</td><td>277.83 ± 0.17</td></tr><tr><td>58593.1</td><td>MUSE</td><td>...</td><td>701 ± 20</td><td>278.3 ± 1.2</td></tr></tbody></table>

**Note.** The true north, -1.75 ± 0.08 deg, and pupil offset, -135.99 ± 0.11 deg, for SPHERE/IRDIS have been adopted from Maire et al. (2016), and the true north for NACO, -0.44 ± 0.10 deg, has been adopted from Cheetham et al. (2019). The listed position angles have been corrected for these offsets, and the uncertainties have been propagated into the error budget.

Soummer et al. (2012) to remove the stellar PSF. We followed mostly the procedure as outlined in Stolker et al. (2020b) but provide a few details here.

While a background subtraction based on the mean of the adjacent data cubes was sufficient to suppress the bright background flux, we applied an additional PCA-based background subtraction to remove striped detector artifacts in several of the frames. We followed the approach by Hunziker et al. (2018) and subtracted 3 (NB4.05) and 2 (M') principal components (PCs) of the background emission from the data. This lowered sufficiently the quasi-static detector noise on visual inspection. Before calibrating the data, we collapsed subsets of 2 (NB4.05) and 67 (M') of the preprocessed frames to end up with a stack of 687 (NB4.05) and 499 (M') frames.

### 3.1.2. Relative Calibration

The flux and position of GQ Lup B were calibrated relative to GQ Lup by injecting negative copies of the unsaturated PSF to minimize the flux at the companion position (see Stolker et al. 2020b, for details). We first retrieved the parameters as a function of the PCs to inspect the dependence. Next, we fixed the number of subtracted PCs to 6 (NB4.05) and 10 (M') and used the Markov chain Monte Carlo (MCMC) ensemble sampler *emcee* (Foreman-Mackey et al. 2013) to determine the best-fit values and statistical errors (assuming Gaussian noise). Finally, we estimated the potential bias and systematic errors (e.g., due to speckle and background residuals) by injecting and retrieving artificial sources.

The photometric and astrometric results are listed in Tables 1 and 2. The separations and positions angles are consistent with each other within the 1–2σ errors. The photometric error budget of NB4.05 also includes a calibration error of 0.03 mag, which is calculated from the unsaturated NB4.05 exposures while the stellar M' flux remained unsaturated in all frames.

### 3.1.3. Absolute Flux Calibration

GQ Lup has a CSD so its brightness at 4–5 µm is affected by IR excess. Therefore, to convert the contrast values of GQ Lup B into apparent magnitudes, we calculated the flux of GQ Lup in the NB4.05 and M' filters by simply fitting the 2MASS *H* and *Ks* and WISE W1 and W2 fluxes with a power-law function in log–log space, that is, $\log f_{\lambda} = a + bx^c$, where $x = \log \lambda$ and with $f_{\lambda}$ in W m⁻² µm⁻¹ and $\lambda$ in µm. We used the Bayesian framework (see Section 4.5 for details) in the species<fnref n="19" /> toolkit (Stolker et al. 2020b) to determine the best-fit parameters and uncertainties: $a = -11.80 \pm 0.07$, $b = -1.81 \pm 0.06$, and $c = 1.38 \pm 0.17$.

The posterior distributions were then propagated into synthetic magnitudes by folding the power-law spectra with the filter profiles. From this, we derived $5.75 \pm 0.04$ mag in NB4.05 and $5.47 \pm 0.06$ mag in M' for GQ Lup, which then gave the apparent magnitudes of the companion. The absolute magnitudes were calculated by adopting the Gaia distance of 154.1 ± 0.7 pc (Gaia Collaboration et al. 2016, 2021). Finally, the magnitudes were converted into fluxes with species by using a flux-calibrated spectrum of Vega (Bohlin 2007) and setting its magnitude to 0.03 mag for all filters. The calibrated photometry is provided in Table 1.

## 3.2. MUSE

The reduction and calibration of the MUSE data were done in a similar way as outlined by Haffert et al. (2020). We used version 2.8.2 of the MUSE pipeline recipes for EsoRex (Weilbacher et al. 2020). The pipeline performs a background subtraction, flat field correction, 2D spectral extraction, wavelength calibration (by using the internal arc lamps), and a telluric correction. The PSF of MUSE is not Nyquist sampled so we centered the images only with pixel precision (i.e., 25 mas) to avoid inaccuracies due to interpolation. The final product is a cleaned and flux-calibrated (by using the spectrum of the standard star) data cube, which was used for the spectral extraction of GQ Lup B.

Some spurious features were visible in the spectrum due to an imperfect correction of the telluric transmission, which may have been caused by the difference in airmass (~0.1–0.3) with the standard star. The same residuals were also seen in the spectrum of GQ Lup so we applied a second correction by normalizing the MUSE spectrum with a PHOENIX model spectrum (Husser et al. 2013), for which we adopted $T_{eff} = 4300$ K and $log g = 3.7$ from Donati et al. (2012), while masking and interpolating regions with emission lines. The derived scaling correction was then applied to the spectrum of GQ Lup B, which effectively removed the remaining telluric residuals.

After image registration and calibration, we removed the stellar halo by subtracting an azimuthally averaged radial profile in steps of 12.5 mas. Additionally, a second order 2D polynomial was fitted around GQ Lup B to subtract remaining, low-frequency

<footnote n="19">¹⁹ https://species.readthedocs.io</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Detection of GQ Lup B with the NACO NB4.05 (left panel) and *M'* (central panel) filters, and with MUSE at Hα (right panel). The images show the mean-combined residuals of the PSF subtraction on a linear scale. The insets display the unsaturated PSF fluxes of GQ Lup on a logarithmic scale. The color scales have been normalized for best visibility. North is up and east is left in all images.</figcaption></figure>

structures in the images. The spectrum of GQ Lup B was then extracted by fitting a shifted and normalized copy of the stellar PSF at the position of the companion. For each wavelength, the flux was estimated as the weighted (based on the PSF template) sum of the pixels in an aperture with a diameter of 150 mas, therefore also accounting for the flux contributions outside the aperture that have a low signal-to-noise ratio (S/N). The noise was estimated from the scatter of the extracted fluxes between 6300 and 7000 Å by sampling random errors. In that regime, the continuum of GQ Lup B lies below the noise floor, and there are limited residuals of GQ Lup.

### 3.3. Archival VLT/SPHERE H-band Imaging

We also reanalyzed the archival SPHERE (Beuzit et al. 2019) H-band imaging data set (ESO program ID: 0101.C-0502(B)) that was published by van Holstein et al. (2021), in order to extract the position of the companion from this 2018 epoch (see orbit fit in Section 4.6). The data were obtained with the dual-beam polarimetric imaging mode (de Boer et al. 2020; van Holstein et al. 2020) of the IRDIS camera (Dohlen et al. 2008) but we only used the Stokes *I* images for the astrometry. These polarimetric observations were done in pupil-tracking mode (see van Holstein et al. 2017), resulting in a field rotation of only 6 deg, but the companion is sufficiently bright to be robustly detected at a separation of ≈0.7".

We reduced the data with IRDAP<fnref n="20" /> (van Holstein et al. 2017, 2020) and measured the relative position with Pyn-Point (Stolker et al. 2019), similar to van Holstein et al. (2021) and the calibration in Section 3.1.2. Given the limited field rotation, we subtracted the PSF with 2 PCs and used artificial sources to obtain the best-fit position and error bars. The centering of the frames was achieved with the satellite spots in the dedicated calibration frames to locate the star behind the coronagraph. There was a ~0.2 pixel difference in the position of the star between the frames obtained at the start and end, which we adopted in the astrometric error budget as the centering precision.

The results from the astrometry are listed in Table 2. The errors are dominated by the accuracy of the stellar position and

the true north correction. Additionally, we cross-checked the astrometric extraction by fitting a 2D Gaussian directly to the companion flux in the preprocessed frames. In that case, we obtained a very similar result with differences of 0.4 mas in separation and 0.01 deg in position angle, which is small (~10%) compared to the total errors.

## 4. Results

### 4.1. Mid-infrared and Optical Detection of GQ Lup B

The mean-combined residuals from the PSF subtraction with PCA of the NACO data sets are shown in the left and central panels of Figure 1. The PSF model was created by projecting each frame of the preprocessed data onto the first 3 PCs (for both filters). GQ Lup B is relatively bright (i.e., ~12 mag) at 4–5 µm so the object is clearly detected west of the central star both in NB4.05 and *M'*. At a separation of 715–720 mas and a position angle of ~278 deg (see Table 2), the separation has decreased by ~5 mas, and the position angle increased by ~1 deg with respect to the last epoch by Ginski et al. (2014) from 2012. The orbital architecture of GQ Lup B will be analyzed in more detail in Section 4.6.

The right panel of Figure 1 shows the MUSE image at Hα after PSF subtraction. Some high-frequency stellar residuals remained present but these noise features impact mainly the blue end (λ < 5500 Å) of the spectrum, while beyond 6300 Å the residuals are relatively small at the separation of GQ Lup B. The extracted spectrum will be presented and analyzed in Section 4.3.

### 4.2. Color and Magnitude Comparison

Color–magnitude diagrams (CMDs) provide an insightful approach for comparing the photometric characteristics of substellar and planetary-mass objects and reveal correlations related to temperature and surface gravity. We show in Figure 2 the absolute *H* magnitude as function of the *H* − *M'* color. This CMD was created with the *species* toolkit (see also Stolker et al. 2020b) and includes a sample of regular field dwarfs and young/low-gravity objects (Dupuy & Liu 2012; Dupuy & Kraus 2013; Liu et al. 2016), directly imaged planets and brown dwarfs (Biller et al. 2010, 2012; Galicher et al. 2011;

<footnote n="20">²⁰ https://irdap.readthedocs.io</footnote>
Ireland et al. 2011; Currie et al. 2012, 2013; Bailey et al. 2013; Bonnefoy et al. 2014; Lacour et al. 2016; Chauvin et al. 2017; Daemgen et al. 2017; Garcia et al. 2017; Milli et al. 2017; Rajan et al. 2017; Stolker et al. 2019, 2020a, 2020b), and synthetic photometry calculated from the AMES-Cond/Dusty isochrones at 3 Myr and model spectra (Chabrier et al. 2000; Allard et al. 2001; Baraffe et al. 2003), that is, hot-start cooling tracks with a clear or cloudy atmosphere. We note that narrowband SPHERE H2 magnitudes are shown for HIP 65426b and HD 206893 B, so ignoring minor color effects. The very red H–M' colors of these two objects are probably caused by an enhanced dust content in their atmospheres (Cheetham et al. 2019; Stolker et al. 2020b) although a scenario in which additional reddening is caused by circumplanetary material could also be considered given their age uncertainties.

For GQ Lup B, we computed synthetic H-band photometry from the VLT/SINFONI spectrum that will be analyzed in Section 4.5. The MKO H-band photometry is $14.02 \pm 0.13$ mag so the H–M' color is $2.04 \pm 0.15$ mag. The error is dominated by the uncertainty on the HST/NICMOS F171M flux from Marois et al. (2007) that was used to calibrate the H-band spectrum. A comparison with the spectral sequence of the high-gravity field dwarfs shows that the absolute H-band brightness of GQ Lup B is consistent with a mid M-type object, somewhat similar to PZ Tel B and HD 1160 B. When comparing with the AMES-Cond/Dusty isochrones, we estimated a photometric mass for GQ Lup B of ~30 $M_J$ at an assumed age of 3 Myr.

The color of GQ Lup B is ~1 mag redder than the synthetic predictions from the atmosphere models and ~1.2 mag redder than the field dwarfs of similar spectral type. For a mid M-type object, the atmosphere is expected to be cloudless, which can indeed be seen from the convergence of the AMES-Cond and AMES-Dusty isochrones in Figure 2 for masses $\gtrsim 20 M_J$. Given the previously detected hydrogen emission lines (Seifahrt et al. 2007; Zhou et al. 2014; Wu et al. 2017b), the red H–M' color is likely caused by thermal emission from dust grains in a protolunar disk, that is, a disk around a planetary- or substellar-mass companion in which satellites might be forming (e.g., Wu et al. 2017b; Pérez et al. 2019).

The color of GQ Lup B is comparable to the young ($\lesssim 10$ Myr), planetary-mass objects GSC 06214 B and ROXs 42 Bb. It has been suggested that GSC 06214 B hosts a dusty accretion disk, as inferred from the Paβ (Bowler et al. 2011) and Hα emission (Zhou et al. 2014), and the red SED colors (Bailey et al. 2013). Hydrogen emission lines have not been identified in the spectrum of ROXs 42 Bb (Bowler et al. 2014) and there is no hint of excess emission from a disk although the SED is reddened by $A_V \approx 1.7–1.9$ mag (Bowler et al. 2014; Currie et al. 2014). In contrast to GQ Lup B, the red colors of these early L-type objects may also point to enhanced cloud densities in their low-gravity atmospheres. Indeed, the AMES-Dusty predictions (i.e., with strongly mixed clouds) are consistent with GSC 06214 B and ROXs 42 Bb within $1–3\sigma$.

Thermal emission from a disk may cause a red H–M' color, but extinction alone could in principle have a similar effect. With the presence of an inclined disk, some extinction is expected to occur, in particular in the surface layer where presumably submicron-sized grains reside. In Section 4.5.2, we will estimate a visual extinction of $A_V = 2.3$ by modeling the optical and NIR spectra of GQ Lup B with an atmospheric model and assuming an interstellar medium (ISM)-like

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Color–magnitude diagram of $M_{H\alpha}$ vs. H–M'. The regular field dwarfs are color coded by M, L, and T spectral types (see discrete colorbar), the young and low-gravity objects are indicated with a gray square, and the directly imaged objects are labeled individually. GQ Lup B is highlighted with a red cross. The blue and orange lines show the synthetic colors computed from isochrones at 3 Myr. The black arrow is a reddening vector for the extinction, $A_V = 2.3$, that is inferred from the SED of GQ Lup B (see Section 4.5.2 for details).</figcaption></figure>

extinction, using the empirical relation from Cardelli et al. (1989). In that case, the reddening of H–M' will be ~0.4 mag (see arrow in Figure 2), which is a factor $\gtrsim 3$ smaller than what has been measured. A much larger reddening is therefore required, but such a scenario seems unlikely given the constraints on the atmospheric parameters and the actual detection of GQ Lup B at optical wavelengths.

### 4.3. Optical Spectral Type and Features

Figure 3 shows the extracted optical spectrum from MUSE, with the main spectral features labeled, compared to a main-sequence M9 spectral template from Kesseli et al. (2017). We estimated the optical spectral type of GQ Lup B by using *PyHammer* (Kesseli et al. 2017), which relies on the comparison of spectral indices with templates, and found a best-fit spectral type between an M8 and M9. As a second method, we used spectral type relations that are specific for young stars and found a spectral type of M7.4 with the TiO 8465 index from Herczeg & Hillenbrand (2014). This discrepancy of about one spectral subtype is expected and has been previously noted when comparing young pre-main-sequence objects to main-sequence stars (Herczeg & Hillenbrand 2014). With both methods we have not accounted for extinction (see Section 4.5.2) but the reddening mainly impacts the broadband SED slope, while the spectral indices that were used for inferring the spectral type are not much affected.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Optical spectrum of GQ Lup B with MUSE. The extracted spectrum (black line) is shown in comparison with a standard M9 V stellar template (blue line). The opacity sources of the identified spectral features are labeled in gray and blue with the latter in case a species is only visible in the stellar template. The inset shows a zoom to the Hα line with the arrows pointing to the stellar residuals (see Section 4.4 for details).</figcaption></figure>

The majority of the differences between the spectrum of GQ Lup B and the M9 template can be explained by the lower surface gravity of GQ Lup B. For example, we identify absorption from alkali metals, K I and Na I, but both doublets are narrower and weaker in GQ Lup B due to the decreased pressure broadening in low-gravity atmospheres (Allers & Liu 2013). The many TiO and VO bands dominate the spectrum and set the pseudo-continuum in both GQ Lup B and the main-sequence template. However, in the regions around 7400 and 8000 Å, it seems that VO might be enhanced in GQ Lup B compared to the template, which is another indicator for a low surface gravity (Martin et al. 1996; McGovern et al. 2004). Finally, we do not detect any absorption from metal hydrides (FeH, CrH, or CaH) in GQ Lup B. This is in contrast to older directly imaged, planetary-mass objects where FeH is clearly detected (e.g., 2MASS 0103 (AB)b; Eriksson et al. 2020), suggesting that FeH is not seen in the spectrum of GQ Lup B because its photosphere is located at low pressures. Metal hydrides have long been used as surface gravity indicators (e.g., Schiavon et al. 1997), and the lack of their features is again complementary evidence for the youth and low-gravity nature of GQ Lup B.

## 4.4. Hydrogen Emission Lines

In addition to the detection of molecular and atomic species, the optical spectrum reveals several emission lines, which are also labeled in Figure 3. Specifically, we detect a prominent Hα line and the Ca II infrared triplet, which are both spectral signatures for accretion (e.g., White & Basri 2003). There was no detection of Hβ but we derived a $1\sigma$ upper limit of $2.7 \times 10^{-19}$ W m$^{-2}$ by injecting artificial sources at a range of position angles while avoiding bright diffraction residuals. Previously, Paβ emission was identified by Seifahrt et al. (2007) in the SINFONI spectrum, but the line characteristics had not been analyzed. The NACO NB4.05 filter is used by the MIRACLES program to detect atmospheric emission at ~4 µm, but this narrowband filter is actually optimized for the detection of Brα emission, which could also be emitted by an accreting planet (Aoyama et al. 2020). With the SED analysis later on in Section 5.1, we find no significant excess flux in NB4.05 that may point to hydrogen line emission so it is consistent with a nondetection of Brα.

We inferred the emission line characteristics of the hydrogen (i.e., Hα and Paβ) and Ca II lines with the species toolkit (Stolker et al. 2020b). First, the continuum flux was estimated by fitting a third-order polynomial to a smoothed version of the

spectrum. This step was not required for the Hα line since the continuum level is consistent with the noise floor. Second, we fitted a Gaussian function to the (continuum-subtracted) emission line by evaluating a Gaussian likelihood function and using uniform priors. The posterior distributions of the three parameters (amplitude, mean, and standard deviation) where sampled with the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019), as implemented in UltraNest (Buchner 2021), and propagated into a line flux, equivalent width, FWHM, and radial velocity. The equivalent width could not be determined for Hα since the continuum of GQ Lup B was not detected at those wavelengths.

The measured and derived properties of the emission lines are listed in Table 3. With a resolution of R ~ 2500 (=120 km s$^{-1}$) in Paβ (see Seifahrt et al. 2007), the line might be (marginally) resolved. The width of the Hα line, on the other hand, is consistent with the instrument resolution of R = 2516 (i.e., 119 km s$^{-1}$; see Figure B.3 in Eriksson et al. 2020), considering that there were some stellar residuals due to local variations in the LSF, which may have led to slight inaccuracies with the spectral extraction of the Hα line. These effects were difficult to correct and caused under/over-subtraction of the stellar light on the blue/red side of the Hα line (see inset in Figure 3) so we made sure that these artificial features were excluded from the fit. There is some spread in the radial velocities (RVs) of the line centers compared to $v_{sys} = -2.8 \pm 0.2$ km s$^{-1}$ that was measured by Schwarz et al. (2016). This probably points to slight inaccuracies in the SINFONI and MUSE wavelength solutions although this is only at the ~10% level of their resolving power. Later on, we will interpret the hydrogen line emission in more detail, but the Ca II lines will not be further analyzed.

## 4.5. Atmospheric and Disk Modeling

### 4.5.1. Data and Model Preparation

To quantify in more detail the atmospheric and disk characteristics of GQ Lup B, we have complemented our optical spectrum and 4–5 µm photometry with the VLT/SINFONI JHK-band spectra from Seifahrt et al. (2007), Subaru/CIAO L' photometry from Marois et al. (2007), and the ALMA Band 6 and 7 upper limits from Wu et al. (2017b) and MacGregor et al. (2017). We have also adopted optical and NIR photometry from Marois et al. (2007) and Wu et al. (2017b) but these have not been included in the analysis given their small constraining power with
<table><caption>**Table 3**
Emission Line Measurements</caption><thead><tr><th>Line</th><th><i>F</i><sub>line</sub><br>(W m<sup>-2</sup>)</th><th><i>L</i><sub>line</sub><br>(L<sub>☉</sub>)</th><th>EW<br>(Å)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>RV<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Hα</td><td>(3.31 ± 0.04) × 10<sup>-18</sup></td><td>(2.38 ± 0.03) × 10<sup>-6</sup></td><td>...</td><td>108 ± 1</td><td>19 ± 1</td></tr><tr><td>Hβ</td><td>&lt;2.7 × 10<sup>-19</sup></td><td>&lt;2.0 × 10<sup>-7</sup></td><td>...</td><td>...</td><td>...</td></tr><tr><td>Paβ</td><td>(1.32 ± 0.01) × 10<sup>-18</sup></td><td>(9.53 ± 0.10) × 10<sup>-7</sup></td><td>-3.72 ± 0.04</td><td>237 ± 3</td><td>15 ± 1</td></tr><tr><td>Ca II λ8498</td><td>(3.35 ± 0.56) × 10<sup>-19</sup></td><td>(2.49 ± 0.42) × 10<sup>-7</sup></td><td>-1.78 ± 0.30</td><td>97 ± 16</td><td>-22 ± 8</td></tr><tr><td>Ca II λ8542</td><td>(4.40 ± 0.16) × 10<sup>-19</sup></td><td>(3.26 ± 0.12) × 10<sup>-7</sup></td><td>-2.29 ± 0.08</td><td>84 ± 3</td><td>-18 ± 1</td></tr><tr><td>Ca II λ8662</td><td>(4.06 ± 0.36) × 10<sup>-19</sup></td><td>(3.01 ± 0.27) × 10<sup>-7</sup></td><td>-1.67 ± 0.15</td><td>110 ± 5</td><td>2 ± 2</td></tr></tbody></table>

Note. The measurements have not been corrected for extinction ($A_V$ = 2.3 mag; see Section 4.5.2). The listed values for Hβ are 1σ upper limits.

respect to the spectra. The SINFONI J-band spectrum has been flux calibrated with the J-band magnitude of 14.90 from McElwain et al. (2007), and the H- and K-band spectra with the HST/ NICMOS F171M and F215N magnitudes from Marois et al. (2007). For the error bars of the spectra, we have adopted the per pixel S/Ns from Seifahrt et al. (2007; i.e., 100 in J and 30 in H and K) and sampled Gaussian errors. Additionally, we obtained a telluric spectrum with the SkyCalc interface (Noll et al. 2012) and scaled the error bars with the reciprocal of the transmission to account for systematic errors by telluric lines.

The combined data were analyzed with the species toolkit by fitting a grid of BT-Settl spectra (Allard et al. 2012), for which we adopted the CIFIST release that used the solar abundances from Caffau et al. (2011). BT-Settl is a 1D radiative-convective equilibrium model, which accounts for nonequilibrium chemistry and cloud physics. Clouds are not expected to form though at the photospheric temperature of GQ Lup B, $T_{eff} \sim 2500$ K, and the carbon content will be mainly locked up in CO molecules. To model the ALMA fluxes, we extended the spectra into the millimeter regime by fitting the fluxes at $\lambda > 50$ µm with a power-law function in log-log space. With the grid of model spectra at hand, we carried out several comparisons with the optical and NIR data.

### 4.5.2. Atmospheric Emission and Extinction

We first analyzed the MUSE and SINFONI spectra by considering atmospheric emission alone. This approach allowed us to constrain the main atmospheric parameters ($T_{eff}$, log g, and $R_p$) while assuming negligible emission from a disk up into the K band. In order to account for potential inaccuracies with the absolute flux calibration of the spectra, we fixed the MUSE and K-band spectra (both were calibrated with HST photometry) while fitting a scaling parameter, a, for the SINFONI J- and H-band spectra (the SINFONI spectra were obtained at separate epochs; see Seifahrt et al. 2007). To account for interstellar and/or local (e.g., due to a disk) reddening, we adopted the empirical relation from Cardelli et al. (1989) and included the visual extinction, $A_V$, as an additional parameter while fixing the reddening to the standard value for the diffuse ISM, $R_V = 3.1$. Extinction in a disk is expected to occur mainly in the surface layer where grains are expected to be small. Hence, adopting an ISM relation might be a reasonable first choice for accounting for extinction even though the true grain properties may deviate. Therefore, we also tested a power-law parameterization (to account for the unknown dust properties) but obtained a similar result, in which the wavelength-dependent extinction appeared comparable to the empirical ISM relation.

The modeling approach that we first tested was a parameter retrieval, which sampled spectra from the interpolated model grid and used Bayesian inference for estimating the uncertainties (i.e., similar as in Mollière et al. 2020 and Kammerer et al. 2021). The residuals showed that the fit was limited by systematic errors in the models and/or data such that the uncertainties on the model parameters were highly underestimated. We tested different kernels for modeling the correlated noise with a Gaussian process, but this only had a marginal effect on the width of the posterior distributions. At the level of precision and spectral resolution of the data, the fit was likely limited by inaccuracies caused by the subgrid sampling, while spectral features are expected to change in a nonlinear manner. A more sophisticated approach is required to account for such interpolation errors (e.g., Czekala et al. 2015). In this work, we opt for a simplified approach by comparing the spectra from the discrete model grid and evaluating the G statistic from Cushing et al. (2008; see their Equation (1)). It is mathematically a weighted $\chi^2$ statistic but is not expected to follow a $\chi^2$ distribution because of the systematic errors. As weights for the fluxes, we use the widths of the wavelength bins as suggested by Cushing et al. (2008).

We selected all available model spectra with $T_{eff}$ in the range of 2000–3500 K with a 50 K sampling up to 2400 and 100 K at larger temperatures. The grid spacing of the surface gravity, log g, is 0.5 dex in the range of 2.5–5.5 dex. Additionally, we tested $A_V$ in the range of 0–5 mag in steps of 0.1 mag. Each model spectrum was first smoothed to the approximate resolving power of the instrument ($R = 3000$ for MUSE and $R = 2500, 4000, 4000$ for SINFONI J, H, and K) and then resampled to the wavelengths of the data. When comparing with the data, we fixed the distance and fitted the radius, $R_p$, as flux scaling such that G is minimized. Since optical and NIR wavelengths trace different parts of the atmosphere, we compare the model grid with the MUSE and SINFONI both separately and combined.

The broadband analysis of the SED did not constrain log g since the impact of this parameter (at the considered $T_{eff}$) on the spectral fluxes is smaller than the typical systematic errors. The optical spectrum is however characterized by several low-gravity features (see Section 4.3), which is further confirmed by the NIR spectra. We calculated $K I_J = 1.05$, where $K I_J$ is a gravity-sensitive spectral index that has been defined by Allers & Liu (2013). The value is small compared to regular field dwarfs with a similar spectral type (see Figure 22 in Allers & Liu 2013), confirming the low surface gravity. We determined, on visual inspection, that the optical and NIR alkali lines are best fitted with log g = 3.5–4.0 when fixing $T_{eff} = 2700$ K. Specifically, the K I doublet near 1.25 µm is best matched with log g = 4.0, but the
<table><caption>**Table 4**
Spectral Analysis of GQ Lup B with Atmospheric Models</caption><thead><tr><th></th><th>T<sub>eff</sub><br>(K)</th><th>log g<sup>a</sup></th><th>R<sub>p</sub><br>(R<sub>J</sub>)</th><th>a<sub>J</sub><sup>b</sup></th><th>a<sub>H</sub><sup>b</sup></th><th>A<sub>V</sub></th><th>T<sub>disk</sub><br>(K)</th><th>R<sub>disk</sub><br>(R<sub>J</sub>)</th><th>ΔG<sup>c</sup></th></tr></thead><tbody><tr><td>Parameter range</td><td>2000–3500</td><td>2.5–5.5</td><td>1–5</td><td>0.1–5</td><td>0.1–5</td><td>0–5</td><td>10–2000</td><td>5–1000</td><td></td></tr><tr><td><em>MUSE spectrum</em></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BT-Settl</td><td>2500</td><td>4.0</td><td>2.85</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0</td></tr><tr><td>BT-Settl + ext.</td><td>2700</td><td>4.0</td><td>4.13</td><td>...</td><td>...</td><td>2.7</td><td>...</td><td>...</td><td>-156</td></tr><tr><td><em>SINFONI spectra</em></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BT-Settl</td><td>2600</td><td>4.0</td><td>3.55</td><td>1.84</td><td>1.19</td><td>...</td><td>...</td><td>...</td><td>0</td></tr><tr><td>BT-Settl + ext.</td><td>2600</td><td>4.0</td><td>3.55</td><td>1.84</td><td>1.19</td><td>0.0</td><td>...</td><td>...</td><td>0</td></tr><tr><td><em>MUSE + SINFONI spectra</em></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BT-Settl</td><td>2350</td><td>4.0</td><td>3.93</td><td>1.68</td><td>1.06</td><td>...</td><td>...</td><td>...</td><td>0</td></tr><tr><td>BT-Settl + ext.</td><td>2700</td><td>4.0</td><td>3.77</td><td>1.32</td><td>1.03</td><td>2.3</td><td>...</td><td>...</td><td>-11587</td></tr><tr><td><em>MUSE + SINFONI spectra + 3–5 µm photometry + ALMA</em></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BT-Settl + ext. + disk<sup>d</sup></td><td>2700</td><td>4.0</td><td>3.77</td><td>1.32</td><td>1.03</td><td>2.3</td><td>461 ± 2</td><td>65 ± 1</td><td>...</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The surface gravity, log g, was estimated from the gravity-sensitive alkali doublets in the MUSE and J-band spectra.

<sup>b</sup> Scaling parameter to account for potential inaccuracies in the absolute flux calibration.

<sup>c</sup> The relative goodness-of-fit is calculated for each combination of spectra separately.

<sup>d</sup> The atmospheric, scaling, and extinction parameters were fixed to the best-fit values of fitting the MUSE and SINFONI spectra.

predicted KI and NaI lines at ~0.77 µm and ~0.82 µm, respectively, are a bit too strong and better matched with log g = 3.5. Similarly, with log g = 4.0, the VO band at 0.74 µm is not sufficiently deep, and it would require a lower surface gravity to match the model spectrum with the observed fluxes at those wavelengths. Given the negligible impact on the overall fit, we set log g = 4.0 in the remaining analysis.

The results from the spectral analysis are presented in Table 4. Several things can be noticed. Fitting the MUSE spectrum alone or combined with the SINFONI spectra yields a better fit when including $A_V$. This is seen from the last column in Table 4, which lists the G statistic relative to the case without extinction. The SINFONI spectra, on the other hand, do not constrain $A_V$. Instead, the J- and H-band scaling factors are relatively large so these parameters may mimic any reddening. Regarding $T_{eff}$, we obtained 2700 K when fitting the MUSE spectrum alone or combined with the SINFONI spectra. The NIR fluxes, however, were under/overestimated when fitting the MUSE spectrum without/with extinction. Indeed, both $R_p$ and $A_V$ are different and expected to be more accurate when using spectra across a broad wavelength range.

When fitting the combined MUSE and SINFONI spectra, the H-band scaling is small and therefore consistent with the HST photometry that was used for the calibration, although the residuals reveal a slope at the red end of the H band (see Figure 5). This differential slope between model and data may point to an issue with the calibration of the pseudo-continuum (see also Figure 5 in Lavigne et al. 2009 with a comparison of JHK spectra from different instruments). The flux calibration of the J-band spectrum is expected to be the least accurate since we had adopted the photometric flux from a ground-based observation while the other spectra were calibrated with HST photometry. When using all spectra and including $A_V$, the required flux scaling for the J-band spectrum is relatively large, $a_J = 1.3$ (see also contours in Figure 4). However, the J-band photometry was computed by McElwain et al. (2007), while their spectra had been obtained at an airmass of 1.8 so a correction factor of 1.3 (i.e., ~0.3 mag) seems reasonable.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Contours of the goodness of fit from the comparison of the MUSE and SINFONI spectra of GQ Lup B with the BT-Settl model grid. The color scale shows the G statistic relative to the best-fit parameters ($T_{eff} = 2700$ K, log g = 4.0, and $A_V = 2.3$ mag), which are indicated with a red cross. The white contours show the best-fit scaling factor that is applied for recalibrating the J-band spectrum.</figcaption></figure>

All together, combining the MUSE and SINFONI spectra is expected to give the most accurate constraints. Therefore, the best-fit parameters from this work are $T_{eff} = 2700$ K, log g = 3.5–4.0, $R_p = 3.8 R_J$, and $A_V = 2.3$ mag, which is shown in comparison with the spectra and some of the available photometric fluxes in Figure 5. We can also compare the derived atmospheric parameters with predictions from evolutionary models. At an age of 3 Myr, the AMES-Dusty isochrones (Chabrier et al. 2000) predict that a mass of $M_p = 32 M_J$ has $T_{eff} \approx 2700$ K, log g = 3.8 (and therefore $R_p = 3.7 R_J$), and $M_H = 8$ mag. While hinging on the uncertain age estimate, this is in agreement with the inferred parameters from the spectra, as well as the absolute H-band brightness (see Figure 2). With the estimated $A_V = 2.3$, the extinction in the H band is ≈0.4 mag so the object would be a bit brighter than the considered model prediction at 3 Myr (see reddening vector in Figure 2). This might be reasonable though given the margins of uncertainties and assumptions with the mass estimate.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of the best-fit model spectrum (black line) with the MUSE (pink line; top panel) and SINFONI spectra (blue line; bottom panel) of GQ Lup B. The BT-Settl spectrum has been smoothed to R = 3000, but the residuals are shown at the actual resolution of each spectrum. For reference, we adopted the magnitudes from Marois et al. (2007) and Wu et al. (2017b) and converted these into fluxes (gray markers) with *species*. The horizontal error bars indicate the FWHM of the filter profiles, and open markers are the respective synthetic photometry from the best-fit spectrum.</figcaption></figure>

While the best-fit model matches well with the overall morphology and spectral features, there are some systematic variations in the residuals. Apart from the discrepancy in the *H* band, the optical absorption bands are somewhat muted compared to the model spectrum. Possibly this is caused by an inaccuracy with the calibration of the pseudo-continuum. Alternatively, the effect could be real if the absorption features are veiled by excess emission that originates from accretion hot spots (e.g., Calvet & Gullbring 1998). It remains to be investigated if such continuum emission is to be expected at ~0.7–0.9 µm. Differences between the MUSE spectrum and the MagAO photometry are likely attributed to the variability of GQ Lup (Δ*R* ≈ 0.4 mag, Δ*I* ≈ 0.3 mag; Broeg et al. 2007) since these fluxes were calibrated relative to the star (see Wu et al. 2017b for details). In Figure 6, we show the SED from optical to MIR wavelengths, together with the best-fit model spectrum. For clarity, here every 50th wavelength point is used. Interestingly, the best-fit model underpredicts the 3–5 µm fluxes by 1.3*σ*, 4.9*σ*, and 7.3*σ* in *L'*, NB4.05, and *M'*, respectively, when considering atmospheric emission alone (see dotted line in Figure 6), which points to an additional luminosity component in the SED at those wavelengths.

4.5.3. *Thermal Emission from a Dusty Disk*

The comparison with the atmospheric model spectra
revealed excess emission at 3–5 µm. While there is some
uncertainty in the absolute calibration of the *J*-band spectrum in
particular, the HST photometry by Marois et al. (2007) is
expected to be most accurate (ignoring possible variability of
GQ Lup B) so we do not expect any spurious reddening in the
*M* band. Instead, the youth of GQ Lup B, the detected hydrogen
emission lines, and the red NIR-*M'* color all point to the
presence of a protolunar disk.

To explore this hypothesis, we model the optical to MIR SED
with a combination of atmospheric emission and a blackbody
component that accounts for thermal emission from a dusty disk.
To do so, we fix all parameters to the best-fit values from fitting
the MUSE and SINFONI spectra and include, in addition to the
spectra, the *L'*, NB4.05, and *M'* fluxes, and the ALMA upper
limits in the fit. We retrieve the effective disk temperature, *T*<sub>disk</sub>,
and radius, *R*<sub>disk</sub>, with the Bayesian framework of the *species*
toolkit, for which we made again use of the nested sampling
implementation of UltraNest (Buchner 2021). We used a
Gaussian likelihood function for the model evaluation while
mapping out the posterior space with 1000 live points and
assuming uniform priors for both parameters (see Table 4).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Spectral energy distribution of GQ Lup B. The various spectral and photometric data points are shown with different colors as indicated in the legend. The NACO NB4.05 and *M'* fluxes are highlighted with a red star. The black solid line is the best-fit model (smoothed to *R* = 500), including the additional blackbody component and the extinction applied. The black dotted and black dashed lines are the spectra of the atmospheric and blackbody component, respectively. The open markers are the synthetic fluxes for all filters, with the same shape and edge color as the empirical fluxes. The lower panel shows the residuals (i.e., data minus model) for the data points that were used in the fit.</figcaption></figure>

The best-fit model spectrum, with the atmospheric and blackbody emission combined, is shown in Figure 6, and the posterior distributions of the disk parameters can be found in Figure 10. The modeled disk emission appears at wavelengths longer than ~3 µm with indeed negligible flux in the K band, as was to be expected from the residuals in Figure 5. The residuals of the *L'* NB4.05, and *M'* fluxes are ≤0.5σ, indicating a good fit with a blackbody spectrum. For NB4.05, the transmission of the filter is optimized for the detection of Brα, but there is no indication of such emission from GQ Lup B. Instead, the fit shows that the NB4.05 flux is consistent with continuum emission alone.

The retrieved parameters for a blackbody model are $T_{\text{disk}} = 461 \pm 2$ K and $R_{\text{disk}} = 65 \pm 1 R_J$ (see bottom row in Table 4). From this, we derive a disk-to-companion luminosity ratio of ≈0.25, which implies that 25% of the atmospheric emission may get reprocessed by the disk if we ignore additional processes (e.g., due to accretion) that may heat the disk internally. At millimeter wavelengths, Wu et al. (2017b) and MacGregor et al. (2017) reported an rms noise level of 39 and 50 µJy beam⁻¹ in Bands 6 (λ = 1.3 mm) and 7 (λ = 870 µm), respectively. The synthetic ALMA Band 6 and 7 fluxes from the fit are $3.8 \times 10^{-24}$ and $1.8 \times 10^{-23}$ W m⁻² µm⁻¹, that is, a factor 18 and 11 below, and therefore consistent with, the ALMA upper limits. Finally, resolving a disk radius of 65 $R_J$ would require a sub-mas resolution so the approximate radius that is probed by our observations is too compact to be resolved with ALMA.

## 4.6. Orbital Architecture

The orbital architecture of GQ Lup B may provide clues about its formation and dynamical history. This is of particular interest because the polarimetric imaging observations by van Holstein et al. (2021) revealed an asymmetric, spiral-like structure in the CSD, which points to a gravitational interaction by the companion. To estimate the orbital elements, we used *orbitize*!<fnref n="21" /> (Blunt et al. 2020) to fit the astrometry from Neuhäuser et al. (2008),

Ginski et al. (2014), and this work (see Table 2). The astrometric data are shown in the inset of Figure 7 with both the R.A. and decl. increasing over time. We also adopted the RV of GQ Lup, $v_{\text{sys}} = -2.8 \pm 0.2$ km s⁻¹ from Schwarz et al. (2016) and the RV of GQ Lup B, $2.0 \pm 0.4$ km s⁻¹, which breaks the 180 deg degeneracy in the longitude of the ascending node. The posterior distributions of the orbit parameters were sampled with a parallel-tempered MCMC algorithm (Foreman-Mackey et al. 2013; Vousden et al. 2016) while marginalizing over the uncertainty on the parallax ($6.49 \pm 0.03$ mas; Gaia Collaboration et al. 2016, 2021) and stellar mass ($1.03 \pm 0.05 M_{\odot}$; MacGregor et al. 2017).

For the parameter estimation, we used 20 temperatures, 500 walkers, and 50,000 steps per walker. The first 40,000 steps were removed as burn-in, and we conservatively selected every 20th step of each walker to exclude correlations between steps. The remaining samples appeared converged on manual inspection, which was confirmed by calculating the integrated autocorrelation time of the chains. The posterior distribution of the orbit fit can be found in Figure 11 of the Appendix. We obtained a constraint on the semimajor axis and eccentricity of $a = 117_{-23}^{+24}$ au and $e = 0.24_{-0.17}^{+0.32}$. The fit favors circular and low-eccentricity orbits although intermediate to high values are not excluded. The orientation of the orbital plane relative to the sky is given by the inclination and the longitude of the ascending node, for which we derived $i = 60_{-9}^{+5}$ deg and $\Omega = 265_{-10}^{+19}$ deg.

Figure 7 shows the orbital constraints in comparison with the $Q_φ$ image from van Holstein et al. (2021) that was obtained with VLT/SPHERE in the H band (i.e., the same data sets that we used for the astrometry in Section 3.3). The orbits were calculated by drawing 150 random samples from the posterior distribution. To enhance the visibility of the spiral arm (south in Figure 7), we used *diskmap*<fnref n="22" /> (Stolker et al. 2016) to apply an $r^2$-scaling that accounts for both the orientation of the midplane, for

<footnote n="21">²¹ https://orbitize.readthedocs.io</footnote>

<footnote n="22">²² https://diskmap.readthedocs.io</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Scattered light image of the circumstellar disk around GQ Lup with sampled orbital configurations of GQ Lup B. The image shows the r²-scaled, polarized flux on a linear color scale (see main text for details). The central part of the image (~0.2/2 in diameter) was obscured by a coronagraph. The inset shows a zoom to the region of the astrometric data, obtained from 2004 to 2019. The orange markers are the four astrometry points from this work. North is up and east is left.</figcaption></figure>

which we adopted an inclination of 60.5 ± 0.5 deg and a position angle of 346 ± 1 deg from MacGregor et al. (2017), and a power-law approximation for the vertical extent of the disk surface. We note that the apparent, polarized flux of GQ Lup B is caused by the subtraction of the unresolved, polarized, stellar component, while the actual companion is not polarized up to ~0.2% (see van Holstein et al. 2021, for more details).

With the detected spiral arm by van Holstein et al. (2021) and the kinematics from MacGregor et al. (2017), we can uniquely determine the orientation of the disk by assuming that the spiral arm has a trailing morphology with respect to the rotation direction of the CSD. Since the red/blueshifted velocities are located in the north/south direction, this implies that the near side of the disk is in the western direction. Knowing the inclination and longitude of the ascending node of the CSD, we constrained the mutual inclination between the orbit and the CSD to 84 ± 9 deg by propagating the posterior distributions from the orbital fit and accounting for the uncertainty on the disk orientation.

# 5. Discussion

## 5.1. Mid-infrared Excess from a Protolunar Disk

Infrared colors have been a powerful probe for the inference of CSDs around young stars for decades (e.g., Lada et al. 2000). Surveys of young, substellar objects in nearby star-forming regions have also revealed that many show excess emission from a dusty disk (Muench et al. 2001; Liu et al. 2003). Jayawardhana et al. (2003a) carried out a systematic study in the L' band and reported that objects around the substellar boundary have inner disk lifetimes comparable to T Tauri stars, indicating that isolated brown dwarfs form through similar mechanisms as solar-mass stars. Similar to brown dwarfs, also young, planetary-mass objects

are known to host disks (e.g., Zapatero Osorio et al. 2007). The thermal excess flux correlates with detected Hα emission, pointing to a common origin of a dusty accretion disk (Liu et al. 2003). For low-mass objects, the accretion rates are typically lower, and the Hα line profiles are more narrow compared to T Tauri stars (Jayawardhana et al. 2003b).

The detection of both MIR excess and a strong Hα line in the SED of GQ Lup B seems in line with expectations from isolated, substellar objects. Yet, the formation mechanisms might be different, and the spiral-arm interaction may affect the accretion and disk characteristics of GQ Lup B. Both the disk of GQ Lup and GQ Lup B may have been truncated and/or perturbed by the dynamical interaction, which could be the origin of the somewhat large extinction, $A_V = 2.3$ mag. While the modeling in Section 4.5 revealed a significant excess emission at NB4.05 and $M'$, the $1\sigma$ deviation from the expected atmospheric emission in L' had in fact already been noted by Marois et al. (2007).

Thanks to the large separation and brightness of GQ Lup B, we were able to precisely measure the 4–5 µm fluxes and estimate the blackbody parameters, while, in comparison, this was much more challenging for the closer-in planet PDS 70 b (Stolker et al. 2020a). Assuming that the emission traces a single blackbody temperature (e.g., the strongly irradiated inner radius of the disk), we were able to constrain the disk parameters to $T_{\text{disk}} = 461 \pm 2$ K and $R_{\text{disk}} = 65 \pm 1$ R<sub>J</sub> (see Section 4.5.3). The three photometric fluxes appear well described by a blackbody spectrum although the true disk structure might be more complex (i.e., it will have a vertical and radial temperature gradient). In that regard, we note that the uncertainties on the disk parameters only reflect the statistical errors from the fit. The peak of the blackbody emission lies at ≈6–7 µm so it would make an appealing target for MIR spectroscopy with the James Webb Space Telescope (JWST) to investigate the disk and dust properties in more detail.

## 5.2. Evidence of Satellite Formation?

The MIR excess that is inferred from the SED is expected to trace the thermal emission from (sub)micron-sized dust grains, which are heated by the radiation field of GQ Lup B and accretion processes within the disk. For a centrally irradiated disk, the hottest and brightest part is the inner radius that is directly illuminated. The effective temperature, $T_{\text{disk}} = 461 \pm 2$ K, is therefore expected to trace the inner radius of the disk, which is much cooler than the effective temperature of the companion atmosphere, $T_{\text{eff}} \approx 2700$ K (see Section 4.5). This speculatively suggests that a large central cavity is present in the disk, which has been created by some mechanism. The estimated disk radius, $R_{\text{disk}} = 65 \pm 1$ R<sub>J</sub>, may provide further evidence for a cavity since it is large compared to $R_p \approx 3.8$ R<sub>J</sub> from the atmosphere. While $R_{\text{disk}}$ might be reasonable first-order estimate, a broader wavelength coverage and more detailed modeling are required to constrain the actual inner radius of the disk.

Several mechanisms can create a central cavity in the disk structure. The sublimation radius of the dust lies at $R_s = \frac{1}{2}\sqrt{Q_R}(T_{\text{eff}}/T_s)^2 R_p$, where $Q_R = Q_{\text{abs}}(T_{\text{eff}})/Q_{\text{abs}}(T_s)$ is the ratio of the absorption efficiencies of the dust, $Q(T)$, for radiation at color temperature, $T$, of the incident and reemitted field (Tuthill et al. 2001). When assuming blackbody grains, $Q_R = 1$, with a sublimation temperature of $T_s = 1500$ K, and adopting $T_{\text{eff}}$ and $R_p$ of GQ Lup B, we obtain $R_s \approx 6$ R<sub>J</sub>. For a relatively cool object such as GQ Lup B, $Q_R$ is not expected to be larger than unity by a factor of a few (i.e., due to inefficient
cooling), when considering different sizes of silicate grains (Monnier & Millan-Gabet 2002). The central cavity that is inferred from the SED analysis can therefore not be explained by sublimation of silicate dust grains, which is also not surprising since $T_{\text{disk}} \ll T_{\text{s}}$. Instead, $T_{\text{disk}}$ lies closer to the freeze-out location of H$_2$O, which would be located a bit further outward. Therefore, the H$_2$O ice line could potentially have provided an efficient growth mechanism for the assembly of satellites. Heller & Pudritz (2015) simulated the evolution of such ice lines in disks around super-Jovian gas planets. For their most massive case, 12 $M_J$, the ice line would (at an early phase) be located at a comparable radius as inferred from the SED, but its location evolves strongly over time. While GQ Lup B is at least a factor 10 more massive than Jupiter, it is interesting to note that the four Galilean moons, which have orbits at 6–27 $R_J$, would all fit within the estimated disk radius, $R_{\text{disk}} \approx 65 R_J$, so the cleared area could be sufficiently large for hosting a multisatellite system.

Attempts with ALMA to detect thermal emission from pebble-sized grains at millimeter wavelengths yielded only nondetections, specifically, $3\sigma$ upper limits of 150 µJy at 880 µm by MacGregor et al. (2017) and 120 µJy at 1.3 mm by Wu et al. (2017a). Wu et al. (2017a) argued that the nondetections for GQ Lup B and other PMCs may imply that their disks could be very compact and optically thick in order to sustain a few megayears of accretion. Alternatively, we propose a scenario in which the reservoir of millimeter-sized grains has actually been depleted by the formation of satellites, and a central cavity may have opened within the inner disk structure. Large grains are expected to drift so they may have provided a continuous influx of solids for the assembly of satellites (Shibaike et al. 2017). Without any pressures traps, pebbles would ultimately get accreted by GQ Lup B, but the disk can sustain gas and small grains on longer timescales. The remaining gas gets accreted onto the companion with $\dot{M} \approx 10^{-6.5} M_J \text{ yr}^{-1}$ (see Section 5.4) while the disk may get replenished each time GQ Lup B crosses the plane of the CSD.

Alternatively, the central disk region may have been truncated by a magnetic field of GQ Lup B if such a magnetic field is sufficiently strong and the accretion rate is not too high (Lovelace et al. 2011). In case of magnetospheric accretion, the H$\alpha$ emission may come from the hot spots on the atmosphere that are fed by accretion funnels from the disk. The magnetic field evolution of brown dwarfs and giant planets has been calculated by Reiners & Christensen (2010). For a 25 $M_J$ object, the predicated magnetic field strength is $B_p = 2$ kG at an age of 3 Myr. Together with the estimated $R_p = 3.8 R_J$ and $\dot{M} \approx 10^{-6.5} M_J \text{ yr}^{-1}$, we calculated a truncation radius of $R_T \approx 600 R_J$ by using the relation from Lovelace et al. (2011; see their Equation (2)). This truncation radius appears large compared to the $R_{\text{disk}}$ that we derived from the SED, but $B_p, M_p$, and $\dot{M}$ are particularly uncertain. With the same parameters, but changing $B_p$ to 40 G, we are in the situation where $R_T \approx R_{\text{disk}}$. Therefore, a weak magnetic field may explain the inferred cavity size if at that radius material is able to couple with the field lines. Magnetospheric accretion is also expected to leave an imprint on the morphology of the hydrogen line profiles and would imply a small filling factor. There is no evidence for such effects in the current data (see Section 5.4 for more details) but a confirmation at higher spectral resolution is required.

## 5.3. System Architecture and Formation Pathway

The orbital analysis in Section 4.6 constrained the mutual inclination between the orbit of GQ Lup B and the CSD to

$84 \pm 9$ deg. Similarly, Wu et al. (2017b) already suggested that these two planes are possibly misaligned, as inferred from the orbital solutions by Ginski et al. (2014) and Schwarz et al. (2016). The authors also pointed out that the CSD is misaligned with the stellar spin axis, which has an inclination of $27.5 \pm 5$ deg (Broeg et al. 2007). This would not be unusual for a T Tauri star with a strong magnetic field, but could in this case also have been induced by a torque from GQ Lup B. While the orbit has a near-polar configuration with respect to the CSD, the misalignment with the spin axis of the star would be different although still significantly misaligned.

A large mutual inclination may have two origins. If the companion has formed from the circumstellar disk, then this points to a rather dynamical history in which, for example, a gravitational interaction with the CSD or a catastrophic event has placed GQ Lup B on a highly misaligned orbit. The second companion that was recently discovered at a projected distance of 2400 au may also have influenced the orbital configuration of GQ Lup B (Alcalá et al. 2020). A scattering event typically results in a high eccentricity (Scharf & Menou 2009; Nagasawa & Ida 2011), which is not in agreement with the possibly low eccentricity that we inferred from the orbital fit (see Figure 11). On the other hand, the large separation and moderate companion-to-star mass ratio may also imply a stellar-like formation pathway for GQ Lup B, which could also be the origin of the substantial mutual inclination (e.g., Jensen & Akeson 2014; Bate 2018).

Although the mutual inclination seems constrained from the orbital fit, some caution is required since modeling a small orbit coverage may lead to a bias in the inferred eccentricity and inclination (Ferrer-Chávez et al. 2021). From the posterior, we derived an orbital period of $\log P/\text{yr} = 3.1 \pm 0.1$ so the astrometric points cover about ~1% of the full orbit. This may not be sufficient to robustly distinguish between an inclined orbit with a small eccentricity and a face-on orbit with a large eccentricity. On the other hand, the difference in RV between GQ Lup A and B that was measured by Schwarz et al. (2016) does exclude a face-on orbit. Higher-precision astrometric measurements, such as with VLTI/GRAVITY, may mitigate potential biases and enable more accurate constraints on the orbital parameters (Gravity Collaboration et al. 2019).

In contrast to the companion's orbit, the inclination of its disk is not well constrained. van Holstein et al. (2021) obtained an upper limit of 0.2% polarization on the scattered light flux from the disk of GQ Lup B. The authors suggested that the disk may have a moderate or low inclination and/or has a low dust surface density (see their Figure 14). If the disk of GQ Lup B is coplanar with its orbit, then the inclination would be $\approx 60$ deg, as inferred with the orbital fit (see Figure 11). That seems indeed a moderate inclination and may explain the polarization nondetection. On the other hand, the simulations by van Holstein et al. (2021) assumed a small inner radius of the disk (i.e., comparable to the dust sublimation radius), while a central cavity would decrease the fractional polarization because the scattered light flux is lower at a larger radius. Therefore, a large inner radius in a strongly inclined disk cannot be excluded.

The extinction provides another constraint on the disk orientation. In Section 4.5.2, we estimated $A_V \approx 2.3$ mag (i.e., an optical depth of $\tau \approx 2$) from the broadband SED, which is larger than what has been inferred for the star (e.g., $0.4 \pm 0.2$ mag; Seperuelo Duarte et al. 2008). In case of GQ Lup B, the extinction may occur locally, for example, in the
disk surface, which would require that the disk is sufficiently inclined. Alternatively, the extinction might be caused by material that gets channeled by the spiral arm from the CSD to the companion. Any dust in this inflow of material could somewhat enshroud both the companion and its disk. A more detailed constraint on the disk inclination would be valuable to understand if gravitational torques from the star could have misaligned the disk (e.g., Batygin 2012)—and therefore also the orbits of potential satellites—with respect to the rotation axis of GQ Lup B. If the companion has a strong magnetic field that is coupled to its disk, then it might also be driving an inner disk warp. In that case, the disk could be misalignment with the orbital plane and may cause a variable extinction. High-precision photometric monitoring of GQ Lup B could potentially reveal such variable changes in the disk structure.

## 5.4. Constraints on the Accretion Processes

Hydrogen emission lines have been commonly used as tracers for accretion. Similar to stars, also brown dwarfs and giant planets radiate their accretion energy in the form of line and continuum emission (e.g., Zhou et al. 2014; Zhu 2015; Aoyama et al. 2018). Located in the Lupus I cloud, GQ Lup B is a few million years old so previously detected hydrogen lines had been linked to the possible presence of disk around the companion. Several authors detected $H\alpha$ emission with ground- and space-based optical photometry (Marois et al. 2007; Zhou et al. 2014; Wu et al. 2017b). The $Pa\beta$ line was found in the SINFONI spectrum by Seifahrt et al. (2007) but interestingly the line was not detected in the J-band spectra by McElwain et al. (2007) and Lavigne et al. (2009), even though the S/N was sufficiently high. This indicates that the accretion rate might be variable, for example, due to density perturbations at the shock front, or possibly GQ Lup B even went through a quiescent state. Associated extinction effects that may occur in the vicinity of the atmosphere (e.g., by the disk and/or infalling material) may therefore be variable as well.

Line fluxes and ratios are a powerful diagnostic for the characterization of the accretion shock on a circumbinary disk (CPD) or planet surface (Aoyama et al. 2020; Hashimoto et al. 2020; Uyama et al. 2021). In Section 4.4, we measured the $H\alpha$ flux and estimated an upper limit for $H\beta$, yielding an extinction-corrected, $1\sigma$ upper limit on the line ratio of $F_{H\beta}/F_{H\alpha} < 0.17$. To quantify the accretion rate, we use the line predictions from Aoyama et al. (2018) that were calculated with a detailed radiative hydrodynamical model, of which the free parameters are the preshock velocity, $v_0$, and hydrogen number density, $n_0$. The line fluxes from the model are converted into $L_{H\alpha}$ by fixing the radius of the shock surface to $R_p = 3.8 R_J$ (see Section 4.5.2) and assuming a filling factor of $f_{\text{fill}} = 1$. We can also map $(v_0, n_0)$ to $M_p$ (i.e., assuming that $v_0$ is equal to the freefall velocity) and $\dot{M}$ with Equations (7)–(9) from Aoyama et al. (2020) such that the accretion rate can be extracted at the $L_{H\alpha}$ of GQ Lup B, which is shown in Figure 8. Later on, we will estimate a filling factor of $f_{\text{fill}} > 0.6$ from the line ratio, which would imply that the actual $\dot{M}$ could be smaller by a factor $\sim 2$.

To break the degeneracy between the shock density and velocity, we consider two constraints on the companion mass. First, in Section 4.4, we showed that the $H\alpha$ line is not resolved at the resolving power of MUSE (i.e., $\approx 120 \text{ km s}^{-1}$ in $H\alpha$). If the gas reaches the planet surface with the freefall velocity, this would point to a low mass of GQ Lup B ($\lesssim 15 M_J$; see

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Constraints on the mass accretion rate from the measured $H\alpha$ luminosity. The white contour shows $L_{H\alpha}$ from GQ Lup B with an extinction correction of $A_V = 2.3$ mag, and the surrounding shaded region indicates $A_V$ in the range 0–3 mag. The horizontally dashed contours are the accretion rates, $\log(\dot{M}/(M_J \text{ yr}^{-1}))$, calculated from the model predictions and GQ Lup B parameters (see main text for further details), the vertically dashed lines are the freefall velocity for different masses of GQ Lup B, and the purple dotted line is the resolving power of MUSE at $H\alpha$.</figcaption></figure>

Figure 8), possibly making GQ Lup B a PMC, with an accretion rate of $\dot{M} \approx 10^{-6.5}-10^{-6} M_J \text{ yr}^{-1}$. However, as noted previously, there was possibly some oversubtraction of the stellar residuals at $H\alpha$, which could somewhat bias the line flux and width. Indeed, we measured a broader line for $Pa\beta$ (FWHM $\approx 237 \text{ km s}^{-1}$), which could imply a higher mass object ($M_p > 40 M_J$). Figure 8 shows that this mass lies at the edge of the grid and indicates an accretion rate that is a bit smaller, $\dot{M} \approx 10^{-6.5} M_J \text{ yr}^{-1}$. As a second approach, we adopted the mass estimate from Section 4.5.2, which we derived by comparing the atmospheric parameters and H-band magnitude with the AMES-Cond isochrone at 3 Myr, yielding $M_p \approx 30 M_J$. Therefore, the inferred mass from the $Pa\beta$ line and the H-band flux are roughly in agreement so we adopt $\dot{M} \approx 10^{-6.5} M_J \text{ yr}^{-1}$ as the estimate for the accretion rate.

With $A_V = 2.3$ mag, we corrected the $H\alpha$ luminosity to $L_{H\alpha} = 1.4 \times 10^{-5} L_\odot$ (see Figure 9), which is very similar to the (extinction-corrected) $L_{H\alpha} = 2.0 \times 10^{-5} L_\odot$ from Zhou et al. (2014) although the authors assumed $A_V = 1.5$ mag, while their flux is a factor 6 smaller. The authors estimated a comparable accretion rate of $\dot{M} = 10^{-6.3} M_J \text{ yr}^{-1}$, which is interesting given the different modeling approaches that were followed. Wu et al. (2017b) reported a smaller $H\alpha$ luminosity, $L_{H\alpha} = 10^{-5.9}-10^{-5.4} L_\odot$, but using $A_V = 0.4$ mag so recalibrating to $A_V = 2.3$ mag would make the luminosity comparable to our finding. The authors estimated $\dot{M} \approx 10^{-9}-10^{-8} M_J \text{ yr}^{-1}$, using an empirical relation for T Tauri stars, which is about 2 orders of magnitude smaller than the value derived in this work. Roughly speaking, the $H\alpha$ luminosity appears somewhat stationary, in contrast to the variable $Pa\beta$ line. Multiepoch observations are required to investigate the accretion processes and variability in more detail, in which case it is recommended to target $H\alpha$ and $Pa\beta$ during the same night to mitigate potential biases with a multiline interpretation.

A comparison of the upper limit on the line ratio, $F_{H\beta}/F_{H\alpha} < 0.17$, with the accretion model yields $v_0 \gtrsim 190 \text{ km s}^{-1}$ and $11.2 \lesssim \log n_0 / \text{cm}^{-3} \lesssim 11.8$. This agrees well with the constraint from the $L_{H\alpha}$ that was shown in Figure 8 (i.e., the white contour at high $v_0$). The $n_0$ estimate from $L_{H\alpha}$ is inversely proportional to $f_{\text{fill}}$ because a smaller $f_{\text{fill}}$ means a more concentrated accretion flow and a higher density at the shock. When $f_{\text{fill}} < 0.6$, the $n_0$ derived from $L_{H\alpha}$ is too dense to
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Directly imaged, planetary and (sub)stellar companions with a measured Hα luminosity and H–M' color. The red markers and teal markers are objects for which the Hα emission is expected to be associated with accretion and chromospheric activity, respectively. Model predictions (see main text for details) are shown with gray dashed lines and marked at 4, 5, 6, 8, 10, 20, and 40 M<sub>J</sub> (from right to left). The sizes of the squared markers are all shown on the same logarithmic mass scale. For GQ Lup B and PDS 70 b, we show also the color–Hα characteristics after dereddening with the A<sub>V</sub> = 2.3 mag and A<sub>V</sub> = 6.1 mag, respectively.</figcaption></figure>

reproduce the low $F_{H\alpha}/F_{H\alpha}$. Instead, the combined constraints from $L_{H\alpha}$ and the line ratio point to a filling factor of $f_{fill} > 0.6$. In contrast, $f_{fill} \sim 0.01$ had been derived for PDS 70 b (Hashimoto et al. 2020) and 2MASS 0103 (AB)b (Eriksson et al. 2020) with the same analysis and also using MUSE data. Modeling of UV continuum for low-mass objects typically yields $f_{fill} \lesssim 0.01$ (e.g., Herczeg & Hillenbrand 2008). Therefore, the large $f_{fill}$ for GQ Lup B suggests a spherical-like accretion rather than concentrated accretion flow such as magnetospheric accretion. Possibly, the last crossing of GQ Lup B with the CSD caused an increased inflow of material, which is being accreted with an approximate spherical geometry. This may also be the origin of the extinction if dust is entailed with the accretion flow.

## 5.5. Hα-Color Characteristics of Directly Imaged Planets and Low-mass Companions

There is an increasing number of directly imaged planets and brown dwarfs with detected hydrogen emission lines. We can therefore start to empirically compare various types of objects to reveal possible trends in line luminosities, which may point to differences and commonalities in their formation mechanisms and accretion geometries. We therefore compiled a list of Hα fluxes of directly imaged objects that also have M' measurements. The H–M' color is a useful diagnostic for the temperature/mass of a companion but is also affected by the presence of a circumplanetary/substellar disk.

Figure 9 shows $L_{H\alpha}$ versus H–M'. Here, we have adopted the available magnitudes and line fluxes from species for PDS 70 b (Hashimoto et al. 2020; Stolker et al. 2020a), GSC 06214 B (Ireland et al. 2011; Bailey et al. 2013; Zhou et al. 2014), PZ Tel B (Biller et al. 2010; Musso Barcucci et al. 2019; Stolker et al. 2020b), HD 142527 B (Biller et al. 2012; Lacour et al. 2016; Cugno et al. 2019), and GQ Lup B (this work). Extinction effects are typically nonnegligible during formation, especially at optical wavelengths, but can also be quite uncertain. We estimated $A_V \approx 2.3$ mag for GQ Lup B (see Section 4.5.2) but assuming that the absolute flux-calibration of the MUSE spectrum is accurate and that the object is not variable. Therefore, we show also the case of $A_V = 0.0$ mag for GQ Lup B in Figure 9. Similarly, for PDS 70 b, we adopted

$A_V = 6.1$ mag from Wang et al. (2021), which was shown to provide a good fit with BT-Settl spectra, but also show the $A_V = 0.0$ mag case for reference. For GSC 06214 B, we adopted $L_{H\alpha} = 9 \times 10^{-6} L_\odot$ from Zhou et al. (2014), which had been corrected for the assumed extinction of $A_V = 0.2$ mag.

With the extinction-corrected line luminosities, differences in the origin of the Hα emission become clear from Figure 9. Late-type field objects are known to be chromospherically active, which produces somewhat low Hα fluxes (e.g., Hawley et al. 1996; Pineda et al. 2017). The detected Hα emission from PZ Tel B ($L_{H\alpha} \sim 10^{-7} L_\odot$) is therefore expected to have a chromospheric origin instead of accretion, which seems also likely given the age of 24 Myr (Musso Barcucci et al. 2019). That interpretation is consistent with the lack of reddening in the M' band (see Stolker et al. 2020b and Figure 9). The fraction of active M and L dwarfs peaks at the M/L transition (e.g., West et al. 2004; Schmidt et al. 2015), so a small part of the Hα flux of GQ Lup B may have a chromospheric origin though. The higher Hα luminosities are instead expected to be associated with accretion. Indeed, these directly imaged objects in Figure 9 are all young (i.e., $\lesssim 10$ Myr) with their host stars showing signs of accretion from a CSD. In the case of GQ Lup B, the retrieved $T_{eff} = 2700$ K and $R_p = 3.8 R_J$ correspond to $\log(L_{bol}/L_\odot) = -2.2$, such that $\log(L_{H\alpha}/L_{bol}) = -2.7$. The latter is about 1 order of magnitude larger than typical luminosity ratios of late M-type field objects (West et al. 2004) while ignoring differences in age or mass. Indeed, Barrado y Navascués & Martín (2003) derived a saturation level for chromospheric activity of $\log(L_{H\alpha}/L_{bol}) = -3.3$ at early M-type T Tauri stars. The detection of the Ca IIλ8662 line (see Section 4.4) further strengthens the interpretation, since this line is almost exclusively found in accreting objects, while λ8498 and λ8542 may also have a chromospheric origin (Muzerolle et al. 2003; Mohanty et al. 2005).

In addition to the empirical data, Figure 9 shows model predictions that are derived by combining the AMES-Dusty isochrones and spectra (Chabrier et al. 2000; Allard et al. 2001) with the hydrogen line predictions from the accretion model by Aoyama et al. (2018). Specifically, we extracted an isochrone at 5 Myr and interpolated $T_{eff}$, $M_p$, and $R_p$. From this, we computed synthetic photometry from the model spectra and, together with considered values for $\dot{M}$, we calculated the preshock velocity and density. The latter were then used to interpolate the grid of Hα fluxes that were also used in Section 5.4. The model predictions in Figure 9 show that lower mass objects are redder due to their lower temperature and emit less Hα (assuming $f_{fill} = 1$) since their radii and freefall velocities are smaller. The accretion rate does in particular impact the Hα flux with approximately a linear scaling between $L_{H\alpha}$ and $\dot{M}$ (see Aoyama et al. 2021 for a detailed analysis of the $\dot{M}-L_{H\alpha}$ relation).

The empirical data seem to follow the expectations from the numerical predictions. Lower-mass objects (see symbol sizes in Figure 9) are typically redder in H–M' although dusty material in the vicinity could also redden a more massive object. In that regard, the color of the M-dwarf companion HD 142527 B is expected to be reddened by a circumsecondary disk (Lacour et al. 2016; Christiaens et al. 2018), similar to the detected excess flux from GQ Lup B. For PDS 70 b, it remains uncertain if the protoplanet is accreting from a CPD, in contrast to PDS 70 c, which has been detected with ALMA (Benisty et al. 2021). A slight, but not significant, excess flux was detected in M' by Stolker et al. (2020a), from which an upper
limit of ~256 K and ~245 $R_J$ was derived for potential emission from a CPD (see also Christiaens et al. 2019). The SED of PDS 70 b appears to be reddened as well and shows muted absorption features. Extinction effects are therefore expected to be nonnegligible although the exact value remains under debate (Cugno et al. 2021; Wang et al. 2021). In general, it is also not clear if the extinction inferred from an SED is affecting the accretion luminosity in a similar way. A better understanding of the extinction properties will be key for deriving accurate line luminosities. Indeed, Figure 9 shows that a correction for extinction makes a large impact. The characteristics of GQ Lup B become similar to GSC 06214 B, and PDS 70 b is consistent with an even higher accretion rate. In fact, the extinction-corrected line luminosities from all considered, accreting, directly imaged objects are comparable within 1 order of magnitude even though they span over 2 orders of magnitude in mass.

Compared to PDS 70 b, GQ Lup B orbits further away from its primary star and is strongly misaligned with respect to the CSD (see Section 4.6). Hydrodynamical simulations have shown that wide-orbit, forming planets might be the origin of spiral structures that are seen in CSDs (see e.g., Dong et al. 2016). Since GQ Lup B appears to drive a spiral arm, it remains to be further investigated whether the companion formed from the CSD or as a binary system by fragmentation from the same molecular cloud. The bottom-up formation timescale might be too long at ~120 au but Stamatellos & Whitworth (2009) showed that large-separation brown dwarfs can form through disk fragmentation at an early phase. As already mentioned, some circumstellar gas may get channeled toward GQ Lup B and its disk, in particular when the companion crosses the CSD. ALMA observations by MacGregor et al. (2017) resolved CO emission at the projected distance of GQ Lup B with a resolution of ~0.3". At higher resolution, ALMA may reveal a possible CO counterpart of the spiral arm, as well as kinematical perturbations in the vicinity of GQ Lup B, providing further insight into the most recent crossing with the CSD.

## 6. Conclusions

We have reported on the optical-to-MIR characterization of the directly imaged, substellar companion GQ Lup B. The optical spectrum is best matched with an M9 spectral type and contains several low-gravity indicators. We detected strong $H\alpha$ emission and derived an upper limit for $H\beta$. Analysis of the full SED, including the JHK-band spectra from Seifahrt et al. (2007), yielded $T_{\text{eff}} \approx 2700$ K, $\log g = 3.5-4.0$, $R_p \approx 3.8 R_J$, and a higher extinction, $A_V \approx 2.3$ mag, than previously estimated. The H-band flux and atmospheric parameters are consistent with $M_p \approx 30 M_J$ at an age of 3 Myr.

The H-M' color is ≈1 mag redder than mid/late field dwarfs due to an overluminosity at MIR wavelengths. A blackbody model provides a good fit to the 3–5 µm fluxes, from which we derived $T_{\text{disk}} = 461 \pm 2$ K and $R_{\text{disk}} = 65 \pm 1 R_J$, and is consistent with the ALMA upper limits. We used 15 yr of astrometric measurements to constrain the mutual inclination between the orbital plane of GQ Lup B and the CSD to $84 \pm 9$ deg, and we tentatively showed that a low-eccentricity orbit is favored.

Our main conclusions are the following:

(i) The MIR excess in the SED of GQ Lup B is caused by thermal emission from small dust grains in a protolunar disk.

(ii) The $H\alpha$ line emission in combination with the MIR excess indicates that GQ Lup B is actively accreting from its disk. In addition, material may get delivered directly from the CSD, which could explain the $f_{\text{fill}} \approx 1$ (i.e., a spherical accretion geometry) that we derived from the $H\beta/H\alpha$ upper limit.

(iii) Analysis of the optical and NIR spectra suggests that the SED is reddened by $A_V \approx 2.3$ mag. The extinction is expected to occur in the vicinity of GQ Lup B, for example, in a moderately inclined disk or by infalling material from the CSD, which somewhat enshrouds GQ Lup B and its disk.

(iv) The disk of GQ Lup B is cool compared to its atmosphere, as well as dust sublimation temperatures. We speculate that the disk is in a transitional stage, in which the inner regions have been cleared by the formation of satellites, while the remaining gas is being accreted onto the atmosphere with $\dot{M} \approx 10^{-6.5} M_J \text{ yr}^{-1}$. A depletion of the pebble reservoir would explain the nondetections with ALMA.

(v) The large mutual inclination between the orbit and CSD could imply a tumultuous dynamical history (e.g., driven by the second, large separation companion in the system), which shares similarities with strong spin-orbit misalignments of hot-Jupiters. Alternatively, GQ Lup B may not have formed from the CSD but collapsed directly from the same parent cloud as GQ Lup.

(vi) Multiline measurements during a single observing night are required to mitigate potential biases due to variability in the accretion processes and to take further advantage of the line-ratio diagnostics.

We have presented the first detailed characterization of a disk around a directly imaged, low-mass companion, which highlights the valuable window for ground-based observations at 4–5 µm. Measurements of a larger sample could reveal trends as function of companion mass and formation environment. At millimeter wavelengths, we predict that a factor ≳10 increase in sensitivity is required to detect the disk component of GQ Lup B. This will be challenging with the current capabilities of ALMA but might be in reach with the Next Generation Very Large Array (Rab et al. 2019; Wu et al. 2020). Finally, the nearing launch of JWST provides the exciting opportunity to reveal the spectral appearance of the protolunar disk at medium resolution and will guide detailed modeling efforts of its structure and dust properties.

We want to thank Andreas Seifahrt and Michael McElwain for sharing their near-infrared spectra of GQ Lup B. T.S. acknowledges the support from the Netherlands Organisation for Scientific Research (NWO) through grant VI.Veni.202.230. This work was performed using the ALICE compute resources provided by Leiden University. S.Y.H. acknowledges the support provided by NASA through the NASA Hubble Fellowship grant \#HST-HF2-51436.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (KU 2849/7-1 and MA 9185/1-1), and from the Swiss National Science Foundation under grant BSSGIO_155816 "PlanetsInTime." S.P.Q. and G.C. thanks the Swiss National Science Foundation for financial support under grant number 200021_169131. Parts of this
work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. J.B. acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through research grants UIDB/04434/2020 and UIDP/04434/2020 and work contract 2020.03379.CEECIND.

*Facilities:* VLT:Antu(NACO), VLT: Yepun(MUSE), VLT: Melipal(SPHERE).

*Software:* species (Stolker et al. 2020b), orbitize! (Blunt et al. 2020), PynPoint (Stolker et al. 2019), IRDAP

(van Holstein et al. 2020), diskmap (Stolker et al. 2016), PyHammer (Kesseli et al. 2017).

## Appendix
### Posterior Distributions

Figures 10 and 11 show the posterior distributions from the SED and orbital analysis, respectively.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Posterior distributions from fitting the MUSE and SINFONI spectra, the photometric fluxes at 3–5 µm, and the ALMA upper limits with a blackbody model while fixing the atmospheric parameters. The 1D marginalized distributions are shown in the main-diagonal panels, and the 2D projections of the posterior samples are shown in the off-diagonal panels. The plot was created with corner.py (Foreman-Mackey 2016). See Section 4.5 for further details on the SED analysis.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Posterior distributions of the orbital elements of GQ Lup B. The 1D marginalized distributions are shown in the main-diagonal panels, and the 2D projections of the posterior samples are shown in the off-diagonal panels. The plot was created with corner.py (Foreman-Mackey 2016). See Section 4.6 for further details on the orbit analysis.</figcaption></figure>

ORCID iDs

Tomas Stolker https://orcid.org/0000-0002-5823-3072

Sebastiaan Y. Haffert https://orcid.org/0000-0001-5130-9153

Aurora Y. Kesseli https://orcid.org/0000-0002-3239-5989

Rob G. van Holstein https://orcid.org/0000-0003-1520-8405

Yuhiko Aoyama https://orcid.org/0000-0003-0568-9225

Jarle Brinchmann https://orcid.org/0000-0003-4359-8797

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

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

Gabriel-Dominique Marleau https://orcid.org/0000-0002-2919-7500

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

Julien Milli https://orcid.org/0000-0001-9325-2511

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

Ignas A. G. Snellen https://orcid.org/0000-0003-1624-3667

Kamen O. Todorov https://orcid.org/0000-0002-9276-8118
References

Alcalá, J. M., Majidi, F. Z., Desidera, S., et al. 2020, A&A, 635, L1

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

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

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

Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948

Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84

Aoyama, Y., Marleau, G.-D., Ikoma, M., & Mordasini, C. 2021, ApJL, 917, L30

Aoyama, Y., Marleau, G.-D., Mordasini, C., & Ikoma, M. 2020, arXiv:2011. 06608

Arsenault, R., Madec, P. Y., Hubin, N., et al. 2008, Proc. SPIE, 7015, 701524

Bacon, R., Accardo, M., Adjali, L., et al. 2010, Proc. SPIE, 7735, 773508

Bailey, V., Hinz, P. M., Currie, T., et al. 2013, ApJ, 767, 31

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Barrado y Navascués, D., & Martín, E. L. 2003, AJ, 126, 2997

Bate, M. R. 2018, MNRAS, 475, 5618

Batygin, K. 2012, Natur, 491, 418

Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJL, 916, L2

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Biller, B., Lacour, S., Juhász, A., & Mann, A. W. 2014, ApJL, 753, L38

Biller, B. A., Liu, M. C., Wahhaj, Z., et al. 2010, ApJL, 720, L82

Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89

Bohlin, R. C. 2007, in ASP Conf. Ser. 364, The Future of Photometric, Spectrophotometric and Polarimetric Standardization, ed. C. Sterken (San Francisco, CA: ASP), 315

Bonnefoy, M., Currie, T., Marleau, G.-D., et al. 2014, A&A, 562, A111

Bowler, B. P. 2016, PASP, 128, 102001

Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63

Bowler, B. P., Liu, M. C., Kraus, A. L., & Mann, A. W. 2014, ApJ, 784, 65

Bowler, B. P., Liu, M. C., Kraus, A. L., Mann, A. W., & Ireland, M. J. 2011, ApJ, 743, 148

Broeg, C., Schmidt, T. O. B., Guenther, E., et al. 2007, A&A, 468, 1039

Buchner, J. 2016, Stat. Comp., 26, 383

Buchner, J. 2019, PASP, 131, 108005

Buchner, J. 2021, JOSS, 6, 3001

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

Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464

Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, L9

Cheetham, A. C., Samland, M., Brems, S. S., et al. 2019, A&A, 622, A80

Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJL, 877, L33

Christiaens, V., Casassus, S., Absil, O., et al. 2018, A&A, 617, A37

Cugno, G., Patapis, P., Stolker, T., et al. 2021, A&A, 653, A12

Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A&A, 622, A156

Currie, T., Burrows, A., Madhusudhan, N., et al. 2013, ApJ, 776, 15

Currie, T., Daemgen, S., Debes, J., et al. 2014, ApJL, 780, L30

Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., & Plavchan, P. 2012, ApJL, 755, L34

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

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

Daemgen, S., Todorov, K., Silva, J., et al. 2017, A&A, 601, A65

Dai, Y., Wilner, D. J., Andrews, S. M., & Ohashi, N. 2010, AJ, 139, 626

de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, Proc. SPIE, 7014, 70143L

Donati, J. F., Gregory, S. G., Alencar, S. H. P., et al. 2012, MNRAS, 425, 2948

Dong, R., Fung, J., & Chiang, E. 2016, ApJ, 826, 75

Dupuy, T. J., & Kraus, A. L. 2013, Sci, 341, 1492

Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19

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

Ferrer-Chávez, R., Wang, J. J., & Blunt, S. 2021, AJ, 161, 241

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Vousden, W., Price-Whelan, A., et al. 2016, corner.py: corner.py v2.0.0, Zenodo, doi:10.21105/joss.00024

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

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJL, 739, L41

Garcia, E. V., Currie, T., Guyon, O., et al. 2017, ApJ, 834, 162

Ginski, C., Schmidt, T. O. B., Mugrauer, M., et al. 2014, MNRAS, 444, 2280

Gravity Collaboration, Lacour, S., Nowak, M., et al. 2019, A&A, 623, L11

Greenstein, J. L. 1969, ApJ, 158, 281

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

Haffert, S. Y., van Holstein, R. G., Ginski, C., et al. 2020, A&A, 640, L12

Hashimoto, J., Aoyama, Y., Konishi, M., et al. 2020, AJ, 159, 222

Hawley, S. L., Gizis, J. E., & Reid, I. N. 1996, AJ, 112, 2799

Heller, R., & Pudritz, R. 2015, ApJ, 806, 181

Herbig, G. H. 1977, ApJ, 214, 747

Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594

Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97

Hunziker, S., Quanz, S. P., Amara, A., & Meyer, M. R. 2018, A&A, 611, A23

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Ireland, M. J., Kraus, A., Martinache, F., Law, N., & Hillenbrand, L. A. 2011, ApJ, 726, L13

Janson, M., Brandner, W., Henning, T., & Zinnecker, H. 2006, A&A, 453, 609

Jayawardhana, R., Ardila, D. R., Stelzer, B., & Haisch, K. E., Jr. 2003a, AJ, 126, 1515

Jayawardhana, R., Mohanty, S., & Basri, G. 2003b, ApJ, 592, 282

Jensen, E. L. N., & Akeson, R. 2014, Natur, 511, 567

Kammerer, R., Lacour, S., Stolker, T., et al. 2021, A&A, 652, A57

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Kesseli, A. Y., West, A. A., Veyette, M., et al. 2017, ApJS, 230, 16

Kraus, A. L., Ireland, M. J., Cieza, L. A., et al. 2003, ApJ, 781, 20

Lacour, S., Biller, B., Cheetham, A., et al. 2016, A&A, 590, A90

Lada, C. J., Muench, A. A., Haisch, Karl E., et al. 2000, AJ, 120, 3162

Lavigne, J.-F., Doyon, R., Lafrenière, D., Marois, C., & Barman, T. 2009, ApJ, 704, 1098

Lazzoni, C., Gratton, R., Alcalá, J. M., et al. 2020, A&A, 635, L11

Lee, J.-M., Heng, K., & Irwin, P. G. J. 2013, ApJ, 778, 97

Lenzen, R., Hartung, M., Brandner, W., et al. 2003, Proc. SPIE, 4841, 944

Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96

Liu, M. C., Najita, J., & Tokunaga, A. T. 2003, ApJ, 585, 372

Long, D. E., Zhang, K., Teague, R., & Bergin, E. A. 2020, ApJL, 895, L46

Lovelace, R. V. E., Covey, K. R., & Lloyd, J. P. 2011, AJ, 141, 51

MacGregor, M. A., Wilner, D. J., Czekala, I., et al. 2017, ApJ, 835, 17

Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJL, 794, L12

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, Proc. SPIE, 9908, 990834

Marois, C., Macintosh, B., & Barman, T. 2007, ApJL, 654, L151

Martin, E. L., Rebolo, R., & Zapatero-Osorio, M. R. 1996, ApJ, 469, 706

McClure, M. K., Manoj, P., Calvet, N., et al. 2012, ApJL, 759, L10

McElwain, M. W., Metchev, S. A., Larkin, J. E., et al. 2007, ApJ, 656, 505

McGovern, M. R., Kirkpatrick, J. D., McLean, I. S., et al. 2004, ApJ, 600, 1020

Milli, J., Hibon, P., Christiaens, V., et al. 2017, A&A, 597, L2

Mohanty, S., Jayawardhana, R., & Basri, G. 2005, ApJ, 626, 498

Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131

Monnier, J. D., & Millan-Gabet, R. 2002, ApJ, 579, 694

Muench, A. A., Alves, J., Lada, C. J., & Lada, E. A. 2001, ApJL, 558, L51

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2

Musso Barcucci, A., Cugno, G., Launhardt, R., et al. 2019, A&A, 631, A84

Muzerolle, J., Hillenbrand, L., Calvet, N., Briceño, C., & Hartmann, L. 2003, ApJ, 592, 266

Nagasawa, M., & Ida, S. 2011, ApJ, 742, 72

Neuhäuser, R., Guenther, E. W., Wuchterl, G., et al. 2005, A&A, 435, L13

Neuhäuser, R., Mugrauer, M., Seifahrt, A., Schmidt, T. O. B., & Vogt, N. 2008, A&A, 484, 281

Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13

Noll, S., Kausch, W., Barden, M., et al. 2012, A&A, 543, A92

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

Oppenheimer, B. R., & Hinkley, S. 2009, ARA&A, 47, 253

Pérez, S., Marino, S., Casassus, S., et al. 2019, MNRAS, 488, 1005

Pineda, J. S., Hallinan, G., & Kao, M. M. 2017, ApJ, 846, 75

Rab, C., Kamp, I., Ginski, C., et al. 2019, A&A, 624, A16

Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, 154, 10

Reiners, A., & Christensen, U. R. 2010, A&A, 522, A13

Rousset, G., Lacombe, F., Puget, P., et al. 2003, Proc. SPIE, 4839, 140

Scharf, C., & Menou, K. 2009, ApJL, 693, L113

Schiavon, R. P., Barbuy, B., & Singh, P. D. 1997, ApJ, 484, 499

Schmidt, S. J., Hawley, S. L., West, A. A., et al. 2015, AJ, 149, 158

Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74

Seifahrt, A., Neuhäuser, R., & Hauschildt, P. H. 2007, A&A, 463, 309

Seperuelo Duarte, E., Alencar, S. H. P., Batalha, C., & Lopes, D. 2008, A&A, 489, 349

Shibaike, Y., Okuzumi, S., Sasaki, T., & Ida, S. 2017, ApJ, 846, 81

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392, 413
Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59

Stolker, T., Dominik, C., Min, M., et al. 2016, A&A, 596, A70

Stolker, T., Marleau, G.-D., Cugno, G., et al. 2020a, A&A, 644, A13

Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020b, A&A, 635, A182

Ströbele, S., La Penna, P., Arsenault, R., et al. 2012, Proc. SPIE, 8447, 844737

Tuthill, P. G., Monnier, J. D., & Danchi, W. C. 2001, Natur, 409, 1012

Uyama, T., Xie, C., Aoyama, Y., et al. 2021, AJ, 162, 214

van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020, A&A, 633, A64

van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, Proc. SPIE, 10400, 1040015

van Holstein, R. G., Stolker, T., Jensen-Clem, R., et al. 2021, A&A, 647, A21

Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72

Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919

Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263

Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148

Weilbacher, P. M., Palsa, R., Streicher, O., et al. 2020, A&A, 641, A28

West, A. A., Hawley, S. L., Walkowicz, L. M., et al. 2004, AJ, 128, 426

White, R. J., & Basri, G. 2003, ApJ, 582, 1109

Wu, Y.-L., Bowler, B. P., Sheehan, P. D., et al. 2020, AJ, 159, 229

Wu, Y.-L., Close, L. M., Eisner, J. A., & Sheehan, P. D. 2017a, AJ, 154, 234

Wu, Y.-L., Close, L. M., Males, J. R., et al. 2015, ApJL, 807, L13

Wu, Y.-L., Sheehan, P. D., Males, J. R., et al. 2017b, ApJ, 836, 223

Zapatero Osorio, M. R., Caballero, J. A., Béjar, V. J. S., et al. 2007, A&A, 472, L9

Zhou, Y., Herczeg, G. J., Kraus, A. L., Metchev, S., & Cruz, K. L. 2014, ApJL, 783, L17

Zhu, Z. 2015, ApJ, 799, 16