© ESO, 2016

1. Introduction

The direct-imaging search for substellar companions around nearby stars has led to an increasing number of discoveries in the vicinity of our Sun. The brown dwarf companion, GJ 758 B (Thalmann et al. 2009), is one of t that stands out from the list. The primary star is a nearby (15.76%; van Leeuwen 2007) solar-type (G9V) star, and the inferred effective temperature (T_eff) of GJ 758 B is among the lowest (~600 K) ever recorded for a directly imaged companion. These peculiarities made this system the subject of two separate studies (Currie et al. 2010; Janson et al. 2011) in addition to its discovery paper.

Previous observations have provided good spectral coverage of the GJ 758 system. Common proper motion of the companion with its parent star was determined through two epochs of Subaru/HiCIAO H-band observations detailed in Thalmann et al. (2009). In their work, they also highlight its very low T_eff (550–640 K) and late spectral type (T9). Currie et al. (2010) published MMT/Clio L′-band measurements of GJ 758 B. These data showed the object to have extremely red colors between near- and mid-infrared (H−L′ = 3.29 ± 0.25). The latest publication on GJ 758 B, by Janson et al. (2011), complete the spectral coverage with measurements in J, H, CH4S, CH4L, K_c, L′ and M_s from Subaru/HiCIAO, Gemini/NIRI and Keck/NIRC2. They confirm again the very low T_eff and late spectral type of the companion and, for the first time, they demonstrate the clear methane absorption in H-band from the NIRI measurements in the CH4S and CH4L filters. In general, all three papers converge towards a similar picture of a low mass brown dwarf (30–40 M_Jup), given the old age of the system (5–9 Gyr). First attempts at an orbit determination for GJ 758 B hinted at a large semi-major axis (30 ≤ a ≤ 90 AU) and high eccentricity (0.4 ≤ e ≤ 0.7). However, even though orbital motion is detected, further astrometric monitoring is needed for accurate orbital parameters to be determined. This is due to only a small fraction of the total orbit having been detected.

In this work we present new near-infrared (near-IR) photometric data obtained with the SPHERE instrument (Spectro-Polarimetric High-contrast Exoplanet REsearch; Beuzit et al. 2008), recently commissioned at the Very Large Telescope (VLT) in Chile. We first present our observations (Sect. 2), and the data reduction and analysis (Sect. 3). These new observations cover the full near-IR range at much higher contrast than previous observations, allowing us to detect a new candidate companion at a closer projected separation than GJ 758 B. We provide photometric measurements of GJ 758 B with improved sampling and resolution, including the very first measurements of the companion flux in Y-band. After revisiting the stellar parameters and age indicators for GJ 758 A (Sect. 4), we perform an updated modeling of the properties of GJ 758 B, making a comparison with empirical objects and atmospheric models (Sect. 5). Finally, we use the new astrometric datapoint to improve the orbit determination (Sect. 6) before concluding with our sensitivity to additional closer-in companions (Sect. 7).

2. Observations

The star GJ 758 was observed as part of the third SPHERE commissioning run in August 2014. The SPHERE planet finder instrument installed at the VLT (Beuzit et al. 2008) is a highly specialized instrument, dedicated to high-contrast imaging and spectroscopy of young giant exoplanets. It is based on the SAXO extreme adaptive optics system (Fusco et al. 2006; Petit et al. 2014; Sauvage et al. 2014), which controls a 41 × 41 actuators deformable mirror, and four control loops (fast visible tip-tilt, high-order, near-infrared differential tip-tilt, and pupil stabilization). The common path optics employ several stress polished toric mirrors (Hugot et al. 2012) to transport the beam to the coronagraphs and scientific instruments. Several types of coronagraphic devices for stellar diffraction-suppression are provided, including apodized pupil Lyot coronagraphs (Soummer 2005) and achromatic four-quadrants phase masks (Boccaletti et al. 2008).

The GJ 758 observations were acquired with one of the three scientific subsystems of SPHERE, the infrared dual-band imager and spectrograph (IRDIS; Dohlen et al. 2008) in its dual-band imaging mode (DBI; Vigan et al. 2010) with four different filter pairs in the Y-, J-, H- and K_s-bands. The spectral characteristics of the filters are provided in Table 1. The observations were performed in pupil-stabilized mode to perform angular differential imaging (ADI; Marois et al. 2006) with an apodized pupil Lyot coronagraph (Soummer 2005) optimized for the H-band (ALC_YJH_S), which uses a coronagraphic mask of diameter 185 mas. The data were acquired on two consecutive nights, 13 and 14 August, 2014, with a total integration time of ~26 min in each filter pair. The IRDIS detector was dithered on a 4 × 4 pattern to reduce the effect of the residual flat-field noise. At each detector dithering position, a data cube of DIT × NDIT = 32 × 3 s was acquired, resulting in a total of 16 data cubes for each observation.

The observing sequence in each of the DBI filters was performed as follows:

• one image of the point spread function (PSF) takenoff-axis (~0.4′′) with the neutral density ND3.5, which reduces the flux by a factor ~3000. The PSF is moved off the coronagraph by applying an offset on the near-IR differential tip-tilt plate. During this observation, the AO visible tip-tilt and high-order loops remain closed to provide a diffraction-limited PSF;

• a “star center” coronagraphic image where four symmetric satellite spots are created by introducing a periodic modulation on the deformable mirror. This data is used in subsequent analysis to determine an accurate position of the star center behind the coronagraph, and hence the center of field rotation;

• the coronagraphic sequence as previously described;

• an additional off-axis PSF to evaluate the variations of the observing conditions between the beginning and end of the sequence.

For commissioning purposes, data from the SPARTA real-time computer of the SAXO extreme AO system (Fusco et al. 2014) were collected at regular intervals in parallel of all the observations. This includes, in particular, images from the differential tip-tilt sensor (DTTS). This sensor removes a minute fraction of the incoming flux in the near-IR arm (at 1.6 μm) to image the PSF just before the coronagraph, and uses it as input for the DTTS loop that maintains the PSF that is locked on the coronagraph, once the observing sequence has started. Every 30 s, the 30 s average of the non-coronagraphic PSF on the DTTS is saved in the SPARTA files, allowing a fine monitoring of the PSF motion and flux variation at the level of the coronagraph.

Standard calibrations for the DBI mode were acquired in the morning as part of the IRDIS calibration plan. Instrumental backgrounds were taken for both the coronagraphic and off-axis exposures with proper DIT values. Detector flat fields were also acquired in each of the DBI filter pairs.

3. Data reduction and analysis

The data were analyzed using two separate pipelines, which are described in this section.

The LAM-ADI pipeline is similar to that described in Vigan et al. (2012), after updates to work with the SPHERE/IRDIS data. The calibrations (backgrounds, flat) were created using the preliminary release (v0.14.0-2) of the SPHERE data reduction and handling (DRH) software (Pavlov et al. 2008). Each of the images in the coronagraphic observing sequences were background subtracted and divided by the flat field in the appropriate DBI filters. Bad pixels were corrected using bad pixel maps created with the DRH by replacing them with the median of neighboring good pixels. Finally, all images were aligned to a common center using the star center data acquired at the beginning of the sequence. For this purpose, the four satellite spots inside the AO control radius were fitted with a 2D Gaussian function using the MPFIT non-linear least squares curve fitting software (Markwardt 2009). The accuracy of the centering using this procedure has been determined to be better than 0.1 pixel (~1.2 mas) for bright stars during the first SPHERE commissioning run in May 2014. For the recentering of the science frames, the shift introduced by the detector dithering procedure was also taken into account, and the 0.06 pixel (0.74 mas) accuracy of the dithering motion stage was included into the astrometric error budget. For each filter pair, the calibration process was applied independently to each of the two wavelengths that were acquired simultaneously with IRDIS, resulting in two separate pre-processed ADI data cubes.

