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

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α; λ₀ = 4.05 μm, Δλ = 0.06 μm ¹⁷ ) and M' (λ₀ = 4.78 μm, Δλ = 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 078 ± 004 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, 085 ± 004, 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.

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 mas)² spaxel⁻¹ 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 032 ± 004 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 064 and a coherence time of 9 ms. The calibration target was however observed at a smaller airmass of 1.07.

3.1.1. Image Processing

3.1.2. Relative Calibration

3.1.3. Absolute Flux Calibration

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 $\mathrm{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 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.

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 ≈07.

We reduced the data with IRDAP ²⁰ (van Holstein et al. 2017, 2020) and measured the relative position with PynPoint (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.

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.

Zoom In Zoom Out Reset image size

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.

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; 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.

Zoom In Zoom Out Reset image size

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 ± 0.13 mag so the H − M' color is 2.04 ± 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 ≳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 (≲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 ≈ 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σ.

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 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 ≳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.

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.

Zoom In Zoom Out Reset image size

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.

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σ upper limit of 2.7 × 10⁻¹⁹ W m⁻² 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⁻¹) 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⁻¹; 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/oversubtraction 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 ± 0.2 km s⁻¹ 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.1. Data and Model Preparation

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, $\mathrm{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 χ² statistic but is not expected to follow a χ² 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, $\mathrm{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 $\mathrm{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 $\mathrm{log}g=3.5\mbox{--}4.0$ when fixing T_eff = 2700 K. Specifically, the K i doublet near 1.25 μm is best matched with $\mathrm{log}g=4.0$, but the predicted K i and Na i lines at ∼0.77 μm and ∼0.82 μm, respectively, are a bit too strong and better matched with $\mathrm{log}g=3.5$. Similarly, with $\mathrm{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 $\mathrm{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.

Table 4. Spectral Analysis of GQ Lup B with Atmospheric Models

	Teff	$\mathrm{log}g$ a	Rp	aJ b	aH b	AV	Tdisk	Rdisk	ΔGc
	(K)		(RJ)				(K)	(RJ)
Parameter range	2000–3500	2.5–5.5	1–5	0.1–5	0.1–5	0–5	10–2000	5–1000
MUSE spectrum
BT-Settl	2500	4.0	2.85	⋯	⋯	⋯	⋯	⋯	0
BT-Settl + ext.	2700	4.0	4.13	⋯	⋯	2.7	⋯	⋯	−156
SINFONI spectra
BT-Settl	2600	4.0	3.55	1.84	1.19	⋯	⋯	⋯	0
BT-Settl + ext.	2600	4.0	3.55	1.84	1.19	0.0	⋯	⋯	0
MUSE + SINFONI spectra
BT-Settl	2350	4.0	3.93	1.68	1.06	⋯	⋯	⋯	0
BT-Settl + ext.	2700	4.0	3.77	1.32	1.03	2.3	⋯	⋯	−11587
MUSE + SINFONI spectra + 3–5 μm photometry + ALMA
BT-Settl + ext. + disk d	2700	4.0	3.77	1.32	1.03	2.3	461 ± 2	65 ± 1	⋯

Notes.

^a The surface gravity, $\mathrm{log}g$, was estimated from the gravity-sensitive alkali doublets in the MUSE and J-band spectra. ^b Scaling parameter to account for potential inaccuracies in the absolute flux calibration. ^c The relative goodness-of-fit is calculated for each combination of spectra separately. ^d The atmospheric, scaling, and extinction parameters were fixed to the best-fit values of fitting the MUSE and SINFONI spectra.

Download table as:  ASCIITypeset image

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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, $\mathrm{log}g=3.5\mbox{--}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 ≈ 2700 K, $\mathrm{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.

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.

Zoom In Zoom Out Reset image size

4.5.3. Thermal Emission from a Dusty Disk

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! ²¹ (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_sys = − 2.8 ± 0.2 km s⁻¹ from Schwarz et al. (2016) and the RV of GQ Lup B, 2.0 ± 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 ± 0.03 mas; Gaia Collaboration et al. 2016, 2021) and stellar mass (1.03 ± 0.05 M_⊙; MacGregor et al. 2017).

Zoom In Zoom Out Reset image size

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 ${\rm{\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 ²² (Stolker et al. 2016) to apply an r²-scaling that accounts for both the orientation of the midplane, for 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.

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σ 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_disk = 461 ± 2 K and R_disk = 65 ± 1 R_J (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.

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_disk = 461 ± 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_eff ≈ 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_disk = 65 ± 1 R_J, may provide further evidence for a cavity since it is large compared to R_p ≈ 3.8 R_J from the atmosphere. While R_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}_{{\rm{s}}}=\tfrac{1}{2}\sqrt{{Q}_{{\rm{R}}}}{({T}_{\mathrm{eff}}/{T}_{{\rm{s}}})}^{2}{R}_{{\rm{p}}}$, where Q_R = Q_abs(T_eff)/Q_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_eff and R_p of GQ Lup B, we obtain R_s ≈ 6 R_J. 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_disk ≪ T_s. Instead, T_disk lies closer to the freeze-out location of H₂O, which would be located a bit further outward. Therefore, the H₂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_disk ≈ 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σ 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}_{{\rm{J}}}\,{\mathrm{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α 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}_{{\rm{J}}}\,{\mathrm{yr}}^{-1}$, we calculated a truncation radius of R_T ≈ 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_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 ≈ R_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.

The orbital analysis in Section 4.6 constrained the mutual inclination between the orbit of GQ Lup B and the CSD to 84 ± 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 ± 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 $\mathrm{log}P/\mathrm{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 ≈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 ≈ 2.3 mag (i.e., an optical depth of τ ≈ 2) from the broadband SED, which is larger than what has been inferred for the star (e.g., 0.4 ± 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.

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α emission with ground- and space-based optical photometry (Marois et al. 2007; Zhou et al. 2014; Wu et al. 2017b). The Paβ 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 circumplanetary 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α flux and estimated an upper limit for Hβ, yielding an extinction-corrected, 1σ upper limit on the line ratio of F_Hβ /F_Hα < 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₀, and hydrogen number density, n₀. The line fluxes from the model are converted into ${L}_{{{\rm{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_fill = 1. We can also map (v₀, n₀) to M_p (i.e., assuming that v₀ 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}_{{{\rm{H}}}_{\alpha }}$ of GQ Lup B, which is shown in Figure 8. Later on, we will estimate a filling factor of f_fill > 0.6 from the line ratio, which would imply that the actual $\dot{M}$ could be smaller by a factor ∼2.

Zoom In Zoom Out Reset image size

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α line is not resolved at the resolving power of MUSE (i.e., ≈120 km s⁻¹ in Hα). If the gas reaches the planet surface with the freefall velocity, this would point to a low mass of GQ Lup B (≲15 M_J; see Figure 8), possibly making GQ Lup B a PMC, with an accretion rate of $\dot{M}\approx {10}^{-6.5}\mbox{--}{10}^{-6}\,{M}_{{\rm{J}}}\,{\mathrm{yr}}^{-1}$. However, as noted previously, there was possibly some oversubtraction of the stellar residuals at Hα, which could somewhat bias the line flux and width. Indeed, we measured a broader line for Paβ (FWHM ≈ 237 km s⁻¹), 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}_{{\rm{J}}}\,{\mathrm{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 ≈ 30 M_J. Therefore, the inferred mass from the Paβ line and the H-band flux are roughly in agreement so we adopt $\dot{M}\approx {10}^{-6.5}\,{M}_{{\rm{J}}}\,{\mathrm{yr}}^{-1}$ as the estimate for the accretion rate.

With A_V = 2.3 mag, we corrected the Hα luminosity to ${L}_{{{\rm{H}}}_{\alpha }}=1.4\times {10}^{-5}\,{L}_{\odot }$ (see Figure 9), which is very similar to the (extinction-corrected) ${L}_{{{\rm{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}_{{\rm{J}}}\,{\mathrm{yr}}^{-1}$, which is interesting given the different modeling approaches that were followed. Wu et al. (2017b) reported a smaller Hα luminosity, ${L}_{{{\rm{H}}}_{\alpha }}={10}^{-5.9}\mbox{--}{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}\mbox{--}{10}^{-8}\,{M}_{{\rm{J}}}\,{\mathrm{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α luminosity appears somewhat stationary, in contrast to the variable Paβ line. Multiepoch observations are required to investigate the accretion processes and variability in more detail, in which case it is recommended to target Hα and Paβ during the same night to mitigate potential biases with a multiline interpretation.

Zoom In Zoom Out Reset image size

A comparison of the upper limit on the line ratio, F_Hβ /F_Hα < 0.17, with the accretion model yields v₀ ≳ 190 km s⁻¹ and $11.2\lesssim \mathrm{log}{n}_{0}/{\mathrm{cm}}^{-3}\lesssim 11.8$. This agrees well with the constraint from the L_Hα that was shown in Figure 8 (i.e., the white contour at high v₀). The n₀ estimate from L_Hα is inversely proportional to f_fill because a smaller f_fill means a more concentrated accretion flow and a higher density at the shock. When f_fill < 0.6, the n₀ derived from L_Hα is too dense to reproduce the low F_Hβ /F_Hα . Instead, the combined constraints from L_Hα and the line ratio point to a filling factor of f_fill > 0.6. In contrast, f_fill ∼ 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 ≲ 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.

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α 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 ≈ 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α = 9 × 10⁻⁶ L_⊙ 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}_{{{\rm{H}}}_{\alpha }}\sim {10}^{-7}$ L_⊙) 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., ≲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 $\mathrm{log}({L}_{\mathrm{bol}}/{L}_{\odot })=-2.2$, such that $\mathrm{log}({L}_{{{\rm{H}}}_{\alpha }}/{L}_{\mathrm{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 $\mathrm{log}({L}_{{{\rm{H}}}_{\alpha }}/{L}_{\mathrm{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α and $\dot{M}$ (see Aoyama et al. 2021 for a detailed analysis of the $\dot{M}$–L_Hα 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 ∼03. 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.