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

HAL
open science

# Gemini Planet Imager Spectroscopy of the Dusty Substellar Companion HD 206893 B

K. Ward-Duong, J. Patience, K. Follette, R. J. de Rosa, J. Rameau, M.

Marley, D. Saumon, E. L. Nielsen, A. Rajan, A. Z. Greenbaum, et al.

► To cite this version:

K. Ward-Duong, J. Patience, K. Follette, R. J. de Rosa, J. Rameau, et al. Gemini Planet Imager Spectroscopy of the Dusty Substellar Companion HD 206893 B. The Astronomical Journal, 2021, 161, 10.3847/1538-3881/abc263. insu-03705332

HAL Id: insu-03705332

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

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>

# Gemini Planet Imager Spectroscopy of the Dusty Substellar Companion HD 206893 B

K. Ward-Duong<sup>1,2</sup>①, J. Patience<sup>2</sup>, K. Follette<sup>3</sup>①, R. J. De Rosa<sup>4,5</sup>①, J. Rameau<sup>6,7</sup>①, M. Marley<sup>8</sup>①, D. Saumon<sup>9</sup>①,
E. L. Nielsen<sup>4</sup>①, A. Rajan<sup>10</sup>①, A. Z. Greenbaum<sup>11</sup>①, J. Lee<sup>12</sup>, J. J. Wang<sup>13,14,40</sup>①, I. Czekala<sup>4,14,41</sup>①, G. Duchêne<sup>6,14</sup>①,
B. Macintosh<sup>4</sup>①, S. Mark Ammons<sup>15</sup>①, V. P. Bailey<sup>16</sup>①, T. Barman<sup>17</sup>①, J. Bulger<sup>18,19</sup>①, C. Chen<sup>10</sup>①, J. Chilcote<sup>4,20</sup>①,
T. Cotten<sup>12</sup>①, R. Doyon<sup>7</sup>, T. M. Esposito<sup>14</sup>①, M. P. Fitzgerald<sup>21</sup>①, B. L. Gerard<sup>22,23</sup>①, S. J. Goodsell<sup>24</sup>①, J. R. Graham<sup>14</sup>,
P. Hibon<sup>5</sup>①, J. Hom<sup>2</sup>①, L.-W. Hung<sup>25</sup>①, P. Ingraham<sup>26</sup>①, P. Kalas<sup>14,27</sup>①, Q. Konopacky<sup>28</sup>①, J. E. Larkin<sup>21</sup>①, J. Maire<sup>28</sup>,
F. Marchis<sup>27</sup>①, C. Marois<sup>23,29</sup>①, S. Metchev<sup>30,31</sup>①, M. A. Millar-Blanchaer<sup>16,42</sup>①, R. Oppenheimer<sup>32</sup>①, D. Palmer<sup>15</sup>①,
M. Perrin<sup>10</sup>①, L. Poyneer<sup>15</sup>, L. Pueyo<sup>10</sup>, F. T. Rantakyrö<sup>33</sup>①, B. Ren<sup>34</sup>①, J.-B. Ruffio<sup>4</sup>①, D. Savransky<sup>35</sup>①,
A. C. Schneider<sup>36,37</sup>①, A. Sivaramakrishnan<sup>10</sup>①, I. Song<sup>12</sup>①, R. Soummer<sup>10</sup>①, M. Tallis<sup>4</sup>, S. Thomas<sup>26</sup>①, J. Kent Wallace<sup>16</sup>①,
S. Wiktorowicz<sup>38</sup>①, and S. Wolff<sup>39</sup>①

## Abstract

We present new near-infrared Gemini Planet Imager (GPI) spectroscopy of HD 206893 B, a substellar companion orbiting within the debris disk of its F5V star. The *J*, *H*, *K1*, and *K2* spectra from GPI demonstrate the extraordinarily red colors of the object, confirming it as the reddest substellar object observed to date. The significant flux increase throughout the infrared presents a challenging atmosphere to model with existing grids. Best-fit values vary from 1200 to 1800 K for effective temperature and from 3.0 to 5.0 for log(g), depending on which individual wavelength band is fit and which model suite is applied. The extreme redness of the companion can be partially reconciled by invoking a high-altitude layer of submicron dust particles, similar to dereddening approaches applied to the peculiar red field L dwarf population. However, reconciling the HD 206893 B spectra with even those of the reddest low-gravity L dwarf spectra still requires the contribution of additional atmospheric

<footnote n="1">¹ Five College Astronomy Department, Amherst College, Amherst, MA 01002, USA; kwduong@stsci.edu</footnote>

<footnote n="2">² School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ 85287, USA</footnote>

<footnote n="3">³ Physics & Astronomy Department, Amherst College, 21 Merrill Science Drive, Amherst, MA 01002, USA</footnote>

<footnote n="4">⁴ Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="5">⁵ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile</footnote>

<footnote n="6">⁶ Univ. Grenoble Alpes/CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="7">⁷ Institut de Recherche sur les Exoplanètes, Département de Physique, Université de Montréal, Montréal QC, H3C 3J7, Canada</footnote>

<footnote n="8">⁸ NASA Ames Research Center, Mountain View, CA 94035, USA</footnote>

<footnote n="9">⁹ Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA</footnote>

<footnote n="10">¹⁰ Space Telescope Science Institute, Baltimore, MD 21218, USA</footnote>

<footnote n="11">¹¹ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="12">¹² Department of Physics and Astronomy, University of Georgia, Athens, GA 30602, USA</footnote>

<footnote n="13">¹³ Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="14">¹⁴ Department of Astronomy, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="15">¹⁵ Lawrence Livermore National Laboratory, Livermore, CA 94551, USA</footnote>

<footnote n="16">¹⁶ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA</footnote>

<footnote n="17">¹⁷ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="18">¹⁸ Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="19">¹⁹ Subaru Telescope, NAOJ, 650 North A'ohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="20">²⁰ Department of Physics, University of Notre Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA</footnote>

<footnote n="21">²¹ Department of Physics & Astronomy, University of California, Los Angeles, CA 90095, USA</footnote>

<footnote n="22">²² Department of Astronomy, UC Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="23">²³ National Research Council of Canada Herzberg, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada</footnote>

<footnote n="24">²⁴ Gemini Observatory, 670 North A'ohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="25">²⁵ Natural Sounds and Night Skies Division, National Park Service, Fort Collins, CO 80525, USA</footnote>

<footnote n="26">²⁶ Large Synoptic Survey Telescope, 950 North Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="27">²⁷ SETI Institute, Carl Sagan Center, 189 Bernardo Avenue, Mountain View, CA 94043, USA</footnote>

<footnote n="28">²⁸ Center for Astrophysics and Space Science, University of California San Diego, La Jolla, CA 92093, USA</footnote>

<footnote n="29">²⁹ University of Victoria, 3800 Finnerty Road, Victoria, BC, V8P 5C2, Canada</footnote>

<footnote n="30">³⁰ Department of Physics and Astronomy, Centre for Planetary Science and Exploration, The University of Western Ontario, London, ON N6A 3K7, Canada</footnote>

<footnote n="31">³¹ Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA</footnote>

<footnote n="32">³² Department of Astrophysics, American Museum of Natural History, New York, NY 10024, USA</footnote>

<footnote n="33">³³ Gemini Observatory, Casilla 603, La Serena, Chile</footnote>

<footnote n="34">³⁴ Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="35">³⁵ Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="36">³⁶ US Naval Observatory, Flagstaff Station, P.O. Box 1149, Flagstaff, AZ 86002, USA</footnote>

<footnote n="37">³⁷ Department of Physics and Astronomy, George Mason University, MS3F3, 4400 University Drive, Fairfax, VA 22030, USA</footnote>

<footnote n="38">³⁸ The Aerospace Corporation, El Segundo, CA 90245, USA</footnote>

<footnote n="39">³⁹ Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote>Received 2019 September 16; revised 2020 October 15; accepted 2020 October 15; published 2020 December 2</footnote>

<footnote n="40">⁴⁰ 51 Pegasi b Fellow.</footnote>

<footnote n="41">⁴¹ NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote n="42">⁴² NASA Hubble Fellow.</footnote>
dust, potentially due to the debris disk environment in which the companion resides. Orbit fitting from 4 yr of astrometric monitoring is consistent with a ~30 yr period, an orbital inclination of 147°, and a semimajor axis of 10 au, well within the estimated disk inner radius of ~50 au. As one of a very few substellar companions imaged interior to a circumstellar disk, the properties of this system offer important dynamical constraints on companion-disk interaction and provide a benchmark for substellar and planetary atmospheric study.

*Unified Astronomy Thesaurus concepts: Substellar companion stars (1648); Brown dwarfs (185); Debris disks (363); Exoplanet atmospheric composition (2021); Orbit determination (1175)*

# 1. Introduction

Low-mass directly imaged companions, ranging from brown dwarfs to extrasolar giant planets, present an ideal laboratory to test formation mechanisms, measure atmospheric properties, and constrain the frequency of companions at large orbital separations. The population of directly imaged giant planets remains small (see review and recent survey results; Bowler 2016; Galicher et al. 2016; Nielsen et al. 2019) and merits careful comparison with higher-mass brown dwarf companions observed over a wide range of system separations. Together with both free-floating planetary-mass objects (e.g., Liu et al. 2013; Gagné et al. 2015; Kellogg et al. 2016; Schneider et al. 2016) and field brown dwarfs (e.g., Burgasser et al. 2006; Kirkpatrick et al. 2011), substellar companions provide more readily characterizable analogs for study and comparison with planetary companions.

Within the brown dwarf population itself, a diversity of spectral features, intrinsic color, and luminosity has been observed. Spectroscopic analysis has shown that objects with the same optical spectral classification may have divergent near-infrared (NIR) colors (Leggett et al. 2003), with younger objects occupying a significantly redder and fainter region of color-magnitude space and systematically cooler effective temperatures (e.g., Filippazzo et al. 2015). The NIR color-magnitude analyses of field brown dwarfs, planetary-mass companions, and isolated low-mass, low-gravity objects have demonstrated distinct sequences between the three populations, corresponding to a diversity in physical properties across spectral type, gravity, and age. Liu et al. (2016) identified similar photometric properties for late M to late L isolated ("free-floating") substellar objects and young bound companions but with indications that the young field population is brighter and/or redder (even for systems within the same moving groups). With the relative paucity of substellar companions, each discovery provides a unique opportunity for comparison across various classes of substellar objects, contributing useful context that improves our understanding of brown dwarf/planet physical processes, formation pathways, and composition.

Bound substellar companion systems that also host a circumstellar disk provide opportunities to constrain and characterize the orbital properties of low-mass companions and disk-companion dynamical interaction. Studies with high-contrast imaging techniques have now imaged systems with companions in the region between the warm inner disk (e.g., HR 3549; Mawet et al. 2015) and the inner edge of a substantial outer debris disk (e.g., HR 2562; Konopacky et al. 2016). While directly imaged companions can explain the observed disk morphology for systems like HR 2562, in other cases, the presence of additional planets is needed to explain complex disk geometries (e.g., HD 95086; Rameau et al. 2016). Even fewer systems have proven amenable to resolved imaging of both the companion and the disk simultaneously. With

marginally resolved 70 µm Herschel PACS imaging, faint disk structure from broadband Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) imaging, and a close substellar companion at 11 au (Milli et al. 2017), the HD 206893 system is among only six systems hosting resolved disks and companions, in addition to β Pictoris (Lagrange et al. 2010), Fomalhaut (Kalas et al. 2008), LkCa 15 (Thalmann et al. 2010; Kraus & Ireland 2012), HD 106906 (Bailey et al. 2014), and PDS 70 (Keppler et al. 2018).

## 1.1. HD 206893 System Properties

Using the SPHERE (Beuzit et al. 2008) instrument on the Very Large Telescope (VLT), the SPHERE High Angular Resolution Debris Disk Survey (SHARDDS) discovered a low-mass companion orbiting the F5V star HD 206893 (Milli et al. 2017). The star was initially selected as a target for the SHARDDS resolved debris disk study with the SPHERE/IRDIS instrument (Dohlen et al. 2008) owing to the high fractional IR luminosity of its disk ($L_{dust}/L_* = 2.3 \times 10^{-4}$; Moór et al. 2006) and characterization of the spectral energy distribution (SED) from Chen et al. (2014). After an initial H-band detection of the companion in 2015, subsequent follow-up in 2016 with VLT/NaCo $L'$ imaging using an annular groove phase mask coronagraph (Mawet et al. 2005) recovered the source at $L'$. With a bright $L' = 13.43_{-0.15}^{+0.17}$ mag relative to the initial measurement of $H = 16.79 \pm 0.06$, indicating an extremely red color, the object's position in the color-magnitude diagram (CMD) is analogous to the L5–L9 field dwarf population. With two imaged epochs, Milli et al. (2017) confirmed the comoving nature of the object using an additional nondetection at the expected position of a background object in archival Hubble Space Telescope/NICMOS data. Without a known association with a moving group, and with literature ages ranging from 0.2 to 2.1 Gyr (Zuckerman & Song 2004; David & Hillenbrand 2015), companion mass estimates range from 24 to 73 $M_{Jup}$ using the COND models (Baraffe et al. 2003). With SPHERE integral field spectrograph (IFS) and IRDIS observations, Delorme et al. (2017) conducted a follow-up study of the system, including revised stellar properties of the primary (estimating solar metallicity and an age from 50 to 700 Myr) and analysis of $R \sim 30$ YJ spectra (0.95–1.64 µm) and 2.11 and 2.25 µm K-band photometry of the companion. Stolker et al. (2020) conducted a photometric study of the companion in the $L'$, NB4.05, and $M'$ filters with VLT/NaCo, further demonstrating its redness at longer wavelengths. These data and analyses substantiated the unusual properties of the companion, establishing it as the reddest of any currently known substellar objects, and provided an estimation of its spectral properties as those of a dusty, intermediate-gravity L dwarf.

In addition to the brown dwarf detection, the debris disk of the system was tentatively detected in the wide-field IRDIS data and also marginally resolved in archival Herschel 70 µm
data. From the Herschel imaging, Milli et al. (2017) applied a joint SED and imaging fit to estimate the parameters of the debris disk, reporting an inner edge of 50 au, position angle of 60° ± 10°, and inclination of 40° ± 10° from face-on.

The orbital properties of HD 206893 B were revisited by Grandjean et al. (2019), who performed a joint fit to radial velocity (RV) monitoring, Hipparcos and Gaia astrometry, and SPHERE and NaCo imaging to infer an orbital period for the companion of 21–33 yr. They derived its dynamical mass to be a potentially planetary mass of $10_{-5}^{+5} M_{Jup}$. However, Grandjean et al. (2019) noted that the fit appears dominated by the Hipparcos–Gaia proper-motion constraints, as the estimated mass from the joint fit is inconsistent at the $2\sigma$ level with the mass derived from the combination of only direct imaging astrometry and RV variation. The authors thus posit that the “B” companion cannot be responsible for the observed 1.6 yr RV drift, and that this variation may correspond instead to an additional ~15 $M_{Jup}$ interior companion with a short 1.6–4 yr period.

