Orbit and Dynamical Mass of the Late-T Dwarf GL 758 B^*

Gl 758 was included in the target sample of the long-duration RV planet search at McDonald Observatory (e.g., Cochran et al. 1997; Endl et al. 2016). The Tull Coudé spectrograph (Tull et al. 1995) was used at the Harlan J. Smith 2.7 m telescope in combination with an I₂ cell in the light path to obtain precise differential RVs. We commenced observations of Gl 758 on 1998 November 4, and have accumulated 118 precise RV measurements over the past 19 years. Beginning in 2009 an exposure meter was used to provide the optimal exposure level and compute precise barycentric corrections. We measure precise RVs from these spectra using our Austral I₂ cell data code (Endl et al. 2000). The exposure times of the GJ 758 spectra range from 365 to 1200 s, primarily controlled by atmospheric seeing conditions to reach the desired signal-to-noise ratio (S/N) of ∼300 per pixel in the I₂ bandpass (500 to 650 nm). The RV data have a total rms of 13 m s⁻¹ and a median uncertainty of 4.6 m s⁻¹. Our measurements are listed in Table 1 and displayed in Figure 1.

Zoom In Zoom Out Reset image size

Gl 758 has been continuously monitored with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) at the Keck I telescope since 2006. This star was initially targeted as an RV standard for FGK stars in the Kepler field due to its proximity on the sky to the Kepler footprint and its stable RVs, but after several years of coverage it became apparent that Gl 758 was undergoing a shallow radial acceleration (G. Marcy 2016, private communication). Altogether 262 spectra were gathered following the standard setup, observing strategy, and procedure for measuring relative RVs implemented by the California Planet Search program (Howard et al. 2010). An iodine cell is mounted in the optical path before the slit entrance to provide a set of stable reference lines (Marcy & Butler 1992), an exposure meter is used to consistently achieve a S/N of about 225 per reduced pixel near 550 nm, and relative RVs are extracted by forward modeling the stellar and iodine spectra convolved with the instrument line spread function (Valenti et al. 1995). The median uncertainty from these measurements is 1.2 m s⁻¹. A secular trend of −2.82 ± 0.03 m s⁻¹ yr⁻¹ is apparent in our Keck RVs (Figure 1). This slope is slightly steeper than the trend from McDonald Observatory, suggesting a recent change in the acceleration. This is readily apparent by considering only the latest HIRES data, which exhibits a slope of −3.08 m s⁻¹ yr⁻¹ since 2010, and −4.15 m s⁻¹ yr⁻¹ since 2013.5. Our HIRES RVs and uncertainties are listed in Table 1.

Observations of Gl 758 have been carried out autonomously at the 2.4 m APF telescope at Lick Observatory since 2013. A total of 250 echelle spectra were gathered as part of APF's automated Doppler search for rocky planets (Fulton et al. 2015) with the Levy Spectrograph, which employs an iodine cell to measure precise relative RVs (Vogt et al. 2014). Each spectrum spans 3740–9700 Å at a resolving power of ≈100,000 with the 1'' decker. RVs are measured by forward-modeling 848 spectral regions, and the resulting variance is adopted as the RV uncertainty. The median instrumental precision of these measurements for Gl 758 is 1.4 m s⁻¹. The APF RVs reveal an acceleration that is significantly steeper than the McDonald and Keck RVs (Figure 1). The linear trend from APF is −4.24 ± 0.07 m s⁻¹ yr⁻¹, indicating a substantial evolution in recent years. Fortunately, these evolving slopes provide curvature that can better constrain the orbit and dynamical mass of the companion. Our APF RVs and uncertainties are listed in Table 1.

We observed Gl 758 with the NIRC2 camera in its narrow mode (102 × 102 field of view) using natural guide star adaptive optics (Wizinowich 2013) at the Keck II telescope on four occasions: UT 2010 May 2, UT 2013 July 3, UT 2016 June 27, and UT 2017 October 10. All observations were taken in the pupil-tracking mode to employ the angular differential imaging (ADI) method (Marois et al. 2006). The star was placed behind the partly transparent 600 mas coronagraph, which has an attenuation of 7.51 ± 0.14 mag at 1.6 μm (Bowler et al. 2015a) and enables precise image registration. The total on-source integration time of our observations was 49, 14, 35, and 30 minutes for our 2010, 2013, 2016, and 2017 epochs, respectively, and the total sky rotation of these sequences was 65°, 13°, 59°, and 48°. Images were corrected for cosmic rays and bad pixels, then dark subtracted and flat fielded. Details about the observations are summarized in Table 2. Further processing of the images is described in Section 4.2.