The ADI data cubes were processed with the LAM-ADI pipeline, using a principal component analysis (PCA) implementation following the Karhunen-Loève image projection (KLIP) approach (Soummer et al. 2012). The number of subtracted modes, minimum and maximum radii for the analysis were varied over a wide range, but the companion was recovered in all analyses. Figure 1 shows the signal of the companion in all of the IRDIS DBI filters. The companion is recovered in all filters with a signal-to-noise ratio (S/N) greater than six, except in the K2 filter where it is only marginally detected with an S/N of ~2.5. As already presented in Janson et al. (2011), the companion displays a clear methane absorption in H-band with a flux about nine times fainter in H3 than in H2. Images were also processed using a combination of spectral differential imaging (SDI; Racine et al. 1999) and ADI to attenuate the speckle noise even more and look for additional fainter candidates. In addition to the detection of GJ 758 B, we report the re-detection of the background star already identified by Janson et al. (2011), and the detection of a new candidate located ~1.1′′ south of the star in all filters except K2. Although not directly detectable in H2 and H3 with only ADI, the candidate was easily identifiable in the SDI+ADI processed image.

The precise astrometry and photometry of the companion and new candidate was estimated using “negative fake companion” subtraction in the pre-processed ADI data cubes (Marois et al. 2010). A rough estimation of the object position and contrast is first performed using a 2D Gaussian fit. Then these initial guesses are used as a starting point for a Levenberg-Marquardt least-squares minimization routine, where the position and contrast of the negative fake companion are varied to minimize the residual noise after ADI processing in a circular aperture of radius λ/D that is centered on the position of the companion. When a minimum is reached, the position and contrast of the fake companion are taken as the optimal values for the astrometry and photometry. We note that this procedure is also applicable for analyses combining SDI and ADI, by minimizing the residuals in an aperture that covers the position of the companion in the first filter, and in the second filter after spatial rescaling. The error bars for the fitting process are then calculated by varying the position and contrast of the fake companion until the variation of the reduced χ² reaches a level of 1σ.

The data were analyzed independently with the LESIA pipeline for a cross-check of the astrometry and photometry. This pipeline uses a similar approach for the pre-processing of the ADI data cubes, but for the speckles subtraction it uses an upgrade of the Template Locally Optimized Combination of Images (TLOCI) algorithm derived from the one presented in Marois et al. (2014). Only ADI is used (no SDI) to avoid issues with the photometry calibration (Maire et al. 2014). Hence, for each dual-band filter sequence, it calibrates the speckle pattern in each individual frame, rotates the frames to align North up, and median-combines all the frames to obtain the final image. To derive the photometry and the astrometry of the detected sources, the pipeline uses the unsaturated PSF of the central star (recorded before and after the coronagraphic sequence) to build a data cube composed of fake companions at the positions of the detected sources that account for the field-of-view rotation in each frame and smearing during exposures. Then, the frames of this data cube are combined using the TLOCI coefficients that were used to obtain the image where the point-source was detected. The resulting frames are aligned in the same way as the science data to obtain an image that gives a model of the off-axis sources in the TLOCI images at the positions of the detections, accounting for TLOCI self-subtraction and distortions. Finally, the sub-pixel position and the flux of the modeled images are adjusted to optimize the subtraction of the model to the real image within a 1.5 λ/D-radius disk centered on the detection (Galicher & Marois 2011). The error bars account for the variations of the stellar flux during the sequence (estimated from the global speckle intensity variations), and the accuracy of the fitting of the companion image models to the real images.

For calibrating the distortion, plate scale, and orientation of the IRDIS images, a field in the outer regions of the 47 Tuc globular cluster was observed in different instrumental configurations (Maire et al. 2016). The 47 Tuc field was selected because it includes a bright star for adaptive optics guiding and was accurately calibrated using Hubble Space Telescope (HST) observations (Bellini et al. 2014). The plate scales for the different DBI filters are summarized in Table 3. Since it was not calibrated in the K12 filter pair during the commissioning, we assumed the same value as for the H23 filter pair. The true North correction measured for this commissioning run is −1.636 ± 0.013deg, and the correction of the orientation also takes into account the zero point orientation of the derotator in pupil-stabilized mode, which was measured to be 135.87 ± 0.03deg.

Relative photometry and astrometry of the companion and the newly detected candidate are reported in Table 4. Both pipelines agree within their respective error bars. The values reported in the table correspond to the average of the results from both pipelines, and the respective error bars have been quadratically added. The final error bars for the photometry include the fitting error detailed above, the variation of the non-coronagraphic PSF measured on the DTTS images (see Sect. 2), and the level of speckle residuals estimated at the same separation as the detections. The astrometric error bars include the fitting error, and the uncertainties on the star center, dithering motion, plate scale, derotator zero point, and true North correction. We note that for astrometry, the reference values are those from the H23 filter pair, which has been the most accurately calibrated.

4. Stellar parameters

A reassessment of stellar parameters of GJ 758 is warranted, when taking their relevance in the derivation of the properties of its substellar companion into consideration, and to explain its peculiar features discussed in Sect. 5. The star GJ 758 is classified as old (age 0.7–8.7 Gyr; Janson et al. 2011), and we revisit here the various age indicators, following the procedures and calibrations described in Desidera et al. (2015), as well as the chemical composition.

4.1. Kinematic parameters

Adopting the trigonometric parallax, the proper motion and error bars by van Leeuwen (2007), and the absolute radial velocity by Nidever et al. (2002) with an error of 0.50 km s^-1, space velocities U,V,W = −21.1 ± 0.2; −14.1 ± 0.5; −3.0 ± 0.2 km s^-1 are obtained. These are very similar to those of the Argus association (U,V,W = −21.5 ± 0.9, −12.2 ± 1.7, −4.6 ± 2.7). Although the BANYAN II on-line tool (Gagné et al. 2014) yields a membership probability of 97.8%, which would correspond to a very young age of 40 Myr, the full version of the BANYAN I tool (Malo et al. 2013), which takes into account both kinematic and photometric information, yields a 100% probability to the hypothesis that GJ 758 is an old field star.

4.2. Abundance analysis

We determined spectroscopic stellar parameters and chemical abundances for GJ 758 by exploiting a high-resolution (R = 42 000), high S/N (S/N = 164 at 5500 Å) ELODIE spectrum¹, which provides a wavelength coverage from 3850 Å to 6800 Å. The spectrum was downloaded from the online ELODIE archive (Moultaka et al. 2004), which provides reduced data products. This investigation aims at chemically tagging our target star, to ascertain whether the abundance pattern is compatible with the Argus association, whose chemical composition has been recently presented by De Silva et al. (2013). Argus reflects a roughly solar chemical composition with [Fe / H] = −0.06 ± 0.05 dex and all [X/Fe] ratios within 0.15 dex from the solar values, with the notable exception of barium (see discussion below).