In this study, we present Gemini Planet Imager (GPI; Macintosh et al. 2008) observations of the HD 206893 system at the J, H, K1, and K2 bands with the goal of characterizing the imaged companion and providing the most comprehensive and highest-resolution spectral coverage to date. In Sections 2 and 3, we detail the GPI observations, data reduction, postprocessing, and analysis techniques. Assessment of the host star properties is described in Section 4. Results are described in Section 5, including companion photometry, location in the CMD, atmospheric characterization, astrometric analysis, and limits on additional companions. We examine potential explanations for the reddening of the system in Section 6 and conclude by comparing HD 206893 B with the known population of directly imaged companions and disk systems in Section 7.

## 2. Observations

The GPI Exoplanet Survey (GPIES; Macintosh et al. 2014) is a dedicated 900 hr direct imaging survey to discover and characterize giant planets within a sample of over 600 nearby young stars using the GPI IFS at the Gemini South observatory. As part of the GPIES campaign, IFS coronagraphic observations of HD 206893 were first obtained in H-band spectroscopic mode on UT 2016 September 22 (program ID GS-2015B-Q-01). The target was included in the campaign sample owing to its proximity, youth, and strong infrared excess, and the companion was noted in the preliminary data reductions. A total of 38 × 60 s IFS data cubes were taken under atmospheric conditions with ~0.0"68 average seeing, close to the median DIMM seeing for Gemini South (~0.0"65). Observations were taken near target transit at a nearly constant airmass, achieving a total of 32.7 of field rotation and enabling angular differential imaging (ADI; Marois et al. 2006a). To obtain wavelength calibration and account for instrumental flexure, argon arc lamp frames were taken immediately preceding the science observations. Additional H-band observations were taken in polarization mode following the spectroscopic imaging sequence but are not included in this analysis.

Follow-up spectroscopy was obtained on UT 2016 October 21 in K1 (74 × 60 s frames; 59.7 rotation) and UT 2016 November 17/18 in K2 (for a 2 night total of 148 × 60 s frames; 29.5 rotation). Seeing conditions were worse than median in all three follow-up data sets, as shown in Table 1.

The object was reobserved in K2 on UT 2017 November 9 with the same observing methodology, achieving additional field rotation under improved atmospheric conditions with lower residual wave-front error. The object was observed once more to obtain a J-band spectrum on UT 2018 September 24. A summary of the instrument modes, observations, and seeing data is provided in Table 1.

## 3. Data Reduction and Analysis

### 3.1. GPI Data Reduction

Calibrated spectral data cubes (x, y, λ) for data sets in each of the four bands were obtained from the raw IFS data using the GPI Data Reduction Pipeline v.1.4.0 (DRP; Perrin et al. 2016). Within the GPI DRP, raw IFS frames are dark-subtracted then interpolated over bad pixels. The effects of instrumental flexure on the positioning of individual microspectra in the IFS frames are accounted for by aligning the contemporaneous arc lamp frames against a library of deep reference argon arcs (Wolff et al. 2014), producing a wavelength calibration used to assemble spectral data cubes from the extracted microspectra (Maire et al. 2014). The resulting data cubes are then flatfielded, interpolated to a common wavelength axis, and corrected for geometric distortion (Konopacky et al. 2014), resulting in spectral data cubes with 37 channels each. The astrometric solution for the GPI has been shown to remain constant within the uncertainties (Konopacky et al. 2014; De Rosa et al. 2015), and we adopt plate scale and position angle values of $14.166 \pm 0.007$ mas pixel$^{-1}$ and $-0.10 \pm 0.13$ E of N.

To align and register the images, the stellar position behind the GPI coronagraphic occulting mask was estimated using the four satellite spots present in each slice of the data cube (Wang et al. 2014). These spots are attenuated replica images of the stellar point-spread function (PSF) introduced by the placement of a two-dimensional amplitude grating at the pupil plane (Marois et al. 2006b; Sivaramakrishnan & Oppenheimer 2006) and enable estimation of the object flux-to-stellar flux ratio (Wang et al. 2014). The adopted spot-to-star flux ratios for each filter and associated apodizer were $1.798 \times 10^{-4}$ for J, $2.035 \times 10^{-4}$ for H, $2.695 \times 10^{-4}$ for K1, and $1.905 \times 10^{-4}$ for K2 (Maire et al. 2014).

Outside of the GPI DRP, the fully calibrated data cubes were further postprocessed to subtract the contribution of the PSF and speckle noise. This was accomplished using three postprocessing techniques: Locally Optimized Combination of Images (LOCI; Lafrenière et al. 2007) for the J-, H-, and K2-band data; classical ADI (cADI; Marois et al. 2006a) for the K1 data, and pyKLIP for all four bands (Wang et al. 2015). In this work, the pyKLIP reductions from the autoreduced GPI pipeline (Wang et al. 2018) are used to measure the achieved contrast and sensitivity to companions (see Section 5.5 and the reduced pyKLIP images in Appendix A), and the LOCI/cADI reductions are used for the companion spectral extraction and astrometry.

An apodized high-pass filter following a Hanning profile (cutoff of 4 equivalent pixels) was applied to each frame to suppress slowly varying spatial frequencies within the data that are not well captured by PSF subtraction algorithms. The small amount of field rotation at H and increased speckle noise at J and K2 prevent a clean extraction of the companion with cADI, so LOCI was used to further process these bands. In order to
<table><caption>**Table 1**
Observation Summary</caption><thead><tr><th>Night (UT)</th><th>Band</th><th>λ/Δλ</th><th>Int. Time (minutes)</th><th>Field Rot. (deg)</th><th>Airmass Range</th><th>Seeing (arcsec)</th></tr></thead><tbody><tr><td>2016-9-22</td><td><i>H</i></td><td>44–49</td><td>37.7</td><td>32.7</td><td>1.05–1.06</td><td>0.68</td></tr><tr><td>2016-10-21</td><td><i>K1</i></td><td>62–70</td><td>73.5</td><td>59.7</td><td>1.05–1.09</td><td>1.17</td></tr><tr><td>2016-11-17*</td><td><i>K2</i></td><td>75–83</td><td>75.5</td><td>14.9</td><td>1.12–1.54</td><td>1.46</td></tr><tr><td>2016-11-18*</td><td><i>K2</i></td><td>75–83</td><td>71.5</td><td>14.6</td><td>1.12–1.40</td><td>1.44</td></tr><tr><td>2017-11-9</td><td><i>K2</i></td><td>75–83</td><td>88</td><td>24.1</td><td>1.08–1.32</td><td>†</td></tr><tr><td>2018-9-24</td><td><i>J</i></td><td>35–39</td><td>89</td><td>68.1</td><td>1.05–1.11</td><td>†</td></tr></tbody></table>

Note. Data sets denoted with an asterisk are shown for reference but not used in the modeling analyses presented in this work due to lower S/N. The † denotes that facility seeing data were unavailable during the observations.

reproduce images optimized to match the spatial distribution of speckles, the best-matched PSFs for each band were then estimated from a library of reference images with the following routine parameters: a PSF subtraction annulus of $d_r = 5$ pixels, optimized region for PSF comparison of $300 \times$ the FWHM, geometry factor of $g = 1$, and separation criteria between regions of $N_δ = 0.75$. For K2, the decreasing signal-to-noise ratio (S/N) of satellite spots near the red edge of the band, as well as the increased brightness of the companion, impacted the fidelity of the spectral extraction, requiring a modified high-pass filter with a cutoff of 10 pixels. (Extractions of the K2 spectrum using various algorithms and filter sizes are presented in Appendix A.) The K1 data benefit from a large field rotation, and simple cADI was used to perform speckle subtraction, subtracting the stellar halo contribution by taking the median of all science frames and subtracting the median from each individual frame.

## 3.2. Astrometry, Spectroscopy, and Covariance

The astrometry of the companion was measured using the collapsed broadband PSF-subtracted data cubes. To extract the spectroscopy of the companion, the process of injecting a negative fake planet was applied (see Marois et al. 2010). A negative template PSF with the estimated position and flux of the companion was inserted into the raw data set, and then the postprocessing pipelines applied above were repeated iteratively to minimize the residuals within a $2 \times 2$ FWHM region centered at the companion position, adjusting the x- and y-positions and flux of the companion in each iteration. This allowed for correction of any differences in algorithm throughput between the LOCI and cADI reductions. To correct for flux offsets between the K1 and K2 data sets in the overlap region of the two bands (from 2.10 to 2.20 µm), the K1 spectrum flux scaling was shifted to minimize the $\chi^2$ value between the two spectra, resulting in a 4% downward shift.

As the raw IFS data include over 36,000 microspectra, the process of spectral extraction and interpolation (from the 16 pixel microspectra to 37 wavelength elements) introduces correlations between the resulting 37 spectral channels in the final data cubes. Applying the methodology of Greco & Brandt (2016), we calculate the covariance between channels to better estimate spectral uncertainties (see Appendix A), an approach that has been applied to other recent analyses of GPI companion spectra (see De Rosa et al. 2016; Johnson-Groh et al. 2017; Rajan et al. 2017). This effect can then be incorporated into the error budget of the individual spectro-photometric points for model fitting, as described in Section 5.3.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Stellar SED and best-fit BT-NextGen model spectrum for the primary star HD 206893 A, with residuals shown in the bottom panel. Best-fit parameters correspond to intermediate age (601 Myr, albeit with large uncertainties of $^{+420}_{-380}$ Myr) and low extinction ($A_V = 0.07^{+0.05}_{-0.04}$).</figcaption></figure>

## 4. Host Star Properties

We derive stellar properties from comparison of the stellar SED to BT-NextGen stellar model atmospheres (Allard et al. 2012) combined with the MIST stellar isochrones (Choi et al. 2016; Dotter 2016). For each stellar model grid at a given mass and age, stellar radius, temperature, and surface gravity are estimated and synthetic photometry is derived. Figure 1 shows the combined Strömgren, Geneva, Tycho, Hipparcos, Two Micron All Sky Survey (2MASS), and Wide-field Infrared Survey Explorer (WISE) photometry for HD 206893 A with the best-fit stellar atmospheric model overplotted. A Markov Chain Monte Carlo (MCMC) approach was applied to determine the best-fitting model, with stellar age, mass, extinction, and distance as free parameters (following De Rosa et al. 2016; Nielsen et al. 2017). Uniform priors in age, mass, and $A_V$ were applied, as well as a Gaussian distance prior centered at the observed Gaia parallax. The metallicity is assumed to be solar, concordant with the findings of Delorme et al. (2017). Figure 2 shows the posterior distribution of stellar properties from the MCMC analysis. Based upon the age posterior from these SED fits, we find that the stellar age is $601^{+420}_{-380}$ Myr with $1\sigma$ confidence. The best-fit stellar parameters are provided in Table 2. The large uncertainty on the stellar age from this approach is commensurate with the difficulty in determining precise age estimates for F5V stars, as noted in the detailed analysis of HD 206893 A by Delorme et al. (2017).

The host star HD 206893 A is not a known member of a stellar association, which would allow for a more precise age determination. We cross-checked the possibility of young moving-group membership using the Bayesian Analysis for
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Posterior distribution of stellar properties of the primary star HD 206893 A. The asymmetric distribution on the stellar age suggests a broad range of ~200 Myr to ~1 Gyr at 1σ confidence.</figcaption></figure>

Nearby Young Associations (BANYAN) Σ tool (Gagné et al. 2018). Given Gaia DR2 values for stellar parallax, proper motion, and RV<sup>43</sup> (Table 2; Gaia Collaboration et al. 2018), BANYAN Σ provides a 61.1% probability of membership in the Argus moving group and a 38.9% probability of field membership.

The existence and age of the Argus association have been the subject of debate, with estimates ranging from ~40 (Torres et al. 2008) to 268 (Bell et al. 2015) Myr. The reality of the association has been called into question (Bell et al. 2015); however, it has recently been reestablished as a field association of age 40–50 Myr by Zuckerman (2019). From the low BANYAN Σ probabilities of both Argus and field membership (<80%–90%), no definitive conclusion can be drawn. However, due to the unique kinematics of Argus and signatures of youth among its member stars, the association is also included in a recent Bayesian analysis of nearby young moving groups by Lee & Song (2019). Calculating the membership probabilities using the models of Lee & Song (2019) yields similar results to BANYAN Σ, with *p* (Argus) = 63%, *p* (Field) = 37%; kinematic comparisons are shown in Figure 3. Given the slight preference for Argus membership from both moving-group assessments, we consider a range of possible ages in estimating companion properties. Specifically, we use the Baraffe et al. (2002) DUSTY and Baraffe et al. (2003) COND models at 50, 100, 120, and 500 Myr to estimate companion mass in Section 5.1, where the upper limit on age is set by the stellar SED analysis. At the same time, we note that the significant mid- and far-IR excess of the system (Chen et al. 2014) also likely points to a younger age.

<table><caption>**Table 2**
Properties of the HD 206893 System</caption><thead><tr><th>System Property</th><th>Value</th><th>Unit</th><th>References</th></tr></thead><tbody><tr><td>Parallax</td><td>26.08 ± 0.53, 24.51 ± 0.06</td><td>mas</td><td>1, 2</td></tr><tr><td>Distance</td><td>38.3 ± 0.8, 40.8 ± 0.1</td><td>pc</td><td>1, 2</td></tr><tr><td>μ<sub>α</sub></td><td>93.67 ± 0.66, 93.78 ± 0.09</td><td>mas yr<sup>-1</sup></td><td>1, 2</td></tr><tr><td>μ<sub>δ</sub></td><td>0.33 ± 0.37, 0.02 ± 0.08</td><td>mas yr<sup>-1</sup></td><td>1, 2</td></tr><tr><td>RV</td><td>−12.45 ± 0.59</td><td>km s<sup>-1</sup></td><td>2</td></tr><tr><td>Age</td><td>~50 (kinematics); 601<sup>+420</sup><sub>-380</sub> (stellar SED)</td><td>Myr</td><td>1, 3, 4</td></tr><tr><td>A<sub>V</sub></td><td>0.07<sup>+0.05</sup><sub>-0.04</sub></td><td>mag</td><td>4</td></tr><tr><td></td><td colspan="3" style="text-align: center;">HD 206893 B</td></tr><tr><td>Spectral type</td><td>F5V</td><td>L4–L8<sup>†</sup></td><td>...</td><td>1/4</td></tr><tr><td>M<sub>J</sub></td><td>2.82 ± 0.06</td><td>15.39 ± 0.04</td><td>mag</td><td>4</td></tr><tr><td>M<sub>H</sub></td><td>2.64 ± 0.06</td><td>13.68 ± 0.04</td><td>mag</td><td>4</td></tr><tr><td>M<sub>Ks</sub></td><td>2.54 ± 0.06</td><td>12.00 ± 0.07</td><td>mag</td><td>4</td></tr><tr><td>M<sub>L'</sub></td><td>2.48</td><td>10.39 ± 0.22</td><td>mag</td><td>1/4</td></tr><tr><td>J</td><td>5.87 ± 0.02</td><td>18.44 ± 0.03</td><td>mag</td><td>4</td></tr><tr><td>H</td><td>5.69 ± 0.03</td><td>16.73 ± 0.03</td><td>mag</td><td>4</td></tr><tr><td>K<sub>s</sub></td><td>5.59 ± 0.02</td><td>15.05 ± 0.07</td><td>mag</td><td>4</td></tr><tr><td>J<sub>MKO</sub></td><td></td><td>18.38 ± 0.03</td><td>mag</td><td>4</td></tr><tr><td>H<sub>MKO</sub></td><td></td><td>16.82 ± 0.03</td><td>mag</td><td>4</td></tr><tr><td>K<sub>MKO</sub></td><td></td><td>15.02 ± 0.07</td><td>mag</td><td>4</td></tr><tr><td>L'</td><td>5.52</td><td>13.43<sup>+0.17</sup><sub>-0.15</sub></td><td>mag</td><td>1/3</td></tr><tr><td>Mass</td><td>1.31 ± 0.03</td><td>12/20/20/40‡</td><td>M<sub>☉</sub>, M<sub>Jup</sub></td><td>4</td></tr><tr><td>T<sub>eff</sub></td><td>6500 ± 100*</td><td>1200–1800<sup>†</sup></td><td>K</td><td>3, 4</td></tr><tr><td>log(g)</td><td>4.45 ± 0.15*</td><td>3.0–5.0<sup>†</sup></td><td>log<sub>10</sub><br>(cm s<sup>-2</sup>)</td><td>3, 4</td></tr></tbody></table>