PSF subtraction for the NIRC2 imaging data was carried out with the Locally Optimized Combination of Images algorithm (LOCI; Lafrenière et al. 2007) using the ADI processing pipeline described in Bowler et al. (2015a). Images were individually corrected for geometric distortions by bilinearly interpolating pixel values to the rectified locations based on the solution from Yelda et al. (2010) for observations taken prior to 2015 April when the Keck II AO system was re-aligned, and the solution from Service et al. (2016) was used for observations taken after pupil realignment. Each frame was then registered by fitting a 2D elliptical Gaussian to the host star located behind the partly transparent coronagraph spot. Two reductions were carried out using aggressive and conservative implementations of LOCI by varying the angular tolerance parameter used to select PSF templates (N_δ). Two point sources are recovered with high significance in the 2010 and 2013 epochs (Gl 758 B and bkg1; Figure 2), and three point sources are recovered in the 2016 and 2017 epochs (Gl 758 B, bkg1, and bkg2; Figure 3) in both of the reductions. We adopt the aggressive implementation for all data sets with LOCI geometric parameters of W = 5, N_A = 300, g = 1, N_δ = 0.5, and dr = 2 following the definitions in Lafrenière et al. (2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The S/N for each point source is calculated using aperture photometry with a 5 pix aperture radius. The sky background is subtracted from the summed flux centered on the source using the mean of 100 sky measurements at the same angular radius but spanning a range of azimuthal angles surrounding (but not overlapping) the object of interest. The standard deviation of these sky values represents the background noise level, and the ratio of these two is used to determine the signal to noise of the detection. Gl 758 B and the two nearby background stars are detected with S/Ns between 10 and 51; full details can be found in Table 2.

PSF subtraction biases astrometry of point sources in processed images as a result of both over-subtraction and self-subtraction. To mitigate these effects, we follow the strategy outlined by Marois et al. (2010) of injecting a negative PSF template into the raw data and iteratively identifying the true position and flux of the sources. Three parameters were optimized using the downhill simplex AMOEBA algorithm (Nelder & Mead 1965)—the separation, position angle, and amplitude of the PSF template—to minimize the resulting rms in a 20 pix aperture radius at the location of the point source in the processed image. Although the host star is visible behind the coronagraph in the science frames, the mask transmission has historically been difficult to characterize in detail and may be non-uniform across the face of the occulting spot. To avoid using a potentially distorted PSF of the host star, we instead utilize unsaturated PSF templates of other stars taken in the same filter (see Table 2).

Results of the negative injection are shown in Figure 4. This procedure successfully removes most of the flux and over-subtracted azimuthal "wings," leaving only slight residual structure that likely originates from an imperfect PSF template, changing atmospheric conditions and Strehl ratios throughout the sequence (which are not taken into account in the modeling) and/or slight blurring of the PSF if substantial rotation occurs during individual exposures—something that preferentially affects sources at wider separations.

Zoom In Zoom Out Reset image size

The astrometric error budget is dominated by measurement errors in the positions of point sources; uncertainty in the residual optical distortion correction; errors in the plate scale and north orientation of the detector; and azimuthal shear caused by sky rotation within individual frames.

Our strategy for estimating point source measurement uncertainties is inspired by the method described in Rajan et al. (2017) and is intended to mimic the inverted PSF template procedure we carried out in Section 4.2. For Gl 758 B, bkg1, and bkg2, we inject a positive PSF template into the raw images at the same separation and amplitude as the objects of interest but a different position angle. We then iteratively inject an inverted PSF template of a different star into the raw data and perform local PSF subtraction with LOCI, then use AMOEBA to identify the optimal position and amplitude that minimizes the local noise. Each injection/recovery provides an estimate of the systematic difference in the separation, P.A., and amplitude of the injected (positive) object compared with what was recovered with the inverted (negative) PSF. This process is repeated ten times at equally spaced position angles for each object, and the average positional differences (${\sigma }_{\bar{\rho },\mathrm{meas}}$ and ${\sigma }_{\bar{\theta },\mathrm{meas}}$) are adopted as estimates of the positional measurement errors.

After correcting images for optical distortion effects, there remain small residual systematic positional uncertainties, σ_d, of about 1 mas which limit the achievable astrometric accuracy across the detector (Yelda et al. 2010; Service et al. 2016). Here we adopt one σ_d term associated with the host star and one for the companion. In addition, the NIRC2 plate scale, s, and its associated uncertainty, σ_s, are taken into account and vary slightly between pre- and post-pupil realignment (9.952 ± 0.002 mas pix⁻¹ from Yelda et al. 2010; 9.971 ± 0.004 mas pix⁻¹ from Service et al. 2016).

The final separation measurement in mas is

$$\begin{eqnarray}&&\rho =s{\bar{\rho }}_{\mathrm{meas}}\,\pm \,s{\bar{\rho }}_{\mathrm{meas}}{\left({\left(\displaystyle \frac{{\sigma }_{s}}{s}\right)}^{2}+{\left(\displaystyle \frac{{\sigma }_{\bar{\rho },\mathrm{tot}}}{{\bar{\rho }}_{\mathrm{meas}}}\right)}^{2}\right)}^{1/2},\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{\bar{\rho },\mathrm{tot}}$ is the combined uncertainty from our injection-recovery exercise and from the imperfect distortion correction:

$$\begin{eqnarray}&&{\sigma }_{\bar{\rho },\mathrm{tot}}^{2}={\sigma }_{\bar{\rho },\mathrm{meas}}^{2}+2{\sigma }_{d}^{2}.\end{eqnarray} \tag{ 2 }$$

The P.A. is determined as follows:

$$\begin{eqnarray}&&\theta ={\bar{\theta }}_{\mathrm{meas}}-{\theta }_{\mathrm{North}}+{\theta }_{\mathrm{shear}}/2,\end{eqnarray} \tag{ 3 }$$

where θ_North is the rotational offset required to align the NIRC2 detector columns with north on the sky: 0.252 ± 0009 for the Yelda et al. (2010) distortion solution, and 0.262 ± 0002 for the Service et al. (2016) distortion solution.¹³ θ_shear is the shear (blurring) per individual frame. To account for this, each frame is de-rotated to the midpoint P.A. of the exposure after PSF subtraction, and prior to coaddition of the sequence.

The error in the P.A., σ_θ, includes the injection-recovery measurement uncertainty (σ_θ,meas), uncertainty in the north alignment (σ_θ,North), residual positional errors after applying the distortion solution (σ_θ,d), and the systematic error from shearing of point sources from sky rotation within each frame (σ_θ,shear):

$$\begin{eqnarray}&&{\sigma }_{\theta }={({\sigma }_{\theta ,\mathrm{meas}}^{2}+{\sigma }_{\theta ,\mathrm{North}}^{2}+{\sigma }_{\theta ,d}^{2}+{\sigma }_{\theta ,\mathrm{shear}}^{2})}^{1/2}.\end{eqnarray} \tag{ 4 }$$

The residual positional distortion errors are about 1 mas, so here we approximate σ_θ,d as ≈1 mas/ρ. The dominant term in the P.A. error budget is the shear per frame, which varies among individual frames and across observation epochs. For this work we conservatively adopt half the average shear for each epoch: 030, 021, 038, and 032 for our 2010, 2013, 2016, and 2017 epochs. Our final astrometry of Gl 758 B and the two background sources is listed in Table 2. Note that the resulting astrometry does not appear to be significantly sensitive to changes in the LOCI parameters used for PSF subtraction based on the same astrometric analysis with N_δ set to 1.5.

Gl 758 B has been observed by many other telescopes and instruments over the past decade (Figure 5). The companion displays clear orbital motion; its separation has contracted from 188 in 2009 to 159 with our latest epoch from NIRC2 in 2017, and has moved by ≈16° in P.A. during that time. Our astrometry of Gl 758 B is broadly consistent with published values, although the separations from Nilsson et al. (2017) are significantly smaller than our measurements and those of Vigan et al. (2016) taken over the same time period. For example, our NIRC2 observations from 2013 were taken within three weeks of the 2013 July 21 data set obtained by Nilsson et al., but these two separation measurements are discrepant at the 4.3-σ level. However, the P.A. measurements from Nilsson et al. are in much better agreement.

Zoom In Zoom Out Reset image size

Before carrying out a detailed joint orbit fit of the RVs and astrometry, we first demonstrate here that the observed acceleration of Gl 758 is consistent with and likely to be caused by the companion Gl 758 B. The minimum mass of an imaged companion needed to produce an observed instantaneous acceleration ${\dot{v}}_{r}$ is

$$\begin{eqnarray}&&M\approx 0.0145{\left(\displaystyle \frac{d}{\mathrm{pc}}\displaystyle \frac{\rho }{{\prime\prime} }\right)}^{2}\left|\displaystyle \frac{{\dot{v}}_{r}}{{\rm{m}}\ {{\rm{s}}}^{-1}\ {\mathrm{yr}}^{-1}}\right|{M}_{\mathrm{Jup}},\end{eqnarray} \tag{ 5 }$$

where d is the distance to the system in pc and ρ is the projected separation in arcseconds (Torres 1999; Liu et al. 2002). Note that the generalized form of this equation includes information about the orbital elements of the system in the form of a multiplicative constant, the minimum value of which (≈2.6) is included here in the prefactor. The measured range of accelerations and angular separations of Gl 758 B implies a corresponding mass range of ≈20–50 M_Jup. Are these reasonable values for Gl 758 B? Evolutionary models from Saumon & Marley (2008) suggest that a brown dwarf with those masses should have effective temperatures between about 430 and 1400 K for ages of 1–6 Gyr; this is in good agreement with the inferred effective temperatures of 600–750 K for Gl 758 B from multi-band imaging and spectroscopy (Vigan et al. 2016; Nilsson et al. 2017). Based on this consistency and the fact that the companion's orbital period must be much longer than the time baseline of the RV observations (≫20 years), we conclude that Gl 758 B is likely the source of the acceleration.

We performed a joint orbit analysis of the three RV data sets and astrometry via a Markov Chain Monte Carlo (MCMC) algorithm. For this analysis we only use our NIRC2 astrometry to avoid systematic errors that may be present in previously published astrometry caused by multiple instruments and PSF subtraction strategies. We used the parallel-tempering (PT) ensemble sampler in emcee v2.1.0 (Foreman-Mackey et al. 2013) which is based on the algorithm described by Earl & Deem (2005). Thirty temperatures were adopted, of which only the coldest chain describes the posterior, together with 100 walkers to sample our 15-parameter model. Six of those parameters describe the orbit: semimajor axis (a), inclination (i), P.A. of the ascending node (Ω), mean longitude (λ_ref) at a reference epoch (t_ref) of 2455197.5 JD, and finally eccentricity (e) and the argument of periastron (ω, for the host star) parameterized as $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$, which avoids the Lucy–Sweeney bias toward non-zero eccentricities and imposes a uniform prior on eccentricity. We assumed a log-flat prior on a, randomly distributed viewing angles for i (i.e., a prior of sini), and uniform priors for the other orbit parameters. The next three parameters used in the fit were the parallax (π) and mass of the host star (M_host) and the mass of the companion (M_comp). We assumed priors of 63.45 ± 0.35 mas on the parallax (van Leeuwen 2007), 0.97 ± 0.02 M_⊙ for M_host (Vigan et al. 2016), and a log-flat prior for M_comp, which is motivated by the broad range of potential masses for the companion spanning ≈10 M_Jup (if the system is younger than expected) to over 100 M_Jup (if the companion is an unresolved binary). The remaining six parameters are the zero points (Δ_zero,McD, Δ_zero,Keck, Δ_zero,APF) and jitters (σ_jit,Keck, σ_jit,APF, σ_jit,McD) for the three RV data sets. The zero points are simply the offsets needed to bring the sets of relative RVs into accord with the orbit model, and the jitter terms account for small random and systematic RV epoch-to-epoch measurement errors from the star and the instrument not captured in the quoted relative RV uncertainties. We assumed uniform priors for the zero points and log-flat priors for the jitters. These are summarized in Table 3, and our complete likelihood function is as follows:

$$\begin{eqnarray}\mathrm{ln}({ \mathcal L }) & = & -0.5\left(\displaystyle \sum _{k=1}^{{N}_{\mathrm{ast}}}{\left(\displaystyle \frac{{\rho }_{k}-\rho ({t}_{k})}{{\sigma }_{\rho ,k}}\right)}^{2}+\displaystyle \sum _{k=1}^{{N}_{\mathrm{ast}}}{\left(\displaystyle \frac{{\theta }_{k}-\theta ({t}_{k})}{{\sigma }_{\theta ,k}}\right)}^{2}\right.\\ & & \left.+\displaystyle \sum _{j=1}^{{N}_{\mathrm{inst}}}\displaystyle \sum _{k=1}^{{N}_{\mathrm{RV}}}\displaystyle \frac{{({\mathrm{RV}}_{\mathrm{rel},k}+{{\rm{\Delta }}}_{\mathrm{zero},j}-\mathrm{RV}({t}_{k}))}^{2}}{{\sigma }_{\mathrm{RV},k}^{2}+{\sigma }_{\mathrm{jit},j}^{2}}\right)\\ & & +\mathrm{ln}(\sin (i))-\mathrm{ln}(a)-\mathrm{ln}({M}_{\sec })\\ & & -0.5{\left(\displaystyle \frac{\pi -63.45\mathrm{mas}}{0.35\mathrm{mas}}\right)}^{2}.\end{eqnarray} \tag{ 6 }$$

Table 3.  MCMC Posteriors for the Orbit of Gl 758B

Property	Mode ± 1σ	Best fit	95.4% c.i.	Prior
Fitted parameters
Companion mass Mcomp (MJup)	${42}_{-7}^{+19}$	55	33, 106	1/M (log-flat)
Host-star mass Mhost (M⊙)	${0.967}_{-0.018}^{+0.022}$	0.969	0.929, 1.009	0.970 ± 0.020 M⊙ (Gaussian)
Parallax (mas)	${63.39}_{-0.32}^{+0.37}$	63.51	62.73, 64.13	63.45 ± 0.35 mas (Gaussian)
Semimajor axis α (mas)	${1340}_{-80}^{+170}$	1350	1210, 1820	$1/\alpha$ (log-flat)
Inclination i (°)	${28}_{-10}^{+12}$	27	10, 49	$\sin (i)$, $0^\circ \lt i\lt 180^\circ$
$\sqrt{e}\sin \omega$	$-{0.05}_{-0.09}^{+0.11}$	−0.06	−0.26, 0.24	uniform
$\sqrt{e}\cos \omega$	$-{0.76}_{-0.03}^{+0.08}$	−0.75	−0.82, −0.47	uniform
Mean longitude at tref = 2455197.5 JD, ${\lambda }_{\mathrm{ref}}$ (°)	72 ± 8	72	52, 89	uniform
PA of the ascending node Ω (°)	175 ± 5	174	163, 183	uniform
McDonald RV zero point (m s−1)	61 ± 19	59	27, 108	uniform
Keck RV zero point (m s−1)	${50}_{-19}^{+20}$	50	17, 98	uniform
APF RV zero point (m s−1)	${44}_{-21}^{+18}$	42	9, 90	uniform
McDonald RV jitter σMcD (m s−1)	${3.1}_{-0.6}^{+0.8}$	3.0	1.6, 4.6	1/σ (log-flat)
Keck RV jitter σKeck (m s−1)	${2.33}_{-0.18}^{+0.17}$	2.32	2.00, 2.70	1/σ (log-flat)
APF RV jitter σAPF (m s−1)	${2.46}_{-0.17}^{+0.18}$	2.44	2.13, 2.84	1/σ (log-flat)
Computed properties
Orbital period P (yr)	${96}_{-9}^{+21}$	97	79, 153	⋯
Semimajor axis a (au)	${21.1}_{-1.3}^{+2.7}$	21.3	18.9, 28.7	⋯
Eccentricity e	${0.58}_{-0.11}^{+0.07}$	0.57	0.26, 0.67	⋯
Argument of periastron ω (°)	184−9+8	184	153, 201	⋯
Time of periastron ${T}_{0}={t}_{\mathrm{ref}}-P\tfrac{\lambda -\omega }{360^\circ }$ (JD)	${2465800}_{-800}^{+2000}$	2466300	2464800, 2470500	⋯
Mass ratio q = Mcomp/Mhost	${0.042}_{-0.008}^{+0.018}$	0.054	0.032, 0.105	⋯

Download table as:  ASCIITypeset image

We are able to fit for the RV jitter because of the numerous independent data points that sample its orbit in each data set. In other words, there are many degrees of freedom in the RV model of the host star. In contrast, with only four epochs of companion astrometry we do not have the same ability to fit for additional astrometric errors. Therefore, we performed an initial orbit fit using the nominal astrometric errors and examined the residuals. The rms of the separation measurements was 4.3 mas about this initial best-fit orbit and χ² = 14.6, while for the P.A.s the rms was 034 with χ² = 3.0. As a point comparison, when we simply fit a line to separation and P.A. as a function of time, we found similar rms values of of 4.7 mas (χ² = 17.5) and 045 (χ² = 4.9). Given that the RVs constrain some of the same orbit parameters that are relevant to the astrometric fit, it is not obvious what number of degrees of freedom is correct to assume here. If we assume two degrees of freedom in each, then p(χ²) = 0.0007 for the separations and p(χ²) = 0.22 for the P.A.s. Ultimately we add 4.3 mas in quadrature to our separation measurements, resulting in effective errors of ≈5 mas at every epoch. We do not add any additional error to our P.A. uncertainties as these are already substantially larger (≈03, or ≈10 mas) and the χ² value is not unreasonable. The source of the 4–5 mas epoch-to-epoch uncertainties in our separation measurements is not known, so it may represent a fundamental floor to astrometry derived from ADI sequences with the NIRC2 coronagraph.

The initial state of the PT-MCMC sampler was determined using a Monte Carlo rejection sampling analysis similar to the method used in Dupuy et al. (2016). First, 2 × 10⁶ randomly distributed orbital periods (10⁴ d < P < 10⁷ d), eccentricities, and times of periastron passage (T₀) were drawn. Using the formalism of Lucy (2014), we computed the corresponding set of a, i, ω, and Ω that best fit the astrometry for each of these trials. For each trial, we computed the χ²_ast of the trial orbit's predicted astrometry and our measured astrometry. To incorporate the RVs, we assumed at this stage that each data set could be represented as a simple linear trend with time. For each orbit trial, we computed the instantaneous slope of the host-star RV at the mean epoch for each RV data set. (In our actual PT-MCMC runs, we fit all the individual relative RV measurements directly.) Because each trial has a and P independent of assumptions about mass, each trial also effectively samples an associated total mass through M_tot ∝ a³/P². We computed the mass ratio that would best bring each orbit's RV slopes into agreement with the measured slopes and then computed the RV zero points needed to bring the measured relative RVs into agreement. Because there are multiple RV slopes to reproduce, the agreement is not perfect for a given orbit trial, and we computed the ${\chi }_{\mathrm{RV}}^{2}$ of the trial orbit's predicted RV slopes and our measured RV slopes. Finally, we computed the ${\chi }_{\mathrm{mass}}^{2}$ of the trial orbit's predicted host-star mass with the estimated mass of 0.97 M_⊙ from Vigan et al. (2016), with an inflated uncertainty of ±0.20 M_⊙ to allow us to perform orbit fits with no mass prior as well. We combined these constraints into ${\chi }_{\mathrm{tot}}^{2}={\chi }_{\mathrm{ast}}^{2}+{\chi }_{\mathrm{RV}}^{2}\,+{\chi }_{\mathrm{mass}}^{2}$, computed rejection probabilities of ${p}_{\mathrm{rej}}=1-\exp (-({\chi }_{\mathrm{tot}}^{2}-\min ({\chi }_{\mathrm{tot}}^{2}))/2)$, and then drew random samples to pass on based on ${p}_{\mathrm{rej}}\gt { \mathcal U }(0,1)$, where ${ \mathcal U }(0,1)$ was a uniformly distributed, randomly drawn number ranging from 0 to 1.

In our PT-MCMC analysis, we experimented with different chain lengths and found that after ∼10⁵ steps our 100-walker chains had clearly stabilized in the mean and rms of the posterior for each of the parameters. We saved every 20th step of our chains and discarded the first 50% of the chain as the burn-in portion, leaving 2.5 × 10⁵ PT-MCMC samples in the cold chain. Table 3 lists information on the posterior distributions of our fitted parameters, as well as parameters that are directly computed from them. To compute the modes of our distributions we binned the posterior and found the bin with the most elements. The 1-σ and 2-σ confidence intervals are computed as the minimum range in that parameter that contains 68.3% and 95.4% of the values, respectively. The quoted best-fit solution is the one with the maximum likelihood, which includes the prior.

Figure 6 displays the companion mass posterior and the most relevant parameter correlations, and Figure 7 shows several of the other marginalized posterior distributions from our fit. As expected, the inclination is highly correlated with companion mass (i.e., M_comp sin(i) is well constrained from the RV orbit). The companion mass posterior extends to high masses (>80 M_Jup) that are likely unphysical assuming that Gl 758 B is a single object. This high-mass tail corresponds to low inclinations (i ≲ 20°) that our astrometry cannot rule out, high eccentricities (e ≳ 0.6), and short periods (P ≲ 100 yr). The companion mass posterior has a sharp lower limit of 30.5 M_Jup at the 4-σ level.¹⁴ Figures 8 and 9 show our orbit solutions relative to our RV and astrometric data. The sky-projected and de-projected solution for Gl 758 is displayed in Figure 10, where orbits are drawn from the MCMC posteriors and are color-coded according to the corresponding companion mass from low mass (pink) to high mass (green).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Vigan et al. (2016) identified a new point source near Gl 758 at a separation of ≈11 based on observations taken with SPHERE in 2014. They found that the photometry of this new object is broadly consistent with the colors of L dwarfs, raising the possibility that this could be a second companion in this system. This object is easily recovered in our 2016 and 2017 data sets at wider separations of 146 and 162. Together with the single-epoch detection from Vigan et al. taken in 2014, this object closely follows the expected motion for a background star (Figure 11) and appears to be unassociated with Gl 758.

Zoom In Zoom Out Reset image size

Despite extensive efforts to characterize this system with follow-up photometry and spectroscopy, only a single bolometric luminosity estimate by Currie et al. (2010) exists in the literature: $\mathrm{log}(L/{L}_{\odot })=-{6.1}_{-0.2}^{+0.3}$ dex. To improve on this value we use use existing near-infrared measurements to constrain the 1–3 μm spectral energy distribution together with atmospheric models for a bolometric correction. We first anchor the 1.0–1.75 μm spectrum of Gl 758 B from Nilsson et al. (2017) by flux calibrating the P1640 observations to the H-band apparent magnitude from Thalmann et al. (2009). To this we add the photometry from Vigan et al. (2016) and Janson et al. (2011) to directly account for the comparably high K-band flux from this object. A solar-metallicity BT-Settl "CIFIST2011" atmospheric model with T_eff = 650 K and $\mathrm{log}g$ = 5.0 dex is used for the long-wavelength bolometric correction (2.5–500 μm) by flux-calibrating the model to the L'-band photometry from Janson et al. (2011). This approach also agrees with the M_S-band upper limit from Janson et al. (2011). The same model is used for the short-wavelength correction (0.1–1.0 μm) by scaling that region to the blue end of the P1640 spectrum. Uncertainties in the spectral measurements, photometry, and flux calibration scale factors for the model and spectrum are all accounted for in a Monte Carlo fashion by integrating under new realizations of the complete 0.1–500 μm spectrum. This procedure yields a bolometric luminosity of $\mathrm{log}(L/{L}_{\odot })$ = −6.07 ± 0.03 dex for Gl 758 B. To assess possible systematic errors, we experimented with alternative atmospheric models from the same grid with effective temperatures of 600 K and 700 K. The results following the same procedure are within 0.02 dex of the value we obtained with the 650 K model, which is smaller than the impact of random measurement errors.

With a measured luminosity, age, and dynamical mass, Gl 758 B offers a rare opportunity to test substellar evolutionary models. For this analysis we begin with the assumption that the age range spans 1–6 Gyr following Vigan et al. (2016), but ultimately re-evaluate this constraint based on recent results from isochrone fitting. We select a variety of publicly available models from the literature for this exercise: the Cond models from Baraffe et al. (2003); three versions of evolutionary models from Saumon & Marley (2008) with no clouds ("SM-NC"), a hybrid prescription for the evolution of clouds at the L/T transition ("SM-Hybrid"), and the retention of thick clouds at all temperatures ("SM-f2"); and the grid from Burrows et al. (1997). All have solar compositions. These models mainly differ in their treatment of atmospheric clouds and molecular opacities, which act as boundary conditions that control the evolution of brown dwarfs and giant planets as these objects radiatively cool over time (see, e.g., Burrows et al. 2001 and Marley & Robinson 2015 for detailed reviews).

Our approach for comparing the models to the observations utilizes a one-tailed hypothesis test. We adopt a null hypothesis in which the posterior probability density function for the dynamical mass of Gl 758 B is statistically consistent with the inferred mass distribution from evolutionary models at some threshold probability; we choose 0.3% (within 3-σ) for this study. In other words, we calculate the probability that random draws from the dynamical mass distribution differ from the inferred model-based mass distribution. If these two values disagree by at least 0.997, then the null hypothesis is rejected and the two distributions are considered to be inconsistent with each other.

For each set of evolutionary models, we randomly draw ages from a uniform distribution (τ = ${ \mathcal U }(1,6)$ Gyr) and luminosities from a log-normal distribution ($\mathrm{log}(L/{L}_{\odot })$ = ${ \mathcal N }(\mu$ = −6.07, σ = 0.03) dex); a visual reference is shown in Figure 12 for the Burrows models. For each {τ, L} pair we find the corresponding mass by finely interpolating the model grid. This Monte Carlo process is repeated of order 10⁶ times to create a distribution of expected masses for each model. The predicted and dynamical mass distributions are then quantitatively compared for consistency.

Zoom In Zoom Out Reset image size

Results from this analysis are listed in Table 4 and illustrated in Figure 13. Based on the input age and luminosity distributions together with our threshold criterion for agreement, only the Burrows models are formally consistent with the dynamical mass distribution, which peaks at 42 M_Jup and has a robust (4-σ) lower limit of 30.5 M_Jup. However, even the formal agreement with the Burrows models is marginal and a slightly lower threshold for consistency would have rejected the null hypothesis; random draws from the dynamical mass distribution result in higher masses 99% of the time. All other models fail our hypothesis test. The SM-f2 models disagree the most with the dynamical mass of Gl 758 B. By retaining clouds to temperatures well below the rainout limit for various grains, this prescription is the most unrealistic for T dwarfs like Gl 758 B (which is expected to be cloud-free) so this result is unsurprising.

Zoom In Zoom Out Reset image size

Table 4.  Predictions from Evolutionary Models

Constraint			Model
Age	Luminosity	Mass
(Gyr)	($\mathrm{log}(L/{L}_{\odot })$)	(MJup)	Cond	SM-NC	SM-Hybrid	SM-f2	Burrows
			Predicted Mass (95.4% Credible Interval)
${ \mathcal U }$(1, 6)	${ \mathcal N }(-6.07,0.03)$	⋯	14–33 MJup	13–33 MJup	13–32 MJup	11–27 MJup	15–38 MJup
			Minimum Compatible Age (99.7% Lower Limit)
⋯	${ \mathcal N }(-6.07,0.03)$	PDF(Mcomp)	>6.2 Gyr	>6.1 Gyr	>6.4 Gyr	>9.2 Gyr	>4.4 Gyr
			Minimum Compatible Luminosity (99.7% Lower Limit)
${ \mathcal U }$(1, 6)	⋯	PDF(Mcomp)	>−5.96 dex	>−5.95 dex	>−5.92 dex	>−5.76 dex	>−6.13 dex

Note. Predictions from evolutionary models for various permutations of input age, luminosity, and mass distributions. Here ${ \mathcal U }$(a, b) refers to a linearly uniform distribution from a to b, ${ \mathcal N }$(μ, σ) is a normal distribution with a mean μ and standard deviation σ, and PDF(M_comp) is our measured probability density function for the dynamical mass of Gl 758 B.

Download table as:  ASCIITypeset image

The general tension between the models and our dynamical mass measurement may point to physics not yet incorporated into current substellar models, which could originate from several sources: interior physics and thermal structure, sources of atmospheric opacity and their evolution with temperature, initial entropy and accretion history of Gl 758 B, or ill-matched metallicities of the models and the companion. On the other hand, it is also possible that the discrepancy originates from the observational side, most likely with the age of the system. For our default analysis we adopted the 1–6 Gyr estimate by Vigan et al. (2016) based on isochrone fitting of the host star (which resulted in younger ages of ≈1–4 Gyr) and activity indicators (which resulted in older ages of ≈3–8 Gyr). Older ages result in higher predicted masses for the same luminosity, so this could also be a natural explanation for the disagreement.

To explore the possibility of an older age or a systematic offset in the luminosities of the models, we perform a series of tests to identify the minimum compatible ages and the minimum compatible luminosities that render the inferred and dynamical mass distributions into agreement. For the former, we begin by randomly drawing masses from interpolated evolutionary model grids following a normal distribution in log-luminosity (${ \mathcal N }(-6.07,0.03)$ dex) and fixing the starting age at 1.0 Gyr. Using the same threshold requirement of 0.3%, we compare the inferred and dynamical mass distributions for consistency. This process is then repeated by increasing the age in 0.1 Gyr steps until the mass distributions agree at the threshold level. The results are summarized in Table 4; the Cond, SM-NC, and SM-Hybrid models all imply similar ages of ≳6 Gyr. The SM-f2 grid is only consistent with the dynamical mass for extremely old ages of ≳9 Gyr, and the Burrows models agree for ages beyond 4.4 Gyr.

A similar process is carried out to identify the minimum luminosity consistent with the measured mass. We randomly draw masses based on a uniform distribution of ages (${ \mathcal U }(1,6)$ Gyr) and a fixed starting luminosity of −6.50 dex, then test the inferred and dynamical mass distributions for consistency. The luminosity is increased in increments of 0.01 dex until agreement is reached. If the 1–6 Gyr age estimate is correct, that would imply that the Cond, SM-NC, SM-Hybrid, and SM-f2 evolutionary models are overluminous by 0.11 dex, 0.12 dex, 0.15 dex, and 0.31 dex, respectively (Table 4). This potential discrepancy is in the opposite sense from results by Dupuy et al. (2009, 2014), who found that substellar cooling models under-predict the luminosities of brown dwarfs with dynamical masses by ≈0.2–0.4 dex, at least at relatively young age of ≈0.5–1 Gyr.

Altogether, the most likely culprit for the disagreement in mass probably resides in the age of Gl 758. Older ages of 6–9 Gyr would readily put the predicted and dynamical distributions in excellent agreement and are indeed suggested from the low activity level, lack of X-ray emission, and slow projected rotational velocity of Gl 758 (Mamajek & Hillenbrand 2008; Thalmann et al. 2009; Vigan et al. 2016). Although there are a wide range of age estimates for the host star from isochrone fitting in the literature, more recent analyses are converging on an older value that agrees better with activity indicators. For example, a recent study by Luck (2017) found an average age of 7.5 Gyr (with a range of 5.3 Gyr about that value) using four sets of isochrones, and Brewer et al. (2016) found an isochronal age of 7.5 Gyr with a range of 4.6–10.4 Gyr using the Yonsei–Yale models.

We searched the residual RVs for closer-in planet candidates using a Lomb–Scargle periodogram after removing the best-fit orbital solution of Gl 758 B. The strongest power from 1–10⁴ days is at a period of 245.5 days, but that potential signal has a false alarm probability of 0.4% and its corresponding velocity semi-amplitude is at the level of the noise in the data, so it is unlikely to be real. No frequencies have powers that exceed our threshold false alarm probability of 0.1% for planet detection. We conclude that there is no convincing evidence of any close-in planet candidates in our data.

Detection limits are quantified using injection-recovery tests as described in Howard & Fulton (2016). Synthetic planets on circular orbits are sequentially injected into the RV residuals by randomly drawing pairs of minimum mass and period surrounding the detection threshold. A periodogram is used to search for planets within each artificial data set with a 1% false alarm probability threshold for recovery and the requirement of a similar period and phase as the injected planet. Results are shown in Figure 14. Gl 758 is devoid of close-in giant planets (≳100 M_⊕) within 10 au, sub-Saturns (≈10–100 M_⊕) within 3 au, and super-Earths (≈2–10 M_⊕) interior to 0.1 au.

Zoom In Zoom Out Reset image size

Given its periastron passage of ${8.9}_{-2.0}^{+3.7}$ au, it is likely that Gl 758 B has impacted the formation efficiency and dynamical stability of closer-in planets in this system. If Gl 758 B formed relatively quickly (≲1 Myr), perhaps from turbulent fragmentation (e.g., Bate 2009) with subsequent migration to its current orbit, this implies that the circumprimary protoplanetary disk would have been truncated between ∼0.2 and 0.35a, or about 4–7 au (Artymowicz & Lubow 1994). The lack of planets outside of this region is not unexpected, but planet formation interior to this region may still have been possible (e.g., Kepler-444; Campante et al. 2015; Dupuy et al. 2016).

However, the lack of planets at small orbital distances from Gl 758 is not particularly unusual compared to the statistical properties of planets orbiting GK dwarfs in general. For example, Cumming et al. (2008) find that about 10% of Sun-like stars host giant planets with minimum masses between 0.3 and 10 M_Jup within ≈3 au. This value increases to about 14% by extrapolating the planet period distribution out to 10 au. Wittenmyer et al. (2016) infer the frequency of Jupiter analogs between 3 and 7 au around solar-type stars to be ${6.2}_{-1.6}^{+2.8}$% (see also Wittenmyer et al. 2011; Zechmeister et al. 2013). Petigura & Howard (2013) measure a completeness-corrected frequency of about 55% for all planets orbiting GK stars between 1 and 12 R_⊕ and orbital periods from 5–100 days (0.05–0.4 au), with lower-mass planets outnumbering gas giants by a factor of ≈50:1 (see also Youdin 2011; Howard et al. 2012; Fressin et al. 2013). While it is possible that the absence of close-in planets could be related to the presence of Gl 758 B, this apparent desert is also broadly consistent with the overall statistical properties of single FGK stars.

Figure 15 shows the mass and separation regimes over which our imaging data and the residuals from our RVs are sensitive. Together this rules out giant planet and brown dwarf companions at close separations as well as massive companions at wide orbital distances. However, there exists a large region beyond about 10 au and at masses less than about 30 M_Jup where additional companions could evade detection, assuming an older age of about 6–9 Gyr. If a giant planet or another brown dwarf resides in this system—which remains possible both below our detection limits or simply at unfavorable viewing geometries—the acceleration we observe would be the superposition from one or more additional companions besides Gl 758 B. This could potentially influence the amplitude and shape of the evolving acceleration and may even reconcile the dynamical mass measurement and the younger age. At this point there are no signs of another companion, but continued RV monitoring and deeper high-contrast imaging would be beneficial to further map the architecture of this system.

Zoom In Zoom Out Reset image size