We carried out a homogeneous and strictly differential analysis for GJ 758 with respect to Argus members published in that previous work by utilizing the same code (MOOG by Sneden 1973, 2014 version), line lists, techniques, and grid of model atmospheres (Kurucz 1993, solar-scaled models and no convective overshooting). Effective temperature (T_eff) and surface gravity (log  g) were derived by imposing excitation and ionisation equilibrium, so that there is no spurious trend of A(Fe) with the excitation potentials of the spectral features and agreement (within 0.05 dex) of iron abundances from Fe i and Fe ii, respectively. Instead, the microturbulence velocity (ξ) was calculated requiring that abundances from Fe i lines show no trend with reduced equivalent widths. We performed equivalent width analysis for Fe, Na, Mg, Al, Si, Ca, Ti, Cr, and Ni, whereas the Ba abundance was inferred via spectral synthesis, including hyperfine structure and isotopic splitting (see De Silva et al. (2013).

Internal (random) uncertainties affecting our derived abundances were computed in the standard way, that is by adding in quadrature errors resulting from the equivalent-width (EW) measurements (or to the best-fit determination in the case of spectral synthesis) and those related to the adopted set of atmospheric parameters (T_eff, log  g, and ξ). The total internal errors for [Fe/H], as well as for [X/Fe] ratios, are given in Table 5 (see De Silva et al. 2013 for further details on the error budget calculation).

We found a metallicity of [Fe / H] = 0.18 ± 0.05, which agrees very well with previous determinations by e.g., Soubiran et al. (2008), Takeda (2007), and Maldonado et al. (2012) and points to super-solar heavy element abundances for this star. The abundances of α-elements Si and Ca, as well as the Fe-peak Cr and Ni, match a solar-scaled pattern, whereas Na, Mg, Al, and Ti (though to a less extent) seem to exhibit a modest enhancement, albeit still consistent with solar abundances within the observational uncertainties. The metallicity distribution as a function of effective temperatures is shown in the left-hand panel of Fig. 2: we report [Fe/H] values for GJ 758 along with stars belonging to Argus (filled circles) and to the open cluster IC 2391 (triangles), deemed to share a common origin with the young association. It is clear from Fig. 2 that GJ 758 stands out from the cluster/association distribution, with its [Fe/H] being roughly ~0.25 dex higher.

Barium deserves a brief, separate discussion. First identified by D’Orazi et al. (2012), and subsequently confirmed by several studies (e.g., Yong et al. 2012; Jacobson & Friel 2013; Mishenina et al. 2013), the Ba abundance shows a decreasing trend with the open cluster’s age. The younger the cluster, the higher its Ba content. The reason for such a peculiar and unique pattern is still a matter for debate: it has been suggested that the efficiency in the production of the s-process elements in low-mass AGB stars is higher than that predicted for standard stellar evolution models, and input physics has still to be revised (D’Orazi et al. 2009; Maiorca et al. 2012). However, subsequent investigations show that the picture might not be that straightforward. The fact that the Ba overabundance is not accompanied by a similar behavior in other s-process elements (e.g., Y, La) makes this explanation unlikely. We refer the reader to (D’Orazi et al. 2012) for a wider discussion of this topic. Regardless of the nature of the super-solar Ba content, [Ba/Fe] ratios range from extremely high values of approximately ~0.6 dex for pre-main sequence clusters, such as e.g., IC 2602 and IC 2391 (D’Orazi & Randich 2009) to solar values, or even lower, for clusters a few Gyr old. De Silva et al. (2013) corroborated this observational evidence and obtained a mean abundance of [Ba / Fe] = 0.53 ± 0.03 (rms = 0.08 dex) for the Argus association and [Ba / Fe] = 0.62 ± 0.02 (rms = 0.07) for IC 2391 (see the right-hand panel of Fig. 2). Conversely, we gathered a [Ba / Fe] = 0.00 ± 0.12 for our star, which implies a difference in the Ba content of more than a factor of 3.5. Thus, in terms of chemical composition, Ba provides us with the strongest observational constraint: GJ 758 cannot have been born from the same molecular cloud as Argus.

4.3. Age indicators

The star, GJ 758, is known to have a low activity level that results from several measurements in the literature: log R_HK = −4.94 (Wright et al. 2004), –5.015 (Isaacson & Fischer 2010); and –5.060 (Duncan et al. 1991; Mamajek & Hillenbrand 2008). The calibration by Mamajek & Hillenbrand (2008) yields values of 5.5–7.7 Gyr for these activity values. The availability of multi-epoch measurements of chromospheric activity spanning several years indicate that this is not the result of a poor sampling of an activity cycle. The X-ray non-detection in the ROSAT All Sky Survey (Voges et al. 1999, 2000) (which would imply log L_X/L_bol< −5.8 and then an age > 3 Gyr), the small projected rotational velocity (0–2 km s^-1), and the small photometric variability (0.008 mag from Hipparcos) further support the low activity level of GJ 758, as expected for a few-Gyr old star.

Lithium is another highly sensitive age indicator for young stars. From the analysis of the spectrum described in Sect. 4.2, the Li 6708 Å resonance line is not detected, confirming the null result by Takeda & Kawanomoto (2005). For stars with GJ 758 colors, detectable amounts of lithium vanish at about the age of the Hyades. Therefore, the lack of lithium allows us to infer a stellar age that is older than 600 Myr.

While stellar members of young moving groups display significant scatter in the age indicators (see Desidera et al. 2011, for the case of Argus), we are not aware of late G-type stars, which are confirmed members of young moving groups, and which have such a low activity level and lack of lithium. The analysis of these indicators therefore converges with the chemical tagging in ruling out Argus membership for GJ 758.

Using the spectroscopic effective temperature and metallicity, and the Hipparcos V magnitude and trigonometric parallax, we derive age and masses from isochrone using the PARAM interface (da Silva et al. 2006)² and the stellar models by Bressan et al. (2012). Limiting possible input values to an age larger than 0.6 Gyr, the result of a lack of lithium, the resulting age would be 2.2 ± 1.4 Gyr and the stellar mass 0.97 ± 0.02 M_⊙.

4.4. Summary

All age indicators suggest that GJ 758 is an old star with lithium providing a tight lower limit at 600 Myr. Chemical tagging derived from a homogeneous comparison of abundances of several elements with those of confirmed members of the Argus association and IC 2391 open cluster also rules out a link between GJ 758 and Argus, with Barium abundance suggesting an age similar to the Sun. Therefore, we conclude that the kinematic parameters of GJ 758 are similar to those of the Argus association by chance, which confirmsing the statistical nature of kinematic ages and the need for independent youth indications to conclusively infer membership in young moving groups (Gagné et al. 2014; Desidera et al. 2015). The young disk kinematics decrease the probability of a star being significantly older than the Sun. The age of the system is likely to be within one to six Gyr, and the most probable value around three Gyr, with isochrone fitting yielding younger values than chromospheric activity. We also confirm the moderate super-solar metallicity of the star.

5. Spectrophotometric analysis

The new SPHERE photometry is complementary to the existing set of photometric data points on the spectral energy distribution (SED) of the companion obtained by Janson et al. (2011). In the following we use the more complete SED to refine the properties of GJ 758 B.

5.1. Fluxes and magnitudes

We retrieved the apparent fluxes that correspond to the SPHERE/IRDIS photometry of the companion by using the contrast ratio listed in Table 4 and by following a three-step process:

• We first built the 0.4–22.1 μm SED of the star from the Tycho B_T, V_T (Hoeg et al. 1997), USNO-B R, and I (Monet et al. 2003), and WISE W3-W4 photometry (Cutri et al. 2013). The 2MASS J, H, K_s, (Cutri et al. 2003) and W1-W2 photometry could not be used because of the saturation of the star (see Janson et al. 2011). The optical photometry was converted to apparent fluxes using the Gemini flux-conversion tool³. We considered the WISE zero points reported in Jarrett et al. (2011) for the infrared part.

• We adjusted a Gaia-COND model (Brott & Hauschildt 2005) with T_eff = 5400 K, log  g = 4.0 dex, and M / H = 0.0 to GJ 758 A fluxes values. This model has atmospheric parameters close to the ones determined from high-resolution spectra of the star (T_eff = 5435 K, log  g = 4.0, M / H = 0.12; Kovtyukh et al. 2004). The Gaia model reproduces the SED of GJ 758 (Fig. 3) well, including the 2MASS K_s band photometry, which appears to be less affected by the saturation.

• We derived the mean stellar flux into the SPHERE/IRDIS passbands using the flux-calibrated Gaia spectrum and the tabulated filter widths reported in Table 1.

The remaining fluxes of GJ 758 B were estimated directly from a flux-calibrated spectrum of Vega (Bohlin 2007), the Keck/NIRC2 and Gemini/NIRI magnitudes of the companions reported in Janson et al. (2011), and corresponding filter transmission curves. The effect of the telluric absorption on the final flux estimates for the companion was simulated using the ESO sky model calculator⁴ (Noll et al. 2012; Jones et al. 2013). We considered two altitudes of targets above the horizon (90 and 30°) to simulate dry and wet conditions. The effect is found to be negligible, compared to the error in the companion photometry. The final estimated fluxes of GJ 758 B, which we consider for the following analysis, are reported in Table 6. The fluxes in the overlapping narrow-band K1 and K_c filters are almost identical. This is an indication that our flux-conversion methods yield consistent results.

5.2. Comparison of GJ 758 B to empirical objects

The Y3/J2, J2/J3, H2/H3, and K1/K2 flux ratios provide a clear detection of water and methane absorptions around 1.15, 1.6, and 2.3 μm in the atmosphere of the brown-dwarf companion. We compared its 1–2.5 μm SED to those of 101 T0–T8 field dwarfs with near-infrared spectra taken from the SpeXPrism library (Burgasser 2014). The mean flux F_k,i and error σ_Fk,i associated with each template spectrum k and filter passband i was estimated and compared to the companion SED f and error σ_f using the G′′ goodness-of-fit indicator defined by Bowler et al. (2010): ${\mathit{G}}_{\mathit{k}}^{\mathrm{\prime \prime }}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{w}}_{\mathit{i}}\frac{\mathrm{(}{\mathit{f}}_{\mathit{i}}\mathrm{-}{\mathit{C}}_{\mathit{k}}^{\mathrm{\prime \prime }}{\mathit{F}}_{\mathit{k,i}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathit{f}}_{\mathit{i}}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{C}}_{\mathit{k}}^{\mathrm{\prime \prime }}{\mathit{\sigma }}_{{\mathit{F}}_{\mathit{k,i}}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(1)where ${\mathit{C}}_{\mathit{k}}^{\mathrm{\prime \prime }}$ is a renormalization factor applied to the template SED k, which minimizes ${\mathit{G}}_{\mathit{k}}^{\mathrm{\prime \prime }}$. w_i is the renormalized FWHM Δλ_i of each filter i following ${\mathit{w}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{\Delta }{\mathit{\lambda }}_{\mathit{i}}}{{\sum }_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{\Delta }{\mathit{\lambda }}_{\mathit{j}}}\mathrm{·}$(2)The indicator enables us to compare SEDs with an inhomogeneous wavelength sampling and with measurement errors on both the templates and the object. We rejected solutions which exceeded the upper limit of the flux into the CH4L passband (Janson et al. 2011). The G′′ indicator is minimized for the T6.5 dwarf 2MASS J22282889-4310262 (Burgasser et al. 2004), which is known to experience wavelength-dependent photometric variability (Buenzli et al. 2012). The comparison is shown in the upper panel of Fig. 4, and we report the G′′ values as a function of spectral type in Fig. 5. When flux-calibrated and scaled to the distance of GJ 758 B (using the parallax of Faherty et al. 2012), the spectrum of 2MASS J22282889-4310262 is over-luminous and a multiplication factor of 0.08 must be applied to fit the companion SED. This indicates that GJ 758 B is most likely later than T6.5. The variation of G′′ with the spectral type also clearly confirms that the companion is later than T5. This agrees with the conclusions of Janson et al. (2011).

In the lower panel of Fig. 4, we show the spectra of standard T8 and T9 dwarfs (Burgasser et al. 2004; Lucas et al. 2010) with measured trigonometric parallaxes and fluxes brought to the distance of the GJ 758 system. The companion SED is midway between the renormalized SED of the T8 and T9 standards. Nevertheless, the templates fail to reproduce the J3, H2, and K2 fluxes simultaneously. The companion also appears to have a luminosity intermediate between these two objects. Its J- and H-band absolute magnitudes agree well with the mean values reported in Dupuy & Kraus (2013) for T8-T8.5 objects.

The causes of the peculiar SED of GJ 758 B are unclear. The companion spectrophotometric properties could be related to a non-solar composition, or a surface gravity that is different to those of standard T8-T9 dwarfs. Both parameters produce opposite effects on 1–5 μm SEDs that are difficult to disentangle (e.g. Leggett et al. 2010). We used the spectra of wide companions to stars with a known age and metallicity to investigate the effect of peculiar atmospheric parameters, making the assumption that these objects share the same composition as their host star.

Ross 458 C (Goldman et al. 2010; Scholz 2010) appears as the only object with an estimated age (150–800 Myr) that is younger than the typical field dwarf ages (≫500 Myr), which has an estimated T_eff (625–755 K Burgasser et al. 2010; Burningham et al. 2011) and near-infrared spectral type (T8.5p) in the same range as that of GJ 758 B (Janson et al. 2011). It is also reported to have a super-solar metallicity (Fe/H = +0.2–0.3; Burgasser et al. 2010), e.g., similar to that of GJ 758 A (+0.2 dex, see Sect. 4). The spectra of both objects are also compared in the lower panel of Fig. 4. The spectrum of Ross 458 C from Burningham et al. (2011) represents the SED of GJ 758 B less well than the T8 standard. Its enhanced flux at K-band suggests that the two companions do not span the same surface gravity and/or metallicity interval.

We considered the opposite case of the peculiar T8 companion to the metal-poor ([Fe / H] = −0.38 ± 0.06 dex) G-type star BD+01 292 (Pinfield et al. 2012) and of the T8 companion to the sdM1.5+WD binary Wolf 1130 ([Fe / H] = −0.64 ± 0.17). Both companions have a suppressed flux at K-band, possibly due to the enhanced collision-induced absorption of H₂ encountered into clear/low-metallicity/higher-pressure atmospheres (Saumon et al. 1994; Borysow et al. 1997). They clearly produce a worse fit to the SED of GJ 758 B than the T8 standard does. In summary, we see an opposite trend for GJ 758 B’s departure from the SED of the standard T8.

The T8.5 companion to the old (3.5–6 Gyr) solar-metallicity star Wolf 940 (Burningham et al. 2009, [Fe / H] = −0.06 ± 0.20) represents the J-band flux better, at the price of a degradation in the fit in the Y band. We do not find a good fit with earlier type companions such as GJ 229 B (T7pec) or Gl 570 D (T7.5) (Geballe et al. 1996, 2001) nor primaries with roughly solar-metallicities (Neves et al. 2014).

We extended the comparison to additional peculiar dwarfs with red near-IR colors but no a priori knowledge of their age and metallicity (e.g. Mace et al. 2013, and references therein). We find that the red T8 dwarfs WISEP J161705.75+180714.0 and WISEP J231336.41-803701.4 (Burgasser et al. 2011) provide the best fit among all other aforementioned objects. They notably represent the Y-band flux well, compared to the other objects. Burgasser et al. (2011) note that the spectral properties of these two objects suggest cool (T_eff = 600 K), low surface gravity (log  g = 4.0), and cloudy atmospheres.

In summary, we cannot find an empirical object with known metallicity and distance that accurately represents all the near-IR narrowband and broadband fluxes of GJ 758 B simultaneously. We estimate a T8 spectral type from this comparison. The analysis is, however, certainly limited by the small amount of spectra of T8–T9 dwarfs with robust constraints on their age and metallicity.

5.3. Comparison of GJ 758 B to atmospheric models

We compared the SED of GJ 758 B to four sets of atmospheric models, i.e. BT-Settl, Exo-REM, Morley+12, and Saumon+12, in order to refine the estimate of log  g, T_eff, and Fe/H, and to understand its peculiar photometry. The models are described in Allard et al. (2013), Baudino et al. (2015), Morley et al. (2012), and Saumon et al. (2012) respectively. The specificities and parameter space of the models are described in more detail in Appendix A and Table A.1. We expect that the use of these different classes of models will allow the best possible approach for the accurate modeling of the atmospheric parameters.

To account for the inhomogeneous sampling of the real SED during the fitting process, we decided to use the goodness-of-fit G_k indicator defined by Cushing et al. (2008). This indicator contains a dilution factor, C_k, similar to the ${\mathit{C}}_{\mathrm{k}}^{\mathrm{\prime \prime }}$ factor defined in Eq. (1). C_k usually equals (R/d)², where d is the distance of the source and R its radius. Given the Hipparcos distance from GJ 758 A, we were able to retrieve the optimal average object radii for each given model.

Confidence levels cannot be derived directly from G. Therefore, we followed the approach of Cushing et al. (2008) to determine the most meaningful fitting solution for each model grid. For each photometric data point of the object, we generated a normal Monte-Carlo (MC) distribution of 10 000 draws with mean values of f_i and standard deviations of σ_i. The G_k values were computed for each of the resulting 10 000 SEDs. For each model of the grid, we computed the fraction of the 10 000 MC simulated SEDs that were best fitted by this given model. This f_MC indicator, ranging from 0 to 1, enabled us to test the significance of any fitting solution. The models with the highest f_MC value represents the most significant solution. However, we note that f_MC is sensitive to the sampling and extent of the model grid. Therefore, despite the criterion being useful for estimating the robustness of a given solution within a grid, it should not be used to evaluate the quality of the solutions found with different grids. We performed a visual inspection of the three solutions with the highest f_MC for each model grid, but only reported the atmospheric parameters and f_MC of the most probable solution in Table 8.

The Monte Carlo method works as long as the errors associated with f_i are uncorrelated. In the case of GJ 758 B, the errors associated with the flux-calibrated SED of the object combine uncorrelated errors. These correspond with companion contrast values that are associated with each filter to a correlated error, which arises for the flux-calibrated spectrum of the star. We accounted for both sources of errors in our MC simulations by multiplying the 10 000 MC SEDs of the companion by 10^{− 0.4 × N(μ = 0,σ_ph)} with N an additional MC normal distribution. From this distribution, 10 000 values were drawn with mean values of 0 and a standard deviation σ_ph equal to the magnitude error on the flux-scaling of the companion spectrum. We took σ_ph = 0.03 mag, which corresponds to the highest photometric error on the SED of the star.

The results of the fits are reported into Table 8 and shown in Fig. 6. The solutions with the highest f_MC always correspond to the solution with the minimum G. The corresponding dilution factors inferred from our MC simulations for the most probable fitting solution (highest f_MC) are shown in Fig. 7. No one model represents the whole SED well, especially the J3 and H2 fluxes. The Morley+12 models provide the best fits according to the G indicator. The three most significant solutions (>65% of the solutions) found with these models correspond to log g = 4.0–4.5 and T_eff = 550–600 K. The cloud-free Saumon+12 models only provides a better fit to the J-band flux. But their poorer representation of the other bands indicate that clouds are still needed in the photosphere to reproduce the GJ 758 B SED. Visually, the BT-SETTL14 models seem to provide a better fit to the H-band flux. The flux drop at M_s band in the BT-SETTL14 models is in better agreement with the upper limit found by Janson et al. (2011). New deeper observations at the M-band of the GJ 758 system could help to further discriminate the models. The BT-SETTL14 models do not provide any meaningful constraints on the log g. A reanalysis with a classical χ² confirms the conclusions. The radii (dilution factors) needed to adjust the surface flux, which were predicted by the models onto the apparent flux of the companion, are unphysical in the case of the BT-SETTL14 models. This may indicate that the T_eff of GJ 758 could be lower than the one corresponding to the best fit. The Exo-REM models fail to represent correctly the SED of the companion, especially in the Y-, H-, and K-bands.

We conclude that the companion has T_eff = 600 ± 100 K from the above analysis. This is in good agreement with Janson et al. (2011), whose analysis relied on the models of Burrows et al. (2006) but extended to colder temperatures (Hubeny & Burrows, in prep.). The Morley+12 parametric model points toward a low surface gravity, in agreement with the hints found in Sect. 5.2. But the non-existant exploration of the effect of the metal-enrichment in these grids of models, which are associated with model uncertainties, certainly biases the analysis. A low-resolution spectrum of the source is needed to determine log  g and metallicity (M/H) with good confidence.

From the T_eff and the derived age for the system (${\mathrm{3}}_{-2}^{\mathrm{+}\mathrm{3}}$ Gyr, see Sect. 4), we estimate a mass of $\mathrm{2}{\mathrm{3}}_{-13}^{\mathrm{+}\mathrm{17}}{\mathit{M}}_{\mathrm{Jup}}$ for GJ 758 B using the BT-SETTL13 grid of models (Allard et al. 2013). This value is in the low range of the masses inferred by Janson et al. (2011), as a direct consequence of our estimated age range for the system that is slightly younger than the one they considered.

Finally, we report in Table 9 the predictions from the Baraffe et al. (2003) and Saumon & Marley (2008) models. The models predict radii corresponding to the estimate T_eff in the range 0.80–1.21 R_Jup, which are marginally consistent with the radii derived from the SED fit with the Exo-REM, Saumon12+, and Morley12+ models. They are still 25 and 11% larger than those infered from the SED fit with the BT-SETTL14-Y dwarfs and BT-SETTL14 models, respectively (see Fig. 7). The difference may arise from the different boundary conditions considered for the evolutionary models and the atmospheric models used for the SED fit. Nevertheless, it is more likely that the T_eff derived from the SED fit is slightly overestimated by the SETTL models (by 100–200 K) which leads to this inconsistency.

5.4. Nature of the new candidate companion

The detection of a new candidate companion around a star with an already known companion is particularly interesting. It is going to become very common with the new generation of high-contrast imagers because of the boost in sensitivity that they provide at smaller angular separations. For the new candidate detected in our IRDIS data, we make use of the large multi-wavelength coverage (Y- to K-band) to perform a photometric analysis.

Although the error bars on the photometry of the new candidate are large (>0.5 mag, see Table 7), we attempt a first-order estimation of its spectral type by comparing its observed flux with the SEDs of stellar and substellar objects. The star data are taken from the IRTF stellar library⁵ (Cushing et al. 2005; Rayner et al. 2009), while the brown dwarf data are taken from the NIRSPEC brown dwarf spectroscopic survey (McLean et al. 2003) and from Leggett et al. (2001). For the comparison, we use the ${\mathit{G}}_{\mathit{k}}^{\mathrm{\prime \prime }}$ indicator as in Sect. 5.2. The results are presented in Fig. 8, where we show the ${\mathit{G}}_{\mathit{k}}^{\mathrm{\prime \prime }}$ values as a function of spectral types (inset), and the SEDs of four objects that can equally well fit the photometry of the new candidate within the error bars.

The ${\mathit{G}}_{\mathit{k}}^{\mathrm{\prime \prime }}$ distribution shows a rather flat minimum in the L3–T1 range, indicating that our candidate could likely be of substellar nature. However, reaching a final conclusion is difficult from our current data because of the significant uncertainties on the photometry. As shown in Fig. 8, the candidate photometry is compatible with mid-L to early T-types, but late-M and early-L (not shown) would also provide decent fit. We note that the low galactic latitude of GJ 758 (+8 deg) significantly increases the probability of background contamination, particularly with late M stars, which are the main source of contamination at high-contrast (e.g., Chauvin et al. 2015).

Other possibilities for the nature of the candidate could include solar system bodies, such as asteroid and transneptunian objects, or extra-galactic objects. However, an asteroid basically reflects the near-IR light from the Sun, resulting in a very flat G2V spectrum that is not compatible with the photometry. In addition, these objects would be characterized by a very large proper motion of several mas to several dozens of mas per second. Our observations, taken over two consecutive nights, completely rule out this possibility. On the other hand, extragalactic sources such as galaxies are another possibility, but they would be resolved by the very fine plate scale of IRDIS (~12.25 mas, see Table 3), even at significant redshifts. The point-like structure of the candidate also rules out this possibility.

In conclusion, we cannot rule out the possibility that the new companion is indeed bound to GJ 758 since its photometry is broadly compatible with L-type objects. A second epoch will be required to clear any possible doubt. Given the high proper motion of the star (~180 mas/yr), a confirmation of the status of this candidate is already possible.

6. Astrometry and orbital properties

6.1. Least-Square Monte Carlo orbital fitting

We used the new IRDIS astrometric measurement to put constraints on the orbital solution of the system. In previous studies by Thalmann et al. (2009) and Janson et al. (2011), it was already shown that the system presents significant orbital motion and Monte Carlo simulations were used to get a first estimate of the orbital elements. In this study we first used a Least-Square Monte-Carlo (LSMC) approach to study the parameter space of possible orbits. For this purpose we created 5 × 10⁶ sets of orbital elements, which were drawn from uniform distributions. These sets of orbital elements were then used as starting points for a least-squares minimization routine. The method is described in detail in Ginski et al. (2013). To limit the parameter space we fixed the total mass of the system to the nominal value of ~1 M_⊙: 0.97 M_⊙ for the star (Takeda et al. 2007) and ~0.03 M_⊙ for the companion at the probable age of the system. In addition, we limited the semi-major axis to values smaller than 63.45′′ (1000 AU at a distance of 15.76 pc). This is assuming that the system is stable long-term against disruption in the galactic disk as described in Close et al. (2003). Given the high age of the system (${\mathrm{3}}_{-2}^{\mathrm{+}\mathrm{3}}$ Gyr, see Sect. 4) and the fact that we still find the companion close to the host star, this assumption seems reasonable.

The results of our simulations are shown in Fig. 9. We do not show the results for the longitude of the ascending node and the argument of the periastron, since they are not well constrained yet by the available astrometry. In Fig. 10 we show the best fitting orbit solution that was recovered by the LSMC orbit fit. The corresponding orbital elements are shown in Table 10, alongside the results recovered from our Markov chain Monte Carlo (MCMC) simulation, which we discuss in the following section.

Since the orbit does not show significant curvature yet, we cannot put an upper limit on the semi-major axis or the eccentricity. However, we find a lower limit of 0.14 for the eccentricity and 21.9 AU (1.39′′) for the semi-major axis. In general, the semi-major axis of possible orbits scales with the eccentricity, as can be seen in Fig. 9a. The minimum values of the semi-major axis and the eccentricity, as well as the general behavior of the well fitting orbit solutions, is consistent with the results presented in Janson et al. (2011), which were derived from simple Monte Carlo simulations.

In Fig. 9b we show the inclination of possible orbital solutions as a function of eccentricity. For close to face-on orbits (inclination close to 0 deg) we can constrain the eccentricity of the orbit to values between 0.47 and 0.55. This range grows continuously larger with increasing orbit inclination. For an inclination of ~50 deg, the full range of recovered eccentricities gives results that are consistent with the astrometric measurements. We can put an upper limit on the inclination of 70.8 deg, i.e., we can exclude edge-on orbit solutions. If we compute a simple median of the recovered orbit inclinations, we get a value of 58.9 ± 18.8 deg. This is, within the given uncertainties, consistent with the interval found in Janson et al. (2011). Inclinations smaller than 40 deg correspond to small semi-major axes, with an upper limit of 77.5 AU (4.92′′) while, for larger values of the inclination, orbit solutions with the full range of recovered semi-major axes are possible.

Finally, in Fig. 9c we show the times of the periastron passage that we recovered from our simulations. The vast majority (87.5%) of our solutions pass the periastron between the years 2000 and 2065. The solutions that show the periastron passage at the time of the observations are generally highly eccentric and have large semi-major axes, which would explain that no curvature of the orbit has yet been observed.

While these solutions fit the orbit very well geometrically, they are, however, very unlikely, given that the companion would spend the vast majority of time at much larger separations from the primary star than where it was discovered. Indeed, if we use the orbital period of roughly 26 000 yr of the best-fitting LSMC orbit, we can estimate that the probability of finding the companion within 30 yr of the periastron passage is only approximately 0.1%. However, the orbits that pass the periastron within the next few decades could have lower eccentricities and semi-major axes. For an eccentricity of around ~0.5 there is a strong peak for the time of the periastron passage in the year 2040. It is thus of great interest to continue an astrometric monitoring of this system, since significant acceleration (i.e., curvature of the orbit) is to be expected, especially for cases with non-extreme eccentricities.

Since it will be very interesting if the system does indeed exhibit a high eccentricity (i.e., for the plausibility of scattering scenarios during its early formation), we examine how reliable the eccentricities of our recovered orbit solutions are. Pearce et al. (2014) study the possibility that an unknown inner (sub)stellar companion could introduce a false-positive eccentricity signal in the relative astrometry between the primary star and the known, directly imaged, companion. This is due to the astrometric displacement of the primary star as it orbits around the common center of mass with the hypothetical inner companion. We use their formalism to calculate the mass and angular distance that would be required for such an inner companion to make the orbit solution for GJ 758 B appear eccentric when, in fact, the real orbit is circular. We do this for all the orbit solutions that fit the astrometric measurements and the results are shown in Fig. 11. We find that a hypothetical inner companion would need a mass between 0.02 M_⊙ and 0.14 M_⊙ with an orbit separation of 0.42′′ (6.6 AU; depending on the mass of the system and the epoch difference of the astrometric observations). For the large majority of our solutions, we can reject such a companion because it would have been discovered in our deep IRDIS images (see Sect. 7 for detection limit estimations). However, because of the old age of the system, we would not have been able to recover inner companions with masses below ~0.05 M_⊙ at the required angular separation.

To exclude the remaining possible solutions, we retrieved archival radial velocity data of GJ 758 A obtained with the ELODIE high-precision fiber-fed echelle spectrograph (Baranne et al. 1979), which coveris a time baseline of 7.8 yr, as well as archival data from the Lick Planet Search program, which covers 13.2 yr (Fischer et al. 2014). These combined data, shown in Fig. 12, covers a total of 16 yr. The data can be used to reject any hypothetical companion on a 17-yr period more massive than 0.02 M_⊙ that has any inclination greater than five degrees. Given the spherical symmetry of the system, this translates into a rejection of 98% of the orbital solutions for a hypothetical 0.02 M_⊙ inner companion and an increasing rejection rate for higher masses. It is thus extremely unlikely that the observed eccentricity is due to an inner companion causing an astrometric signal.

6.2. Markov chain Monte Carlo orbital fitting

The use of the MCMC technique to fit orbits of companions, either detected by radial velocity or by direct imaging, has become very popular in recent years. In relation to imaged planetary or substellar companions, it was for instance successfully applied to β Pictoris b (Chauvin et al. 2012; Nielsen et al. 2014; Macintosh et al. 2014), Fomalhaut b (Kalas et al. 2013; Beust et al. 2014), and to the four-planet system of HR 8799 (Pueyo et al. 2015).

MCMC is particularly well-suited for imaged companions for which the observational follow-up usually covers only a small part of the whole orbit (because of large orbital periods). To fit GJ 758 B’s orbit, we first used the code already used to fit β Pictoris b’s (Chauvin et al. 2012) and Fomalhaut b’s (Beust et al. 2014) orbits. But, given the number of solutions at very large eccentricities that were hard to reach (the best-fit solution has eccentricity ~0.93; see Table 10), we moved to the use of another code that we have developed recently, which is based on the use of universal Keplerian variables with the Metropolis-Hastings algorithm and using Gibbs sampling as a convergence test. Universal variable formulation (Danby & Burkardt 1983; Burkardt & Danby 1983; Danby 1987) is an elegant way to provide a unique and continuous description of the Keplerian motion and is valid for any kind of orbit, whether bound or unbound. The details of this code will be presented in Beust et al. (in prep.). This code can handle both bound and unbound orbits, and is therefore not limited to elliptic orbits. It is thus well suited for very eccentric orbits.

Finding very eccentric orbital solutions should indeed not be surprising. As was shown by Pearce et al. (2015), whenever astrometric orbits are followed over small orbital arcs, more or less arbitrarily eccentric solutions can be found depending on the unknown values of the z-coordinate and velocity along the line of sight. By and large, this situation applies here. At least from a mathemetical point of view, unbound solutions should be valid as well. This motivated us to use the universal variable code.

Ten chains were run in parallel until the Gelman-Rubin parameters and (Ford 2006) repeatedly reached convergence criteria for all parameters, i.e., and . This occurred after 5.2 × 10⁸ steps. At this point, a sample of 10⁶ orbital solutions is taken from the chains as representative for the posterior distribution of orbits. The orbital parameters considered are the periastron q, the eccentricity e, the inclination i with respect to the sky plane, the longitude of ascending node Ω (counted from north), the argument of periastron ω, and the time for periastron passage t_p. Here, we consider the periastron instead of the semi-major axis, as the periastron assumes a continuous distribution from elliptical to hyperbolic orbits. The priors on those elements are assumed as being uniform for Ω, ω, e, and t_p, logarithmic for q and ∝sini for i. Combined with a uniform prior for Ω, the latter choice ensures a uniform probability distribution over the sphere for the direction of the orbital angular momentum vector. We emphasize that this choice of prior, especially concerning the eccentricity, is not dictated by physical likelihood considerations, but rather by mathematical constraints on the sole basis on the available astrometric data. While a linear eccentricity prior between 0 and 1 can be realistic, clearly unbound orbits appear unprobable, given the age of the star. The probability of witnessing an ejection or a flyby right now is indeed very low. But, as very eccentric solutions appeared to be compatible with the astrometric data, we wanted to allow the MCMC code to explore the unbound regime to estimate the actual constraints on the data and to avoid the introduction of artificial cut-offs.

The resulting posterior distribution is shown in Fig. 13, where probability histograms for individual elements are displayed as well as density maps for all possible pairs of parameters. The red bars that appear on the histogram plots, as well as the black stars in the bidimensional maps, correspond to the best χ² solution that was derived using a least-square Levenberg-Marquardt fitting scheme before launching MCMC. This solution has a reduced χ² = 0.419, but more than 80% of the solutions in our posterior sample have reduced χ² < 1.5. Peak values, confidence intervals, as well as details about the best χ² solution, are given in Table 10. Figure 13 shows that the eccentricity distributions extend beyond e = 1, so that we have both bound and unbound solutions in our sample. The upper limit at e = 2 in the eccentricity distribution is not physical. This threshold was fixed at the beginning of the simulation to save computing time.

The plots involving Ω and ω appear twofold, with similar patterns saparated by ± 180°. This is a direct consequence of the degeneracy of the projected astrometric motion (Beust et al. 2014). To each solution with (Ω, ω) values, we find a corresponding twin solution with the same other orbital elements, but with (Ω + π,ω + π). Both generate the same projected orbital motion.

Our initial comment on the result is that the orbit is clearly eccentric. However, despite the presence of unbound solutions in our sample, and although the best χ² solution appears very close to e = 1, most solutions have moderate eccentricities ≲0.7, with a peak around e = 0.5. In our sample, 68% of orbits are bound. This is enough to confirm that GJ 758 B is very probably a bound companion to GJ 758, as an unbound orbit would mean an ongoing flyby or a very recent ejection. As both configurations can be regarded as improbable (though not impossible), finding more than two-thirds of bound solutions in our sample is a very strong indication of a bound orbit. If 68% can be taken as the minimum probability to have a bound orbit, then the actual probability is, in fact, much higher.

As noted above, this value is very probably far below the actual probability, since an unbound orbit would mean an ongoing flyby or a very recent ejection. This is a very improbable configuration given the age of the star. We regard 68% as the minimum probability to have a bound orbit without any physical consideration about the likelihood of unbound configurations. It is thus sufficiently high to allow us to claim that GJ 758 B is actually a bound companion to GJ 758 A. Based on the ratio between the timescale of an ejection event and the age of the star, the actual probability o observing one today should not exceed ~10^-6.

The periastron lies in the range 10–40 AU for about 70% of solutions, so that this must be regarded as the most probable range, with a clear probability peak at q = 20 AU. Indeed, solutions with higher q values correspond mostly to unbound solutions, and must therefore be considered as less probable.

All orbital solutions have inclination i well below 90°, compatible with a prograde motion as seen from the Earth. The inclination peaks around 60°, while the longitude of ascending node Ω exhibits a clear peak at ~−40°. This shows that the orbital plane of GJ 758 B is rather well constrained. Conversely, the argument of periastron ω is very badly constrained. This is a direct consequence of the eccentricity dispersion of the solutions, as can be seen from the (ω,e) density maps (Fig. 13). The periastron passage is, however, better constrained with the next occurrence of a peak in 2039.

6.3. Conclusion on astrometry

To better compare LSMC and MCMC results, we created matching bidimensional density maps, which are restricted to bound orbits, and which consider the semi-major axis a and the eccentricity e, as well as the inclination i. For the MCMC results we restricted ourselves to bound solutions only. To match the MCMC results closely, we cut off the LSMC results at semi-major axes smaller than 10′′. The results are shown in Fig. 14. We first show posterior MCMC distribution for the semi-major axis and the orbital period (left plots). Although both histograms exhibit tails towards large values (and thus approaching unbound orbits), we see that clear peaks appears around ~40 AU and ~200 yr. These values must be taken as the most probable ones. Then we compare the 2D maps that were generated by MCMC and LSMC (middle and right plots). An immediate comparison shows that both approaches agree very well and derive the same well-fitting range of bound orbital solutions. The differences in the density of solutions are likely caused by the difference of prior distributions that were used as input for both methods. While we considered uniform distributions in a and i for LSMC, we used distributions that were uniform in log  q and ∝sini for i. We did this because the aim of the methods is somewhat different. With LSMC, we aim to find the full possible range of geometrically well-fitting orbits, as well as the best-fitting orbit in terms of ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$. With MCMC, on the other hand, the goal is to find the correct posterior probability distribution, given our prior knowledge of the system. This includes knowing that shorter period orbits are more likely given where we find the companion in our observation epochs, as well as the statistical likelihood of orbit inclinations.

We emphasise that the general results obtained by both methods agree very well. In particular, with LSMC we find an upper limit of the inclination of bound orbits of 70.8 deg, while the upper limit of the 1σ confidence interval recovered by LSMC is 70.3 deg. Both methods also find strong peaks in the time of the periastron passage between the years 2039 and 2040. Furthermore, the MCMC fit recovered 68% of bound orbit solutions. Given this high likelihood of a bound orbit and the high chance of a periastron passage within the next few decades, the GJ 758 system remains an interesting target for continued orbital monitoring. Orbit curvature will likely be discovered in this timeframe, allowing for a much better constrained determination of the companion’s orbit.

7. Sensitivity to additional companions

To conclude our analysis, it is interesting to look at our sensitivity to additional massive companions in the GJ 758 system. We calculated detection limits in the different DBI filters following an ADI analysis with KLIP. The limits were estimated by measuring the standard deviation of the residuals in annuli of width 1 λ/D at increasing angular separation, divided by the maximum of the off-axis PSF of the star in the same filter. To properly take the effect of self-subtraction that was induced on the detection limits into account, the algorithm throughput was estimated by injecting fake companions, regularly spaced from 0.1′′ to 2.0′′, into the pre-processed data cubes. They were injected at a level ten times higher than the noise residuals in the final images at the separation of each of the companions. This process was repeated ten times with different orientations of the pattern of fake companions to average out possible variations of the throughput as a function of the position in the field. The throughput at each separation was then calculated to be the mean throughput over the ten measurements.

Taking into account the throughput in the different filters, the final detection limits are plotted in the right-hand panel of Fig. 15. Because of the small amount of field rotation in all observing sequences (≲8^o), the throughput of the analysis decreases significantly towards small angular separations, resulting in a sharp deterioration of the detection limits. It is only at separations larger than 0.21′′, 0.25′′, 0.33′′, and 0.43′′ that the field rotates by more than λ/D over the course of the complete sequence in Y-, J-, H-, and K-band respectively. The sensitivity is, nonetheless, improved compared with previous measurements below ~0.5′′.

The right-hand panel of Fig. 15 shows the conversion of these detection limits into physical units of projected separation and physical mass, using the known distance of the star (15.76 pc) and the AMES-COND evolutionary models (Baraffe et al. 2003) that were calculated in the IRDIS DBI filters. The two sets of curves represent the limit for the two extremes of the system’s age range, 1–6 Gyr (see Sect. 4). For the youngest part of the age range, our observations clearly probe the low-mass brown-dwarf regime down to 4–5 AU, and even the planetary-mass regime beyond 20–30 AU. However, if we assume an older age for the system, i.e., more in line with the different age indicators, only massive brown dwarfs could be detected.

8. Summary and conclusions

Our new study of GJ 758 offers an improved overview of this interesting system. The brown-dwarf companion is redetected and we confirm some of its already known properties using our finer spectral sampling from SPHERE/IRDIS observations. In particular, we recover a low T_eff = 600 ± 100 K from a comparison of four different sets of models, which are in good agreement with previous studies (Thalmann et al. 2009; Currie et al. 2010; Janson et al. 2011). There are, however, some interesting peculiarities that are worth discussing.

In comparison to empirical objects, GJ 758 B appears as a very interesting object because we cannot find objects with known ages and metallicity to match all of its observed fluxes. We estimate a T8 spectral type for this object, but this estimation is limited by the small amount of spectra of T8–T9 dwarfs with robust constraints on their ages and metallicity. One of the most likely explanations for the peculiar SED of this brown dwarf is the super-solar metallicity of the primary ([Fe / H] = 0.18 ± 0.05; see Sect. 4.2) and of the companion, if we assume they share the same metal enrichment. This is supported by the fact that the T8pec companion to the metal-poor dwarf BD+01 292 (Pinfield et al. 2012) shows a K-band spectral deviation opposite to that of GJ 758 B. Unfortunately, the current lack of such companions precludes us from performing a meaningful comparison. Similarly, a comparison with synthetic grids of models is constrained by the limited extension of most grids toward non-solar metallicities and low T_eff values. None of the four types of models that we tested was able to accurately reproduce the photometry of GJ 758 B in all filters. The J3 and H2 fluxes are especially difficult to reproduce and cannot be readily explained with any of the models. The BT-SETTL14 models (Allard et al. 2013), with an enrichment of 0.3 dex in α elements with respect to solar composition, do provide a better fit of the J3 flux, but at the expense of much smaller dilution factors, which correspond to unphysical values of the companion radius (~0.7 R_Jup). Even so, the high H2 flux is not reproduced. This could actually indicate that the T_eff of GJ 758 B might be even lower than the one of the best fit but we are limited by the absence of metal-enriched models at very low T_eff. As a result, our analysis confirms the low T_eff of the object, corresponding to a mass of $\mathrm{2}{\mathrm{3}}_{-13}^{\mathrm{+}\mathrm{17}}$M_Jup in the considered age range (using the Baraffe et al. 2003 evolutionary tracks), but we cannot infer any precise value for log g and M/H. The study of this object would strongly benefit from low-resolution spectroscopy, e.g., with IRDIS long-slit spectroscopy mode (Vigan et al. 2008).

The new astrometry confirms the picture of a very eccentric companion. As reported by Janson et al. (2011), the curvature of the orbit is still not detectable yet. The ~0.28′′ motion of the companion along its orbit roughly appears as a straight line from the previous measurements taken in 2010. Our LSMC and MCMC simulations favour an eccentric – but bound – orbit, with a high likelihood of e ≃ 0.5. In particular, no orbital solution shows an eccentricity lower than 0.14, which is consistent with previously reported results and is between our two approaches ([0.133–1.78] 95% confidence interval from MCMC). In addition, we have ruled out the possibility that the observed eccentricity is just caused by a massive closer-in companion that would create a false positive eccentricity by astrometric displacement of the primary. Indeed, although not extremely accurate, our radial velocity (RV) data reject a companion that is more massive then 0.02 M_⊙ on periods shorter than 17 yr and an inclination that is larger than 5 deg. Finally, our new IRDIS observations reject the possibility of an additional companion that is more massive than ~30–40 M_Jup (for ages of 1–3 Gyr) above 4–5 AU.

In the light of our constraints on the orbit and mass of GJ 758 B, it is interesting to look into the formation of this object. While our study does not bring enough new material in favor or against planet-like formation scenarios (core accretion vs. gravitational instability, migration, etc.), we can instead focus on stellar-like formation scenarios. Several past studies have argued for a universal companion mass function (CMF) for stellar and substellar companions. In particular, in their in-depth anaysis focused on brown dwarfs around solar-type stars, Metchev & Hillenbrand (2009) find tentative evidence for such a universal CMF and predict a peak in semi-major axes for brown dwarfs at ~30 AU. From this point of view, it makes sense to compare the properties of the GJ 758 system to the properties of solar-type multiple systems. (Raghavan et al. 2010) have published the most complete study on this topic to date. While we cannot compare quantitatively our results on a single object with their global multiplicity analysis, it is interesting to note a few qualitative facts. Firstly, with a most likely period of 312.48 yr (from MCMC, see Table 10), GJ 758 B falls exactly at the peak of their period distribution (293.57 yr; Raghavan et al. 2010, Fig. 13). Secondly, with a most likely eccentricity of 0.525, GJ 758 B falls in the bulk of their eccentricity-period distributions. And thirdly, with a ratio q = M₂/M₁ = 0.023 ± 0.013, the GJ 758 systems falls at the very edge of the study of Raghavan et al. (2010), but in their q-period plot, the system is not completely isolated and could be part of the tail of the distribution. While these facts do not prove that the system was formed in a stellar way, they certainly support this possibility qualitatively.

In conclusion, GJ 758 B remains a very interesting object that warrants i) deeper observations to look for additional companions in the system; ii) spectroscopic observations to better constrain its physical properties, and iii) astrometric monitoring to get tighter constraints on its eccentric orbit.