**Note.** Photometry is provided in the 2MASS system, with MKO also provided for the companion. Companion photometry has been synthesized by integrating the GPI spectra over the relevant bandpasses accounting for differences in the 2MASS and MKO transmission profiles, as described in Section 5.1. ‡Estimated masses for the companion in $M_{Jup}$ are provided assuming four distinct system ages of 50, 100, 120, and 500 Myr. Asterisks indicate stellar parameters derived from FEROS spectroscopy by Delorme et al. (2017). The † symbols indicate that the large reported ranges in spectral type, effective temperature, and surface gravity reflect the difficulty in matching the observed companion spectrum to L dwarf standard templates and atmospheric models; see Sections 5.2 and 5.3. References: (1) Milli et al. (2017); (2) Gaia DR2: Gaia Collaboration et al. (2018); (3) Delorme et al. (2017); (4) this work.

# 5. Results

## 5.1. GPI Images and Integrated Photometry

Figure 4 shows the final median-collapsed GPI images in the *J*, *H*, *K1*, and *K2* bands from the LOCI and cADI reductions. The images show significant orbital motion of the companion between the earliest 2016 *H* epoch and the latest 2018 *J* epoch, as well as the higher-S/N detections at redder wavelengths, consistent with the upward slope of the observed spectrum.

Broadband photometry for HD 206893 B was estimated from the GPI spectral data cubes by interpolating both 2MASS and Maunakea Observatory (MKO) transmission filter profiles (Tokunaga et al. 2002) to the GPI wavelength scaling. The *J*, *H*, and *Ks* (or *K*) filter profiles were convolved with a Gaussian filter of FWHM corresponding to each GPI band, then interpolated to the corresponding GPI wavelength grid. The GPI spectra (combining both *K1* and *K2* as described in Section 3.2 to cover the full *K* band) were then integrated over the *J*, *H*, and *Ks* or *K* bandpasses to derive the magnitudes in

<footnote n="43"><sup>43</sup> The RV drift of −0.1 ± 0.5 km s<sup>−1</sup> identified in Grandjean et al. (2019) is within the Gaia DR2 uncertainty (−12.45 ± 0.59 km s<sup>−1</sup>) and likely does not affect this kinematic estimate.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 3. The XYZUVW position of HD 206893 (blue point) compared to moving-group models from Lee & Song (2019).</figcaption></figure>

Table 2, accounting for zero-point corrections. The derived photometry allows for positioning HD 206893 B on NIR CMDs of substellar objects, as shown in Figure 5. The derived photometry from the GPI is consistent within the uncertainties with the VLT/SPHERE magnitudes reported in Delorme et al. (2017) and Milli et al. (2017). The position of HD 206893 B on the NIR CMDs demonstrates its extraordinary red color, distinct among even the reddest planetary-mass companions, free-floating objects, and low-gravity and/or otherwise peculiar field objects.

## 5.2. Composite Spectrum and Spectral Classification

The resulting combined spectrum is shown in Figure 6 with each of the extractions from the reduced, postprocessed *J*, *H*, *K*1, and *K*2 GPI data cubes. The VLT spectra and photometry from Milli et al. (2017) and Delorme et al. (2017) are also shown, demonstrating good agreement between GPI data and SPHERE. The GPI *K*2 spectral slope beyond 2.2 µm was noted to be unusually steep in the first epoch of observations and motivated a deeper second-epoch set of observations. The spectral morphology persisted between the two epochs of *K*2 data taken 2 yr apart, and the steepness appears to be dependent upon the spectral extraction as described in Appendix A. For the remaining analysis in this work, only the extraction from the higher-S/N *K*2 spectrum is used.

To determine the spectral properties of HD 206893 B in comparison with known brown dwarfs and atmospheric models, the combined spectra from SPHERE and GPI spanning *Y* to *K*2 were analyzed using the SpeX Prism Library Analysis Toolkit (SPLAT; Burgasser & the SPLAT Development Team 2017). Spectral typing was performed within SPLAT with a $\chi^2$ minimization approach, classifying the object by template across the full ensemble of L0–T5 spectra in the SpeX Prism Library. Based on the exceptional redness of the spectrum, the resulting matches within the spectral library correspond to those of the coolest L-type objects, at $L6 \pm 1$ subclasses. While no single spectrum provides a close morphological match across the full *Y*–*K*2 range, the three closest matches within the spectral library were late-type, low-gravity objects, all with reduced $\chi^2$ values of $\sim$9: WISE J174102.78–464225.5 (>L5γ; Faherty et al. 2016), 2MASS J11193254–1137466 (L7γ, likely TW Hya member; Kellogg et al. 2015), and 2MASS J11472421–2040204 (L7γ, likely TW Hya member; Kellogg et al. 2016). When the SPLAT library comparison is extended to earlier spectral types, the closest-matching objects correspond to M9 pre-main-sequence stars with extreme levels of extinction (e.g.,

2MASS J03444520+3201197 in IC 348 and 2MASS J04325026+2422115 in Taurus, with reduced $\chi^2$ values of 7–8); given the distance and low extinction to HD 206893 and its companion absolute magnitude, these matches are not likely to share similar physical properties to HD 206893 B, and goodness of fit is driven by the redness of the object alone.

The large observed spread in NIR colors observed among L dwarfs of the same optical spectral type (e.g., Leggett et al. 2003; Faherty et al. 2013) has motivated efforts to produce a systematic means of spectral typing, including qualitative comparison on a band-by-band basis to account for variations in NIR spectral morphology for objects of the same optical spectral classification. Following the work of Cruz et al. (2018), which aimed to reconcile differences in NIR spectral features for the L dwarf population, we compare the HD 206893 B spectrum to an ensemble of L dwarf field and low-gravity standards using the Ultracool Typing Kit (UTK;<fnref n="44" /> Abrahams & Cruz 2017). The UTK analysis is shown in Figure 7 and involves normalizing each *J*, *H*, or *K* spectrum individually and comparing the normalized spectra with the spectral templates. The comparison demonstrates the similarities between HD 206893 B and typical L dwarfs. The spectra are qualitatively similar to the later field and young L populations, commensurate with the SPLAT typing analysis. However, the steep slope of the blue side of the *H* band and the peaked morphology more closely match the low-gravity population than field objects, potentially indicating youth (although challenges exist in using the triangular shape to reliably distinguish dusty atmospheres from low-gravity ones; Allers & Liu 2013). Figure 8 shows the three closest spectral matches from SPLAT, including 10 Myr old candidate TW Hya members, in addition to the unusually red AB Doradus member WISEP J004701.06+680352.1 (Gizis et al. 2012). Spectra have been normalized per band over the same normalization ranges used in UTK, extending a similar analysis to later low-gravity L dwarfs. These late-type comparisons show the greatest morphological similarities to HD 206893 B, particularly across the *J* and *H* bands, when normalized on a band-by-band basis. As such, we adopt a conservative range of $L6 \pm 2$ for the companion spectral type, noting stronger similarities to young low-gravity objects in this spectral range than to the late-type field population seen in Figure 7.

<footnote n="44">⁴⁴ https://github.com/BDNYC/UltracoolTypingKit</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The GPI PSF-subtracted images of HD 206893 B in the *J* (upper left), *H* (upper right), *K*1 (lower left), and *K*2 (lower right) bands from the LOCI and cADI reductions. Each individual image is median-collapsed and scaled separately to minimize the contributions of residual speckle noise. A noticeable change in the companion position angle can be seen between the earliest 2016 epoch (upper right) and the latest 2018 epoch (upper left).</figcaption></figure>

5.3. *Brown Dwarf Atmospheric Modeling*

The combined spectra and photometry of HD 206893 B from GPI, SPHERE, and NaCo, spanning 0.95–3.8 µm, were compared to four grids of theoretical atmosphere models in order to derive estimates of effective temperature, surface gravity, and metallicity. The model grids used were (1) a subset of the Sonora models, a new model grid designed for substellar and young giant planet atmospheres (M. S. Marley et al. 2020, in preparation; see also the currently available solar-metallicity cloud-free grid<fnref n="45" /> from Marley et al. 2018); (2) the Cloudy-AE<fnref n="46" /> models, originally developed for the HR 8799 planets (Madhusudhan et al. 2011); (3) the DRIFT-PHOENIX<fnref n="47" /> atmosphere models (Helling et al. 2008; Witte et al. 2009, 2011); and (4) the BT-Settl (2015)<fnref n="48" /> grid (Allard et al. 2012; Allard 2014). For each grid, individual models spanning

