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

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 × 10⁻⁴; 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}^{{\prime} }={13.43}_{-0.15}^{+0.17}$ mag relative to the initial measurement of H = 16.79 ± 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 ∼ 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}_{-4}^{+5}{M}_{\mathrm{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σ 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.

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 flat-fielded, 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 ± 0.007 mas pixel⁻¹ and −010 ± 013 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 × 10⁻⁴ for J, 2.035 × 10⁻⁴ for H, 2.695 × 10⁻⁴ for K1, and 1.905 × 10⁻⁴ 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 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× 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.

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 × 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 χ² 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 spectrophotometric points for model fitting, as described in Section 5.3.

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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

The resulting combined spectrum is shown in Figure 6 with each of the extractions from the reduced, postprocessed J, H, K1, and K2 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 K2 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 K2 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 K2 spectrum is used.

Zoom In Zoom Out Reset image size

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 K2 were analyzed using the SpeX Prism Library Analysis Toolkit (SPLAT; Burgasser & the SPLAT Development Team 2017). Spectral typing was performed within SPLAT with a χ² 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 ± 1 subclasses. While no single spectrum provides a close morphological match across the full Y–K2 range, the three closest matches within the spectral library were late-type, low-gravity objects, all with reduced χ² values of ∼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 χ² 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; ⁴⁴ 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 dissimmilarities 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 simliar 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 ± 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ⁴⁵ from Marley et al. 2018); (2) the Cloudy-AE ⁴⁶ models, originally developed for the HR 8799 planets (Madhusudhan et al. 2011); (3) the DRIFT-PHOENIX ⁴⁷ atmosphere models (Helling et al. 2008; Witte et al. 2009, 2011); and (4) the BT-Settl (2015) ⁴⁸ grid (Allard et al. 2012; Allard 2014). For each grid, individual models spanning ranges of T_eff 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_Jup in steps of 0.01. The resulting χ² 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 χ² approach are also shown on a 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.

Zoom In Zoom Out Reset image size

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_Jup. 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. Accounting 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 broad range of surface gravities and radii depending on the band in question.

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\buildrel{\circ}\over{.} 6}_{-6\buildrel{\circ}\over{.} 6}^{+13\buildrel{\circ}\over{.} 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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&\cos ({i}_{m})=\cos ({i}_{c})\cos ({i}_{d})+\sin ({i}_{c})\sin ({i}_{d})\cos ({{\rm{\Omega }}}_{c}-{{\rm{\Omega }}}_{d}),\end{eqnarray} \tag{ 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 visual orbit are not sufficient to exclude moderately misaligned configurations (i_m  ∼ 20°).

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 periapse 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 ∼50 au by Milli et al. (2017).

Zoom In Zoom Out Reset image size

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⁴ 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 ∼20 au and in the brown dwarf regime down to ∼10–20 M_Jup at ∼10 au. For comparison, the inner additional companion signal predicted from RV variation by Grandjean et al. (2019) suggests an ∼15 M_Jup object at an orbital separation of 1.4–2.6 au. The requisite contrast to detect such an object (<10⁻⁴ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the existence of circumstellar material in the HD 206893 system, dust extinction may have an effect on 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 circumsecondary 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.21 × 10²¹ 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 × 10⁻⁴ g cm⁻² for A_V  = 10. For comparison, visual extinction for the edge-on AU Mic debris disk system with higher fractional infrared luminosity (3.9 × 10⁻⁴; 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⁻³) yields a particle column density of 2 × 10⁸ grains cm⁻², suggesting that a debris disk with A_V  = 10 is unrealistically optically thick (τ ∼ 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 argues against a significant reddening effect due to interstellar extinction, leaving dust in the companion atmosphere as the most viable scenario.

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 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_v (up to A_v  = 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 be 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₃), iron (Fe), forsterite (Mg₂SiO₄), and corundum (Al₂O₃) 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₂. As all three species produced similar extinction profiles and grain parameters, we focus here on the results from the forsterite grain fitting.

Zoom In Zoom Out Reset image size

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) ⁴⁹ 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 × 10⁸ particles cm⁻², 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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

However, while the blue slope of K1 is roughly similar between the L3 and GPI spectra, the decline seen at the red end of K2 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 K2. 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.

The high fractional infrared luminosity of HD 206893 A (L_dust/L_* = 2.3 × 10⁻⁴; Moór et al. 2006) and the existence 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 × 10⁻³, companion at ∼56 au; Chen et al. 2012; Rameau et al. 2016), HR 2562 (L_dust/L_* = 1.1 × 10⁻⁴, companion at 20 au; Moór et al. 2006; Konopacky et al. 2016), and, most recently, HD 193571 (L_disk/L_* = 2.3 × 10⁻⁵, 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 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. ⁵⁰ 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 from 850 μm JCMT/SCUBA observations (Holland et al. 2017). Adopting a conservative limit of 1.1 × 10⁻⁶ M_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⁻⁶ M_⊕ yr⁻¹ for gas-rich disks (8 × 10⁻⁵ M_Moon yr⁻¹; 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⁻⁴ M_disk 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 × 10⁻⁶ M_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⁻⁴ M_disk orbit^–1, which implies a mass accretion rate of 4 × 10⁻¹² M_Moon yr⁻¹. 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⁻¹³ g m⁻² s⁻¹ of 0.5–100 μm material (Poppe 2016), corresponding to a much lower accretion rate of 6 × 10⁻¹⁶ M_Moon yr⁻¹ (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.

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 the results from a standard χ² 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 χ² minimization.

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 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 ≤ [M/H] ≤ 2. The subset of solar-metallicity cloudy models tested here spans T_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_sed. This parameter sets the balance between grain sedimentation and upward motion from turbulent mixing. Small values of f_sed correspond to slow sedimentation, smaller grains, and vertically extended clouds. The value of f_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_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; 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_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 χ² 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.36R_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 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_Jup).

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₂ 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_eff of 1600 K, a low surface gravity (log(g) = 3.5), and a radius of 0.96 R_Jup. 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.