---

Acknowledgments

A.V., M.B., G.C., G.S., J.-L.B., and D.M. acknowledge support in France from the French National Research Agency (ANR) through project grant ANR10-BLANC0504-01, the CNRS-D2P PICS grant, and the Programmes Nationaux de Planétologie et de Physique Stellaire (PNP & PNPS). J.-L.Ba.’s Ph.D is funded by the LabEx Exploration Spatiale des Environnements Planétaires (ESEP) #2011-LABX-030. V.D. is partially supported by the Australian Research Council. V.D., S.D., A.L.M., R.G., and D.M. acknowledge support from the Progetti Premiali funding scheme of the Italian Ministry of Education, University, and Research. E.B. and J.H. are supported by the Swiss National Science Foundation (SNSF). A.Z. acknowledges support from the Millennium Science Initiative (Chilean Ministry of Economy), through grant Nucleus RC130007. The authors warmly thank A. Bellini and J. Anderson for kindly providing the catalog positions of the stars in the 47 Tuc field before their publication. We are very grateful to D. Kirkpatrick, A. Burgasser, D. Pinfield, B. Burningham, and G. Mace for sending us the spectra of benchmark T–Y objects. We thank C. Morley, D. Saumon, and F. Allard for free online access to their atmospheric model grids, and F. Bouchy for his help with the ELODIE data. This research has benefitted from the SpeX Prism Spectral Libraries, maintained by Adam Burgasser at http://pono.ucsd.edu/~adam/browndwarfs/spexprism. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland), and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012), and grant number 312430 for FP7 (2013–2016).

References

Appendix A: Description of the atmospheric models

The specificity of the models in the range of T_eff suitable for GJ 758 B have not yet all been described in the literature. So it is important to make a description of the most relevant hypothesis in the models and differences between the models in this paper. The parameter space of the models is summarized in Table A.1.

The BT-Settl model couples a cloud model to a 1D radiative transfer code PHOENIX (Allard et al. 1994; Hauschildt et al. 1997). The model considers the formation of a cloud deck, which is composed of up to 55 grain species. The grain size and density, the abundances of chemicals in the gas phase, including the effect of element depletion induced by the grain formation, is computed layer per layer through the photosphere, following a comparison of the timescales for nucleation, condensation, gravitational settling or sedimentation, and mixing. Once rained out below the photosphere, the grain opacities are not accounted for in the radiative transfer. Nevertheless, these grains can still interact chemically with the gas phase. These models can predict the flux at the surface of a given object that is only defined by log  g, T_eff, and [Fe/H]. They account for the non-equilibrium chemistry of CO, CH₄, N₂, NH₃, and CO₂. The models predict the formation of a secondary (resurgent) cloud layer into the photosphere made up of Na₂S and MnS then of KCl, NaCl, and some ZnS that lies above the rained-out primary cloud layer that is located below the photosphere and that originally remains in the atmosphere of L and early-T dwarfs. Here we used the 2014 releases of the models (hereafter BT-SETTL14), which include revised alkali cloud opacities and the latest CIA opacities (Abel et al. 2012). A specific grid was computed for the project to cover the Y-dwarf temperature domain (hereafter BT-SETTL14-Y), in addition to the already existing grid that covers a broader interval of T_eff and that considers α-element enrichment (hereafter BT-SETTL14).

The 1D Exo-REM models (Baudino et al. 2015) propose a simplified approach of substellar atmospheres. They predict the

equilibrium-temperature profile and mixing-ratio profiles of the most important gases (H2-He collision-induced absorption, H₂O, CO, CH₄, NH₃, VO, TiO, Na, and K). The absorption by iron and silicate cloud particles is added above the expected condensation levels with a fixed-scale height and a given optical depth at a given reference wavelength. For the purpose of the GJ 758 B study, two grids of models – NC and T3 – were computed. The NC models consider photospheres with no cloud opacity. The T3 models consider a photosphere with 30 μm Iron (Fe) and Forsterite (Mg₂SiO₄) grains and an optical depth of reference τ_cloud = 3. The grains are located between the condensation level and a 100 times lower pressure. They have scale heights equal to the gas-scale height, and optical depths of 3 and 0.45 at 1.2 μm, respectively, for Fe and Mg₂SiO₄.

Similarly to the BT-Settl models, the Morley+12 models account for the possible resurgence of clouds in late-T dwarf atmospheres. These 1D models build on the cloud model of Ackerman & Marley (2001). The cloud content (and opacity) is determined by a balance between the upward transport by turbulent mixing with the sedimentation. The Ackerman & Marley (2001) models do not compute the microphysics associated in the clouds, but instead leave as free parameters the vertical eddy diffusion coefficient K_zz and a sedimentation efficiency parameter f_sed. A higher f_sed corresponds to thinner (rained-out) clouds. In the case of the Morley+12 models, non-equilibrium chemistry is not included, so that only models with K_zz = 0 are available and they only enable an exploration of the f_sed. As a second difference to the BT-Settl, the Morley+12 models do not account for chemical reactions between the condensed species and the gas phase. Finally, we added the models of Saumon et al. (2012) to this grid, which are also based on Ackerman & Marley (2001) models, to explore the case of an extreme sedimentation efficiency (cloud-free atmospheres).

All Tables

All Figures