ranges of *T*<sub>eff</sub> and log(g) space were binned to the resolution of the GPI spectra and interpolated over the GPI wavelength axes. All models used were of solar metallicity, with the exception of DRIFT-PHOENIX, which included both solar and supersolar [*M*/*H*] values. Object radius was treated as a fit parameter, with each model scaled to match the object SED using the Gaia DR2 distance and a range of plausible object radii. The explored parameter ranges for the model grids are summarized in Table 3, scaling the object radius from 0.5 to 3.5 *R*<sub>Jup</sub> in steps of 0.01. The resulting χ<sup>2</sup> fit statistics were calculated accounting for covariances in spectral channels, following the approach outlined in De Rosa et al. (2016), with zero-valued covariance assigned to the SPHERE IFS *Y* and *J*/*H* data not covered by the GPI. Best-fit models were determined separately for each band and the full spectrum, allowing us to investigate disparities between the best-fit models between wavelength regimes. The resulting best-fit models from the Sonora, DRIFT-PHOENIX, Cloudy-AE, and BT-Settl grids are shown in Figure 9 and incorporate covariances. Results from a standard covariance-free χ<sup>2</sup> approach are also shown on a

<footnote n="45">⁴⁵ https://zenodo.org/record/1309035#/.W2JIYC2ZNjs</footnote>

<footnote n="46">⁴⁶ https://www.astro.princeton.edu/~burrows/8799/8799.html</footnote>

<footnote n="47">⁴⁷ http://svo2.cab.inta-csic.es/theory/newov</footnote>

<footnote n="48">⁴⁸ https://phoenix.ens-lyon.fr/Grids/BT-Settl/CIFIST2011_2015/</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The NIR CMDs of L- and T-type brown dwarfs and known planetary-mass objects in the 2MASS photometric system. Data are compiled from the following resources: field brown dwarfs with Gaia parallaxes from Smart et al. (2019; orange circles); the Database of Ultracool Parallaxes (http://www.as.utexas.edu/~tdupuy/plx/Database_of_Ultracool_Parallaxes.html; gray triangles), maintained by T. Dupuy, originally in Dupuy & Liu (2012), Dupuy & Kraus (2013), and Liu et al. (2016; shown separately with cyan diamonds); and substellar objects considered analogous to planets (blue squares) from Faherty et al. (2016). Objects with low-gravity indicators are shown in magenta, demonstrating the known red offset of these objects from the field sequence. The left panel shows the *M*<sub>J</sub> vs. *J* − *Ks* CMD, with HD 206893 B indicated by a red star and benchmark planetary-mass objects labeled (yellow circles). The right panel shows the *M*<sub>H</sub> vs. *H* − *Ks* CMD. At *H* − *Ks* = 1.68, HD 206893 B occupies a far redder space in the CMD than any other known substellar object. The additional points for HD 206893 B correspond to *A*<sub>V</sub> values of 10–3.4 (left to right), adopted to match a similar range of extinction values estimated for the dusty, younger CT Cha system, and *A*<sub>V</sub> = 10 (green star) indicates the extreme values needed to reconcile the object position with the red edge of the substellar sequence.</figcaption></figure>

band-by-band basis with faint dashed lines in Figure 9, the
results of which are summarized in Appendix B. While the
standard fits without covariance estimation sometimes more
closely visually match the GPI spectral data points, the derived
physical parameters are in better agreement across bands when
accounting for spectral covariance.

Given the extraordinary redness of the companion, none of
the models provides a qualitatively good fit to the full
spectrum, and reduced χ² statistic values range from
~1.2 to 8. The best fits to the Sonora and Cloudy-AE60
models all lie at the edge of the available grid parameters
(Table 3), suggesting that the best fits likely reside outside of
the range of the current model grids. As noted in Delorme et al.
(2017), the unusual SED of the object and its extreme spectral
features present significant discrepancies relative to existing
model atmosphere grids, necessitating additional sources of
reddening. Challenges in atmospheric model fitting of L
dwarfs, particularly in the reproduction of the NIR spectral
slope and H-band morphology of young/low-gravity objects,
have been noted in previous studies (e.g., Manjavacas et al.
2014) and attributed to the treatment of dust in current
theoretical atmosphere models. We explore this possibility
further in Section 6.2 and note that determining the cause of
mismatches between observations and current models may
involve considering a wide range of particle properties, overall
dust content, and multispecies gas opacities, tasks well suited
for retrieval-based methods.

When we fit the full spectrum, the best-fit results from each
model grid show a range of physical parameters; temperatures
range from 1200 to 1700 K, surface gravities from log(g) of

3.5 to 5.0, and object radii from 0.96 to 1.55 R<sub>Jup</sub>. However, the derived properties vary significantly when comparing results across fits to individual bands. The derived log(g) values from the model fits and the inferred radii necessary to scale the flux are inconsistent, a discrepancy noted for previous fits of exoplanet spectra and photometry (e.g., HR 8799 bcd, β Pic b, 51 Eri b, and others; Marois et al. 2008; Marley et al. 2012; Morzinski et al. 2015; Rajan et al. 2017). A detailed assessment of the full-spectrum fitting using each model grid, as well as brief summaries of the model grid properties, are provided in Appendix B.

Following the methodology applied for spectral typing the
object, we performed fits to the four GPI bands individually for
each set of model grids, in addition to the full spectral fit. These
fits are shown in the four subpanel rows of Figure 9, and each
fit accounts for spectral covariance within that band. Account-
ing for spectral covariance, the best-fit models agree within
100 K across all four bands for the Sonora and Cloudy-AE60
models and across JHK2 for DRIFT-PHOENIX and BT-Settl.
Departures in K1 for the DRIFT-PHOENIX and BT-Settl
model grids may be attributed to the higher level of spectral
covariance in the K1 band. Similarly, the high level of spectral
correlation observed in the H band leads to a mismatch
between the best-fit models and data points for all four suites of
models, an effect observed in previous IFS studies (e.g., Rajan
et al. 2017). For comparison, the best-fit models adopting a
standard χ² approach are shown by faint dashed lines in each
panel. The estimated log(g) values are consistent across bands
only for the Sonora and Cloudy-AE60 models, favoring lower-
gravity fits. Each of the other three model suites provides a
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Combined GPI J, H, K1, and K2 spectra of HD 206893 B, overlaid with the spectra and photometry from Milli et al. (2017) and Delorme et al. (2017), which show good agreement between the SPHERE IFS and GPI spectroscopy. IRDIS K1/K2 photometry has been synthesized from the GPI spectra (white stars) to enable a comparison with the measured SPHERE photometry (white diamonds).</figcaption></figure>

broad range of surface gravities and radii depending on the band in question.

## 5.4. Astrometric Analysis

The four new astrometric epochs from the GPI observations presented here were combined with five additional astrometric measurements from VLT/SPHERE and VLT/NaCo (Delorme et al. 2017; Milli et al. 2017; Grandjean et al. 2019) to explore the range of potential orbital parameters for HD 206893 B. The full set of astrometric measurements used in this analysis is provided in Table 4. Following the astrometric analyses described in De Rosa et al. (2015) and Rameau et al. (2016), a parallel-tempered Bayesian MCMC approach using the emcee package (Foreman-Mackey et al. 2013) was used to simultaneously fit eight companion orbital parameters, namely: semimajor axis (a), inclination (i), eccentricity (e), sums and differences of longitude of ascending node and argument of periastron (Ω ± ω), epoch of periastron passage (τ), orbital period, and system parallax. We adopt the following standard prior distributions on the orbital parameters: log uniform in a, uniform in e and cos (i), and Gaussian for the system mass and parallax. We initialized 512 walkers at 16 different temperatures; the lowest temperature samples the posterior distribution, while the highest temperature samples the prior. We advanced each chain for 1 million steps, saving every 100 steps. The first half of the chain was discarded as a “burn-in” period to ensure the final posterior distributions were not a function of the initial positions of the walkers.

The resulting posterior distributions for a subset of selected orbital parameters are shown in Figure 10. The 1σ confidence intervals from the posterior distributions suggest a tightly constrained semimajor axis of $10.4_{-1.7}^{+1.8}$ au, orbital inclination of $145.6_{-6.6}^{+13.8}$, and apoapsis distance of $11.7_{-0.6}^{+2.7}$ au and place less restrictive constraints on the eccentricity ($0.23_{-0.16}^{+0.13}$) with a corresponding period range of $29.1_{-6.7}^{+8.1}$ yr. The narrower posterior distribution on the apoapsis distance is consistent with observing the companion near apastron, with the

uncertainty on the semimajor axis dominated by a less-stringent constraint on the radius of periapsis. This can also be seen in the covariance of semimajor axis and eccentricity posterior distributions in Figure 10, which favor smaller values of a for more eccentric ($e > 0.2$) orbits.

We investigated the coplanarity of the orbit of the companion with the spatially resolved debris disk by comparing the inclination and position angle of the visual orbit on the sky to those estimated for the debris disk by Milli et al. (2017). The disk inclination and position angles are defined over more restrictive ranges than for a visual orbit; the disk is not measured to orbit in a particular direction around the star, and a 180° ambiguity exists in the position angle measurement. Under the assumption that the disk and orbital plane of the companion are nonorthogonal, we adopt an inclination for the disk of $i_d = 140° ± 10°$ and a position angle of $Ω_d = 60° ± 10°$. We drew random variates from two Gaussian distributions describing these parameters for each sample of the posterior distribution of the visual orbit, assuming that the measured geometry of the disk and the fit to the visual orbit are uncorrelated. We measure the mutual inclination, $i_m$, between the disk and orbit using the following relation:

$$ \cos(i_m) = \cos(i_c)\cos(i_d) + \sin(i_c)\sin(i_d)\cos(Ω_c - Ω_d), \quad (1) $$

where the c and d subscripts refer to the companion and disk, respectively. The resulting distribution of $i_m$ and correlations with the other orbital elements are shown in Figure 10. The position angle of nodes for the orbit of the companion (Ω) is plotted from 0° to 180°, the same restricted range as the position angle of the disk. An equally good fit to the visual orbit can be found with $ω + π$ and $Ω + π$; however, these orbits would be significantly misaligned with the assumed disk position angle of $60° ± 10°$. The distribution of mutual inclinations for two randomly oriented planes is also shown for reference, and $i_m$ appears significantly shifted toward more coplanar solutions relative to this prior distribution. The current measurement precision for both the disk geometry and the
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison of the combined GPI+SPHERE spectra (black) to a suite of L dwarf standard templates of varying gravity (red) using the UTK (Abrahams & Cruz 2017). The low-gravity (γ) standards from Cruz et al. (2018) include only spectral types up to L4. Each band is independently normalized using the band-specific normalization ranges described in Cruz et al. (2018): 0.87–1.39 µm for zJ, 1.41–1.89 µm for H, and 1.91–2.39 µm for K. Morphologically, the HD 206893 B spectrum appears most similar to late, low-gravity objects.</figcaption></figure>

visual orbit are not sufficient to exclude moderately misaligned configurations ($i_m \sim 20^\circ$).

A representative sample of 250 orbits drawn randomly from the posterior distribution is shown in Figure 11. The plotted orbit colors indicate the mutual inclination value for the companion orbit relative to the disk geometry inferred from modeling in Milli et al. (2017), favoring more coplanar configurations. Misaligned configurations with greater values of $i_m$ are generally more eccentric, corresponding to smaller semimajor axis and periape distance, which can be seen from the posterior covariances in Figure 10. The long tail of the posterior distribution for apoapse distance, $r_{apo}$, extends to large values coincident with the estimated inner edge of the disk at $\sim 50$ au by Milli et al. (2017).

## 5.5. Sensitivity to Additional Companions

From the pyKLIP reductions described in Section 3 and shown in Figure 18, we calculate contrast curves for each GPI band. These curves are shown in Figure 12; they are calculated by estimating noise levels from the PSF-subtracted images and correcting for instrument throughput, calculated by injecting and recovering false planets with L-type spectra (see Wang et al. 2018). The sensitivity to companions is then estimated from the contrast curves following the methodology of Nielsen et al. (2019). In brief, a synthetic population of $10^4$ companions is drawn from a grid sampling mass and semimajor axis values, and each companion is randomly assigned orbital parameters

(inclination, eccentricity, argument of periastron, and epoch of periastron passage). Projected separation and magnitude difference are estimated and compared to the contrast curve to determine detection probability. Variation in contrast by epoch is accounted for by stepping the companion orbit forward in time in order to generate a new estimation for companion separation and contrast that can be compared to the observed contrast curve. The recovery rate for companions at a given mass and semimajor axis is then translated into a completeness percentage. The results are shown in Figure 13, where the mass of HD 206893 B corresponds to an assumed system age of 250 Myr. It is possible to exclude additional companions in the GPI images in the planetary regime down to 5 $M_{Jup}$ at $\sim 20$ au and in the brown dwarf regime down to $\sim 10–20 M_{Jup}$ at $\sim 10$ au. For comparison, the inner additional companion signal predicted from RV variation by Grandjean et al. (2019) suggests an $\sim 15 M_{Jup}$ object at an orbital separation of 1.4–2.6 au. The requisite contrast to detect such an object ($<10^{-4}$ in the NIR at a separation of 50–60 mas) is beyond the sensitivity limits of current extreme adaptive optics direct imaging instruments and can be compared to the contrast curves and sensitivity of our data shown in Figures 12 and 13.

# 6. Discussion

## 6.1. Companion Redness and System Extinction Properties

Given the existence of circumstellar material in the HD 206893 system, dust extinction may have an effect on
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the combined GPI+SPHERE spectra (black) to young, low-gravity L5–L7 objects, including the closest substellar spectral matches in the SpeX prism library. Each panel has been normalized on a band-by-band basis to the peak of *J*, *H*, or *K* (top to bottom). Normalization at shorter wavelengths shows good morphological matches to each of *J* and *H* and emphasizes the shallowness of the water band at 1.4 µm.</figcaption></figure>

<table><caption>**Table 3**
Parameter Ranges in the Model Grid Search</caption><thead><tr><th>Model Grid</th><th><i>T</i><sub>eff</sub> (K)</th><th>log[g (cgs)]</th><th>Notes</th></tr></thead><tbody><tr><td>Sonora</td><td>1100–1600</td><td>3.75–5.0</td><td>Solar metallicity, <i>f</i><sub>sed</sub> = 1</td></tr><tr><td>Cloudy-AE60</td><td>600–1700</td><td>3.5–5.0</td><td>Solar metallicity, fixed particle size of 60 µm</td></tr><tr><td>DRIFT-PHOENIX</td><td>1000–2000</td><td>3.0–5.0</td><td>Metallicity includes solar, supersolar [M/H] = 0.3</td></tr><tr><td>BT-Settl</td><td>1200–2050</td><td>2.5–5.5</td><td>Solar metallicity</td></tr></tbody></table>

**Note.** For all models tested, the object radius was scaled by factors ranging from 0.5 to 3.5 $R_{\text{Jup}}$ in the model fitting, in steps of 0.01.

the measured companion magnitudes; however, the nature of reddening in the system is currently unknown. The unusual redness of the companion could be explained by a range of reddening sources, namely, conventional interstellar reddening, extinction due to the debris disk, the presence of circumstellar disk material around the companion itself, or dust within the atmosphere of the companion itself.

We adopted reddening values of $A_V = 3.4$ and 5.2, in conjunction with the extinction laws of Rieke & Lebofsky (1985), to estimate the potential shift in color–magnitude space due to extinction caused by dust in the vicinity of the companion. These values for $A_V$ were chosen from best-fit extinction estimates for the substellar companion in the CT Cha disk system (Schmidt et al. 2008; Wu et al. 2015). As CT Cha is considerably younger (~2 Myr) than

HD 206893, this is likely an overestimate of potential extinction caused by the far less optically thick HD 206893 debris disk. The resulting dereddened CMD positions are shown in the right panel of Figure 5, with an additional value of $A_V = 10$ shown for comparison to illustrate the extreme and likely unphysical amount of extinction required to meet the red edge of the field L dwarf sequence.

To determine which extinction estimates correspond to realistic values for a companion viewed edge-on within a debris disk, we compare dust column density and optical depth estimates to typical circumstellar disk values. Using standard relations for visual extinction and gas column density for the interstellar medium ($N_H$ (cm$^{-2}$) = $2.21 \times 10^{21} A_V$ (mag); Güver & Özel 2009) and translating this relation into a dust column density using the canonical 100:1 gas-to-dust ratio yields a dust mass column density of $3.7 \times 10^{-4} \text{ g cm}^{-2}$ for $A_V = 10$. For comparison, visual extinction for the edge-on AU Mic debris disk system with higher fractional infrared luminosity ($3.9 \times 10^{-4}$; Matthews et al. 2015) has been estimated as $A_V = 0.5$ (Yee & Jensen 2010). Assuming small grains and a typical silicate density (0.5 µm and 3.27 g cm$^{-3}$) yields a particle column density of $2 \times 10^8 \text{ grains cm}^{-2}$, suggesting that a debris disk with $A_V = 10$ is unrealistically optically thick ($\tau \sim 1.7$). We thus conclude that the companion’s extreme redness is not likely to be caused by the debris disk, nor is this level of extreme optical depth likely in a circumsecondary disk. Further effects of reddening due to the disk environment are explored in detail in Section 6.3. The $A_V$ estimate of 0.03 from stellar SED modeling in Figure 2 also
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Best-fit models from a grid comparison of the GPI spectra (black points) and SPHERE+NaCo spectra and photometry (gray points) with four theoretical atmosphere models (Sonora in purple; M. S. Marley et al. 2020, in preparation; Cloudy-AE60 in blue; Madhusudhan et al. 2011; DRIFT-PHOENIX in dark green; Helling et al. 2008; Witte et al. 2009, 2011; and BT-Settl in light green; Allard et al. 2012). The top panel shows the best fits to the full spectral coverage from the combined SPHERE, GPI, and NaCo data. The integrated L' flux is shown for each set of models with squares. The four rows of subpanels show the results from fitting the four GPI bands separately. Corresponding $\chi^2$ best-fit model parameters with (solid lines) and without (faint dashed lines) incorporation of covariances are shown in each panel.</figcaption></figure>

argues against a significant reddening effect due to interstellar extinction, leaving dust in the companion atmosphere as the most viable scenario.

## 6.2. Atmospheric Reddening from Dusty Aerosols

The known L dwarf population exhibits a wide range of atmospheric color in the NIR, spanning ~1.5 mag in $J - K_s$ (Faherty et al. 2013). In particular, young L dwarfs associated with moving groups have redder NIR colors than field counterparts of the same spectral type, an effect attributed to lower surface gravities that retain clouds at high altitude (Marley et al. 2012). From a large survey of 420 ultracool dwarfs, Kellogg et al. (2017) found similar (~2%) fractions of

unusually red objects in both younger (<200 Myr) and older (≥200 Myr) populations, an age delineation corresponding to the halting of gravitational contraction within ultracool evolutionary models (e.g., Burrows et al. 2001; Baraffe et al. 2015) and roughly corresponding to the demarcation between very low and intermediate gravity indicators from Allers & Liu (2013). While the age and peculiarity of ultracool objects are difficult to precisely quantify, a significant population of reportedly older red L dwarfs exists (see Looper et al. 2008). Inclination angle has also been invoked as a potential explanation for discrepant red colors, as substellar objects viewed equator-on have redder infrared colors for the same
<table><caption>**Table 4**
Astrometry for HD 206893 B</caption><thead><tr><th>Epoch</th><th>Separation (mas)</th><th>Position Angle (deg)</th><th>Instrument, Band</th><th>Ref.</th></tr></thead><tbody><tr><td>2015-10-5</td><td>270.4 ± 2.6</td><td>69.95 ± 0.55</td><td>SPHERE/IRDIS, H</td><td>1</td></tr><tr><td>2016-8-8</td><td>269.0 ± 10.4</td><td>61.6 ± 1.9</td><td>VLT/NaCo, L'</td><td>1</td></tr><tr><td>2016-9-16</td><td>265 ± 2</td><td>62.25 ± 0.11</td><td>SPHERE/IRDIS, K1/K2</td><td>2</td></tr><tr><td>2016-9-22</td><td>267.6 ± 2.9</td><td>62.72 ± 0.62</td><td>GPI, H</td><td>3</td></tr><tr><td>2016-10-21</td><td>265.0 ± 2.7</td><td>61.33 ± 0.64</td><td>GPI, K1</td><td>3</td></tr><tr><td>2017-7-14</td><td>260.3 ± 2</td><td>54.2 ± 0.4</td><td>SPHERE/IRDIS, H</td><td>4</td></tr><tr><td>2017-11-9</td><td>256.9 ± 1.1</td><td>51.01 ± 0.35</td><td>GPI, K2</td><td>3</td></tr><tr><td>2018-6-20</td><td>249.11 ± 1.6</td><td>45.5 ± 0.37</td><td>SPHERE/IRDIS, H2/H3</td><td>4</td></tr><tr><td>2018-9-24</td><td>251.7 ± 5.4</td><td>42.6 ± 1.6</td><td>GPI, J</td><td>3</td></tr></tbody></table>

**References.** (1) Milli et al. (2017), (2) Delorme et al. (2017), (3) this work, (4) Grandjean et al. (2019).

spectral type than those viewed pole-on (Kirkpatrick et al. 2010; Vos et al. 2017).

The reasons for the persistence of redder colors in moderate-age, high-gravity objects remain uncertain. However, the presence of upper-atmosphere dusty aerosols of submicron-sized grains has been put forth as a potential explanation. In this framework, small, cool, dusty aerosols high in the atmosphere scatter (but do not emit) thermal flux. The sensitivity of scattering efficiency to wavelength for such particles is small compared to the wavelength of interest, thus producing a reddening. Marocco et al. (2014) and Hiranaka et al. (2016) examined such extinction effects in peculiar red L dwarf atmospheres by assuming a population of submicron grains in the upper atmosphere. Marocco et al. (2014) adopted a dereddening approach with wavelength-dependent corrections based on standard interstellar extinction laws (Cardelli et al. 1989; Fitzpatrick 1999) and showed that applying dereddening using common dust compositions could make peculiar red objects appear similar to the field L standard population. The HD 206893 B analysis by Delorme et al. (2017) identified the spectral type as L5–L7, concordant with its extremely red *Js* − K1 color, and they applied a similar method, generating extinction profiles with a range of *A*<sub>v</sub> (up to *A*<sub>v</sub> = 10) and particle sizes (from 0.05 to 1 µm). However, corresponding objects in the same region of the NIR CMD exhibit significantly deeper water absorption features at 1.4 µm and much bluer slopes than those evident in the spectrophotometry of HD 206893 B, regardless of youth. By comparing the dereddened spectrum of HD 206893 B to standard objects using a Cushing G analysis, Delorme et al. (2017) determined the closest matches to those of field or young L3.5 dwarfs, reproducing the slope of the SPHERE spectrophotometry.

In Hiranaka et al. (2016), extinction profiles generated from Mie scattering models with various grain sizes, size distributions, and compositions were fit to extinction profile estimates derived by dividing the spectra of unusually red L dwarfs by those of typical field L dwarfs. In this work, we apply the same methodology and compare the full NIR spectrum of HD 206893 B to known L dwarfs in order to investigate potential small dust grain populations that may be responsible for the observed reddening.

Figure 14 shows the full spectrum of HD 206893 B divided by L5, L6, and L7 spectral standards, selected as objects with median *J*–*Ks* colors for their spectral classes from the field gravity template sample in Cruz et al. (2018). Dividing the normalized spectrum of the standard by the normalized

companion spectrum yields a visualization of the observed reddening, shown as the flux ratio of the two objects. The reddening profile is then fit by grain populations with various scattering properties, which may be used, in turn, to deredden the companion spectrum. Initial treatments of dust in very low mass stellar atmospheres recognized that the condensation temperatures of species such as enstatite (MgSiO<sub>3</sub>), iron (Fe), forsterite (Mg<sub>2</sub>SiO<sub>4</sub>), and corundum (Al<sub>2</sub>O<sub>3</sub>) occurred in the photospheres of cool M dwarfs (e.g., Tsuji et al. 1996), and many additional grain compositions with cool condensation temperatures have since been incorporated into atmospheric models treating dust (e.g., ~30 various species of dust; Allard et al. 2001). As in Hiranaka et al. (2016), we use Mie scattering models and the refractive indices of various species to determine extinction coefficients. We adopt a standard power-law distribution $n(a) = a^{-3.5}$ with grain sizes from 0.1 to 1 µm and use the Mie scattering code LX-MIE (Kitzmann & Heng 2018) to calculate extinction curves for particles of forsterite, corundum, and TiO<sub>2</sub>. As all three species produced similar extinction profiles and grain parameters, we focus here on the results from the forsterite grain fitting.

Given the previous Delorme et al. (2017) analysis demonstrating the closest spectral approximation of HD 206893 B to a reddened atmosphere of spectral type L3 (consistent with both field and younger AB Dor member L3 objects), we perform extinction fitting with an extinction curve derived from the flux ratio of HD 206893 B to a low-gravity L3γ object (2MASSW J2208136+292121; Kirkpatrick et al. 2000)<sup>49</sup> and an unusually red field object with NIR spectral type L6.5pec (2MASSW J2244316+204343; Dahn et al. 2002; McLean et al. 2003; Looper et al. 2008). Fits to each of the extinction curves were performed using a Bayesian MCMC approach (the *emcee* package; Foreman-Mackey et al. 2013) to estimate particle column density, mean grain radius, and opacity scaling. The best-fit curves from the MCMC analysis are shown in Figure 15. For the earlier, very low gravity L3 spectral type, the corresponding posterior distributions are shown in Figure 16, with a best-fit column density of $2.8 \times 10^8$ particles cm<sup>-2</sup>, a mean particle radius of 0.27 µm, and a constant vertical offset term of $C = -1.9$, which corresponds to a gray atmospheric opacity scaling term between the object and the field spectrum. For the peculiar L6.5, the best-fit values have slightly lower

<footnote n="49"><sup>49</sup> We note that the L3γ object J2208 was identified as a β Pic candidate by Gagné et al. (2014), while an updated parallax and CMD analysis by Liu et al. (2016) show that it is marginally more similar to the field CMD sequence than that of very low gravity objects, making its youth determination slightly uncertain; however, the latter authors concluded that it remains a promising candidate member of the young moving group.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Posterior distributions for semimajor axis, eccentricity, inclination, longitude of the ascending node, mutual inclination of the companion orbit and circumstellar disk (with respect to the inferred disk geometry from Milli et al. 2017), and radius of apoapsis. The overplotted gray dashed line represents the prior distribution in mutual inclination, showing that the posterior distribution slightly favors more coplanar over more misaligned configurations.</figcaption></figure>

column density (2.5 × 10⁸ particles cm⁻²) and smaller mean particle radius (0.20 µm), with a constant vertical offset term of $C = -0.43$.

In Figure 17, we show the result of applying the best-fit extinction curve for forsterite to effectively “deredden” the full spectrum of HD 206893 B and compare it with both the low-gravity L3 object and the unusually red field L6. The dereddened spectra closely correspond to L6.5 over the *J* band and similarly to both L3 and L6.5 over the *H* band, albeit with a weaker water absorption feature at 1.4 µm, and more closely replicates the peaky *H*-band shape of the lower-gravity L3 than the broader morphology of L6.5. The shallowness of the water band is consistent with the presence of significant dust in the object atmosphere; for late M and L dwarfs, Leggett et al. (2001) noted that below $T_{\text{eff}} < 2500$ K, the presence of dust heats the atmosphere in the line-forming region, and the water features may become broader and shallower, depending upon the dust properties (e.g., metallicity). A detailed analysis of

such feedback effects between the posited dust layer and the atmosphere is beyond the scope of this investigation.

However, while the blue slope of *K*1 is roughly similar between the L3 and GPI spectra, the decline seen at the red end of *K*2 deviates sharply from the L3 object spectrum. In comparison, the dereddened spectrum assuming HD 206893 B is similar to an unusually red late field object (L6.5pec) provides a more similar match to the general morphology of the *K* band, with a less significant departure from the template object at the longest wavelengths in *K*2. As a slight red slope departure appears to persist after the dereddening application of a small particle extinction profile, if physical, its origin may be attributed to factors beyond a high-altitude aerosol layer.

## 6.3. A Possible Atmospheric Reddening Scenario for HD 206893 B

The high fractional infrared luminosity of HD 206893 A ($L_{\text{dust}}/L_* = 2.3 \times 10^{-4}$; Moór et al. 2006) and the existence
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Subsample of 250 randomly sampled MCMC orbits to HD 206893 B based on four epochs of GPI astrometric monitoring (blue diamonds), in addition to five epochs of VLT/SPHERE (black circles) and VLT/NaCo (red square) astrometry, covering a total time baseline from 2015 to 2018 (astrometric points and errors shown in inset). The darker orbits correspond to orbits with lower values of mutual inclinations between the companion orbit and the disk plane, with the mutual inclination distribution in Figure 10 peaking at ~20°.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** The 5σ contrast curves calculated for the pyKLIP-reduced GPI images in each band of observation, assuming L-type object spectra.</figcaption></figure>

of a companion interior to its debris disk make its architecture similar to companion-disk systems such as HD 95086 ($L_{dust}/L_* = 1.4 \times 10^{-3}$, companion at ~56 au; Chen et al. 2012; Rameau et al. 2016), HR 2562 ($L_{dust}/L_* = 1.1 \times 10^{-4}$, companion at 20 au; Moór et al. 2006; Konopacky et al. 2016), and, most recently, HD 193571 ($L_{disk}/L_* = 2.3 \times 10^{-5}$, companion at 11 au; Musso Barcucci et al. 2019). Of these systems, all of which have notably red companions, HD 206893 B remains the reddest, despite the fact that the system does not have the highest fractional infrared luminosity. As described in Section 6.2, additional dust contributions (e.g., in the form of high-altitude aerosols) are necessary to reconcile

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Combined sensitivity of all epochs and bands of GPI observations, expressed in terms of companion masses as a function of orbital separation. PyKLIP contrasts from Figure 12 are converted into physical parameters following the methodology of Nielsen et al. (2019) using an assumed age of 250 Myr and the CIFIST2011 BT-Settl atmosphere models (Caffau et al. 2011; Allard 2014; Baraffe et al. 2015) combined with the COND evolutionary sequences (Baraffe et al. 2003). The “divot” seen at 10 au and ~12 $M_{Jup}$ corresponds to the closer inner working angles accessible in only the shorter-wavelength observations. The color bar and scale show corresponding completeness levels, with the approximate mass and semimajor axis of HD 206893 B for reference (red point).</figcaption></figure>

the spectrum of HD 206893 B with that of field and low-gravity L dwarfs. Here we examine the plausibility of accretion from the debris disk as a potential source of dusty material in the HD 206893 B atmosphere.

As a simple approximation, the amount of dust required to redden the companion atmosphere can be compared to the debris disk properties derived from previous observations and modeling. We assume that any accreted dust remains solid and is not fully vaporized during the accretion process.<fnref n="50" /> From the derived extinction required to reconcile HD 206893 B's spectrum to that of a VLG L3 object in Section 6.2, the derived column density of ~3 × 10⁸ cm⁻² can be used to estimate the excess dust contribution in an aerosol layer. Assuming that the object radius is 1 $R_{Jup}$ and that high-altitude aerosols consist of forsterite particles of size ~0.3 µm and density 3.27 g cm⁻³ suggests the presence of ~9 × 10⁻¹⁰ $M_{Moon}$ of dust within the atmosphere. Previous analyses of the HD 206893 B debris disk fit the SED with both single- and two-component dust temperature models: the two-component dust temperature model included a hot (499 K) dust component at 0.8 au consisting of 1.1 × 10⁻⁶ $M_{Moon}$ of material and a colder (48 K) dust belt at 261 au of 5.6 × 10⁻¹¹ $M_{Moon}$ of material (Chen et al. 2014). The single-component models consisted of a blackbody of 54 K dust at 43 au, with an estimated mass of 0.030 $M_⊕$ (2.4 $M_{Moon}$) derived

<footnote n="50">⁵⁰ We note that the high temperature of an accretion shock similar to those seen in T Tauris could destroy dust, while the ablation and/or vaporization of incident meteoroids is often observed in the solar system (e.g., on the Earth and Mars; Hunten et al. 1980; Hartwick et al. 2019, respectively) and thought to seed high-altitude clouds. Also see the impact of comet Shoemaker-Levy 9, which produced debris aerosols in the Jovian atmosphere with particle sizes of 0.13–0.3 µm (West et al. 1995).</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Comparison of the full SPHERE+GPI spectrum of HD 206893 B with three late field L dwarfs (L5: 2MASS J06244595-4521548; L6: 2MASSI J1010148-040649; L7: 2MASSI J0825196+211552). The right panel plots the flux ratios of the standard spectra divided by that of HD 206893 B. The spike at ~1.55 µm is an artifact corresponding to the gap between the SPHERE IFS band edge and the onset of GPI H-band coverage.</figcaption></figure>

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

from 850 µm JCMT/SCUBA observations (Holland et al. 2017). Adopting a conservative limit of $1.1 \times 10^{-6} M_{\text{Moon}}$ of hot dust in the inner portion of the disk, we can roughly estimate the relative timescales of dust accretion onto the substellar companion.

Accretion rates of solids onto young planets have been estimated at ~$10^{-6} M_{\oplus} \text{yr}^{-1}$ for gas-rich disks ($8 \times 10^{-5} M_{\text{Moon}} \text{yr}^{-1}$; e.g., Alibert et al. 2018). Assuming this high solid accretion rate, the time required to collect enough material to redden HD 206893 B would be significantly less than 1 yr. However, the HD 206893 debris disk is significantly more evolved, with a much lower dust surface density than younger gas-rich disks. We therefore examine two limiting cases for potential accretion rates: (1) accretion of solids onto a giant planet based on orbital parameters and consistent with estimates for younger protoplanetary disks ($10^{-4} M_{\text{disk}} \text{orbit}^{-1}$; Paardekooper 2007), and (2) estimates of interplanetary dust flux at Jupiter derived from in situ measurements from Pioneer 10 and the New Horizons Student Dust Counter (Poppe 2016).

In the high accretion rate (protoplanetary) scenario, we assume an orbital period of ~30 yr, as estimated in Section 5.4; a dust mass of $1.1 \times 10^{-6} M_{\text{Moon}}$ in the vicinity of the companion (corresponding to the hot dust estimate from Chen et al. 2014); and a solid accretion rate of $10^{-4} M_{\text{disk}} \text{orbit}^{-1}$, which implies a mass accretion rate of $4 \times 10^{-12} M_{\text{Moon}} \text{yr}^{-1}$. Using the derived column density for the companion, this suggests that only a very short period of time (~260 yr) would be required to accrete sufficient dust from the environment and redden the atmosphere to the extent currently observed. In contrast, in the low-accretion solar system-like scenario, typical values of incident dust flux at Jupiter are on the order of $10^{-13} \text{ g m}^{-2} \text{ s}^{-1}$ of 0.5–100 µm material (Poppe 2016), corresponding to a much lower accretion rate of $6 \times 10^{-16} M_{\text{Moon}} \text{yr}^{-1}$ (for a Jovian-sized body) and requiring ~1 Myr to accumulate the estimated dust content in HD 206893 B.

The steady state reached by the atmosphere depends upon the lifetime of the accreted dust, which in turn depends upon the particle size-dependent fall speed and the strength of atmospheric eddy mixing. In lower-gravity atmospheres, the fall speed is lower and the effect of a given strength of eddy mixing is stronger, all else being equal, favoring longer dust

lifetimes. A complete analysis, accounting for the coupled problems of radiative heating of the atmosphere by the accreted dust and the dust sedimentation, will be considered in the future.

Either of the two accretion scenarios is plausible, considering the age of the system and potential replenishment of dust within the disk, and a significantly higher dust estimate at the ~10 au separation of HD 206893 B than the conservative estimate adopted here would increase the interception of grains from the environment into the atmosphere. If the companion indeed has low surface gravity, sustaining a small dust grain population at a high altitude in the atmosphere may be possible for extended periods, particularly if debris disk dust production replenishes the aerosol layer. In comparison to the red colors of the young/low-gravity substellar population, which are already postulated to result from the presence of thick high-altitude clouds, the even redder color of HD 206893 B could potentially be attributed to an additional source of reddening from the dusty disk environment in which it resides.

## 7. Conclusions

We present GPI spectroscopy of the substellar companion HD 206893 B obtained at the J, H, K1, and K2 bands. Consistent with the extraordinary red nature at $H - L'$ initially noted by Milli et al. (2017), the broader spectral coverage of the GPI further supports its exceptionally red nature. The overlapping wavelength regimes between the GPI and SPHERE observations shows excellent agreement with the SPHERE YJH observations presented by Delorme et al. (2017). The addition of full H and K1/K2 coverage made possible with GPI spectroscopy suggests that the companion may have low surface gravity and that its spectrum appears morphologically most similar to that of the young low-gravity late L dwarf population.

The GPI photometry of HD 206893 B consistently demonstrates its extraordinary position in color-magnitude space, with $H - Ks = 1.68 \pm 0.08$, significantly redder than the previously reddest substellar object, 2MASS J22362452 +4751425 b ($2\text{M}2236\text{b}$, $H - K = 1.26 \pm 0.18$; Bowler et al. 2017). From a comparison of the GPI spectra to brown dwarf spectral libraries, we find that the closest-matching
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Best-fit forsterite grain aerosol models (red) with the corresponding parameters from the MCMC analysis shown in Figure 16. The green curve represents the observed extinction for HD 206893 B when compared with a VLG L3 object, 2MASSW J2208136+292121 (left panel), and the blue curve represents that for a red field L6.5pec, 2MASSW J2244316+204343 (right panel).</figcaption></figure>

spectrum is that of a late-type, possibly low-gravity L dwarf, enabling comparison with other extremely red/peculiar objects like 2M2236b ("late L pec"), 2M1207b (M8.5–L4; Patience et al. 2010), and PSO J318.5–22 (L7; Liu et al. 2013).

We apply a dereddening approach akin to that of Hiranaka et al. (2016) to determine whether a small grain, submicron aerosol layer could reconcile the observed spectrum with that of field and low-gravity L dwarfs. We find that for reasonable column densities and grain properties, both a low-gravity L3 and “peculiar” red field L6.5 provide good matches to HD 206893 B, with the overall spectrum matching that of the later-type object but the spectral shape of features in, e.g., the *H* band more closely approximated by that of the lower-gravity object.

The emergent spectrum of HD 206893 B proves challenging to fit when conducting comparisons with a suite of various atmospheric model grids, owing to its enhanced luminosity at the *K* band relative to shorter wavelengths and a slightly unusual *K*2 morphology. Conducting a model grid fit to each of the GPI bands separately provides different companion parameters as compared to fitting the full *Y* to *L* spectral coverage of SPHERE, GPI, and NaCo. Each of the four model grids provided more internally consistent effective temperatures for the individual bands (ranging from 1400 to 1800 K, with an average across all bands and grids of ~1540 K), albeit over the full range of log(*g*) = 3.5–5.0, making the surface gravity of the object ambiguous.

We provide uncertainties on the stellar age ranging from 40 to 600 Myr, adopting 250 Myr for completeness analyses, but highlight that the high infrared luminosity of the disk and the nonnegligible (61%–63%) likelihood of membership in the Argus association point to a younger age for the system. As the mass of the companion depends on the assumed age, the companion mass ranges from 12 to 40 *M*<sub>Jup</sub> for ages in the 50–500 Myr range of COND models. The peaky morphology of the *H*-band spectra, good fit of the dereddened spectrum to that of a late-type L dwarf, morphological similarity of the spectrum to known young moving-group late L dwarfs, and potential lower dynamical mass of 10 *M*<sub>Jup</sub> estimated by Grandjean et al. (2019) all also point toward a self-consistent scenario for HD 206893 B being significantly younger and

lower-mass than initial age estimates and its luminosity have implied.

We have combined the five epochs of VLT astrometry presented in Milli et al. (2017), Delorme et al. (2017), and Grandjean et al. (2019) with four new multiband GPI epochs spanning 4 yr of orbital coverage. We estimate the orbital period to be 29.1<sup>+8.1</sup><sub>-6.7</sub> yr, with a probable semimajor axis of 10.4<sup>+1.8</sup><sub>-1.7</sub> au and orbital inclination of 145°.6<sup>+13.8</sup><sub>-6.6</sub>. From previous estimates of the debris disk inclination and inner gap radius of ~50 au (Milli et al. 2017), these data are consistent with the companion being well within the inner edge of the dust emission (as previously suggested by Delorme et al. 2017). The potential significant gap between the ~10 au best-fit semimajor axis and the 50 au inner disk radius has motivated the search for additional companions within the gap. However, no additional companions are detected in our images. The GPI data are 25%–75% complete to 5 *M*<sub>Jup</sub> at orbital separations of 20–30 au, respectively (assuming an older system age of 250 Myr). If the system is significantly younger, this may suggest that any additional companion responsible for carving the disk edge would be sub-Jovian. Comparing our MCMC fits to the visual orbit with the initial disk geometry estimates from Milli et al. (2017) favors less misaligned and more coplanar orbital configurations; however, the measurement precision on disk and orbital elements cannot currently exclude moderately misaligned (*i*<sub>m</sub> ~ 20°) configurations. As recent work has shown that *i*<sub>m</sub> correlates with the orbital period of companions in circumbinary systems (Czekala et al. 2019), with tighter binaries exhibiting greater coplanarity, this system provides a useful laboratory for disk–companion dynamics. Milli et al. (2017) offered the caveat that the marginally resolved nature of the debris disk from 70 µm Herschel observations does not tightly constrain the disk geometry, necessitating future disk observations at higher spatial resolution and longer-term orbital monitoring to explore the dynamical interplay between the companion and disk environment.

As the reddest substellar companion identified to date, HD 206893 B adds to the growing population of remarkably red planetary-mass and free-floating objects. With its exceptional color and close orbital separation, HD 206893 B is an important comparison for similarly red free-floating objects (e.g., PSO 318); red, wide substellar objects like 2M2236b; and
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Posterior distributions for the best-fit extinction profiles in Figure 15, with the left panel corresponding to the fit to the assumed L3 extinction profile and the right panel corresponding to the assumed L6.5pec extinction profile. Distributions shown include grain column density in $10^8$ particles cm$^{-2}$ (N), mean grain size in cm (a), vertical offset due to gray opacity (C), and an error tolerance parameter (ln(s)). The covariance of the mean grain size and column density reflect the inverse proportionality of the parameters, with similar extinction properties for small grain sizes.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Dereddened spectra of HD 206893 B shown with spectra of a VLG L3 object, 2MASSW J2208136+292121 (top panel), and a red field L6.5pec, 2MASSW J2244316+204343 (bottom panel). The dereddening approach is capable of producing realistic *J*, *H*, and *K*1-band spectral features in either case, albeit with too-shallow H$_2$O absorption, and is well approximated by the peaky *H*-band morphology of the lower surface gravity L3 object. However, the *K*2 band deviates significantly from that of an earlier L dwarf and is more similar to that of a late red field object with additional extinction.</figcaption></figure>

planetary-mass companions posited to have high photospheric dust content like HD 95086 b (De Rosa et al. 2016). The redness of the companion and its extremely dusty nature also make it well suited to exploring the effects of clouds, metallicity, and disequilibrium chemistry at cool temperatures in atmospheric models. Interpreting these observations, and future observations extended further into the infrared (e.g., with the James Webb Space Telescope), will likely be subject to how dusty clouds, collisionally induced absorption, and atmospheric chemistry are treated in models (Baudino et al. 2017), and unusual systems such as HD 206893 B will provide important laboratories for these processes. Furthermore, if the system is indeed young, this companion would present an exceptional addition to the growing population of directly imaged giant exoplanets. If the system is ~250 Myr, HD 206893 B represents only the second detection of a brown dwarf orbiting within the inner gap of its host debris disk (after HR 2562; Konopacky et al. 2016) and is otherwise one of only six substellar companion systems with both resolved disks and companions. This makes it a key testing ground to explore the dynamical interaction between disks, companions, and their host stars. The dustiness of the system also merits a search for polarized signal from the companion (e.g., as observed in companion systems like CT Cha and HD 142527; Rodigas et al. 2014; Ginski et al. 2018). With multiwavelength spectral coverage, detailed modeling, and continued dynamical monitoring, constraining the physical properties of this system will provide critical context for our understanding of the atmospheric and evolutionary histories of substellar objects.

We thank the anonymous referee for an exceptionally thoughtful and comprehensive report that greatly benefited this manuscript. K.W.D. thanks Alan Jackson for helpful conversations on debris disk dust properties, and Brian Svoboda, Sarah
Betti, Mike Petersen, Daniella Bardalez-Gagliuffi, Sarah Logsdon, and Emily Martin for helpful discussions and feedback. This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministerio da Ciencia, Tecnologia e Inovacao (Brazil). This research has made use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC02-05CH11231. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. ACI-1548562. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work was supported by NSF grants AST-1411868 (K.W.D., A.R., J.P., B.M., E.L.N., and K.B.F.), AST-141378 (G.D.), and AST-1518332 (R.D.R., J.J.W., T.M.E., J.R.G., and P.G.K.). Some authors were supported by NASA grants NNX14AJ80G (E.L. N., S.C.B., B.M., F.M., and M.P.), NNX15AC89G and NNX15AD95G (B.M., J.E.W., T.M.E., R.J.D.R., G.D., J.R. G., and P.G.K.), NN15AB52I (D.S.), and NNX16AD44G (K. M.M.). K.W.D. was supported by an NRAO Student Observing Support Award (SOSPA3-007). D.S.'s contribution was supported by NASA grant NNH15AZ59I. J.R., R.D., and D. L. acknowledge support from the Fonds de Recherche du Quebec. J.R.M.'s work was performed in part under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. Support for M.M.B.'s work was provided by NASA through Hubble Fellowship grant No. 51378.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.J. W. is supported by the Heising-Simons Foundation 51 Pegasi b postdoctoral fellowship. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. This work benefited from NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate.

# Appendix A
## Pipeline Postprocessing and Covariance Estimation

### A.1. pyKLIP Reduction

The images from the GPI autoreduced pyKLIP pipeline (Wang et al. 2018) are shown in Figure 18, with reduction parameters of nine annuli, four subsections per annulus, and an exclusion criterion corresponding to 1 pixel of movement,

using the first 10 KL modes to generate the reference PSF. The images show the close separation of the companion from the coronagraphic mask.

### A.2. K2 Spectral Extraction

As noted in Section 5.2, the red end of the GPI K2 spectrum shows a steep spectral slope beyond 2.2 µm. This slope was observed in both epochs of K2 data taken 1 yr apart and appears to be affected predominantly by the very low S/N in the GPI satellite spots at the red end of the band, which are critical to calibrate the measured flux in individual frames. The spectral slope is also impacted by the excessive brightness of the companion in these frames, the flux of which was significantly altered by the narrow high-pass filter. Therefore, the shape of the extracted spectrum varies significantly depending upon decisions made in postprocessing and PSF subtraction. We show the results of the 2016 CADI extraction in Figure 19, keeping in mind the strong bias of this spectrum due to the very low S/N of the satellite spots, compared with the 2017 CADI and LOCI extractions with a high-pass filter of 4 equivalent pixels. Reducing the 2016 data set with both LOCI and a broader high-pass filter resulted in the same spectrum as with cADI and a narrow filter (4 pixels). The 2017 epoch of data has a higher S/N owing to significantly more field rotation and longer integration, and because LOCI performed with a moderately large high-pass filter (10 pixels) recovered more flux from the extended wings of the PSF in the redder channels, we use the LOCI extraction throughout the analysis in this study.

### A.3. Covariance Estimation

Figure 20 shows the spectral covariance matrices derived from measuring interpixel correlation within the final PSF-subtracted LOCI and cADI images. Each covariance matrix was calculated independently for each GPI band by use of a parallel-tempered MCMC to fit the 18 parameters of spectral correlation in the IFS data cubes as outlined in Greco & Brandt (2016) and applied in De Rosa et al. (2016). Each MCMC was run with 128 walkers at 16 different temperatures and for 5000 steps saving every 10th step, burning in the first 1000 steps in the chains. Plotted are the spectral correlations as a function of wavelength channel for the final data cube in each GPI band, where the off-diagonal elements correspond to correlated noise terms. For each spectrum, the high- and low-frequency noise components, corresponding to read/background noise and speckle noise, respectively, are extracted and introduced into the error budget of the final spectra separately. Table 5 provides the fitted amplitudes of the correlated noise terms with respect to angular separation ($A_ρ$) and wavelength due to interpolation or crosstalk ($A_λ$). Detailed results from applying the spectral covariance to the model fitting of the spectra are provided in Section 5.3 and Appendix B.

# Appendix B
## Detailed Results for Model Atmospheric Grids

### B.1. Modeling Results with and without Covariances

Section 5.3 describes the modeling results implementing spectral covariance as summarized in Figure 20 and Table 5, and here we provide additional details on each of the model grids explored in the fitting of the HD 206893 B spectrum and
<figure><img src="image_22.png" /><figcaption>**Figure 18.** Autoreduced and pipeline-processed PSF-subtracted images of HD 206893 B in the *J* (upper left), *H* (upper right), *K*1 (lower left), and *K*2 (lower right) bands. Images shown are from the pyKLIP reductions, oriented north up and east left, with parameters of nine annuli, four subsections, 1 pixel movement criterion, and using the first 10 KL modes. Each individual image is scaled separately to minimize the contributions of residual speckle noise. A noticeable change in the companion position angle can be seen between the earliest 2016 epoch (upper right) and the latest 2018 epoch (upper left).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Comparisons of postprocessing approaches for the GPI K2-band data, demonstrating the recovery of flux at longer wavelengths. The 2017 November 9 LOCI extraction (black line) with a broader high-pass filter is used in this study.</figcaption></figure>

the results from a standard $\chi^2$ minimization. Table 6

summarizes the properties of the best-fit models to the full spectrum and all four bands, demonstrating the wide range of derived physical properties dependent upon the wavelength range, choice of models, and whether or not covariances were implemented. As noted in Section 5.3, the best-fit models accounting for spectral covariances do not always pass through the spectral data points in cases where the spectra are highly correlated; however, the resulting derived parameters across all four bands using covariances agree more closely than the wider range of temperatures, surface gravities, and radii estimated from a standard $\chi^2$ minimization.

## B.2. Model Descriptions

*Sonora models.* We conducted model comparisons with a subset of cloudy, solar-metallicity models from the upcoming Sonora 2020 model grid (M. S. Marley et al. 2020, in preparation). The Sonora models are applicable to brown dwarfs and giant exoplanets and incorporate revised solar abundances (Lodders 2010); rainout chemistry; updated
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Spectral correlation matrices calculated from the LOCI and cADI PSF-subtracted data cubes, with the correlation between pixel values at different wavelength slices shown as intensity. The greatest spectral covariance can be seen at the *H* and *K1* bands, with broad off-diagonal terms corresponding to a high correlation between adjacent wavelength channels.</figcaption></figure>

opacities for H₂, CH₄, and alkali species; and equilibrium and disequilibrium chemistry, and they span both cloudless and cloudy models with temperatures of 200–2400 K and metallicities of $-0.5 \le [M/H] \le 2$. The subset of solar-metallicity cloudy models tested here spans $T_{\text{eff}} = 1100–1600$ K and log (g) = 3.75–5.0. The cloudy Sonora models use the cloud model of Ackerman & Marley (2001), which is parameterized by the grain sedimentation efficiency, $f_{\text{sed}}$. This parameter sets the balance between grain sedimentation and upward motion from turbulent mixing. Small values of $f_{\text{sed}}$ correspond to slow sedimentation, smaller grains, and vertically extended clouds. The value of $f_{\text{sed}}$ was set equal to 1 (moderate particle settling efficiency) in the tested model subset, corresponding to a vertically extended thick cloud layer. Fitting over the full wavelength range favors a very low temperature (1200 K) object with high surface gravity and a physically plausible radius of ~1.36 $R_{\text{Jup}}$. The results from fitting individual bands favor slightly higher temperatures and generally lower surface gravities, suggesting that the low temperature over the full wavelength range is driven by both spectral correlation and increasing redness in the relative band-to-band flux. These models more closely approximate the emergent flux across shorter wavelengths and the shallow water absorption features between *J* and *H*, but they underestimate the flux at the longest wavelengths.

**DRIFT-PHOENIX models.** We also compared the spectra of HD 206893 B with the DRIFT-PHOENIX atmosphere models, which incorporate small grains in their mineral cloud physics prescription. The DRIFT-PHOENIX models are built upon the PHOENIX radiative transfer atmosphere models (Hauschildt 1992) but also incorporate microparticle physics through the DRIFT code (Dehn et al. 2007; Helling et al. 2008;

<caption>**Table 5**</caption>

Measured Amplitudes of the Correlated Noise Terms for Each GPI Data Cube ($A_\rho$ and $A_\lambda$), with the Relevant Correlation Lengths Measured for Each Band ($\sigma_\rho$, $\sigma_\lambda$)

<table><thead><tr><th>Band</th><th>ρ (mas)</th><th>A<sub>ρ</sub></th><th>A<sub>λ</sub></th></tr></thead><tbody><tr><td rowspan="7">J (σ<sub>ρ</sub> = 0.51, σ<sub>λ</sub> = 4.05)</td><td>200</td><td>0.70</td><td>0.25</td></tr><tr><td>225</td><td>0.80</td><td>0.14</td></tr><tr><td>250</td><td>0.57</td><td>0.33</td></tr><tr><td>265</td><td>0.51</td><td>0.38</td></tr><tr><td>280</td><td>0.75</td><td>0.14</td></tr><tr><td>300</td><td>0.89</td><td>0.02</td></tr><tr><td>325</td><td>0.89</td><td>0.01</td></tr><tr><td>350</td><td>0.86</td><td>0.01</td></tr><tr><td rowspan="7">H (σ<sub>ρ</sub> = 1.14, σ<sub>λ</sub> = 0.07)</td><td>200</td><td>0.59</td><td>0.37</td></tr><tr><td>225</td><td>0.60</td><td>0.34</td></tr><tr><td>250</td><td>0.56</td><td>0.39</td></tr><tr><td>265</td><td>0.70</td><td>0.23</td></tr><tr><td>280</td><td>0.68</td><td>0.25</td></tr><tr><td>300</td><td>0.55</td><td>0.37</td></tr><tr><td>325</td><td>0.37</td><td>0.51</td></tr><tr><td>350</td><td>0.04</td><td>0.81</td></tr><tr><td rowspan="7">K1 (σ<sub>ρ</sub> = 0.37, σ<sub>λ</sub> = 0.13)</td><td>200</td><td>0.02</td><td>0.93</td></tr><tr><td>225</td><td>0.73</td><td>0.18</td></tr><tr><td>250</td><td>0.73</td><td>0.14</td></tr><tr><td>265</td><td>0.63</td><td>0.22</td></tr><tr><td>280</td><td>0.66</td><td>0.16</td></tr><tr><td>300</td><td>0.74</td><td>0.01</td></tr><tr><td>325</td><td>0.64</td><td>0.05</td></tr><tr><td>350</td><td>0.63</td><td>0.01</td></tr><tr><td rowspan="7">K2 (σ<sub>ρ</sub> = 0.31, σ<sub>λ</sub> = 0.01)</td><td>200</td><td>0.83</td><td>0.10</td></tr><tr><td>225</td><td>0.82</td><td>0.10</td></tr><tr><td>250</td><td>0.79</td><td>0.10</td></tr><tr><td>265</td><td>0.78</td><td>0.10</td></tr><tr><td>280</td><td>0.77</td><td>0.08</td></tr><tr><td>300</td><td>0.65</td><td>0.17</td></tr><tr><td>325</td><td>0.54</td><td>0.22</td></tr><tr><td>350</td><td>0.51</td><td>0.26</td></tr></tbody></table>

**Note.** The separations within the data cube are selected to span the range of companion separations at various epochs of observation.

Witte et al. 2009), which includes the motion of particles within forming clouds and treats a variety of particle sizes and compositions. The DRIFT-PHOENIX models tested the ranges $T_{\text{eff}} = 1000–3000$ K and log(g) = 3.0–6.0 and metallicities of solar [M/H] = 0 or +0.3 (supersolar). For model comparisons with DRIFT-PHOENIX, solar-metallicity models generally produced better fits than supersolar, concordant with the metallicity, estimate for HD 206893 A from Delorme et al. (2017); only two individual GPI bands were best fit with [M/H] = +0.3 in either the standard $\chi^2$ or covariance fitting. The best model fit to the full spectrum favors a moderate temperature (1400 K) and surface gravity (log(g) = 4.5) and a plausible object radius (1.36$R_{\text{Jup}}$). The DRIFT-PHOENIX best-fit model most closely reproduces the general rise in flux from the bluest to reddest wavelengths in this study but underestimates the flux over *K*, in addition to shallower water absorption and overestimation of the *L'* photometry.

**Cloudy-AE60** (Madhusudhan et al. 2011) models. Given the red nature of the companion, we fit its spectra with the publicly available cloudy atmosphere models generated by Madhusudhan et al. (2011) for the HR 8799 planets, which are also anomalously red relative to the field population. The full model
<table><caption>**Table 6**
Best-fit Model Properties for Each of the Four Suites of Models, Both Using a Standard $\chi^2$ Minimization (Visualized in Figure 9) and when Incorporating Covariance Matrices</caption><thead><tr><th rowspan="2">Model Grid</th><th rowspan="2">Range</th><th colspan="4">Standard χ² Fitting</th><th colspan="4">Fitting Incorporating Covariance</th></tr><tr><th>T<sub>Eff</sub> (K)</th><th>log(g)</th><th>Radius (R<sub>Jup</sub>)</th><th>χ²</th><th>T<sub>Eff</sub> (K)</th><th>log(g)</th><th>Radius (R<sub>Jup</sub>)</th><th>χ²</th></tr></thead><tbody><tr><td rowspan="5">Sonora</td><td>Full</td><td>1100</td><td>5.0</td><td>1.68</td><td>4.4</td><td>1200</td><td>5.0</td><td>1.36</td><td>3.5</td></tr><tr><td>J</td><td>1200</td><td>5.0</td><td>1.19</td><td>0.3</td><td>1400</td><td>5.0</td><td>0.82</td><td>1.4</td></tr><tr><td>H</td><td>1300</td><td>3.7</td><td>1.06</td><td>0.9</td><td>1400</td><td>3.7</td><td>1.13</td><td>1.4</td></tr><tr><td>K1</td><td>1400</td><td>3.7</td><td>1.28</td><td>0.2</td><td>1400</td><td>3.7</td><td>1.28</td><td>2.7</td></tr><tr><td>K2</td><td>1400</td><td>3.7</td><td>1.28</td><td>0.2</td><td>1500</td><td>3.7</td><td>1.18</td><td>0.7</td></tr><tr><td rowspan="5">Cloudy-AE60</td><td>Full</td><td>1200</td><td>3.5</td><td>1.67</td><td>4.4</td><td>1700</td><td>3.5</td><td>1.55</td><td>8.6</td></tr><tr><td>J</td><td>1700</td><td>3.5</td><td>1.31</td><td>1.8</td><td>1700</td><td>3.5</td><td>1.25</td><td>5.2</td></tr><tr><td>H</td><td>1400</td><td>3.5</td><td>2.56</td><td>2.0</td><td>1700</td><td>3.5</td><td>1.77</td><td>2.4</td></tr><tr><td>K1</td><td>1700</td><td>3.5</td><td>2.68</td><td>0.5</td><td>1700</td><td>3.5</td><td>2.82</td><td>7.3</td></tr><tr><td>K2</td><td>1700</td><td>3.5</td><td>2.77</td><td>0.2</td><td>1700</td><td>3.5</td><td>2.80</td><td>1.1</td></tr><tr><td rowspan="5">DRIFT-PHOENIX</td><td>Full</td><td>1400</td><td>4.5</td><td>1.68</td><td>0.8</td><td>1400</td><td>4.5</td><td>1.36</td><td>1.2</td></tr><tr><td>J</td><td>1400</td><td>3.0</td><td>1.23</td><td>0.1</td><td>1500</td><td>4.5</td><td>1.33</td><td>0.9</td></tr><tr><td>H</td><td>1100*</td><td>5.0</td><td>3.01</td><td>0.2</td><td>1400</td><td>3.0</td><td>1.34</td><td>0.9</td></tr><tr><td>K1</td><td>1800*</td><td>3.0</td><td>0.81</td><td>0.2</td><td>1800*</td><td>3.5</td><td>0.80</td><td>2.5</td></tr><tr><td>K2</td><td>2000</td><td>5.0</td><td>0.71</td><td>0.2</td><td>1300*</td><td>5.0</td><td>1.63</td><td>0.6</td></tr><tr><td rowspan="5">BT-Settl</td><td>Full</td><td>1600</td><td>3.5</td><td>0.94</td><td>1.6</td><td>1600</td><td>3.5</td><td>0.96</td><td>1.2</td></tr><tr><td>J</td><td>1600</td><td>3.5</td><td>0.91</td><td>0.1</td><td>1600</td><td>3.5</td><td>0.88</td><td>0.9</td></tr><tr><td>H</td><td>1600</td><td>3.5</td><td>0.87</td><td>0.3</td><td>1500</td><td>3.0</td><td>1.12</td><td>1.1</td></tr><tr><td>K1</td><td>1800</td><td>4.5</td><td>0.76</td><td>0.1</td><td>1700</td><td>4.0</td><td>0.80</td><td>2.2</td></tr><tr><td>K2</td><td>1700</td><td>4.5</td><td>0.83</td><td>0.2</td><td>1500</td><td>3.5</td><td>1.28</td><td>0.7</td></tr></tbody></table>

Note. Asterisks indicate best-fit models from DRIFT-PHOENIX with supersolar metallicity ([M/H] = +0.3).

grid includes both cloudy and cloud-free atmospheres with a range of forsterite cloud thicknesses and optical depths. The Madhusudhan et al. (2011) modal grain sizes for the cloudy atmospheres are in the 60 µm size range (although similar results can be achieved within these models for ~1 µm Fe grains at 1% supersaturation, demonstrating some degeneracies in opacity that depend upon composition and grain size). Due to the red nature of HD 206893 B, we focus on the Cloudy-AE grid, which has the greatest vertical extent in cloud depth; this generally serves to suppress emergent flux at shorter wavelengths and produce higher fluxes at longer wavelengths for the same effective temperature. However, we find that the best-fit model incorporating covariances corresponds to low surface gravity (log(g) = 3.5) and the hottest temperature (1700 K) of the four grids but is remarkably consistent across the four separate band-to-band fits. The full spectral fit corresponds to a model that is significantly bluer than that of HD 206893 B, with the largest estimated corresponding radius (1.55 R<sub>Jup</sub>).

**BT-Settl models.** The wide range of potential masses and temperatures for the companion estimated by the other model grids motivates further spectral comparison with the BT-Settl (2015) models (Allard et al. 2012; Allard 2014), which cover larger ranges in these values. BT-Settl was designed to model the atmospheres of very low mass stars, brown dwarfs, and planets using the PHOENIX radiative transfer atmosphere models (Hauschildt 1992) in conjunction with updated molecular line lists for water, methane, ammonia, and CO<sub>2</sub> and revised solar abundances derived from the solar radiative hydrodynamical simulations from Asplund et al. (2009). The BT-Settl models include a cloud treatment involving supersaturation, which accounts for cloud growth and mixing, and aim to reproduce the stellar M–L and brown dwarf L–T transitions. To conduct a comparison with the HD 206893 B spectra, and given the likely solar metallicity of the host star,

we use the solar-metallicity BT-Settl (2015) models with grid parameters ranging from 1200 to 2050 K and surface gravities ranging from log(g) = 2.5 to 5.5. The resulting best-fit full spectral model shown in Figure 9 corresponds to a T<sub>eff</sub> of 1600 K, a low surface gravity (log(g) = 3.5), and a radius of 0.96 R<sub>Jup</sub>. Along with the best-fit DRIFT-PHOENIX model, this fit has the lowest reduced χ² incorporating spectral covariance. This model reproduces flux over the YJ bands well and most closely approximates the L' photometry of all four model suites, but it is insufficiently red and departs significantly from the K spectra.

## ORCID iDs

K. Ward-Duong https://orcid.org/0000-0002-4479-8291

K. Follette https://orcid.org/0000-0002-7821-0695

R. J. De Rosa https://orcid.org/0000-0002-4918-0247

J. Rameau https://orcid.org/0000-0003-0029-0258

M. Marley https://orcid.org/0000-0002-5251-2943

D. Saumon https://orcid.org/0000-0001-6800-3505

E. L. Nielsen https://orcid.org/0000-0001-6975-9056

A. Rajan https://orcid.org/0000-0002-9246-5467

A. Z. Greenbaum https://orcid.org/0000-0002-7162-8036

J. J. Wang https://orcid.org/0000-0003-0774-6502

I. Czekala https://orcid.org/0000-0002-1483-8811

G. Duchêne https://orcid.org/0000-0002-5092-6464

B. Macintosh https://orcid.org/0000-0003-1212-7538

S. Mark Ammons https://orcid.org/0000-0001-5172-7902

V. P. Bailey https://orcid.org/0000-0002-5407-2806

T. Barman https://orcid.org/0000-0002-7129-3002

J. Bulger https://orcid.org/0000-0003-4641-2003

C. Chen https://orcid.org/0000-0002-8382-0447

J. Chilcote https://orcid.org/0000-0001-6305-7272

T. Cotten https://orcid.org/0000-0003-0156-3019
T. M. Esposito https://orcid.org/0000-0002-0792-3719

M. P. Fitzgerald https://orcid.org/0000-0002-0176-8973

B. L. Gerard https://orcid.org/0000-0003-3978-9195

S. J. Goodsell https://orcid.org/0000-0002-4144-5116

P. Hibon https://orcid.org/0000-0003-3726-5494

J. Hom https://orcid.org/0000-0001-9994-2142

L.-W. Hung https://orcid.org/0000-0003-1498-6088

P. Ingraham https://orcid.org/0000-0003-3715-8138

P. Kalas https://orcid.org/0000-0002-6221-5360

Q. Konopacky https://orcid.org/0000-0002-9936-6285

J. E. Larkin https://orcid.org/0000-0001-7687-3965

F. Marchis https://orcid.org/0000-0001-7016-7277

C. Marois https://orcid.org/0000-0002-4164-4182

S. Metchev https://orcid.org/0000-0003-3050-8203

M. A. Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

R. Oppenheimer https://orcid.org/0000-0001-7130-7681

D. Palmer https://orcid.org/0000-0001-7128-0802

M. Perrin https://orcid.org/0000-0002-3191-8151

F. T. Rantakyrö https://orcid.org/0000-0002-9667-2244

B. Ren https://orcid.org/0000-0003-1698-9696

J.-B. Ruffio https://orcid.org/0000-0003-2233-4821

D. Savransky https://orcid.org/0000-0002-8711-7206

A. C. Schneider https://orcid.org/0000-0002-6294-5937

A. Sivaramakrishnan https://orcid.org/0000-0003-1251-4124

I. Song https://orcid.org/0000-0002-5815-7372

R. Soummer https://orcid.org/0000-0003-2753-2819

S. Thomas https://orcid.org/0000-0002-9121-3436

J. Kent Wallace https://orcid.org/0000-0001-5299-6899

S. Wiktorowicz https://orcid.org/0000-0003-4483-5037

S. Wolff https://orcid.org/0000-0002-9977-8255

## References

Abrahams, E., & Cruz, K. 2017, BDNYC/UltracoolTypingKit: UltracoolTypingKit v1.0.0, Zenodo, doi:10.5281/zenodo.840582

Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872

Alibert, Y., Venturini, J., Helled, R., et al. 2018, NatAs, 2, 873

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

Allard, F. 2014, in IAU Symp. 299, Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham (Cambridge: Cambridge Univ. Press), 271

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

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Bailey, V., Meshkat, T., Reiter, M., et al. 2014, ApJL, 780, L4

Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563

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

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Baudino, J.-L., Mollière, P., Venot, O., et al. 2017, ApJ, 850, 150

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, Proc. SPIE, 7014, 701418

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

Bowler, B. P., Liu, M. C., Mawet, D., et al. 2017, AJ, 153, 18

Burgasser, A. J. & the SPLAT Development Team 2017, arXiv:1707.00062

Burgasser, A. J., Geballe, T. R., Leggett, S. K., Kirkpatrick, J. D., & Golimowski, D. A. 2006, ApJ, 637, 1067

Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, RvMP, 73, 719

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

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

Chen, C. H., Mittal, T., Kuchner, M., et al. 2014, ApJS, 211, 25

Chen, C. H., Pecaut, M., Mamajek, E. E., Su, K. Y. L., & Bitner, M. 2012, ApJ, 756, 133

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Cruz, K. L., Núñez, A., Burgasser, A. J., et al. 2018, AJ, 155, 34

Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ, 883, 22

Dahn, C. C., Harris, H. C., Vrba, F. J., et al. 2002, AJ, 124, 1170

David, T. J., & Hillenbrand, L. A. 2015, ApJ, 804, 146

De Rosa, R. J., Nielsen, E. L., Blunt, S. C., et al. 2015, ApJL, 814, L3

De Rosa, R. J., Rameau, J., Patience, J., et al. 2016, ApJ, 824, 121

Dehn, M., Helling, C., Woitke, P., & Hauschildt, P. 2007, in IAU Symp. 239, Convection in Astrophysics, ed. F. Kupka, I. Roxburgh, & K. L. Chan (Cambridge: Cambridge Univ. Press), 227

Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79

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

Dotter, A. 2016, ApJS, 222, 8

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

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

Faherty, J. K., Rice, E. L., Cruz, K. L., Mamajek, E. E., & Núñez, A. 2013, AJ, 145, 2

Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016, ApJS, 225, 10

Filippazzo, J. C., Rice, E. L., Faherty, J., et al. 2015, ApJ, 810, 158

Fitzpatrick, E. L. 1999, PASP, 111, 63

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

Gagné, J., Faherty, J. K., Cruz, K. L., et al. 2015, ApJS, 219, 33

Gagné, J., Lafrenière, D., Doyon, R., Malo, L., & Artigau, É 2014, ApJ, 783, 121

Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, ApJ, 856, 23

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Galicher, R., Marois, C., Macintosh, B., et al. 2016, A&A, 594, A63

Ginski, C., Benisty, M., van Holstein, R. G., et al. 2018, A&A, 616, A79

Gizis, J. E., Faherty, J. K., Liu, M. C., et al. 2012, AJ, 144, 94

Grandjean, A., Lagrange, A. M., Beust, H., et al. 2019, A&A, 627, L9

Greco, J. P., & Brandt, T. D. 2016, ApJ, 833, 134

Güver, T., & Özel, F. 2009, MNRAS, 400, 2050

Hartwick, V. L., Toon, O. B., & Heavens, N. G. 2019, NatGe, 12, 516

Hauschildt, P. H. 1992, JQSR, 47, 433

Helling, C., Woitke, P., & Thi, W.-F. 2008, A&A, 485, 547

Hiranaka, K., Cruz, K. L., Douglas, S. T., Marley, M. S., & Baldassare, V. F. 2016, ApJ, 830, 96

Holland, W. S., Matthews, B. C., Kennedy, G. M., et al. 2017, MNRAS, 470, 3606

Hunten, D. M., Turco, R. P., & Toon, O. B. 1980, JAtS, 37, 1342

Johnson-Groh, M., Marois, C., De Rosa, R. J., et al. 2017, AJ, 153, 190

Kalas, P., Graham, J. R., Chiang, E., et al. 2008, Sci, 322, 1345

Kellogg, K., Metchev, S., Gagné, J., & Faherty, J. 2016, ApJL, 821, L15

Kellogg, K., Metchev, S., Geißler, K., et al. 2015, AJ, 150, 182

Kellogg, K., Metchev, S., Miles-Páez, P. A., & Tannock, M. E. 2017, AJ, 154, 112

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

Kirkpatrick, J. D., Cushing, M. C., Gelino, C. R., et al. 2011, ApJS, 197, 19

Kirkpatrick, J. D., Looper, D. L., Burgasser, A. J., et al. 2010, ApJS, 190, 100

Kirkpatrick, J. D., Reid, I. N., Liebert, J., et al. 2000, AJ, 120, 447

Kitzmann, D., & Heng, K. 2018, MNRAS, 475, 94

Konopacky, Q. M., Rameau, J., Duchêne, G., et al. 2016, ApJL, 829, L4

Konopacky, Q. M., Thomas, S. J., Macintosh, B. A., et al. 2014, Proc. SPIE, 9147, 914784

Kraus, A. L., & Ireland, M. J. 2012, ApJ, 745, 5

Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É 2007, ApJ, 660, 770

Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Sci, 329, 57

Lee, J., & Song, I. 2019, MNRAS, 486, 3434

Leggett, S. K., Allard, F., Geballe, T. R., Hauschildt, P. H., & Schweitzer, A. 2001, ApJ, 548, 908

Leggett, S. K., Golimowski, D. A., Fan, X., Geballe, T. R., & Knapp, G. R. 2003, in 12th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun: The Future of Cool-Star Astrophysics, ed. A. Brown, G. M. Harper, & T. R. Ayres (Boulder, CO: Univ. Colorado), 120

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

Liu, M. C., Magnier, E. A., Deacon, N. R., et al. 2013, ApJL, 777, L20

Lodders, K. 2010, ASSP, 16, 379

Looper, D. L., Kirkpatrick, J. D., Cutri, R. M., et al. 2008, ApJ, 686, 528

Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, PNAS, 111, 12661

Macintosh, B. A., Graham, J. R., Palmer, D. W., et al. 2008, Proc. SPIE, 7015, 701518

Madhusudhan, N., Burrows, A., & Currie, T. 2011, ApJ, 737, 34

Maire, J., Ingraham, P. J., De Rosa, R. J., et al. 2014, Proc. SPIE, 9147, 914785

Manjavacas, E., Bonnefoy, M., Schlieder, J. E., et al. 2014, A&A, 564, A55
Marley, M., Saumon, D., Morley, C., & Fortney, J. 2018, Sonora 2018: Cloud-free, Solar Composition, Solar C/O Substellar Atmosphere Models and Spectra, v1.0, Zenodo, doi:10.5281/zenodo.1309035

Marley, M. S., Saumon, D., Cushing, M., et al. 2012, ApJ, 754, 135

Marocco, F., Day-Jones, A. C., Lucas, P. W., et al. 2014, MNRAS, 439, 372

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006a, ApJ, 641, 556

Marois, C., Lafrenière, D., Macintosh, B., & Doyon, R. 2006b, ApJ, 647, 612

Marois, C., Macintosh, B., Barman, T., et al. 2008, Sci, 322, 1348

Marois, C., Macintosh, B., & Véran, J.-P. 2010, Proc. SPIE, 7736, 77361J

Matthews, B. C., Kennedy, G., Sibthorpe, B., et al. 2015, ApJ, 811, 100

Mawet, D., David, T., Bottom, M., et al. 2015, ApJ, 811, 103

Mawet, D., Riaud, P., Absil, O., & Surdej, J. 2005, ApJ, 633, 1191

McLean, I. S., McGovern, M. R., Burgasser, A. J., et al. 2003, ApJ, 596, 561

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

Moór, A., Ábrahám, P., Derekas, A., et al. 2006, ApJ, 644, 525

Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108

Musso Barcucci, A., Launhardt, R., Kennedy, G. M., et al. 2019, A&A, 627, A77

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

Nielsen, E. L., Rosa, R. J. D., Rameau, J., et al. 2017, AJ, 154, 218

Paardekooper, S. J. 2007, A&A, 462, 355

Patience, J., King, R. R., de Rosa, R. J., & Marois, C. 2010, A&A, 517, A76

Perrin, M. D., Ingraham, P., Follette, K. B., et al. 2016, Proc. SPIE, 9908, 990837

Poppe, A. R. 2016, Icar, 264, 369

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

Rameau, J., Nielsen, E. L., De Rosa, R. J., et al. 2016, ApJL, 822, L29

Rieke, G. H., & Lebofsky, M. J. 1985, ApJ, 288, 618

Rodigas, T. J., Follette, K. B., Weinberger, A., Close, L., & Hines, D. C. 2014, ApJL, 791, L37

Schmidt, T. O. B., Neuhäuser, R., Seifahrt, A., et al. 2008, A&A, 491, 311

Schneider, A. C., Windsor, J., Cushing, M. C., Kirkpatrick, J. D., & Wright, E. L. 2016, ApJL, 822, L1

Sivaramakrishnan, A., & Oppenheimer, B. R. 2006, ApJ, 647, 620

Smart, R. L., Marocco, F., Sarro, L. M., et al. 2019, MNRAS, 485, 4423

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

Thalmann, C., Grady, C. A., Goto, M., et al. 2010, ApJL, 718, L87

Tokunaga, A. T., Simons, D. A., & Vacca, W. D. 2002, PASP, 114, 180

Torres, C. A. O., Quast, G. R., Melo, C. H. F., & Sterzik, M. F. 2008, in Handbook of Star Forming Regions, ed. B. Reipurth, Vol. II (San Francisco, CA: ASP), 757

Tsuji, T., Ohnaka, K., & Aoki, W. 1996, A&A, 305, L1

Vos, J. M., Allers, K. N., & Biller, B. A. 2017, ApJ, 842, 78

Wang, J. J., Perrin, M. D., Savransky, D., et al. 2018, JATIS, 4, 018002

Wang, J. J., Rajan, A., Graham, J. R., et al. 2014, Proc. SPIE, 9147, 914755

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, v2.0, Astrophysics Source Code Library, ascl:1506.001

West, R. A., Karkoschka, E., Friedson, A. J., et al. 1995, Sci, 267, 1296

Witte, S., Helling, C., Barman, T., Heidrich, N., & Hauschildt, P. H. 2011, A&A, 529, A44

Witte, S., Helling, C., & Hauschildt, P. H. 2009, A&A, 506, 1367

Wolff, S. G., Perrin, M. D., Maire, J., et al. 2014, Proc. SPIE, 9147, 91477H

Wu, Y.-L., Close, L. M., Males, J. R., et al. 2015, ApJ, 801, 4

Yee, J. C., & Jensen, E. L. N. 2010, ApJ, 711, 303

Zuckerman, B. 2019, ApJ, 870, 27

Zuckerman, B., & Song, I. 2004, ApJ, 603, 738