Uniform Forward-Modeling Analysis of Ultracool Dwarfs. I. Methodology and Benchmarking
Zhoujian Zhang, Michael C. Liu, Mark S. Marley, Michael R. Line, William M. J. Best
DD RAFT VERSION N OVEMBER
26, 2020Typeset using L A TEX preprint style in AASTeX61
Uniform Forward-Modeling Analysis of Ultracool Dwarfs.I. Methodology and Benchmarking Z HOUJIAN Z HANG ( 张 周 健 ),
1, 2 M ICHAEL
C. L IU , M ARK
S. M
ARLEY , M ICHAEL
R. L
INE , AND W ILLIAM
M. J. B
EST Institute for Astronomy, University of Hawaii at Manoa, Honolulu, HI 96822, USA Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract80HQTR19D0030 with the National Aeronautics and Space Administration. NASA Ames Research Center, Mail Stop 245-3, Moffett Field, CA 94035, USA School of Earth & Space Exploration, Arizona State University, Tempe AZ 85287, USA Department of Astronomy, University of Texas at Austin, Austin, Texas 78712, USA (Received July 31, 2020)
Submitted to ApJABSTRACTWe present a forward-modeling framework using the Bayesian inference tool Starfish and cloudlessSonora-Bobcat model atmospheres to analyze low-resolution ( R ≈ − . − . µ m)spectra of T dwarfs. Our approach infers effective temperatures, surface gravities, metallicities, radii, andmasses, and by accounting for uncertainties from model interpolation and correlated residuals due to in-strumental effects and modeling systematics, produces more realistic parameter posteriors than traditional( χ -based) spectral-fitting analyses. We validate our framework by fitting the model atmospheres them-selves and finding negligible offsets between derived and input parameters. We apply our methodology tothree well-known benchmark late-T dwarfs, HD 3651B, GJ 570D, and Ross 458C, using both solar andnon-solar metallicity atmospheric models. We also derive these benchmarks’ physical properties using theirbolometric luminosities, their primary stars’ ages and metallicities, and Sonora-Bobcat evolutionary mod-els. Assuming the evolutionary-based parameters are more robust, we find our atmospheric-based, forward-modeling analysis produces two outcomes. For HD 3615B and GJ 570D, spectral fits provide accurate T eff and R but underestimated log g (by ≈ . Z (by ≈ .
35 dex), likely due to the systematics frommodeling the potassium line profiles. For Ross 458C, spectral fits provide accurate log g and Z but overes-timated T eff (by ≈
120 K) and underestimated R (by ≈ . × ), likely because our model atmospheres lackclouds, reduced vertical temperature gradient, or disequilibrium processes. Finally, the spectroscopicallyinferred masses of these benchmarks are all considerably underestimated. INTRODUCTIONThe past two decades have witnessed a wealth of exoplanet discoveries. Characterizing such a large anddiverse population is essential to understanding their formation pathways, evolutionary stages, and habit-ability. Emission (e.g., Line et al. 2016; Samland et al. 2017a), transmission (e.g., Deming et al. 2013;Kreidberg et al. 2014), and reflection (e.g., Parmentier et al. 2016; MacDonald et al. 2018) spectra and a r X i v : . [ a s t r o - ph . S R ] N ov photometry are all useful for characterizing the thermal structure, cloud properties, and composition of exo-planet atmospheres, which are crucial in shaping the appearance and evolution of exoplanets (e.g., Burrowset al. 2001; Fortney et al. 2008b; Marley & Robinson 2015). However, robust spectroscopic analyses ofexoplanets are often hampered by low signal-to-noise (S/N) data due to the planets’ faintness and smallplanet-to-star area ratios, as well as by data with fairly narrow and/or non-contiguous wavelength coverage.Fortunately, non-irradiated, self-luminous brown dwarfs ( ≈ −
70 M
Jup ; Spiegel et al. 2011; Dupuy & Liu2017) overlap in temperature and surface gravity to imaged and transiting planets and are readily viable forhigh-quality emission spectroscopy. Investigating the properties of brown dwarf atmospheres can thereforeprovide insight into the similar atmospheric chemical and physical processes operating in exoplanets.Studying atmospheres of brown dwarfs has commonly relied on forward modeling (e.g., Saumon et al.2000; Geballe et al. 2001; Burgasser et al. 2006; Cushing et al. 2008, 2011; Stephens et al. 2009; Liuet al. 2011; Bonnefoy et al. 2014, 2018; Manjavacas et al. 2016), i.e., comparing the observed spectrawith a grid of theoretical model spectra, a.k.a. “grid models”, which are pre-computed with a few freeparameters (e.g., effective temperature, surface gravity, and metallicity) and self-consistent assumptions(including radiative-convective equilibrium, [dis]equilibrium chemistry, and simplified cloud properties).Such models are essential for characterizing the properties of brown dwarfs and exoplanets, developing ourunderstanding of ultracool atmospheres, and planning ground- and space-based observations.Traditional forward-modeling analysis compares the observed spectrum with the model spectrum at eachgrid point and then typically derives best-fit parameters with a least-squares approach (e.g., Cushing et al.2008; Rice et al. 2010; Zhang et al. 2020). This process usually assumes flux measurements follow Gaussianstatistics and the (data − model) residuals between different wavelengths are independent. These assump-tions lead to a diagonal covariance matrix when evaluating free parameters, with inverse weights (usuallyassumed as squared flux uncertainties) placed along the diagonal axis. In addition, since model grids arecoarsely created, interpolation of models is required to evaluate parameters in between grid points.However, this traditional method has two limitations. First, model interpolation is an important sourceof parameter uncertainties, but this is usually ignored. A typical approach is to use linear interpolation todetermine the model spectrum for physical parameters in between model grid points, which returns a singlespectrum with no uncertainty. This can result in artificially small uncertainties for inferred physical param-eters and may bias parameter posteriors toward the model grid points (e.g., Cottaar et al. 2014; Czekalaet al. 2015; Zhang et al. 2020). Second, residuals from the data-model comparison are correlated betweendifferent wavelengths, given that (1) spectrographs are designed to over-sample the instrumental line spreadfunction, so fluxes recorded on adjacent pixels are correlated, and (2) the number of free parameters for gridmodels is usually too small to comprehensively describe real data (as many other free parameters that de-scribe, e.g., cloud sedimentation efficiency and vertical mixing of chemical species, are often pre-assumedwithin the grid models). Failure to account for interpolation uncertainties and correlated residuals can causethe resulting derived parameters to have underestimated errors, which in practice are often masked simplyby adopting a fraction of the model grid spacing for the final parameter uncertainties (e.g., Leggett et al.2007; Cushing et al. 2008; Stephens et al. 2009; Zhang et al. 2020).Czekala et al. (2015; C15) have developed a new Bayesian framework for forward-modeling analysisof stellar spectroscopy, dubbed Starfish , which can mitigate limitations of traditional methods. To syn-thesize spectra in between grid points, Starfish is equipped with a “spectral emulator”, which generates a https://github.com/iancze/Starfish . probability distribution of spectra at the desired grid parameters and propagates the associated interpolationuncertainties into the resulting inferred parameters. Also, when evaluating physical parameters, Starfishconstructs a covariance matrix with off-diagonal components to account for correlated residuals causedby the instrumental effect and model imperfections. Therefore, Starfish derives larger but more realisticparameter uncertainties than the traditional spectral-fitting technique.In this work, we apply Starfish to the spectroscopic analysis of three benchmark brown dwarfs, HD 3651B(T7.5; Mugrauer et al. 2006), GJ 570D (T7.5; Burgasser et al. 2000), and Ross 458C (T8; Goldman et al.2010), using the cloudless Sonora-Bobcat model grids (Marley et al. 2017; Marley et al., in prep) with bothsolar and non-solar metallicities. These 3 objects are common proper-motion companions to stars at wideorbital separations. Hence, their distances, metallicities, and ages can be independently obtained from theirprimary hosts. Being late-T dwarfs, these objects are believed to be free of optically thick clouds in theatmospheric extent probed by near-infrared wavelengths (although optically thin clouds might exist; e.g.,Morley et al. 2012). Therefore, they constitute an optimal sample for testing cloud-free models, which serveas the starting point for the more complex cloudy models representative of directly imaged exoplanets (e.g.,Marois et al. 2008; Bowler et al. 2017; Chauvin et al. 2017).We start with a brief review of the three benchmarks and observations (Section 2) and the cloudlessSonora-Bobcat models (Section 3). We then construct and validate our forward-modeling framework usingthe Sonora-Bobcat atmospheric models, infer physical properties of our sample, and assess systematics ofthe model atmospheres (Section 4). We also derive the objects’ physical properties using their ages andmetallicities and the Sonora-Bobcat evolutionary models (Section 5), and we compare these properties withatmospheric-based spectral fits to test the cloudless model predictions (Section 6). Finally, we provide asummary and an outlook of future work (Section 7). BENCHMARK COMPANIONS AND OBSERVATIONS2.1.
HD 3651B
HD 3651B (T7.5) is the first directly-imaged brown dwarf companion to an exoplanet host, found inde-pendently by Mugrauer et al. (2006) and Luhman et al. (2007). The companion is at 43 (cid:48)(cid:48) (479 au) projectedseparation to HD 3651A (K0V), which has a
Gaia
DR2 distance of 11 . ± .
007 pc (Bailer-Jones et al.2018) and hosts a Saturn-mass planet ( M sin i = . Jup ) on a close ( a = . e = . R ≈ × ), and here we adopt Z = . − . −
13 Gyr based ona variety of stellar parameters and model sets. Gyrochronology suggests an age of 6 . − . R (cid:48) HK or the X-ray emission index log R X ≡ log L X / L bol using various activity-age relations(e.g., Donahue 1993; Soderblom et al. 1991; Gaidos 1998; Mamajek & Hillenbrand 2008). We re-computeages from these measured activity indices using the Mamajek & Hillenbrand (2008) R (cid:48) HK -age and R X -agerelations, leading to an age of 2 . − . draws from each age range assuming a uniformdistribution. We then combine these draws, compute their 16th and 84th percentiles, and adopt an age of4 . − . GJ 570D
GJ 570D (T7.5) was found by Burgasser et al. (2000) in a hierarchical quadruple system with GJ 570A(K4V) and a M1.5V + M3V close binary GJ 570BC (Duquennoy & Mayor 1988). This system is located at5 . ± .
003 pc (Bailer-Jones et al. 2018). The projected separations of A–BC, B–C, and A–D componentsare 24 . (cid:48)(cid:48) (145 au), 0 . (cid:48)(cid:48) (0.9 au), and 258 . (cid:48)(cid:48) (1520 au), respectively.Table 2 summarizes the metallicity and age of GJ 570A (also see Liu et al. 2007). We adopt a metallicityof Z = − .
05 to 0 .
05 dex based on high-resolution spectroscopy ( R ≈ × ). The age of GJ 570A hasbeen estimated as 1 − . − . R (cid:48) HK -age and R X -age relations, leading to an age of 0 . − . . − . Ross 458C
Ross 458C (T8) was found by Goldman et al. (2010) and Scholz (2010) as a 102 (cid:48)(cid:48) (1175 au) companion tothe M0.5Ve + M7Ve binary Ross 458AB (separated by 6 au; Beuzit et al. 2004), located at 11 . ± .
020 pc (Bailer-Jones et al. 2018). Montes et al. (2001) noted the space motion of Ross 458AB is close to theHyades cluster, but the object’s Hyades membership was later rejected by Burgasser et al. (2010), whofound a significant difference between the updated space motions of Ross 458AB and the Hyades. UsingBANYAN Σ (version 1.0; Gagné et al. 2018) and the Gaia
DR2 proper motion, parallax, and radial velocity,we find Ross 458AB is a field dwarf with a 99 .
9% probability.The metallicity and age of the Ross 458 system were summarized by Burgasser et al. (2010), which webriefly describe here. Based on V - and K -band photometry, Burgasser et al. (2010) derived a photometricmetallicity of 0 . − . . ± .
08 dex by Gaidos & Mann (2014) using the moderate-resolution ( R ≈ Z = . − .
33 dex. The fast rotation of Ross 458A (1 . − < α emission and a very bright X-ray luminosity, suggests an age of < . > .
15 Gyr. We therefore follow Burgasser et al. (2010) andadopt an age of 0 . − . Observations
The near-infrared spectra of the three benchmarks were all observed using the SpeX spectrograph (Rayneret al. 2003) mounted on the NASA Infrared Telescope Facility (IRTF). We obtain spectra of HD 3651B andGJ 570D from the SpeX Prism Library (Burgasser 2014) as observed by Burgasser (2007) and Burgasseret al. (2004), respectively . We observed Ross 458C using IRTF/SpeX in prism mode on 2015 July 07 (UT).We took 30 exposures with 120 seconds each for the target in a ABBA pattern and contemporaneously We note the
Gaia
DR2 astrometry of Ross 458AB is consistent with
Gaia ’s single-star model and unlikely affected by theorbital motion of the close binary, given its renormalized unit weight error (RUWE = .
09) is smaller than 1.4 (the quality cutsuggested by Lindegren 2018) and its astrometric excess noise is 0 mas. The IRTF/SpeX spectrum of HD 3651B was also observed by Liu et al. (2007). observed a nearby A0V standard star HD 116960 for telluric correction. We reduced the data using version4.1 of the Spextool software package (Cushing et al. 2004).Spectra of all these objects were observed with the SpeX 0 . (cid:48)(cid:48) × (cid:48)(cid:48) slit and thus have same spectral reso-lution ( R ≈ − H MKO magnitudes and the WFCAM H -band filter response (Hewett et al. 2006) and the zero-point flux (Lawrence et al. 2007). Astrometry andphotometry of the three benchmarks are listed in Table 3. THE CLOUDLESS SONORA-BOBCAT MODELS3.1.
Model Description
Throughout this work, we use the cloudless Sonora-Bobcat models (Marley et al. 2017; Marley et al., inprep). The models are 1D radiative-convective equilibrium atmosphere models that iteratively solve for aself-consistent atmospheric structure given a specified gravity and effective temperature. The model atmo-sphere code was originally developed for modeling Titan’s atmosphere (McKay et al. 1989) and was subse-quently modified by Marley & McKay (1999) to study the atmosphere of Uranus. It has since been appliedto the study of brown dwarfs (e.g., Marley et al. 1996; Burrows et al. 1997; Cushing et al. 2008; Saumon& Marley 2008), and solar and extrasolar giant planets (e.g., Marley et al. 1999, 2002, 2012; Fortney et al.2005, 2008a, 2013). The modeling framework is described in Marley & Robinson (2015). Chemical equi-librium is computed by interpolation of pre-computed tables using a modified version of the NASA CEAGibbs minimization code (see Gordon & McBride 1994), computed accounting for rainout based upon priorthermochemical models of substellar atmospheres (e.g., Lodders 1999; Visscher et al. 2010; Moses et al.2013). After a thermal structure solution in chemical and radiative-convective equilibrium has been found,the profile can be post-processed to yield model spectra at arbitrary spectral resolution. The resulting mod-els span a parameter space of 200 − T eff ), 3 . − . g ), and − .
5, 0, and 0 . Z ).3.2. Synthetic Spectra
In order to understand the effects of physical parameters on the resulting spectra, we compare three sets ofsynthetic spectra with varying T eff from 600 K to 1200 K, log g from 3 . . Z from − . + . . (cid:48)(cid:48) slit of the SpeX prism (Rayner et al. 2003). We thennormalize all model spectra using their peak J -band flux to highlight the variations in spectral shape.To compare the spectral variations among different grid parameters, we differentiate the synthetic spectra(which have been already smoothed and normalized by their peak J -band fluxes) with respect to the physicalparameters, by computing the average flux increase at each wavelength relative to a spectrum at { T eff =
900 K , log g = . , Z = } , when T eff is increased by 10 K, log g by 0 . Z by 0 . Oand CH to become broader and deeper, due to the increasing column densities of absorbers (e.g., Zahnle& Marley 2014). This can also been seen from Figure 2, where the spectral derivative in terms of T eff ispositive across almost the entire near-infrared wavelength range. In general, the cooling of atmospheresleads to a series of flux peaks at Y , J , H , and K bands, which symbolize the brown dwarf evolution fromearlier to later spectral types.As shown in the middle panel of Figure 1, changing surface gravity causes strong spectral variations at K band and the Y -band peak relative to the J -band peak. The K -band flux is suppressed at high gravity, dueto enhanced absorption from collisions between H molecules and the H and He atoms (e.g., Lenzuni et al.1991). In Y band, the spectral variation due to gravity is not significant for log g (cid:38) .
5, in agreement with thespectra of benchmark substellar companions whose surface gravities are independently derived from agesof their primary stars (e.g., HD 3651B and GJ 570D; Liu et al. 2007; Leggett et al. 2007). Figure 1 showsthat low-gravity objects with log g ≈ . Y − J colors than higher-gravityobjects (log g (cid:38) . g = . (cid:46) . −
70 M
Jup , which are muchyounger than current search efforts, or (2) a mass of 1 − Jup for a field age of 1 −
10 Gyr, which aretoo faint to be directly imaged with existing instruments. More discoveries of young, low-mass ultracooldwarfs are thereby needed to validate such blue Y − J colors at very low surface gravities as predicted bythe cloudless Sonora-Bobcat models.As shown in the lower panel of Figure 1, lower metallicity has qualitatively similar effects as higher sur-face gravity (except for the blue wing of Y band). A metal-poor atmosphere has relatively more abundantH and He atoms, and is less opaque, with the photosphere residing at lower-altitude, higher-pressure lay-ers. These lead to enhanced collision-induced H absorption (CIA) relative to other molecular absorption,thereby suppressing the flux in K band. While such K -band flux suppression also occurs for objects withhigh surface gravities, the amount of the suppression has a different dependence on log g and Z (e.g., Bur-rows et al. 2006; Leggett et al. 2007; Liu et al. 2007; Burningham et al. 2009). Quantitatively, our Figure 2shows that the log g should be increased/decreased by a factor of ≈ Z should bedecreased/increased, in order to achieve comparable K -band flux changes in the cloudless Sonora-Bobcatmodels. In addition, higher metallicity results in the suppressed flux at the blue wing of Y band and at ≈ . − . µ m, primarily because the enhancement of Na and K abundances leads to larger opacity inthe optical resonance lines of these species, whose pressure-broadened red wings extend to near-infrared(e.g., Tsuji et al. 1999; Burrows & Volobuyev 2003; Allard et al. 2007a). ATMOSPHERIC MODEL ANALYSISHere we use the cloudless Sonora-Bobcat models and the Starfish methodology to study atmosphericproperties of HD 3651B, GJ 570D, and Ross 458C. We first construct and validate our forward-modelingframework using Sonora-Bobcat models (Section 4.1). Then we present results of our atmospheric modelanalysis and compare them to those derived from the traditional forward-modeling approach and previousspectroscopic analyses (Section 4.2). Finally, we use the Starfish hyper-parameters to quantify the system-atic difference between data and models (Section 4.3).4.1.
Forward-Modeling Analysis with Starfish
Spectral Emulator Construction
We start our forward-modeling analysis by constructing Starfish’s spectral emulator, which can gener-ate the cloudless Sonora-Bobcat model spectra with an arbitrary set of grid parameters. (1) Followingthe Appendix of C15 (also see Habib et al. 2007 and Heitmann et al. 2009), we first trim the cloudlessSonora-Bobcat models to wavelengths of 0 . − . µ m over the parameter space relevant to late-T dwarfs:[600 , T eff , [3 . , .
5] dex in log g , and [ − . , + .
5] dex in Z , encompassing 286 grid points. Wethen pick 198 grid points where the Sonora-Bobcat models cover the same { T eff , log g } values among thethree metallicity points Z = { + . , , − . } . The selected 198 grid points are used to train the spectral em-ulator and are denoted as { θ (cid:63) } train , where θ (cid:63) = { T eff , log g , Z } . Model spectra from these training grid points, f λ (cid:0) { θ (cid:63) } train (cid:1) , have four dimensions (i.e., wavelength, T eff , log g , and Z ) with a spacing of 50 K and 100 K in T eff , 0 .
25 dex and 0 . g , and 0 . Z . The remaining 88 ( = − { θ (cid:63) } test , and their corresponding model spectra, f λ (cid:0) { θ (cid:63) } test (cid:1) , will be used to validate the performance ofour spectral emulator (Section 4.1.2).(2) We downgrade the resolution of the model spectra using a Gaussian kernel corresponding to the 0 . (cid:48)(cid:48) slit of SpeX. Our convolution process accounts for the wavelength-dependent spectral resolution of theSpeX prism mode (Rayner et al. 2003). We therefore end up with convolved models at both the training( { θ (cid:63) } train ) and the testing ( { θ (cid:63) } test ) grid points that have the same spectral resolution as the data.(3) We perform principal component analysis (PCA) of the model spectra from the training grid points f λ (cid:0) { θ (cid:63) } train (cid:1) . We standardize f λ (cid:0) { θ (cid:63) } train (cid:1) by subtracting their mean spectrum ξ µ and then dividing by theirstandard deviation spectrum ξ σ . Then we decompose the standardized model spectra into a few principaleigenspectra (e.g., Ivezi´c et al. 2014). We obtain 9 eigenspectra ξ k ( k = , , . . . , ) with a PCA thresholdvalue of 0 .
99 (Figure 3), meaning that the 99% summative variance of the standardized model spectra can bemodeled with these their corresponding eigenspectra. The model spectrum at each grid point can be therebyapproximated using the linear combination of ξ µ , ξ σ , and 9 eigenspectra ξ k (Figure 3; also see Equation 23of C15), f λ (cid:0) { θ (cid:63) } train (cid:1) = f rec λ (cid:0) { θ (cid:63) } train (cid:1) + (cid:15) emu (cid:0) { θ (cid:63) } train (cid:1) ≡ ξ µ + ξ σ (cid:88) k = w k (cid:0) { θ (cid:63) } train (cid:1) ξ k + (cid:15) emu (cid:0) { θ (cid:63) } train (cid:1) (1)Here f rec λ (cid:0) { θ (cid:63) } train (cid:1) is the reconstructed model spectra at the training grid points, and (cid:15) emu is the reconstruc-tion error, i.e., the flux difference between the original and reconstructed models. The standardized weightsof eigenspectra, w k (cid:0) { θ (cid:63) } train (cid:1) , is a function of grid parameters.(4) We train a Gaussian process on the standardized weights w k (cid:0) { θ (cid:63) } train (cid:1) for the eigenspectra at differentgrid parameters. For a given set of parameters θ (cid:63) = { T eff , log g , Z } , this Gaussian process provides aprobability distribution of w k ( θ (cid:63) ), thereby leading to a distribution of the interpolated spectra (Equation 1)as well as the propagation of the interpolation uncertainties into the resulting inferred parameters.Starfish uses a squared exponential kernel to describe the Gaussian process covariance matrix for thestandardized weights of each eigenspectrum. This kernel is characterized by four hyper-parameters, φ int , k = { a int , k , (cid:96) T eff , (cid:96) log g , (cid:96) Z } . The parameter a int , k describes the covariance along the diagonal axis, and (cid:96) T eff , (cid:96) log g , (cid:96) Z describe the length scales in grid parameters over which the covariance decreases in an exponential manner.We follow the suggestions of C15 and assume a Γ function prior for each length scale, so that each Γ function reaches its peak probability at a value ( ≈
100 K in T eff , ≈ . g , and ≈ . Z )which is one or two times the grid spacing of each parameter.When the PCA process decomposes grid models into 9 principal components, the flux information fromthe remaining less important components is truncated, leading to residuals (cid:15) emu between the original traininggrid models and the reconstructed ones (Equation 1). C15 accounts for this reconstruction error by including The T eff spacing is 50 K for [600 , , g spacing is 0 .
25 dex for [3 . , .
5] dex and0 . . , .
5] dex. an additional hyper-parameter λ ξ into the likelihood function of the Gaussian process. We assume a Γ function prior for λ ξ and compute its rate and shape following Equation 43 −
44 of C15.(5) We determine the posteriors of all 37 hyper-parameters (4 φ int , k parameters for each of 9 eigenspectraand one λ ξ ) using the Markov chain Monte Carlo (MCMC) algorithm emcee (Foreman-Mackey et al. 2013)with 148 walkers. We use the medians of the hyper-parameter posteriors to construct the spectral emulator.These 37 hyper-parameters constitute the core of the spectral emulator, allowing computation of the prob-ability distribution of the 9 eigenspectra’ standardized weights w k ( θ (cid:63) ) at a given set of grid parameters θ (cid:63) (Equation 47 − of C15). Combining ξ µ , ξ σ , ξ k , and the distribution of their standardized weights, Equa-tion 1 allows us to compute the distribution of interpolated model spectra (Equation 52 of C15). Starfish alsouses these hyper-parameters to define an “emulator covariance matrix” K E , which describes the covarianceof the interpolated spectral fluxes between different wavelengths (Equation 59 of C15). We include K E intothe final covariance matrix for the subsequent model evaluation, so that the interpolation uncertainties arepropagated into the inferred physical parameters (Section 4.1.5).4.1.2. Model Reconstruction
Before inferring physical parameters of objects, we verify that our spectral emulator can reconstruct boththe training and the testing Sonora-Bobcat models within an acceptable level. At each grid point, we gener-ate a probability distribution of spectra using the spectral emulator and then compare them with the originalmodel spectrum. Representative examples of the reconstructed training and testing models are shown inFigures 4 and 5, respectively.We note the spectra generated from our spectral emulator are nearly identical to the Sonora-Bobcat modelsover all grid points, confirming that our emulator is a robust model interpolator. At several grid points andin a certain wavelength ranges (primarily near the flux peaks in
Y JHK bands), there are noticeable fluxdifferences between the original and reconstructed model spectra. These differences mostly constitutes <
6% of the local flux and < .
5% of the peak J -band flux in the original models, much smaller than thesystematic difference between the cloudless Sonora-Bobcat models and late-T dwarf spectra ( ≈ − J -band flux; Section 4.3). In Section 4.1.6, we quantify the systematics in the derived physicalparameters due to imperfect model reconstruction.4.1.3. Physical Parameters and Priors
There are six physical parameters to be determined: effective temperature T eff , logarithmic surface grav-ity log g , metallicity Z , radial velocity v r , projected rotational velocity v sin i , and logarithmic solid anglelog Ω = log ( R / d ) , where i , R , and d is the inclination of the rotation axis, the radius, and the distance,respectively. For a given set of these parameters, we use our spectral emulator to generate a probability dis-tribution of Sonora-Bobcat model spectra at { T eff , log g , Z } , convolve the interpolated spectra with kernelsthat apply the Doppler shift by radial velocity v r and the broadening by rotation v sin i (e.g., Gray 2008), andthen scale the spectra using log Ω . The resulting model spectra are thereby ready to be compared with datato infer physical parameters.We assume uniform priors for T eff , log g , and Z within the parameter space where our spectral emulatoris constructed, i.e., [600 , T eff , [3 . , .
5] dex in log g , and [ − . , + .
5] dex in Z . We assume auniform prior of ( −∞ , + ∞ ) for both v r and log Ω . The component “ ( λ ξ Φ T Φ) ” in Equation 47 −
48 of C15 should be “ ( λ ξ Φ T Φ) − ”, though the Starfish code itself is free of thistypo. We assume a uniform prior of [0 , v max ] for v sin i and determine the v max value using an object’s oblate-ness. We follow Zhang et al. (2020) and express the rotationally induced oblateness as f = Cv / gd √ Ω (e.g., Barnes & Fortney 2003), where d is the distance, v rot is the equatorial rotational velocity, and C = . n = . n = . f crit = .
385 (James 1964). Altogether, we derive thefollowing constraints for v sin i : 0 (cid:54) v sin i (cid:54) v rot (cid:54) Ω / (cid:18) f crit gd C (cid:19) / ⇒ v sin i ∈ (cid:34) , v max ≡ Ω / (cid:18) f crit gd C (cid:19) / (cid:35) (2)For the three benchmarks, we follow Equation 2 and compute v max using their d , log g , and Ω in each stepof the spectral-fitting process. 4.1.4. Covariance Hyper-Parameters and Priors
In addition to physical parameters, Starfish uses three hyper-parameters { a N , a G , (cid:96) } to describe the co-variance matrix for subsequent model evaluation (Section 4.1.5). The hyper-parameter a N characterizes a“noise covariance matrix” K N , which describes the covariance in (data − model) residuals between wave-length pixels i and j , namely K Ni j = a N δ i j σ i σ j , where δ i j is the Kronecker delta and σ i , σ j is observed fluxuncertainty in wavelength pixels i , j , respectively. This diagonal K N assumes the residuals from differentwavelengths are independent and is usually adopted as the full covariance matrix by traditional forward-modeling analyses (Section 4.2.1). We assume a Gaussian prior on a N with mean of 1 .
0, standard deviationof 0 .
25, and a range of [0 , + ∞ ) .The other two hyper-parameters a G and (cid:96) characterize a “global covariance matrix” K G , which describesthe covariance in residuals between two wavelength pixels using a customized Matérn kernel. The hyper-parameter a G describes the covariance along the diagonal axis of the matrix and (cid:96) describes the wavelengthscale over which the covariance between two wavelengths decreases in nearly an exponential manner (seeAppendix A for details). This K G matrix has off-diagonal components and can account for the importantcorrelation in residuals caused by (1) the over-sampled instrumental line spread function, and (2) systematicsin model atmospheres. We assume a uniform prior of ( −∞ , + ∞ ) in log a G and assume a uniform prior on (cid:96) based on the prism-mode line spread function (Appendix A).4.1.5. Model Evaluation
There are nine free parameters { T eff , log g , Z , v r , v sin i , log Ω , a N , a G , (cid:96) } to be determined by ourforward-modeling framework. For a given spectrum, we first generate model spectra M using the spectralemulator with physical parameters incorporated. We then compare M with the data D to compute theresidual spectra, R = D − M ( T eff , log g , Z , v r , v sin i , log Ω) (3) C15 originally used the letter “ b ” for the hyper-parameter “ a N ” that we describe here. C S and likelihood function L S , C S = K E ( φ int , k , λ ξ ) + K N ( a N ) + K G ( a G , (cid:96) ) L S = − (cid:0) R T C − R + ln | C S | (cid:1) − n π ) (4)where n is number of wavelength pixels in the spectra. The C S matrix has three components: the emulatorcovariance matrix K E (Section 4.1.1), the noise covariance matrix K N , and the global covariance matrix K G (Section 4.1.4). The K E matrix is already determined from the spectral emulator construction and thereforeis fixed for a given slit width of the data. The K N and K G matrices (described by the three hyper-parameters)are determined as part of the analysis process.We use emcee (Foreman-Mackey et al. 2013) to fit our 1 . − . µ m spectra with 48 walkers. Duringthe fitting the process, we remove the first 5000 iterations as the burn-in phase and then iteratively esti-mate the chains’ average autocorrelation length using both the Goodman & Weare (2010) method and therevised version suggested by Daniel Foreman-Mackey . We terminate the fitting process once the numberof iteration exceeds 50 times the autocorrelation length of all parameters (e.g., Zhang et al. 2020). TheMCMC analysis for all objects converges with 1 × iterations, and we use the resulting chains to derivethe parameter posteriors for each object.4.1.6. Systematics of Inferred Physical Parameters
Here we validate our forward-modeling framework by investigating the systematics of the inferred physi-cal parameters caused by the Starfish machinery and data reduction.As described in Section 4.1.2, Starfish’s spectral emulator decomposes the Sonora-Bobcat grid models intoeigenspectra, leading to small flux differences between the original model spectra and those reconstructed bythe spectral emulator (e.g., Figures 4 and 5). While such reconstruction error is generally insignificant overthe entire grid, it may introduce systematic errors to the inferred physical parameters. These “emulator-induced systematics” are not fully accounted for when evaluating model parameters. Here we quantifythem by running our forward-modeling analysis directly on the original Sonora-Bobcat model spectra andcomparing the resulting parameter posteriors to the input model properties.We use all 286 Sonora-Bobcat model spectra with their spectral resolution downgraded to the 0 . (cid:48)(cid:48) slit.We treat them the same as the on-sky data and use the aforementioned spectral emulator, parameter pri-ors, covariance matrix, and likelihood function to infer the properties of each model spectrum. We adoptzero flux uncertainties for the original model spectra and thus do not include the noise covariance ma-trix K N . We then compare the resulting parameters to the input parameters of each model spectrum { T grideff , log g grid , Z grid , v r ≡ , v sin i ≡ , log Ω ≡ } to estimate the emulator-induced systematics.Figure 6 summarizes each physical parameter over all 286 grid points. We find the median and 1 σ errorsof the “derived-minus-input” parameter bias of + + − K in T eff , 0 . + . − . dex in log g , − . ± .
09 dex in Z , + + − km s − in v r , + + − km s − in v sin i , and − . + . − . dex in log Ω . The median parameter biasof each parameter (except v sin i ) is insignificant from zero, illustrating that our forward-modeling analysisis robust. The derived v sin i values are notably larger than zero at many grid points, and this occurs because v sin i is defined to be non-negative and the rotational broadening cannot be well-constrained by spectra with https://emcee.readthedocs.io/en/stable/tutorials/autocorr/ . When constructing the spectral emulator, Starfish uses the hyper-parameter λ ξ to account for the reconstruction error (Sec-tion 4.1.1). However, λ ξ only characterizes an average level of such error among all the training models across the entire0 . − . µ m wavelengths (also see Equation 33 of C15). Therefore, this single, constant λ ξ parameter is not enough to absorball reconstruction errors that occur at specific grid points and/or in short wavelength ranges, leading to occasionally noticeableflux differences between the original and reconstructed model spectra as seen in Figures 4 and 5. σ T eff , emu =
20 K σ log g , emu = . σ Z , emu = .
12 dex σ v r , emu =
10 km s − σ v sin i , emu =
40 km s − σ log Ω , emu = .
05 dex (5)These values are chosen conservatively to encompass the derived-minus-input values of parameters overmost of the grid points (Figure 6).In addition to these errors, the objects’ radial velocities v r has another source of systematic error due tothe uncertainty in the wavelength calibration of the SpeX data. This uncertainty is 5 . −
70 km s − over the 1 . − . µ m wavelength range. We thus adopt an additional error of σ v r , wcal =
180 km s − for theinferred radial velocity. The total systematic error of v r is ( σ v r , emu + σ v r , wcal ) / ≈
180 km s − .We incorporate these systematic errors for each object by modifying the MCMC chains obtained from thespectral-fitting process (Section 4.1.5). For each of the six physical parameters { T eff , log g , Z , v r , v sin i , log Ω } ,we draw errors from a normal distribution centered at zero with a standard deviation being the adopted sys-tematic error, and we make the number of draws the same as the number of MCMC samples. We thenadd these errors to the corresponding chain values to produce the final parameter posteriors. For the fourparameters { T eff , log g , Z , v r } , if any of their chain values are outside the range of [600 , T eff ,[3 . , .
5] for log g , [ − . , + .
5] for Z , and [0 , v max ] for v sin i ( v max is computed from Equation 2), we forcethose values to be at the lower or upper boundary. The resulting medians of all the final parameter posteriorsare nearly unchanged from those without the inclusion of these emulator-induced systematics.4.2. Results
Figure 7 presents the resulting parameter posteriors of our three benchmarks and also compares the ob-served data with the fitted Sonora-Bobcat model spectra, which are interpolated at the physical parametersdrawn from the MCMC samples (with systematic errors incorporated). The fitted Sonora-Bobcat modelsmatch the observations well, but large residuals are found in J and especially H bands, suggesting somephysical and chemical processes of late-T dwarf atmospheres are inadequately modeled or missing in thecloudless Sonora-Bobcat models. In J band ( ≈ . − . µ m), the fitted models over-predict spectra of thethree benchmarks. This discrepancy plausibly arises from clouds (e.g., Morley et al. 2012) and/or reduc-tions in the vertical temperature gradient (e.g., Tremblin et al. 2015), which can redistribute fluxes fromshorter to longer wavelengths and thereby produce spectra that could better match the observations. In H band ( ≈ . − . µ m), the fitted models under-predict the emergent flux. This discrepancy might resultfrom the spectral-fitting process itself aiming to balance the residuals from over-estimated fluxes in J band.Also, it might be related to the disequilibrium abundance of NH , which would be less than the amountcomputed by the equilibrium chemistry of Sonora-Bobcat models and thereby could weaken absorption in2the blue wing of H band to better match the data (e.g., Fegley & Lodders 1994; Saumon et al. 2006; Zahnle& Marley 2014). A detailed discussion of the cloudless Sonora-Bobcat models is in Section 6.We further use the fitted parameter posteriors and the objects’ parallaxes (from their primary stars) tocompute their radii ( R ) and masses ( M ). For each object, we draw parallaxes from a normal distributioncorresponding to the measured parallax and uncertainty, truncated within ( , + ∞ ) . We make the number ofdraws the same as their MCMC samples ( =
48 walkers × iterations) and combine these parallaxes andthe log Ω posterior to compute the R posterior. We then use log g and R to determine the M posterior.We summarize the inferred physical parameters and covariance hyper-parameters of the three benchmarksin Table 4 and 5, respectively. Our derived uncertainties in T eff , log g , and Z are close to 1 / − / R ≈ − v r of the three benchmark companions range from +
200 km s − to +
500 km s − , which are not consistentwith their primary stars’ Gaia
DR2 radial velocities and are even comparable with the local Galactic escapespeed ( ≈
540 km s − ; Smith et al. 2007; Piffl et al. 2014). However, the inferred v r are mostly insignificantfrom zero ( < σ ) and well within one resolution pixel (1200 − − ) of the prism spectra. Also,the inferred v sin i of these objects are only at 1 σ significance. Finally, we note that while both v r and v sin i are coupled with modeling systematics and cannot be well-constrained from the data, their posteriorshave at most weak correlations with (and thereby negligible impact on) the posteriors of the other physicalparameters. 4.2.1. Comparison with Traditional Forward-Modeling Approach
We have also conducted a forward-modeling analysis for the three benchmarks following the traditionalapproach (e.g., Cushing et al. 2008; Stephens et al. 2009; Rice et al. 2010; Zhang et al. 2020), where we uselinear interpolation to generate the model spectrum in between grid points and adopt a diagonal covariancematrix (defined by observed flux uncertainties) for model evaluation.We use the cloudless Sonora-Bobcat model atmospheres over their entire parameter space of [200 , T eff , [3 . , .
5] dex in log g , and [ − . , + .
5] dex in Z , amounting to 1014 grid points. The grid spacingis 25 K, 50 K, and 100 K in T eff , 0 .
25 dex and 0 . g , and 0 . Z . Then we determine the sixphysical parameters { T eff , log g , Z , v r , v sin i , log Ω } with same priors as our Starfish analysis (Section 4.1).At a given set of { T eff , log g , Z } , we use linear interpolation of the flux (conducted in logarithmic units for T eff ) to generate a single model spectrum. Compared to the spectral emulator of Starfish, linear interpolationprovides the exact model spectrum at each grid point but does not provide interpolation uncertainties forspectra in between grid points. We apply { v r , v sin i , log Ω } values to the interpolated spectra with the sameapproach as in Starfish and then compare data to models to compute the residual spectra R (Equation 3). We The T eff spacing is 25 K for [200 , , , g spacing varieswith metallicity. At Z = − .
5, the log g spacing is 0 . . , .
0] dex and 0 .
25 dex for [3 . , .
5] dex. At Z =
0, the log g spacing is 0 .
25 dex. At Z = + . g spacing is 0 . . , .
5] dex and 0 .
25 dex among the rest of points. C T and likelihood function L T , C T = K N ( a N ≡ ) L T = − (cid:0) R T C − R + ln | C T | (cid:1) − n π ) (6)Here we construct the covariance matrix C T with solely the diagonal noise covariance matrix and with a N fixed at 1. We do not include other covariance hyper-parameters as adopted by Starfish.We use emcee to fit the objects’ 1 . − . µ m spectra with 24 walkers, remove the first 5000 iterations asthe burn-in phase, and terminate the fitting process with 3 × iterations so these MCMC chains converge.We incorporate the systematic error of 180 km s − into the inferred radial velocity (caused by the uncertaintyin the wavelength calibration of the SpeX data), following the same approach as in our Starfish-basedanalysis (Section 4.1.6). We compute the objects’ radii and masses using their parallaxes and log g andlog Ω posteriors, and summarize all inferred parameters and their uncertainties in Table 4.As shown in Figure 7, the inferred parameters and fitted models derived from our Starfish-based andthe traditional methods are well consistent. However, the traditional approach produces artificially smallparameter errors, smaller than those of Starfish by factors of 2 −
15 in { T eff , log g , Z , R , M } . We note thelarger parameter uncertainties derived from Starfish are not due to the incorporation of the aforementionedemulator-induced systematic errors (which are not applied to the traditional approach; Section 4.1.6), butrather because Starfish propagates the model interpolation uncertainties and the correlated residuals amongadjacent wavelengths to the inferred parameters (which are ignored by the traditional approach). Therefore,our Starfish-based forward-modeling analysis produces more realistic error estimates.4.2.2. Comparison with Previous Spectroscopic Analyses
Table 4 compares our inferred physical parameters of the three benchmarks to those from previous spec-troscopic analyses, which used grid models or the retrieval method . For HD 3651B and GJ 570D, our fitted T eff are consistent with literature values, especially those using grid models with cloudless and chemical-equilibrium atmospheres. However, our log g and Z are smaller by ≈ . ≈ .
35 dex, respectively.Surface gravity and metallicity have similar effect on the spectral morphology, as either a high (low) log g or a low (high) Z leads to the same suppressed (enhanced) K -band flux in late-T dwarf spectra (Section 3.2).As a consequence, these two parameters are degenerate in the spectroscopic analysis (e.g., Burgasser et al.2006; Leggett et al. 2007; Liu et al. 2007). As shown in Table 4, most past studies of HD 3651B andGJ 570D used grid models with a fixed metallicity at solar abundance or that of their primary stars, andtheir fitted log g are therefore anchored at a specific Z . In comparison, we determine log g and Z simul-taneously from the spectral fitting process and therefore our results represent a comprehensive evaluationof the model parameter space. According to the subsequent evolutionary model analysis (Section 5), ourfitted log g and Z of these two objects are underestimated, indicating shortcomings of the adopted modelatmospheres (see Section 6). Moreover, our analysis implies that using the cloudless Sonora-Bobcat modelsto analyze late-T dwarf spectra whose metallicities are not known in priori can lead to inaccurate log g and Z .For Ross 458C, our fitted log g and Z match the literature values and our T eff are consistent with thoseusing cloudless, chemical-equilibrium model atmospheres. However, we note that analyses using cloudy The compilation of literature values for HD 3651B and GJ 570D is mostly from Line et al. (2017), to which we have addeda few new studies and more details of the grid models used in literature. ≈
100 K cooler T eff , which is consistent with our evolutionary model analysis (Section 5). While the specificcondensation opacities of some cloud models are perhaps uncertain (iron and silicate clouds versus sulfideclouds; see Morley et al. 2012), it is clear that cloudless model atmospheres cannot fully interpret the spectraof Ross 458C unless disequilibrium processes and reduced vertical temperature gradient are also considered(e.g., Tremblin et al. 2015).4.3. Assessment of Modeling Systematics using Starfish Hyper-Parameters
We can assess the modeling systematics using the derived covariance hyper-parameters of the three bench-marks listed in Table 5. As described in Section 4.1, the full covariance matrix of Starfish consists of threecomponents, the emulator covariance matrix K E , the noise covariance matrix K N , and the global covariancematrix K G . While K E has been already pre-determined during the spectral emulator construction, K N and K G are computed as part of our spectral-fitting process. The K N matrix has one hyper-parameter a N andaccounts for measurement uncertainties. The K G matrix is characterized by two hyper-parameters, (cid:96) and a G ,with the latter describing the amplitude of K G ’s normalization (Equation A1). The (cid:96) parameter describesthe auto-correlation wavelength scale of the residual and such correlated residuals are caused by (1) theover-sampled instrumental line spread function (LSF) and (2) the modeling systematics. If the latter effectdominates, then we can use these hyper-parameters to examine the performance of the models.As shown in Table 5, the fitted (cid:96) of all the three benchmarks are 2700 − − and exceed theexpected range of 820 − − from LSF (Appendix A). Therefore, the correlated residual of oursample are caused by the modeling systematics rather than the instrumental effect. As a consequence, theglobal covariance matrix K G mainly describes the model imperfections, and we can use its determinant toquantify such systematics.The hyper-parameter a G is a proxy for K G ’s normalization. Ideally, if the models match the observationswithin the measurement uncertainties, then the need for including K G into the spectral fitting is fairly lowand thus the fitted a G should be very small. A larger fitted a G indicates measurement uncertainties cannotsolely explain the difference between data and models, and thus Starfish needs K G with higher values toaccount for this data-model discrepancy. The meaning of a G can be understood from the definition of thenoise covariance K N (Section 4.1.4), which is a diagonal matrix with values of a N σ i , where σ i is the fluxuncertainty at the i -th wavelength pixel. Similarly, the global covariance K G is a band matrix, with a G onits main diagonal and with much smaller non-zero values along diagonals above or below (Equation A1).Following K N , we can express a G ≡ a N σ m and regard σ m as an equivalent flux that characterizes “modeluncertainties”. We can then normalize σ m of the object by its observed peak J -band flux and define (cid:15) J = (cid:112) a G / a N max (cid:0) f obs , J (cid:1) (7)The (cid:15) J parameter quantifies the modeling systematics, with higher values suggesting more significant data-model discrepancy. We compute (cid:15) J values for the three benchmarks (Table 5) and find that the systematicdifference between the cloudless Sonora-Bobcat models and spectra of the three benchmarks comprises2% −
4% of the observed peak J -band fluxes, equivalent to a S/N of 50 −
25. This implies that the modeluncertainties can exceed measurement uncertainties when fitting the cloudless Sonora-Bobcat models tolow-resolution spectra with high S/N ( >
50 in J band) and increasing the S/N of data does not promiseenhanced precision of fitted physical parameters.5 EVOLUTIONARY MODEL ANALYSISWe also derive physical properties of HD 3651B, GJ 570D, and Ross 458C using cloudless Sonora-Bobcatevolutionary models. We extend the Saumon et al. (2006) method from solar- to multi-metallicity evolution-ary models. Specifically, we compute { T eff , log g , Z } at which (1) the companions’ interpolated bolometricluminosities from the evolutionary models match the measured luminosities from their near-infrared spectra,their primary stars’ distances, and the atmospheric model-based ratio between near-infrared and bolometricfluxes, and (2) the companions’ metallicity and the interpolated age from evolutionary models match thoseof their primary stars (Section 2). We implement this method in a Bayesian fashion following Zhang et al.(2020).For each benchmark companion, we first integrate its SpeX spectrum in a wavelength range of 1 . − . µ m and use its primary star’s Gaia
DR2 distance to derive a near-infrared luminosity, namely L . − . µ m ,with uncertainties in spectra and distance propagated accordingly. We then execute a MCMC routine withfree parameters being T eff , log g , and Z . Given a set of these parameters, we use linear interpolation ofthe cloudless Sonora-Bobcat models to construct the corresponding synthetic spectrum, from which wecompute a ratio of the integrated flux in 1 . − . µ m to that in 0 . − µ m. We note such ratio is 0 − . T eff = −
600 K, 0 . − . T eff = − . − . − { T eff , log g , Z } grid spacing of < . < .
04, and < .
01 in each of the above T eff ranges, respectively. We therefore adopt a conservative uncertainty of 0 .
04. We apply this ratio tothe measured near-infrared luminosity to derive the companion’s atmospheric-based bolometric luminosity L bol , atm . The uncertainties in L . − . µ m and the model-based flux ratio are both propagated accordingly.For each set of { T eff , log g , Z } in the MCMC chains, we also interpolate the Sonora-Bobcat evolutionarymodels (in logarithmic scales for T eff and age) and compute the companion’s evolutionary-based bolometricluminosity L bol , evo and age t evo . Altogether, we derive { T eff , log g , Z } using the following likelihood function, L = p ( L bol , atm | L bol , evo ) × p ( Z (cid:63) | Z ) × p ( t (cid:63) | t evo ) × p ( T eff ) × p ( log g ) × p ( Z ) (8)where Z (cid:63) and t (cid:63) stand for the metallicity and age of primary stars summarized in Section 2. We assume thepriors, p ( T eff ) , p ( log g ) , and p ( Z ) are all uniform distributions within the parameter space of the Sonora-Bobcat models, i.e., [200 , T eff , [3 . , .
5] dex for log g , and [ − . , + .
5] dex for Z . To compute p ( L bol , atm | L bol , evo ) , we assume the bolometric luminosity follows a normal distribution corresponding to thevalue and uncertainty of L bol , atm . To compute p ( Z (cid:63) | Z ) and p ( t (cid:63) | t evo ) , we assume the metallicity and age ofthe companion follow the uniform distributions constrained by those of their primary stars.We use emcee to derive { T eff , log g , Z } values of each companion with 30 walkers, removing the first1000 iterations as the burn-in phase and terminating with 1 . × iterations so that our MCMC analysisconverges. We use the companions’ resulting { T eff , log g , Z } posteriors and the interpolated cloudlessSonora-Bobcat evolutionary models to derive their radii, masses, and bolometric luminosities. We presentthe posteriors of their evolutionary model parameters in Figure 8 and summarize their values in Table 6. BENCHMARKINGThe three benchmarks now have the same set of physical parameters { T eff , log g , Z , R , M } derived fromboth atmospheric models (Section 4) and evolutionary models (Section 5). Ideally, the results from thesetwo approaches would be consistent, but they are not. The evolutionary model parameters are expected tobe more robust (though not totally immune from the systematics; e.g., Dupuy et al. 2009a, 2014), given thewell-noted shortcomings of the atmospheric models (e.g., Cushing et al. 2008; Stephens et al. 2009; Leggett6et al. 2017). As a consequence, comparing physical parameters derived from these two model sets, a.k.a. the“benchmarking” process (e.g., Pinfield et al. 2006; Liu et al. 2008), can calibrate our inferred atmosphericmodel parameters and shed light on modeling systematics.To compare the results from atmospheric and evolutionary models, we subtract the values of each ob-ject’s evolutionary-based MCMC chains from those of its atmospheric-based chains . Distributions of theresulting parameter differences are shown in Figure 9, with values summarized in Table 6.HD 3651B has the highest evolutionary-based T eff (809 + − K) and log g (5 . ± .
06 dex) among thethree benchmarks. Its T eff and R derived from the atmospheric and evolutionary models are consis-tent within uncertainties, but its spectroscopically inferred log g and Z are lower by 1 . ± .
29 dex and0 . + . − . dex, respectively. The inaccurate atmospheric log g also leads to the object’s very small inferredmass (2 . + . − . M Jup ).GJ 570D has a slightly lower evolutionary-based T eff (786 ±
20 K) and log g (5 . ± .
13 dex) thanHD 3651B. The spectroscopically inferred log g and Z are underestimated by 1 . ± .
28 dex and0 . ± .
14 dex, respectively. In addition, our atmospheric model analysis slightly overestimates this ob-ject’s T eff by 42 ±
32 K and underestimates R by 0 . ± .
07 R
Jup (or a factor of 1 . + . − . ). The inaccurateatmospheric log g and R together lead to its small inferred mass (2 . + . − . M Jup ).Ross 458C has the lowest evolutionary-based T eff (682 + − K) and log g (4 . ± .
16 dex). This object hasconsistent log g and Z values from the two sets of models. However, its spectroscopically inferred T eff isoverestimated by 122 + − K and R is underestimated by 0 . ± .
07 R
Jup (or a factor of 1 . + . − . ), whichlead to its very small inferred mass (2 . + . − . M Jup ).To summarize, our atmospheric model analysis are in line with the evolutionary model analysis, but someparameters of the former appear to be inaccurate. We find the parameter difference between these twosets of models exhibits two outcomes. At T eff ≈ −
810 K and log g ≈ . − . T eff and R , but underestimate log g and Z by ≈ . − . ≈ . − . T eff ≈
700 K and a lower log g ≈ . g and Z , but overestimate T eff by ≈
120 Kand underestimate R by ≈ . Jup (or a factor of ≈ . { T eff , log g , Z , R } and its primary star’s distance. If the atmospheric models are perfect, then theseevolutionary-based model spectra should exactly match the observed spectrophotometry of each compan-ion. In other words, inconsistencies between these two sets of spectra can inform the specific modelingsystematics.Figure 10 compares the observed and evolutionary-based model spectra of three benchmarks at 1 . − . µ m, and Figure 11 extends such comparisons to 0 . − . µ m . For HD 3651B and GJ 570D, the Given that the number of MCMC samples of our atmospheric model analysis (48 walkers × iterations) is 16 times largerthan that of our evolutionary model analysis (30 walkers × iterations), we simply repeat the chains of the latter 16 timesand then subtract the two to generate the posteriors. We also generated posteriors of parameter differences using 1 /
16 of theatmospheric-based MCMC chains and obtained consistent results. We use the spectral emulator to generate 0 . − . µ m spectra and use the linear interpolation for 2 . − . µ m, over whichthe spectral emulator is not defined. In 0 . − . µ m, spectra of the three companions generated from the linear interpolation areconsistent with those generated from our spectral emulator within uncertainties. JHK bands, but produce notably fainter fluxes in Y band. This islikely related to the potassium resonance doublet at 0 . µ m (Allard et al. 1999; Burrows et al. 2000), whosepressure-broadened wings can extend to Y and J bands. The cloudless Sonora-Bobcat models adopt the Kline shape theory by Allard et al. (2007b) and these authors have recently produced improved calculations ofthe K − H potential (Allard et al. 2016). As demonstrated by Phillips et al. (2020), model atmospheres withthe Allard et al. (2016) prescription of K profiles can better match the observed Y -band spectra of GJ 570Dthan those of Allard et al. (2007b). In addition, the evolutionary-based model spectra of HD 3651B andGJ 570D have slightly brighter fluxes in W .
5] bands, suggesting the disequilibrium abundance ofCO in their atmospheres (e.g., Zahnle & Marley 2014).The data-model mismatch of HD 3651B and GJ 570D might explain why their { log g , Z } derived fromour forward-modeling analysis are underestimated. Since the evolutionary-based model spectra have under-predicted Y -band fluxes, the fitting process favors a much lower log g in order to boost the Y -band emission(at a given T eff and Z ) to match the data (e.g., Figure 1). The decreased log g as such further causes a lower Z , given the log g − Z degeneracy (Section 4.2.2), to ensure the data and models remain consistent in otherbands (e.g., K band).For Ross 458C, the evolutionary-based model spectra have significantly brighter fluxes in Y J bands andfainter fluxes in K band, suggesting the models should probably include clouds to properly interpret theobject’s spectra (e.g., Morley et al. 2012). Alternatively, Tremblin et al. (2015) proposed that the cloudlessdisequilibrium models with reduced vertical temperature gradient (as compared to the adiabatic thermalstructure) can predict redder near-infrared spectra than the chemical equilibrium models and better matchthe observations. We note the evolutionary-based models of Ross 458C have fattener fluxes in the bluewing of H band, likely due to the disequilibrium chemistry of NH which is not included in the cloudlessSonora-Bobcat models. Convective mixing can disturb the chemical equilibrium between N and NH inthe upper atmosphere, reducing the amount of NH3 and thereby weakening absorption in the blue wing ofH band (e.g., Fegley & Lodders 1994; Saumon et al. 2006, 2012; Cushing et al. 2011; Zahnle & Marley2014; Tremblin et al. 2015). Also, the brighter W .
5] fluxes as predicted by model spectra are likelydue to the disequilibrium mixing of CO.Again, we can use the data-model mismatch of Ross 458C to explain why its T eff and R derived fromour forward-modeling analysis are over- and under-estimated, respectively. The evolutionary-based modelspectra have significantly bluer near-infrared colors (e.g., Y − K , J − K , J − H ) than the data. Therefore,the spectral fitting is forced to choose models with redder colors to match the data. Among all three gridparameters, T eff is the primary modulator of near-infrared colors for late-T dwarfs. Consequently, a higher T eff than the evolutionary-based value is favored during the fitting process, which leads to a smaller R topreserve the integrated flux of spectra. SUMMARY AND FUTURE WORKWe have constructed a forward-modeling framework to analyze low-resolution ( R ≈ − . − . µ m) spectra of T dwarfs. We use the Bayesian inference tool Starfish (Czekalaet al. 2015) and state-of-art, cloudless Sonora-Bobcat model atmospheres within a parameter space of T eff = − g = . − . Z = − .
5, 0, + . χ -based) forward-modeling analyses. First, Starfish’s spectral emulator generates a probability distribution of interpolatedspectra at any model grid location, with the associated interpolation uncertainties then propagated into theresulting inferred parameters. Second, Starfish constructs a covariance matrix with off-diagonal compo-8nents to account for correlated residuals caused by the instrumental effect and modeling systematics. Bothof these two aspects are not accounted for in traditional analyses.In our forward-modeling framework, we verified that the spectral emulator can reproduce the originalSonora-Bobcat models, with only slight flux differences that are smaller than the systematics of Sonora-Bobcat models. To further validate our approach, we treated the original model spectra as if they wereon-sky data and used Starfish to infer the properties of each model spectrum. We found negligible offsetsbetween derived and input physical parameters and thereby confirmed the emulator-induced flux differencesfrom the original models do not bias our results.We then applied our forward-modeling framework to three benchmark late-T dwarfs, HD 3651B,GJ 570D, and Ross 458C, which are wide-orbit companions to main-sequence stars. We derived theseobjects’ effective temperatures ( T eff ), surface gravities (log g ), metallicities ( Z ), radii ( R ), and masses ( M ).Our fitted atmospheric models generally match the observed spectra of the three benchmarks. However,residuals are seen in J and H bands, likely arising from clouds and/or reductions of the vertical temperaturegradient in the atmospheres, as neither of these two effects are incorporated in the cloudless Sonora-Bobcatmodels. The data-model discrepancy in the blue wing of H band might also be related to the disequilib-rium abundance of NH , which would be less than the amount assumed by the equilibrium chemistry ofSonora-Bobcat.Our derived physical parameters are consistent with those derived from the traditional spectral-fittingapproach and from previous spectroscopic studies. The parameter uncertainties from the traditional methodare implausibly small, while the error estimates from our Starfish-based analysis are a factor of 2 −
15 largerin { T eff , log g , Z , R , M } and more realistic. Our resulting uncertainties in T eff , log g , and Z are typicallyabout 1 / − / − J -band fluxes, equivalent to a S/N of 50 −
25. Consequently, when usingthe cloudless Sonora-Bobcat models, increasing the S/N of observed spectra beyond 50 in J band will notlead to enhanced precision of fitted physical parameters, given that the model uncertainties will exceed themeasurement uncertainties in the forward-modeling analysis.We also used the Sonora-Bobcat evolutionary models to derive these benchmarks’ physical properties,based on their bolometric luminosities, and their primary stars’ metallicities and ages. As a result, thesebenchmarks have { T eff , log g , Z , R , M } derived from both the atmospheric models (via our forward-modeling analysis) and the evolutionary models. We assumed evolutionary-based parameters are morerobust and then used them to test the accuracy of spectral fits. We find the parameter difference between at-mospheric and evolutionary models exhibits two outcomes. At T eff ≈ −
810 K and log g ≈ . − . T eff and R , but underestimate log g and Z by ≈ . − . ≈ . − . T eff ≈
700 K and a lowerlog g ≈ . g and Z , but overestimate T eff by ≈
120 K and underestimate R by ≈ . Jup (or a factor of ≈ . { T eff , log g , Z , R } of the three benchmarks and thencompared with their observed spectrophotometry. For HD 3651B and GJ 570D, the evolutionary-basedmodel spectra have fainter fluxes in Y band and slightly brighter fluxes in W .
5] bands, suggesting9the modeling systematics mainly come from the uncertainties of potassium line profiles and the lack ofdisequilibrium abundance of CO. For Ross 458C, the evolutionary-based model spectra have brighter fluxesin Y , J , W .
5] bands and fainter fluxes in the blue wing of H band, K , W
1, and [3 .
6] bands, suggestingthe modeling systematics mainly come from the lack of clouds, reductions of vertical temperature gradient,and disequilibrium chemistry of CO/CH and N /NH .In a companion paper, we apply our forward-modeling analysis to a much larger sample of late-T dwarfspectra and investigate their population properties. Unlike benchmark companions, metallicities and agesof most field dwarfs are unknown, but we can examine their spectral-fitting residuals as a function of wave-length and atmospheric properties to test model atmospheres. Also, thanks to the flexibility of Starfish, ourforward-modeling framework can be conveniently extended to earlier/later spectral types, wider/narrowerwavelength coverages, and different spectral resolutions, leading to a stronger understanding of the preci-sion and accuracy of model assumption.Finally, we emphasize the value of benchmark systems, including wide-orbit companions (e.g., Pinfieldet al. 2006; Deacon et al. 2014; Zhang et al. 2020), members of nearby associations (e.g., Liu et al. 2016;Faherty et al. 2016; Gagné et al. 2018; Zhang et al. 2018), and substellar binaries with dynamical masses(e.g., Dupuy et al. 2009b; Dupuy & Liu 2017). As exemplified by HD 3651B, GJ 570D, and Ross 458C,each benchmark can have its physical properties determined independently from spectral fitting and thus canbe used to validate atmospheric models in a specific part of the parameter space. Continued discoveries andanalyses of these systems with diverse ages, compositions, and masses will therefore test and help improvemodel atmospheres over an unprecedentedly large extent of physical parameters.We thank Mark Phillips, Didier Saumon, Caroline Morley, Eugene Magnier, Paul Mollière, Joe Za-lesky, Trent Dupuy, Ehsan Gharib-Nezhad for insightful discussions and comments on this work. Wethank Ian Czekala, Michael Gully-Santiago, and Miles Lucas for helpful discussions about Starfish.This work benefited from the 2017–2019 Exoplanet Summer Program in the Other Worlds Laboratory(OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Founda-tion. M.C.L. acknowledges National Science Foundation (NSF) grant AST-1518339. The advanced com-puting resources from the University of Hawaii Information Technology Services – Cyberinfrastructureare greatly acknowledged, and Z. Z. thanks the technical support received from Curt Dodds. This workpresents results from the European Space Agency (ESA) space mission Gaia . Gaia data are being pro-cessed by the
Gaia
Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is pro-vided by national institutions, in particular the institutions participating in the
Gaia
MultiLateral Agree-ment (MLA). The
Gaia
Gaia
Facilities:
IRTF (SpeX)
Software:
Starfish (Czekala et al. 2015; Czekala et al. 2018), emcee (Foreman-Mackey et al. 2013),Spextool (Cushing et al. 2004), TOPCAT (Taylor 2005), Astropy (Astropy Collaboration et al. 2013, 2018),0IPtyhon (Pérez & Granger 2007), Numpy (Oliphant 2006), Scipy (Jones et al. 2001), Matplotlib (Hunter2007).1APPENDIX A. NOTES ON THE GLOBAL COVARIANCE HYPER-PARAMETERSIn this Appendix, we describe the meaning of the global covariance hyper-parameters in Starfish and theirimplications for the systematics in model assumption. As described in Section 4.1.4, Starfish uses the globalcovariance matrix K G to characterize the correlation in residuals among adjacent pixels, as caused by (1)the over-sampled instrumental line spread function and (2) the systematics in model atmospheres. Thiscovariance matrix has two hyper-parameters { a G , (cid:96) } , with each matrix element expressed as follows (alsosee Equation 9 −
12 of C15) : K Gi j = w i j a G (cid:32) + √ r i j (cid:96) (cid:33) exp (cid:32) − √ r i j (cid:96) (cid:33) where r i j = c λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i − λ j (cid:12)(cid:12)(cid:12)(cid:12) w i j = (cid:2) + cos (cid:0) π r i j / r (cid:1)(cid:3) / r i j (cid:54) r r i j > r r = (cid:96) (A1)Here c is the speed of light and i , j ∈ { , , , . . . , n } , where n is number of wavelength pixels in the spectra.Consequently, K G is a band matrix, with a G on its main diagonal (where i = j ) and with much smallernon-zero values located along diagonals above or below followed by a truncation at r i j = r . The hyper-parameter a G is therefore a proxy of K G ’s normalization or determinant. As described in Section 4.3, thefitted a G values from the forward-modeling analysis inform uncertainties of model atmospheres.The other hyper-parameter (cid:96) describes the wavelength scale over which K G decreases exponentially alonganti-diagonal directions. (cid:96) has units of km s − (same as the r i j in Equation A1) and describes such auto-correlation wavelength in an equivalent velocity space. In order to better understand the meaning of (cid:96) values as inferred from our forward-modeling analysis, here we derive a conversion between (cid:96) and theauto-correlation wavelength.We first convert the covariance matrix K G into a correlation matrix R G , namely R Gi j ≡ K
Gi j / a G , which onlydepends on (cid:96) . Then we compute the autocorrelation as a function of wavelength offset ∆ λ , R G (∆ λ ; (cid:96) ) = n (cid:88) i = R Gi j ( (cid:96) ) such that λ i − λ j = ∆ λ , for j ∈ { , . . . , n } (A2)Using the typical wavelength grid from the spectra in our late-T dwarf sample, we construct K G usingEquation A1 and compute R G (∆ λ ; (cid:96) ) using Equation A2 for a set of different (cid:96) values and then derive theautocorrelation wavelength τ for each (cid:96) from the half width at half maximum (HWHM) of R G (∆ λ ; (cid:96) ) . Wefinally fit a line to the computed τ and their corresponding (cid:96) values and obtain τ = . × (cid:96) − + . The expression of K Gi j in Equation A1 is adopted by Starfish v0.2 (used in this work) and is slightly different from the originaldefinition of C15 (their Equation 9 −
12) where r i j = c (cid:12)(cid:12)(cid:12)(cid:12) λ i − λ j λ i + λ j (cid:12)(cid:12)(cid:12)(cid:12) r = (cid:96) These differences should have negligible impact for the inference of the objects’ physical properties. (cid:96) of each late-T dwarf from our forward-modelinganalysis into an auto-correction wavelength.Equation A3 can also help determine the (cid:96) values that correspond to the prism mode line spread function.Based on the wavelength-dependent spectral resolution of SpeX’s prism mode (Rayner et al. 2003), wederive that the autocorrelation wavelength of the 0 . (cid:48)(cid:48) slit from 1 . − . µ m spans 50 −
102 Å (with amedian of 83 Å), which is equivalent to 2 − (cid:96) = − − for the 0 . (cid:48)(cid:48) slit.In our forward-modeling analysis, we set the prior of (cid:96) based on the prism-mode line spread function. Weassume a uniform prior of [820 , ×
5] in (cid:96) for spectra taken in the 0 . (cid:48)(cid:48) slit and here we multiply 5 to themaximum (cid:96) as expected from the line spread function, in order to capture the correlation in residuals causedby modeling systematics. Therefore, comparing the resulting (cid:96) values of our late-T dwarfs to the rangeexpected from the instrumental line spread function, we can examine whether the systematics in modelatmospheres is playing a crucial role in interpreting the residuals (see Section 4.3).3REFERENCES Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556,872, doi:
Aguilera-Gómez, C., Ramírez, I., & Chanamé, J. 2018,A&A, 614, A55,doi:
Allard, F., Allard, N. F., Homeier, D., et al. 2007a,A&A, 474, L21,doi:
Allard, F., & Freytag, B. 2010, Highlights ofAstronomy, 15, 756,doi:
Allard, F., Hauschildt, P. H., Alexander, D. R.,Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357,doi:
Allard, F., Homeier, D., & Freytag, B. 2011,Astronomical Society of the Pacific ConferenceSeries, Vol. 448, Model Atmospheres From VeryLow Mass Stars to Brown Dwarfs, ed.C. Johns-Krull, M. K. Browning, & A. A. West, 91Allard, N. F., Royer, A., Kielkopf, J. F., & Feautrier, N.1999, PhRvA, 60, 1021,doi:
Allard, N. F., Spiegelman, F., & Kielkopf, J. F. 2007b,A&A, 465, 1085,doi: —. 2016, A&A, 589, A21,doi:
Allende Prieto, C., Barklem, P. S., Lambert, D. L., &Cunha, K. 2004, A&A, 420, 183,doi:
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi:
Astropy Collaboration, Price-Whelan, A. M., Sip˝ocz,B. M., et al. 2018, AJ, 156, 123,doi:
Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M.,Mantelet, G., & Andrae, R. 2018, AJ, 156, 58,doi:
Baliunas, S., Sokoloff, D., & Soon, W. 1996, ApJL,457, L99, doi:
Barnes, J. W., & Fortney, J. J. 2003, ApJ, 588, 545,doi:
Barnes, S. A. 2003, ApJ, 586, 464,doi: —. 2007, ApJ, 669, 1167, doi:
Bertelli, G., Girardi, L., Marigo, P., & Nasi, E. 2008,A&A, 484, 815,doi:
Bertelli, G., Nasi, E., Girardi, L., & Marigo, P. 2009,A&A, 508, 355,doi:
Best, W. M. J., Liu, M. C., Magnier, E. A., et al. 2015,ApJ, 814, 118,doi:
Beuzit, J. L., Ségransan, D., Forveille, T., et al. 2004,A&A, 425, 997,doi:
Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016,A&A, 585, A5,doi:
Bonnefoy, M., Chauvin, G., Lagrange, A. M., et al.2014, A&A, 562, A127,doi:
Bonnefoy, M., Perraut, K., Lagrange, A. M., et al.2018, A&A, 618, A63,doi:
Bowler, B. P., Liu, M. C., Mawet, D., et al. 2017, AJ,153, 18, doi:
Bressan, A., Marigo, P., Girardi, L., et al. 2012,MNRAS, 427, 127,doi:
Burgasser, A. J. 2007, ApJ, 658, 617,doi:
Burgasser, A. J. 2014, in Astronomical Society of IndiaConference Series, Vol. 11, Astronomical Society ofIndia Conference Series, 7–16Burgasser, A. J., Burrows, A., & Kirkpatrick, J. D.2006, ApJ, 639, 1095, doi:
Burgasser, A. J., McElwain, M. W., Kirkpatrick, J. D.,et al. 2004, AJ, 127, 2856, doi:
Burgasser, A. J., Tinney, C. G., Cushing, M. C., et al.2008, ApJL, 689, L53, doi:
Burgasser, A. J., Kirkpatrick, J. D., Cutri, R. M., et al.2000, ApJL, 531, L57, doi:
Burgasser, A. J., Simcoe, R. A., Bochanski, J. J., et al.2010, ApJ, 725, 1405,doi:
Burningham, B., Marley, M. S., Line, M. R., et al.2017, MNRAS, 470, 1177,doi:
Burningham, B., Pinfield, D. J., Leggett, S. K., et al.2009, MNRAS, 395, 1237,doi: Burningham, B., Pinfield, D. J., Lucas, P. W., et al.2010, MNRAS, 406, 1885,doi:
Burningham, B., Leggett, S. K., Homeier, D., et al.2011, MNRAS, 414, 3590,doi:
Burningham, B., Cardoso, C. V., Smith, L., et al. 2013,MNRAS, 433, 457,doi:
Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert,J. 2001, Reviews of Modern Physics, 73, 719,doi:
Burrows, A., & Liebert, J. 1993, Reviews of ModernPhysics, 65, 301,doi:
Burrows, A., Marley, M. S., & Sharp, C. M. 2000,ApJ, 531, 438, doi:
Burrows, A., Sudarsky, D., & Hubeny, I. 2006, ApJ,640, 1063, doi:
Burrows, A., Sudarsky, D., & Lunine, J. I. 2003, ApJ,596, 587, doi:
Burrows, A., & Volobuyev, M. 2003, ApJ, 583, 985,doi:
Burrows, A., Marley, M., Hubbard, W. B., et al. 1997,ApJ, 491, 856, doi:
Casagrande, L., Flynn, C., Portinari, L., Girardi, L., &Jimenez, R. 2007, MNRAS, 382, 1516,doi:
Casagrande, L., Schönrich, R., Asplund, M., et al.2011, A&A, 530, A138,doi:
Chandrasekhar, S. 1939, An introduction to the studyof stellar structureChauvin, G., Desidera, S., Lagrange, A. M., et al.2017, A&A, 605, L9,doi:
Chiu, K., Fan, X., Leggett, S. K., et al. 2006, AJ, 131,2722, doi:
Cottaar, M., Covey, K. R., Meyer, M. R., et al. 2014,ApJ, 794, 125,doi:
Cross, N. J. G., Collins, R. S., Mann, R. G., et al. 2012,A&A, 548, A119,doi:
Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al.2014, AJ, 147, 113,doi:
Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004,PASP, 116, 362, doi:
Cushing, M. C., Marley, M. S., Saumon, D., et al.2008, ApJ, 678, 1372, doi:
Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al.2011, ApJ, 743, 50,doi:
Cutri, R. M., & et al. 2014, VizieR Online DataCatalog, 2328Czekala, I., Andrews, S. M., Mandel, K. S., Hogg,D. W., & Green, G. M. 2015, ApJ, 812, 128,doi:
Czekala, I., gully, Gullikson, K., et al. 2018,iancze/Starfish: ca. Czekala et al. 2015 release w/Zenodo, doi: . https://doi.org/10.5281/zenodo.2221006 Deacon, N. R., Liu, M. C., Magnier, E. A., et al. 2014,ApJ, 792, 119,doi:
Del Burgo, C., Martín, E. L., Zapatero Osorio, M. R.,& Hauschildt, P. H. 2009, A&A, 501, 1059,doi:
Delgado Mena, E., Israelian, G., González Hernández,J. I., et al. 2010, ApJ, 725, 2349,doi:
Demarque, P., Guenther, D. B., Li, L. H., Mazumdar,A., & Straka, C. W. 2008, Ap&SS, 316, 31,doi:
Demarque, P., Woo, J.-H., Kim, Y.-C., & Yi, S. K.2004, ApJS, 155, 667, doi:
Deming, D., Wilkins, A., McCullough, P., et al. 2013,ApJ, 774, 95,doi:
Donahue, R. A. 1993, PhD thesis, New Mexico StateUniversity, University Park.Dupuy, T. J., & Kraus, A. L. 2013, Science, 341, 1492,doi:
Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19,doi: —. 2017, ApJS, 231, 15,doi:
Dupuy, T. J., Liu, M. C., Bowler, B. P., et al. 2010,ApJ, 721, 1725,doi:
Dupuy, T. J., Liu, M. C., & Ireland, M. J. 2009a, ApJ,692, 729,doi: —. 2009b, ApJ, 699, 168,doi: —. 2014, ApJ, 790, 133,doi: Duquennoy, A., & Mayor, M. 1988, A&A, 200, 135Edge, A., Sutherland, W., & Viking Team. 2016,VizieR Online Data Catalog, II/343Faherty, J. K., Burgasser, A. J., Walter, F. M., et al.2012, ApJ, 752, 56,doi:
Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016,ApJS, 225, 10,doi:
Fegley, Bruce, J., & Lodders, K. 1994, Icarus, 110,117, doi:
Feltzing, S., & Gustafsson, B. 1998, A&AS, 129, 237,doi:
Fernandes, J. M., Vaz, A. I. F., & Vicente, L. N. 2011,A&A, 532, A20,doi:
Fischer, D. A., Butler, R. P., Marcy, G. W., Vogt, S. S.,& Henry, G. W. 2003, ApJ, 590, 1081,doi:
Foreman-Mackey, D., Hogg, D. W., Lang, D., &Goodman, J. 2013, PASP, 125, 306,doi:
Fortney, J. J., Lodders, K., Marley, M. S., & Freedman,R. S. 2008a, ApJ, 678, 1419,doi:
Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D.,& Freedman, R. 2005, ApJL, 627, L69,doi:
Fortney, J. J., Marley, M. S., Saumon, D., & Lodders,K. 2008b, ApJ, 683, 1104, doi:
Fortney, J. J., Mordasini, C., Nettelmann, N., et al.2013, ApJ, 775, 80,doi:
Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, ApJ,856, 23, doi:
Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al.2016, A&A, 595, A1,doi:
Gaia Collaboration, Brown, A. G. A., Vallenari, A.,et al. 2018, A&A, 616, A1,doi:
Gaidos, E., & Mann, A. W. 2014, ApJ, 791, 54,doi:
Gaidos, E. J. 1998, PASP, 110, 1259,doi:
Geballe, T. R., Saumon, D., Golimowski, D. A., et al.2009, ApJ, 695, 844,doi:
Geballe, T. R., Saumon, D., Leggett, S. K., et al. 2001,ApJ, 556, 373, doi:
Gelino, C. R., Kirkpatrick, J. D., Cushing, M. C., et al.2011, AJ, 142, 57,doi:
Ghezzi, L., Cunha, K., Smith, V. V., et al. 2010, ApJ,720, 1290,doi:
Goldman, B., Marsat, S., Henning, T., Clemens, C., &Greiner, J. 2010, MNRAS, 405, 1140,doi:
Goodman, J., & Weare, J. 2010, Communications inApplied Mathematics and Computational Science, 5,65, doi:
Gordon, S., & McBride, B. J. 1994, NASA ReferencePubl. 1311Gray, D. F. 2008, The Observation and Analysis ofStellar PhotospheresGray, R. O., Corbally, C. J., Garrison, R. F.,McFadden, M. T., & Robinson, P. E. 2003, AJ, 126,2048, doi:
Habib, S., Heitmann, K., Higdon, D., Nakhleh, C., &Williams, B. 2007, PhRvD, 76, 083503,doi:
Heitmann, K., Higdon, D., White, M., et al. 2009, ApJ,705, 156,doi:
Hempelmann, A., Schmitt, J. H. M. M., Schultz, M.,Ruediger, G., & Stepien, K. 1995, A&A, 294, 515Henry, T. J., Soderblom, D. R., Donahue, R. A., &Baliunas, S. L. 1996, AJ, 111, 439,doi:
Hewett, P. C., Warren, S. J., Leggett, S. K., & Hodgkin,S. T. 2006, MNRAS, 367, 454,doi:
Hunter, J. D. 2007, Computing in Science &Engineering, 9, 90,doi:
Ivezi´c, Ž., Connelly, A. J., Vand erPlas, J. T., & Gray,A. 2014, Statistics, Data Mining, and MachineLearning in AstronomyJames, R. A. 1964, ApJ, 140, 552,doi:
Jarrett, T. H., Cohen, M., Masci, F., et al. 2011, ApJ,735, 112,doi:
Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy:Open source scientific tools for Python.
Kiraga, M., & Stepien, K. 2007, AcA, 57, 149. https://arxiv.org/abs/0707.2577 Kirkpatrick, J. D., Cushing, M. C., Gelino, C. R., et al.2011, ApJS, 197, 19,doi:
Kirkpatrick, J. D., Martin, E. C., Smart, R. L., et al.2019, ApJS, 240, 19,doi:
Knapp, G. R., Leggett, S. K., Fan, X., et al. 2004, AJ,127, 3553, doi:
Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014,Nature, 505, 69, doi:
Lawrence, A., Warren, S. J., Almaini, O., et al. 2007,MNRAS, 379, 1599,doi: —. 2012, VizieR Online Data Catalog, II/314Leggett, S. K., Marley, M. S., Freedman, R., et al.2007, ApJ, 667, 537, doi:
Leggett, S. K., Tremblin, P., Esplin, T. L., Luhman,K. L., & Morley, C. V. 2017, ApJ, 842, 118,doi:
Leggett, S. K., Golimowski, D. A., Fan, X., et al. 2002,ApJ, 564, 452, doi:
Leggett, S. K., Cushing, M. C., Saumon, D., et al.2009, ApJ, 695, 1517,doi:
Leggett, S. K., Burningham, B., Saumon, D., et al.2010, ApJ, 710, 1627,doi:
Leggett, S. K., Saumon, D., Marley, M. S., et al. 2012,ApJ, 748, 74,doi:
Lenzuni, P., Chernoff, D. F., & Salpeter, E. E. 1991,ApJS, 76, 759, doi:
Ligi, R., Creevey, O., Mourard, D., et al. 2016, A&A,586, A94,doi:
Lindegren, L. 2018, Gaia Technical Note:GAIA-C3-TN-LU-LL-124-01Line, M. R., Teske, J., Burningham, B., Fortney, J. J.,& Marley, M. S. 2015, ApJ, 807, 183,doi:
Line, M. R., Stevenson, K. B., Bean, J., et al. 2016, AJ,152, 203,doi:
Line, M. R., Marley, M. S., Liu, M. C., et al. 2017,ApJ, 848, 83,doi:
Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ,833, 96, doi:
Liu, M. C., Dupuy, T. J., & Ireland, M. J. 2008, ApJ,689, 436, doi:
Liu, M. C., Leggett, S. K., & Chiu, K. 2007, ApJ, 660,1507, doi:
Liu, M. C., Delorme, P., Dupuy, T. J., et al. 2011, ApJ,740, 108,doi:
Lodders, K. 1999, ApJ, 519, 793,doi:
Lucas, P. W., Tinney, C. G., Burningham, B., et al.2010, MNRAS, 408, L56,doi:
Luhman, K. L., Patten, B. M., Marengo, M., et al.2007, ApJ, 654, 570, doi:
MacDonald, R. J., Marley, M. S., Fortney, J. J., &Lewis, N. K. 2018, ApJ, 858, 69,doi:
Mace, G. N., Kirkpatrick, J. D., Cushing, M. C., et al.2013, ApJS, 205, 6,doi:
Mainzer, A., Cushing, M. C., Skrutskie, M., et al.2011, ApJ, 726, 30,doi:
Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687,1264, doi:
Manjavacas, E., Goldman, B., Alcalá, J. M., et al.2016, MNRAS, 455, 1341,doi:
Marley, M. S., Gelino, C., Stephens, D., Lunine, J. I.,& Freedman, R. 1999, ApJ, 513, 879,doi:
Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268,doi:
Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53,279, doi:
Marley, M. S., Saumon, D., Cushing, M., et al. 2012,ApJ, 754, 135,doi:
Marley, M. S., Saumon, D., Fortney, J. J., et al. 2017,in American Astronomical Society MeetingAbstracts, Vol. 230, American Astronomical SocietyMeeting Abstracts
Marley, M. S., Seager, S., Saumon, D., et al. 2002,ApJ, 568, 335, doi:
Marley, M. S., Harman, C., Hammel, H. B., et al.2020, arXiv e-prints, arXiv:2007.10549. https://arxiv.org/abs/2007.10549 Marois, C., Macintosh, B., Barman, T., et al. 2008,Science, 322, 1348,doi:
McKay, C. P., Pollack, J. B., & Courtin, R. 1989,Icarus, 80, 23,doi:
Mishenina, T. V., Soubiran, C., Kovtyukh, V. V.,Katsova, M. M., & Livshits, M. A. 2012, A&A, 547,A106, doi:
Mollière, P., van Boekel, R., Dullemond, C., Henning,T., & Mordasini, C. 2015, ApJ, 813, 47,doi:
Montes, D., López-Santiago, J., Gálvez, M. C., et al.2001, MNRAS, 328, 45,doi:
Morel, P. 1997, A&AS, 124, 597,doi:
Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012,ApJ, 756, 172,doi:
Moses, J. I., Madhusudhan, N., Visscher, C., &Freedman, R. S. 2013, ApJ, 763, 25,doi:
Mugrauer, M., Seifahrt, A., Neuhäuser, R., & Mazeh,T. 2006, MNRAS, 373, L31,doi:
Oliphant, T. 2006, NumPy: A guide to NumPy, USA:Trelgol Publishing.
Parmentier, V., Fortney, J. J., Showman, A. P., Morley,C., & Marley, M. S. 2016, ApJ, 828, 22,doi:
Pérez, F., & Granger, B. E. 2007, Computing inScience and Engineering, 9, 21,doi:
Petigura, E. A., & Marcy, G. W. 2011, ApJ, 735, 41,doi:
Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020,A&A, 637, A38,doi:
Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F.2004, ApJ, 612, 168, doi: —. 2006, ApJ, 642, 797, doi:
Pietrinferni, A., Cassisi, S., Salaris, M., Percival, S., &Ferguson, J. W. 2009, ApJ, 697, 275,doi:
Piffl, T., Scannapieco, C., Binney, J., et al. 2014, A&A,562, A91,doi:
Pinfield, D. J., Jones, H. R. A., Lucas, P. W., et al.2006, MNRAS, 368, 1281,doi:
Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., &Ventura, P. 2003, A&A, 397, 147,doi:
Prugniel, P., Vauglin, I., & Koleva, M. 2011, A&A,531, A165,doi:
Ramírez, I., Allende Prieto, C., & Lambert, D. L. 2013,ApJ, 764, 78,doi:
Rayner, J. T., Toomey, D. W., Onaka, P. M., et al.2003, PASP, 115, 362, doi:
Rice, E. L., Barman, T., Mclean, I. S., Prato, L., &Kirkpatrick, J. D. 2010, ApJS, 186, 63,doi:
Rocha-Pinto, H. J., Flynn, C., Scalo, J., et al. 2004,A&A, 423, 517,doi:
Salasnich, B., Girardi, L., Weiss, A., & Chiosi, C.2000, A&A, 361, 1023. https://arxiv.org/abs/astro-ph/0007388
Samland, M., Mollière, P., Bonnefoy, M., et al. 2017a,A&A, 603, A57,doi: —. 2017b, A&A, 603, A57,doi:
Santos, N. C., Israelian, G., & Mayor, M. 2004, A&A,415, 1153,doi:
Santos, N. C., Israelian, G., Mayor, M., et al. 2005,A&A, 437, 1127,doi:
Saumon, D., Geballe, T. R., Leggett, S. K., et al. 2000,ApJ, 541, 374, doi:
Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327,doi:
Saumon, D., Marley, M. S., Abel, M., Frommhold, L.,& Freedman, R. S. 2012, ApJ, 750, 74,doi:
Saumon, D., Marley, M. S., Cushing, M. C., et al.2006, ApJ, 647, 552, doi:
Saumon, D., Marley, M. S., Lodders, K., & Freedman,R. S. 2003, in IAU Symposium, Vol. 211, BrownDwarfs, ed. E. Martín, 345Scholz, R.-D. 2010, A&A, 515, A92,doi: Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007,MNRAS, 379, 755,doi:
Soderblom, D. R., Duncan, D. K., & Johnson, D. R. H.1991, ApJ, 375, 722, doi:
Sorahana, S., & Yamamura, I. 2012, ApJ, 760, 151,doi:
Spiegel, D. S., Burrows, A., & Milsom, J. A. 2011,ApJ, 727, 57,doi:
Stephens, D. C., Leggett, S. K., Cushing, M. C., et al.2009, ApJ, 702, 154,doi:
Takeda, G., Ford, E. B., Sills, A., et al. 2007, ApJS,168, 297, doi:
Taylor, M. B. 2005, in Astronomical Society of thePacific Conference Series, Vol. 347, AstronomicalData Analysis Software and Systems XIV, ed.P. Shopbell, M. Britton, & R. Ebert, 29Testi, L. 2009, A&A, 503, 639,doi:
Thompson, M. A., Kirkpatrick, J. D., Mace, G. N.,et al. 2013, PASP, 125, 809,doi:
Thorén, P., & Feltzing, S. 2000, A&A, 363, 692. https://arxiv.org/abs/astro-ph/0010067
Tinney, C. G., Burgasser, A. J., & Kirkpatrick, J. D.2003, AJ, 126, 975, doi:
Tinney, C. G., Kirkpatrick, J. D., Faherty, J. K., et al.2018, ApJS, 236, 28,doi:
Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015,ApJL, 804, L17,doi:
Tsuji, T. 2002, ApJ, 575, 264,doi: —. 2005, ApJ, 621, 1033, doi:
Tsuji, T., Ohnaka, K., & Aoki, W. 1999, ApJ, 520,L119, doi:
UKIDSS Consortium. 2012, VizieR Online DataCatalog, II/316Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141,doi:
Vican, L. 2012, AJ, 143, 135,doi:
Visscher, C., Lodders, K., & Fegley, Bruce, J. 2010,ApJ, 716, 1060,doi:
Wright, J. T., Marcy, G. W., Butler, R. P., & Vogt, S. S.2004, ApJS, 152, 261, doi:
Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41,doi:
Zalesky, J. A., Line, M. R., Schneider, A. C., &Patience, J. 2019, ApJ, 877, 24,doi:
Zhang, Z., Liu, M. C., Best, W. M. J., et al. 2018, ApJ,858, 41, doi:
Zhang, Z., Liu, M. C., Hermes, J. J., et al. 2020, ApJ,891, 171, doi: T eff = 600K T eff = 900K, log g = 4.5, Z = 0 T eff = 1200K . log g = 3.5 T eff = 900K, log g = 4.5, Z = 0log g = 5.5 Wavelength ( m) Z = 0.5 T eff = 900K, log g = 4.5, Z = 0 Z = + 0.5 N o r m a li z e d F l u x ( r e l a t i v e t o J ) Figure 1.
Comparisons of cloudless Sonora-Bobcat model spectra with varying effective temperature T eff , surfacegravity log g , and metallicity Z . The upper panel compares spectra with T = g = . Z =
0. The middle panel compares spectra with log g = . . . T eff =
900 K and Z =
0. The lower panel compares spectra with Z = + . − . T eff =
900 K and log g = .
5. All model spectra have been smoothedto match the wavelength-dependent resolution of the 0 . (cid:48)(cid:48) slit of SpeX prism ( R ≈ − J -band fluxes. Wavelength ( m) F l u x V a r i a t i o n p e r P a r a m e t e r S t e p T eff (per 10 K)log g (per 0.1 dex) Z (per 0.1 dex) Y J H K
Figure 2.
The average flux increase of model spectra relative to a spectrum at { T eff =
900 K , log g = . , Z = } ,when T eff is increased by 10 K (magenta), log g is increased by 0 . Z is increased by 0 . J -band fluxes (as those plotted in Figure 1). As aguide to the location of prominent spectral features, the upper and lower boundary of the background grey shadowcorrespond to a spectrum with { T eff =
900 K , log g = . , Z = } , scaled with a positive and an equally negativefactor. Figure 3.
Left: The mean spectrum ξ µ , standard deviation spectrum ξ σ , and 9 eigenspectra ξ k of the 198 cloudlessSonora-Bobcat models f λ (cid:0) { θ (cid:63) } train (cid:1) used to train the spectral emulator (Section 4.1.1). These spectra are scaledwith the same scale factor and then offset by constants. Right: Comparisons between the original Sonora-Bobcatmodel spectra (offset by a constant; black) and the linear combination of ξ µ , ξ σ , and ξ k (orange) at { T eff , log g , Z } = {
600 K , . , − . } (top), {
900 K , . , . } (middle), and { , . , + . } (bottom). We scale the flux ofeach Sonora-Bobcat model spectrum and use the same scale factor for the corresponding reconstructed spectra. Thereconstructed spectra by PCA components well-approximate the original models. T eff /log g / Z = 600/3.5/ 0.5residualSonora training modelemulator reconstruction 600/4.5/0.0residual 600/5.5/ + 0.5residual0.4+0.40246810 . S c a l e d F l u x Wavelength ( m)
Figure 4.
Each panel compares a Sonora-Bobcat model spectrum (black) used to define the spectral emulator (a.k.a.,the training models; Section 4.1.1) with a spectra distribution generated from the emulator (blue) at a given grid point(with { T eff , log g , Z } values listed at the upper corner). We normalize the flux of each Sonora-Bobcat model spectrumto the same J -band peak and then use the same scale factor for the reconstructed ones from the spectral emulator.Residuals between the original and reconstructed models (black) are shown at the bottom of each panel (with thesame vertical scale as the spectra shown at the top), overlaid with the 1 σ , 2 σ , and 3 σ dispersions (blue backgroundshadows) of the spectra distribution from the emulator. Our emulator can generate spectra nearly identical to theSonora-Bobcat models. Small residuals can be found near the Y JHK flux peaks, which mostly constitute < .
5% ofthe peak J -band flux of the original models, smaller than the modeling systematics of the cloudless Sonora-Bobcatmodels (Section 4.3). T eff /log g / Z = 600/5.25/ 0.5residualSonora testing modelemulator reconstruction 600/3.75/0.0residual 600/4.75/ + 0.5residual0.4+0.40246810 . S c a l e d F l u x Wavelength ( m)
Figure 5.
Comparisons between the Sonora-Bobcat model spectra (black) not used to define the spectral emulator(a.k.a. the testing models) and the spectra distribution generated from the spectral emulator (blue), with the sameformat as Figure 4. The testing models can be well reconstructed even though they are not used to construct thespectral emulator.
600 700 800 900 1000 1100 1200 T eff l o g g Effective Temperature ( T eff ) log g T eff Z d e r i v e d i n p u t ( K )
600 700 800 900 1000 1100 1200 T eff l o g g Surface Gravity (log g ) log g T eff Z d e r i v e d i n p u t ( d e x )
600 700 800 900 1000 1100 1200 T eff l o g g Metallicity ( Z ) log g T eff Z d e r i v e d i n p u t ( d e x )
600 700 800 900 1000 1100 1200 T eff l o g g Radial Velocity ( v r ) log g T eff Z d e r i v e d i n p u t ( k m s )
600 700 800 900 1000 1100 1200 T eff l o g g Projected Rotational Velocity ( v sin i ) log g T eff Z d e r i v e d i n p u t ( k m s )
600 700 800 900 1000 1100 1200 T eff l o g g Solid Angle (log ) log g T eff Z d e r i v e d i n p u t ( d e x ) Figure 6.
Test of the systematic errors (arising from the spectral emulator) derived by fitting the origi-nal Sonora-Bobcat models using the Starfish analysis. Each panel shows one of the six physical parameters { T eff , log g , Z , v r , v sin i , log Ω } with all the 286 Sonora-Bobcat grid points, where the small squares located atthe { T eff , log g } grid points correspond to solar metallicity ( Z = Z = − .
5; lower left) and super-solar metallicity ( Z = + .
5; upper right). The layout of gridpoints is also indicated by the compass shown at the upper right in each panel. We use large grey squares to mark thetraining grid points used to construct the spectral emulator and do not mark the remaining testing grid points (Sec-tion 4.1.1). We use three different colors to map the derived-minus-input bias (caused by the spectral emulator) of thecorresponding parameter, with boundary values defined by our adopted systematic errors (Equation 5). We note thatthe parameter bias over most grid points are well within our adopted systematics (light blue), and such bias has nosignificant difference between the training and testing grid points. Figure 7.
Forward-modeling results of HD 3651B, GJ 570D, and Ross 458C, using the Starfish (blue) and thetraditional approach (purple). Left: The upper panel for each object shows the observed spectrum (black) and themedian Sonora-Bobcat model spectra of those interpolated at parameters drawn from the MCMC chains based on theStarfish (blue) and traditional (purple) methods. The object’s name, spectral type, the slit width and J -band S/N ofits spectrum, and inferred physical parameters are in the upper right corner. The middle and lower panel for eachobject shows the residual of each method (data − model; black), with the blue and purple shadows being 1 σ and 2 σ dispersions of 5 × draws from the Starfish and traditional covariance matrix, respectively. Right: Posteriors of thesix physical parameters { T eff , log g , Z , v r , v sin i , log Ω } derived from the Starfish-based forward-modeling analysis(blue). We overlay the median values and uncertainties from the traditional method (purple), shown as vertical linesand shadows in the 1-D histograms and as circles and error bars in the 2-D histograms. We use grey vertical andhorizontal lines to mark the { T eff , log g , Z } grids points of the cloudless Sonora-Bobcat models. Figure 8.
Posterior distributions of evolutionary model parameters for HD 3651B (violet), GJ 570D (light blue), andRoss 458C (orange). The median and 1 σ uncertainties of each parameter are labeled over each 1-D histogram, withcolors indicating the corresponding objects. The parameters plotted here are effective temperature ( T eff ), logarithmicsurface gravity (log g ), metallicity ( Z ), radius ( R ), and mass ( M ). The companions’ log g and Z posteriors are close touniform distributions because we assume uniform distributions for their primary stars’ age and metallicity. Figure 9.
Posterior distributions of the atmospheric − evolutionary model parameter differences for HD 3651B (violet),GJ 570D (light blue), and Ross 458C (orange) as described in Section 6. We generate these posteriors by subtractingthe evolutionary-based MCMC chains from the atmospheric-based chains. The median and 1 σ uncertainties of theparameter differences are labeled over each 1-D histogram, with colors indicating the corresponding objects. Thepositions of these posteriors suggest the accuracy of model assumption has two outcomes. At T eff ≈ −
810 Kand log g ≈ . − . T eff and R butunderestimated log g (by ≈ . − . Z (by ≈ . − . T eff ≈
700 K and log g ≈ . g and Z but overestimated T eff (by ≈
120 K) and underestimated R (by ≈ . Jup or ≈ . × ). . HD 3651B IRTF/SpeX (0.5 )Fitted Atmospheric Model SpectraEvolutionary-based Model Spectra . GJ 570D
Wavelength ( m) . Ross 458C F l u x f ( × e r g / s / c m / Å ) Figure 10.
The atmospheric model spectra (red) and their 1 σ uncertainties (orange shadow) generated and scaled atthe evolutionary-based { T eff , log g , Z , R } of three benchmark companions, HD 3651B (top), GJ 570D (middle), andRoss 458C (bottom). We also show the observed spectra (black) and fitted atmospheric model spectra (blue) as inFigure 7. . IRTF/SpeX (0.5 )Fitted Atmospheric Model SpectraEvolutionary-based Model SpectraHD 3651B Y MKO J MKO H MKO K MKO W W . GJ 570D Y MKO J MKO H MKO K MKO W W Wavelength ( m) . Ross 458C Y MKO J MKO H MKO K MKO W W F l u x f ( × e r g / s / c m / Å ) Figure 11.
The observed spectra (black), the fitted atmospheric model spectra (blue), and the evolutionary-basedmodel spectra (red) of three benchmark companions, using the same format as Figure 10. Thinner black lines are usedto plot low-S/N regions of the observed spectra. Black squares show the observed fluxes derived from photometry,and red and blue circles show the synthetic fluxes from evolutionary-based and fitted atmospheric model spectra,respectively. We compute these fluxes using the filter responses and zero-point fluxes from Hewett et al. (2006) andLawrence et al. (2007) for MKO photometry, Jarrett et al. (2011) for
WISE , and the IRAC Instrument Handbook for
Spitzer /IRAC. Table 1.
Metallicity and Age of HD 3651A
Value a Reference Notes b Metallicity (dex) 0 . ± .
06 7 Spectroscopy ( R ≈ . × )0 . ± .
04 10 Spectroscopy ( R ≈ . × )0 .
12 19 Spectroscopy ( R ≈ . × )0 .
16 21 Spectroscopy ( R ≈ . × )0 . ± .
05 23 Spectroscopy ( R ≈ . × )0 . ± .
07 24 Spectroscopy ( R ≈ . × )0 . ± .
04 27 Spectroscopy ( R ≈ . × ) Adopted Metallicity (dex) 0 . − . Isochrone Age (Gyr) 3 − . (Yonsei-Yale) isochrones (8)5 .
13 14 Padova isochrones (5) > . .
059 20 CESAM isochrones (4)6 . ± . . ± .
52 26 PARSEC (Padova and Trieste Stellar Evolutionary Code) isochrones (22)10 . + . − .
27 Y (Yonsei-Yale) isochrones (8)age range 3 − . ± .
99 13 Rotation period 44 days (3)7 . ± . . − . . − . R (cid:48) HK = − .
991 (3)8 . ± . R (cid:48) HK = − .
09 (6,15)3 . ± . R (cid:48) HK = − .
85 (9)7 . ± . R (cid:48) HK = − .
02 (11)2 . ± . R X ≡ log ( L X / L bol ) = − .
69 (2)age range 2 . − . Adopted Age (Gyr) 4 . − . a The gyrochronology and stellar-activity ages from this work are computed based on the rotation-age, R (cid:48) HK -age, and R X -age relationsfrom Mamajek & Hillenbrand (2008). b Notes contain (1) the method and spectral resolution used to compute the metallicity, and (2) the models, rotation periods, stellar-activityindices, and related references (indicated by the number in parenthesese) used to compute the object’s ages.
References —(1) This work, (2) Hempelmann et al. (1995), (3) Baliunas et al. (1996), (4) Morel (1997), (5) Salasnich et al. (2000), (6)Gray et al. (2003), (7) Allende Prieto et al. (2004), (8) Demarque et al. (2004), (9) Rocha-Pinto et al. (2004), (10) Santos et al. (2004),(11) Wright et al. (2004), (12) Valenti & Fischer (2005), (13) Barnes (2007), (14) Casagrande et al. (2007), (15) Liu et al. (2007), (16)Takeda et al. (2007), (17) Demarque et al. (2008), (18) Mamajek & Hillenbrand (2008), (19) Delgado Mena et al. (2010), (20) Fernandeset al. (2011), (21) Petigura & Marcy (2011), (22) Bressan et al. (2012), (23) Mishenina et al. (2012), (24) Ramírez et al. (2013), (25)Bonfanti et al. (2016), (26) Ligi et al. (2016), (27) Aguilera-Gómez et al. (2018) Table 2.
Metallicity and Age of GJ 570A
Value a Reference Notes b Metallicity (dex) 0 . ± .
09 4 Spectroscopy ( R ≈ . × )0 .
04 5 Spectroscopy ( R ≈ . − . × )0 . ± .
10 9 Spectroscopy ( R ≈ . × )0 . ± .
03 19 Spectroscopy ( R ≈ . × ) − . ± .
05 21 Spectroscopy ( R ≈ . × ) − . ± .
17 23 Spectroscopy ( R ≈ . − . × )0 . ± .
04 24 Spectroscopy ( R ≈ . × ) Adopted Metallicity (dex) − . − . Isochrone Age (Gyr) 3 . + . − .
10 Y (Yonsei-Yale) isochrones (7) < . + −
19 Y (Yonsei-Yale) isochrones (7)6 . + . − .
20 BASTI isochrones (8,11,18)6 . + . − .
20 Padova isochrones (15,17)7 . + . − .
24 Y (Yonsei-Yale) isochrones (7)age range 1 − . ± .
56 12 Rotation period 44 . . ± .
24 1 Rotation period 44 . . − . . ± . R (cid:48) HK = − .
75 (2)0 . ± . R (cid:48) HK = − .
49 (3)1 . ± . R (cid:48) HK = − .
63 (22)1 . ± . R X ≡ log ( L X / L bol ) = − .
38 (6)1 . ± . R X = − .
27 (13)age range 0 . − . Adopted Age (Gyr) 1 . − . a The gyrochronology and stellar-activity ages from this work are computed based on the rotation-age, R (cid:48) HK -age, and R X -age relations from Mamajek & Hillenbrand (2008). b Notes contain (1) the method and spectral resolution used to compute the metallicity, and (2) the models, rotationperiods, stellar-activity indices, and related references (indicated by the number in parenthesese) used to computethe object’s ages.
References —(1) This work, (2) Soderblom et al. (1991), (3) Henry et al. (1996), (4) Feltzing & Gustafsson (1998), (5)Thorén & Feltzing (2000), (6) Pizzolato et al. (2003), (7) Demarque et al. (2004), (8) Pietrinferni et al. (2004), (9)Santos et al. (2005), (10) Valenti & Fischer (2005), (11) Pietrinferni et al. (2006), (12) Barnes (2007), (13) Liu et al.(2007), (14) Takeda et al. (2007), (15) Bertelli et al. (2008), (16) Demarque et al. (2008), (17) Bertelli et al. (2009),(18) Pietrinferni et al. (2009), (19) Ghezzi et al. (2010), (20) Casagrande et al. (2011), (21) Prugniel et al. (2011),(22) Vican (2012), (23) Line et al. (2015), (24) Aguilera-Gómez et al. (2018) Table 3.
Astrometry and Photometry of HD 3651B, GJ 570D, and Ross 458C
HD 3651B GJ 570D Ross 458CSpectral Type T7.5 T7.5 T8R.A. (hh:mm:ss.ss; J2000) 00:39:18.91 14:57:15.84 13:00:41.64Decl. (hh:mm:ss.ss; J2000) +21:15:16.8 -21:22:08.1 +12:21:14.6 µ α cos δ (mas yr − ) − . ± .
11 1031 . ± . − . ± . µ δ (mas yr − ) − . ± . − . ± . − . ± . . ± .
06 170 . ± .
09 86 . ± . Y MKO (mag) 17 . ± .
06 15 . ± .
10 17 . ± . J MKO (mag) 16 . ± .
03 14 . ± .
05 16 . ± . H MKO (mag) 16 . ± .
04 15 . ± .
05 17 . ± . K MKO (mag) 16 . ± .
05 15 . ± .
05 16 . ± . W . ± .
04 16 . ± . W . ± .
02 13 . ± . .
6] (mag) 15 . ± .
04 13 . ± .
02 15 . ± . .
5] (mag) 13 . ± .
02 12 . ± .
02 13 . ± . Spitzer /IRAC References 10 31 28
References —(1) Leggett et al. (2002), (2) Tinney et al. (2003), (3) Knapp et al. (2004), (4) Chiu et al. (2006), (5) Luhman et al.(2007), (6) Burgasser et al. (2008), (7) Leggett et al. (2009), (8) Burningham et al. (2010), (9) Lucas et al. (2010), (10) Leggett et al.(2010), (11) Gelino et al. (2011), (12) Kirkpatrick et al. (2011), (13) Cross et al. (2012), (14) Dupuy & Liu (2012), (15) Fahertyet al. (2012), (16) Lawrence et al. (2012), (17) Leggett et al. (2012), (18) UKIDSS Consortium (2012), (19) Burningham et al.(2013), (20) Dupuy & Kraus (2013), (21) Mace et al. (2013), (22) Thompson et al. (2013), (23) Cushing et al. (2014), (24) Cutri &et al. (2014), (25) Best et al. (2015), (26) Edge et al. (2016), (27) Gaia Collaboration et al. (2016), (28) Leggett et al. (2017), (29)Gaia Collaboration et al. (2018), (30) Tinney et al. (2018), (31) Kirkpatrick et al. (2019), (32) Best et al. (submitted) T a b l e . S p ec t r o s c op i ca ll y I n f e rr e d P hy s i ca l P r op e r ti e s o f HD B , G J D , a nd R o ss C F itt e d P a r a m e t e r s a D e r i v e d P a r a m e t e r s b M od e l s c O b j ec t λλ T e ff l og g Z d v r v s i n i l og Ω R M G r i d s C l oud l e ss ? C h e m . E q . ? R e f e r e n ce s ( µ m )( K )( d e x )( d e x )( k m s − )( k m s − )( d e x )( R J up )( M J up ) HD B . − . + ( + ) − ( − ) . + . ( + . ) − . ( − . ) − . + . ( + . ) − . ( − . ) + ( + ) − ( + ) + ( + ) − ( + ) − . + . ( + . ) − . ( + . ) . + . ( + . ) − . ( + . ) . + . ( + . ) − . ( + . ) , Y e s Y e s T h i s W o r k ( S t a r fi s h ) . − . + − . + . − . − . + . − . + ( + ) − ( + ) + − − . + . − . . + . − . . + . − . , Y e s Y e s T h i s W o r k ( T r a d iti on a l ) . − . ± . ± . . ± . (cid:63) ·················· e
10 0 . − . − . − . ≈ . ··············· , , Y e s Y e s
11 1 . − . ≈ . − . (cid:63) ········· . − . ··· , Y e s Y e s
20 1 . − . . (cid:63) ········· . ··· Y e s
20 0 . − . + − . ± . . ± . ········· . ± . ··· , Y e s Y e s
27 1 . − . + − . + . − . . + . − . ········· . + . − . ······ f G J D . − . + ( + ) − ( − ) . + . ( + . ) − . ( − . ) − . + . ( + . ) − . ( − . ) + ( + ) − ( + ) + ( + ) − ( + ) − . + . ( + . ) − . ( + . ) . + . ( + . ) − . ( + . ) . + . ( + . ) − . ( + . ) , Y e s Y e s T h i s W o r k ( S t a r fi s h ) . − . + − . + . − . − . + . − . + ( + ) − ( + ) + − − . + . − . . + . − . . + . − . , Y e s Y e s T h i s W o r k ( T r a d iti on a l ) . − . − . (cid:63) ·················· e . − . − . − . (cid:63) ··············· Y e s ,
14 1 . − . ± . ± . (cid:63) ··· ± ········· Y e s Y e s
13 0 . − . . (cid:63) ··············· Y e s Y e s
15 1 . − . ≈ . − . (cid:63) ········· . − . ··· , Y e s Y e s
20 1 . − . . (cid:63) ········· . ··· Y e s
20 1 . − . . (cid:63) ··············· , Y e s
22 0 . − . + − . ± . . ± . ········· . ± . ··· , Y e s Y e s
27 1 . − . + − . + . − . − . + . − . ········· . + . − . ······ f R o ss C . − . + ( + ) − ( − ) . + . ( + . ) − . ( − . ) . + . ( + . ) − . ( − . ) + ( + ) − ( + ) + ( + ) − ( + ) − . + . ( + . ) − . ( + . ) . + . ( + . ) − . ( + . ) . + . ( + . ) − . ( + . ) , Y e s Y e s T h i s W o r k ( S t a r fi s h ) . − . + − . + . − . . + . − . + ( + ) − ( + ) + − − . + . − . . + . − . . + . − . , Y e s Y e s T h i s W o r k ( T r a d iti on a l ) . − . + − . + . − . ≈ . ··············· Y e s Y e s
17 0 . − . + − ≈ . (cid:63) ··············· Y e s
17 1 . − . − . (cid:63) ··············· Y e s Y e s
19 1 . − . − . . (cid:63) ··············· Y e s Y e s T ab l e c on ti nu e donn ex t pag e T a b l e ( c on ti nu e d ) F itt e d P a r a m e t e r s a D e r i v e d P a r a m e t e r s b M od e l s c O b j ec t λλ T e ff l og g Z d v r v s i n i l og Ω R M G r i d s C l oud l e ss ? C h e m . E q . ? R e f e r e n ce s ( µ m )( K )( d e x )( d e x )( k m s − )( k m s − )( d e x )( R J up )( M J up ) . − . . (cid:63) ··············· Y e s
19 1 . − . − . − . (cid:63) ··············· Y e s
19 1 . − . − . . (cid:63) ··············· Y e s
19 1 . − . − . − . (cid:63) ········· . − . ··· , Y e s Y e s
20 1 . − . . (cid:63) ········· . ··· Y e s
20 0 . − . . (cid:63) ··············· Y e s
21 0 . − . . . (cid:63) ··············· Y e s a P a r a m e t e r s d e r i v e d fr o m t h i s w o r k a r e s ho w n a s m e d i a n a nd1 σ un ce r t a i n ti e s , w it h s y s t e m a ti ce rr o r s ( S ec ti on4 . . ) a l r ea dy i n c o r po r a t e d . V a l u e s i n s i d e p a r e n t h e s e s a r e f o r m a l s p ec t r a l - fi tti ngun ce r t a i n ti e s fr o m S t a r fi s hb e f o r ea pp l y i ng a ny s y s t ea m ti ce rr o r s . b R a d ii ( R ) i s d e r i v e d fr o m t h e ob j ec t s ’ p a r a ll a x e s a nd fi tt e d l og Ω a nd m a ss ( M ) i s d e r v e d fr o m R a nd t h e fi tt e d l og g . A g e ( t ) i s d e r i v e dby i n t e r po l a ti ng t h e e vo l u ti on a r y m od e l s u s i ng t h e fi tt e d { T e ff , l og g , Z } v a l u e s . c G r i d m od e l s a dop t e dbyp r e v i ou s f o r w a r d - m od e li ng a n a l y s e s a nd w h e t h e r t h e s e m od e l s a r e g e n e r a t e d w it h c l oud l e ss a nd / o r c h e m i ca l e qu ili b r i u m a ss u m p ti on s ( on l y s ho w n i f " Y e s " ) . d T h e " (cid:63) " s y m bo l s m a r k m e t a lli c iti e s t h a t a r ea ss u m e d i n s t ea do f fi tt e ddu r i ng t h e s p ec t r a l fi tti ngp r o ce ss . e P a r a m e t e r s a r e d e r i v e d fr o m t h ee m p i r i ca l a n a l y s i s o f s p ec t r a li nd i ce s by B u r g a ss e r e t a l . ( ) . f P a r a m e t e r s a r e d e r i v e d fr o m r e t r i e v a l a n a l y s i s . R e f ere n ce s — ( ) A ll a r d e t a l . ( ) , ( ) A c k e r m a n & M a r l e y ( ) , ( ) M a r l e y e t a l . ( ) , ( ) T s u ji ( ) , ( ) B u rr o w s e t a l . ( ) , ( ) S a u m on e t a l . ( ) , ( ) T s u ji ( ) , ( ) B u r g a ss e r e t a l . ( ) , ( ) S a u m on e t a l . ( ) , ( ) B u r g a ss e r( ) , ( ) L e gg e tt e t a l . ( ) , ( ) S a u m on & M a r l e y ( ) , ( ) D e l B u r go e t a l . ( ) , ( ) G e b a ll ee t a l . ( ) , ( ) T e s ti ( ) , ( ) A ll a r d & F r e y t a g ( ) , ( ) B u r g a ss e r e t a l . ( ) , ( ) A ll a r d e t a l . ( ) , ( ) B u r n i ngh a m e t a l . ( ) , ( ) L i u e t a l . ( ) , ( ) M o r l e y e t a l . ( ) , ( ) S o r a h a n a & Y a m a m u r a ( ) , ( ) M o lli è r ee t a l . ( ) , ( ) T r e m b li n e t a l . ( ) , ( ) L i n ee t a l . ( ) , ( ) M a r l e y e t a l . ( ) , ( ) S a m l a nd e t a l . ( ) , ( ) M a r l e y e t a l . ( i np r e p ) Table 5.
Spectroscopically Inferred Covariance Hyper-Parameters and (cid:15) J of HD 3651B,GJ 570D, and Ross 458C Object a N (cid:96) log a G (cid:15) J (km s − ) (dex)HD 3651B 1 . + . − . + − − . + . − . . + . − . GJ 570D 1 . + . − . + − − . + . − . . + . − . Ross 458C 0 . + . − . + − − . + . − . . + . − . T a b l e . A t m o s ph e r i ca nd E vo l u ti on a r y M od e l P a r a m e t e r s o f B e n c h m a r k C o m p a n i on s HD B G J D R o ss C P a r a m e t e r A t m o s ph e r i c E vo l u ti on a r y A t m . − E vo . A t m o s ph e r i c E vo l u ti on a r y A t m . − E vo . A t m o s ph e r i c E vo l u ti on a r y A t m . − E vo . P r i m a r y S t a r A g e ( t (cid:63) ; G y r) . − . . − . . − . P r i m a r y S t a r M e t a lli c it y ( Z (cid:63) ; d e x ) . − . − . − . . − . E ff ec ti v e T e m p e r a t u r e ( T e ff ; K ) + − + − + − + − + − + − + − + − + − S u rf ace G r a v it y ( l og g ; c g s ) . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . M e t a lli c it y ( Z ; d e x ) − . + . − . . + . − . − . + . − . − . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . R a d i u s ( R ; R J up ) . + . − . . + . − . . + . − . . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . M a ss ( M ; M J up ) . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . B o l o m e t r i c L u m i no s it y ( l og ( L bo l / L (cid:12) ) ; d e x ) − . + . − . − . + . − . . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − ..