© F. Kiefer et al. 2019

1 Introduction

According to the classical convention, brown dwarfs (BD) are substellar objects whose masses are too small to maintain hydrostatic equilibrium thanks to hydrogen-based nuclear reactions, while massive enough to ignite Deuterium nuclear reactions in the core, at least for a few million years. Following this definition, the BD domain is framed within the mass range of 13–80 M_J. These boundaries may vary according to intrinsic stellar properties, such as metallicity (Chabrier & Baraffe 1997; Spiegel et al. 2011). Defining these limits is the subject of much debate (in e.g. Saumon et al. 1996; Chabrier & Baraffe 2000; Luhman et al. 2007; Luhman 2012; Chabrier et al. 2014) from which it has been proposed that the existence of nuclear reactions in the core of a substellar body is not the crucial parameter to define its nature.

The observation of objects with masses as low as 5 M_J in young stellar clusters is strong evidence that molecular cloud fragmentation is not limited in mass and can form objects in the BD and giant planet mass regime (see e.g. De Marchi et al. 2010). This is well reproduced by star formation simulations (Chabrier 2003; Luhman 2012; Lee & Hennebelle 2018). Moreover, the mass distribution of widely separated binaries extends well within the BD domain (see e.g. Burgasser et al. 2007 and reference therein). On the other hand, the scarcity of detections of BD companions with orbital periods shorter than 10 yr, the so-called brown dwarf desert (Halbwachs et al. 2000; Grether & Lineweaver 2006), is followed by an increase of detection frequency at masses lower than 10 M_J (Marcy & Butler 2000; Udry et al. 2002). This shows that giant planets and substellar objects that were formed like stars overlap on afew tens of Jupiter masses.

Planet-formation pathways such as disc instability or core accretion could in principle allow the formation of bodies up to 40 M_J within protoplanetary discs (Pollack et al. 1996; Boss 1997; Ida & Lin 2004; Alibert et al. 2005; Mordasini et al. 2009). Knowing the extent of the tail of the distribution of giant planets within the BD domain could therefore help to constrain the planet formation models. This tail is yet undetermined because the statistics of detections of substellar companions in the 5–40 M_J are still poor, though the observational efforts made in recent years have led to abundant detections of BD companions with diverse instrumental methods (Sozzetti & Desidera 2010; Sahlmann et al. 2011; Díaz et al. 2012; Ranc et al. 2015; Wilson et al. 2016).

A difficulty arises due to the scarcity of BD companions detected at short orbital periods. For some reason, substellar companions do not exist at close distance from a more massive primary star. This implies that the mass–period distribution of BD companions to sun-like stars is affected by several possible perturbing effects, such as tidal interactions, magnetic braking, and tidal dissipation (Guillot et al. 2014). This strongly biases the determination of the real mass distribution of giant planets and very low-mass stars.

It is therefore necessary to constrain the minimum orbital period above which this effect becomes negligible. Mixing the results from several surveys performed with diverse detection methods, Ma & Ge (2014) proposed a restricted BD desert enclosed within P < 100 d and 30 < M < 60 M_J, with a mass separation between star-like and planet-like BDs at 43 M_J. More-recent microlensing detections (Ranc et al. 2015) and the results of RV and astrometry (Wilson et al. 2016) added to already published detections tend to confirm the framing of the desert at periods lower than 100 days. The use of detections arising from several diversely biased or incomplete surveys is perilous however. To our knowledge, there exists no fully complete non-biased statistical sample of detected BD companions.

It would be most valuable to achieve a survey of BD companions that is non-biased, or at least for which the selection function of the followed-up sample is well known and allows meaningful statistics to be derived for the BD population. Some of the most problematic issues with gathering detections from multiple surveys, apart from instrumental bias, stem from the diverse interests of the observers. Observations are usually stopped as soon as the followed-up target is no longer of interest regarding the given study. Typically, on one hand, sources with a companion that is not within the planetary mass domain, beyond about 20 M_J, are not followed up and not always published. On the other, orbits and mass ratio of obvious stellar binaries are easily characterised and published. It follows that BDs within the BD desert and especially at periods larger than 1 yr are undersampled.

The volume-limited FGK star survey program searching for giant planets with the SOPHIE spectrograph installed at the Observatoire de Haute-Provence (Bouchy et al. 2009; Hébrard et al. 2016) offers a well-constrained framework for characterising the statistics of BD companions around solar-like stars. The target sample includes about 2350 sources among all 2950 known FGK stars of the northern sky (δ > +00 : 00 : 00) in the neighbourhood of the Sun below 60 pc, and in the main sequence (± 2 mag), with +0.35 < B–V < +1 (Dalal et al., in prep.). To this date, around 2050 sources have been observed at least three epochs each, with a target signal-to-noise ratio (S/N) per spectrum of at least 50.

The reflex motion due to BDs within the desert leads to RV amplitudes larger than 100 m s⁻¹. Since the SOPHIE spectrograph is able to detect RV signals as low as a few m s⁻¹ (Courcol et al. 2015) on a time baseline of 13 yr, BD companions can easily be detected around nearby bright stars. Therefore, we expect to eventually reach almost 100% completion of RV-detected BD candidates with orbital periods of less than 10 000 days around FGK stars in this volume-limited sample, whose visual magnitude is brighter than 11.

In the continuation of the work of Díaz et al. (2012) and Wilson et al. (2016) which published several new objects in the BD desert with SOPHIE, we present here the latest results of this radial velocity (RV) survey on 54 solar-like sources with spectral types ranging from K5 to F5. With RV only, we report 12 BD candidates with M sini within 15–90 M_J, among which 8 have not yet been published.

We conservatively extend the BD domain above 80 M_J in order to include objects in the grey zone of 80–90 M_J, separating M-dwarfs from BDs. We believe this to be justified for three main reasons. Firstly, there is always an uncertainty (up to a few M_J) on the M sini derived with RV. Secondly, the mass limit for hydrogen burning is not a strict one, and may vary according to, for example, metallicity from 83 to 75 M_J within M/H ~ [−1;0] (Chabrier & Baraffe 1997). Finally, extending towards low-mass M-dwarfs allows the tail of the BD mass distribution on the stellar side to be explored.

Although velocimetry is an efficient means to detect stellar companions, either stellar, sub-stellar, or planetary, it also comes with a drawback. The inclination of the system being unknown, the derivation of orbital parameters of the star can only lead to determination of the companion mass up to a factor depending on inclination. We present in this work exact mass derivations using astrometry with HIPPARCOS and Gaia. In particular, with the intermediate data of Gaia being yet unpublished, we developed the GASTON method to make use of Gaia released data to constrain the inclination of the systems studied here.

In Sect. 2 we present the target selection. In Sect. 3 we review the observations performed and the targets observed. In Sect. 4, the spectroscopic analysis of the SOPHIE observations is discussed, including the result of Keplerian fitting to the RV variations. In Sects. 5 and 6 we study the astrometric measurements made with HIPPARCOS and Gaia. In Sect. 7 we review the seven discovered BDs. Finally, in Sect. 8 we discuss the implication of the presented results on the BD desert localisation. We conclude in Sect. 9.

2 Target selection

The goal of the program in which this study takes place is to complete a meaningful unbiased statistical sample of companions detected within and about the BD mass regime, and with periods of up to 10 yr. Extracting BD candidates from a sample of stars, the selection function of which is well controlled, gives us the opportunity to constrain the location of the BD desert in terms of period and mass.

In the framework of the volume-limited FGK stars survey program for searching giant planets with the SOPHIE spectrograph (Bouchy et al. 2009; Hébrard et al. 2016), observers have collected RVs for many massive objects, including companions with M sini > 15M_J, on a time-span larger than 10 yr. This could allow the determination of the orbit of BD companions with periods as large as 10 yr.

In order to gather the largest possible number of BDs in the BD desert, and to be able to compare the BD population to the low-mass star population, we especially focused on sources with companion masses in the broad M sini range of 20–150 M_J. This range includes the whole BD regime, from the upper end of the giant planet domain, but also extends up to the late M-dwarf domain. We thus continued the RV monitoring of the sources that present any sign of a companion within 20–150 M_J with SOPHIE, along with the giant planet candidates below 20 M_J. Interested only in massive companions producing RV signals with large amplitudes, we aimed for a S/N per spectrum of at least 30.

Table B.1 summarises the basic information on the 54 targets covered by the present study. It includes only sources for which we gathered more than six RV data points and for which a meaningful Keplerian solution of the RV variations, or a lower mass limit beyond 150M_J, could be derived. We excluded SB2 sources from this publication.

3 Observations

The observations were performed with the SOPHIE spectrograph, fibre-fed from the Cassegrain focus of the 1.93 m telescope at the Haute-Provence Observatory (OHP, France). SOPHIE is installed in a temperature-stabilised environment and the dispersive elements are kept at constant pressure in order to provide high-precision radial velocities (Perruchot et al. 2008). The 39 spectral orders of SOPHIE cover the visible range between 3872 and 6943 Å. The spectra were collected in high-resolution mode, which leads to a resolving power of ~75 000 at 550 nm. During exposition of the spectrograph to the stellar photons in the science fibre, the instrument is also exposed to the background sky in a second fibre allowing subtraction of scattered light contamination in the science spectrum. The exposure time was varied to reach a S/N of at least 30 per resolving element under varying weather conditions.

Radial velocities are derived by the standard data-reduction pipeline (Bouchy et al. 2009), including spectrum extraction, telluric line removal, sky spectrum removal, CTI correction, cross-correlation function (CCF) computation, and barycentric Earth RV correction. In the reduction software, the CCFs are fitted by Gaussians to calculate the radial velocities of the sources (Baranne et al. 1996; Pepe et al. 2002). Moreover, the bisector spans (BIS) and full-width-at-half-maximum (FWHM) of each CCF are computed following Queloz et al. (2001). We did not correct the seasonal RVzero point variation from standard star variation (see e.g. Courcol et al. 2015), which is negligible (~m s⁻¹) compared to the expected velocity variation amplitude (~km s⁻¹).

In June2011, the SOPHIE spectrograph hexagonal fibres were installed, greatly improving the precision of the RV (Perruchot et al. 2011; Bouchy et al. 2013). Additionally a shift up to about 50 m s⁻¹ in the measured velocities was observed on standard stars (Bouchy et al. (2013). We therefore separated the data in June 2011 (JD 2 455 731.5). Before that date, the data are referred to as SOPHIE data, and after that date, they are referred to as SOPHIE+ data. Moreover, a systematic noise of 5 m s⁻¹ was quadratically added to the measured RV uncertainty of the SOPHIE data before June 2011 (Hébrard et al. 2016).

Additional non-SOPHIE data were found in the literature, with occasional, previously published orbits. These are summarised in Table B.2. We make use of these additional data to maximise the precision on the derived companion mass and period, and present relevant refinements of the already published companion parameters. Some of the public data were found in the SB9 catalogue¹ (Pourbaix et al. 2004).

Data points with less than half of the median S/N and large uncertainty on the RV measurements were treated as outliers and discarded. The number of points given in Table B.1 takes this into account, with an average of 18 SOPHIE spectra per star. Adding the other published data, the average number of RV points per source rises up to 27, with a minimum of 8 RV points per star and a maximum of 103. The RV coverage of all the stars spans between 475 days and 47 yr, with a median at 8 yr.

4 Spectroscopic analysis

4.1 Stellar parameters

The stellar parameters, effective temperature, surface gravity, microturbulence, and metallicity were derived using the spectroscopic analysis methods described in Santos et al. (2013), and references therein; see also Sousa et al. (2018) for more recent updates. The method makes use of the equivalent widths of a list of Fe I and Fe II lines in the SOPHIE spectra (see Table B.3), and assuming local thermodynamical equilibrium (LTE). The software used for the parameter derivation is the 2014 version of the MOOG software (Sneden 1973) with one-dimensional Kurucz model atmospheres. All derived parameters are given in Table B.3. We estimated the stellar mass of the primary star using the Torres et al. (2010) empirical relation. The log (g) values were corrected in order to be calibrated on log(g) derived using asteroseismology, following Eq. (4) in Mortier et al. (2014): (1)

The host stars presented in this paper are of type K5 to F8 (from the SIMBAD catalogue) on the main sequence with metallicities [Fe/H] ranging from − 0.3 to +0.3 dex. The stars HD 24505, HD 109157, and HD 204613 that were reported as (sub)giants in Simbad are rather located in the dwarf regime according to the present derivation. In particular, the spectral type of HD 204613 is reported in Simbad with the spectral type of a giant CH-star, G1IIIa:CH1.5 according to the analysis of photographic spectrogram done by Keenan & McNeil (1989). Interestingly, the photometry and colorimetry of this star tends to be more compatible with a dwarf (Ginestet et al. 2000). In agreement with the most recent published analysis of spectra of this source done with MOOG by Karinkuzhi & Goswami (2015), the present derivation leads to a G1V-IV, with an effective temperature of 5870 K, a surface gravity of 4.1, and metallicity of − 0.3 dex.

We added the average activity indicator to the table, which was calculated using all spectra of each target and the SOPHIE reduction software (Boisse et al. 2011). The uncertainties are estimated from the standard deviation of the mean, and an error of 0.1 dex was quadratically added to account for typical uncertainty of in SOPHIE spectra (Boisse et al. 2011). Our targets show a medium stellar activity level in general, with 17 targets having , a classical limit for separating active from weakly active stars (Santos et al. 2000).

Among the 39 more-active sources, 9 can be considered as highly active with . Nevertheless, the amplitudes of the derived RVs are all larger than a few hundred m s⁻¹, while activity is expected to influence RV measurements only at the scale of a few tens of m s⁻¹ (Campbell et al. 1991; Saar & Donahue 1997; Saar et al. 1998; Santos et al. 2000; Boisse et al. 2011). All the detections presented in this paper are those of true companions but the magnetic jitter of the most active stars will add scatter and imply larger uncertainty on the measurement of orbital parameters for the hosted companions.

4.2 Search for activity and binarity indicators in the CCF

We calculated the FWHM, the bissector span (BIS), and the S/N variations for all sources. This allows us to check whether the RV variations are polluted by the light of the secondary. This typically happens if the mass ratio q is greater than 0.6 (Halbwachs et al. 2014; Santerne et al. 2015). Any spectrum for which the FWHM or the bissector span showed anomalously large variation uncorrelated with S/N variations was systematically verified for a secondary peak. In the sample presented here, we selected only targets for which no spectra showed obvious secondary peaks.

To verify the absence of weaker secondary peak pollution, we used two indicators, namely the variation of the FWHM and of the BIS (Santerne et al. 2014, 2015). For each indicator, we performed a χ² -test of the “no-variation” null hypothesis, and calculated the Pearson correlation coefficient R of FWHM or BIS with RVs. SOPHIE and SOPHIE+ datasets were considered separately since the instrument update could have introduced changes in the reduction and the quality of the spectra (Díaz et al. 2016). Among all datasets, eight with less than four spectra were not analysed regarding these diagnostics, since any variations would be of little significance.

Initially, the errorbars of the FWHM and the bissector were calculated using the correspondence with RV errors σ_FWHM ~ (2−4) × σ_RV and σ_BIS ~ 2 × σ_BIS proposed in Santerne et al. (2014, 2015). These multiplication factors can be refined here, comparing the scatter of the FWHM or BIS to the median RV uncertainty for every source. The median factors found lead to

Fig. 1 P-value of the χ2 -test vs. the Pearson correlation coefficient, R, for both FWHM (orange) and BIS (blue) indicators. The radii of the symbols are linearly scaled with the number of points considered. The dashed lines show the limits of significance for R (±0.5) and for the p-value at 3σ (<0.01).

(2)

For the BISwe confirm the result of Santerne et al. (2015) but for the FWHM, we found that the multiplication factor stands higher. These corrected factors were used eventually to calculate the χ² test and the Pearson correlation coefficient that are summarised in Table B.4 and presented in Fig. 1.

Only ten systems show an FWHM or BIS dispersion that is significant with p-values lower than the 3σ limit. There is a single source however that also shows a strong correlation coefficient of FWHM with RV variations, R(FWHM, RV) = −0.82, as can be seen also in Fig. 2. This system is HD 77712, a K-type star with a medium activity level at . On the other hand, it presents no significant variations of the bissector. This is similar to the case of a triple system studied in Sect. 2.9 of Santerne et al. (2015) with the pollution of the CCF from a weak secondary peak always present at fixed RV. A possible explanation is therefore that HD 77712 is a triple system with a long-period binary whose secondary is polluting the CCF and a shorter-period binary with a dark companion. The RV amplitude and M sini we report for this system are likely underestimated.

4.3 Keplerian orbit fitting

The yorbit software (Ségransan et al. 2011) was used to calculate the solution, which uses a genetic algorithm to refine initial parameters for a Levenberg-Marquardt optimisation. This leads to priors for an Markov chain Monte Carlo (MCMC) estimation of error barsfollowing Díaz et al. (2014, 2016). The MCMC was applied on 1000 iterations. The varied parameters are the period P, the RV amplitude K, the eccentricity e, the angle of periastron ω, the periastron passage time T₀, and the offsets γ_S and γ_S, + for SOPHIEand SOPHIE+ datasets, respectively. Specific additional offsets are used for supplementary datasets as indicated in Table B.2. When the RV error bars of the previously published data are not given, they are uniformly fixed to the unbiased standard deviation of the residuals as soon as a good orbital solution is found. Following Anderson et al. (2012), eccentricities compatible with zero at 2σ were subsequently fixed to zero and the solution recalculated with these new constraints. In this case, T_p indicates the epoch of the transit – if the system were to be edge-on.

The final parameter values given in the results section are the median of the MCMC distribution and the symmetric error bars calculated by the standard deviation of the MCMC distribution. The error bars defined by the confidence interval (CI) at 68.3% around the value with the best likelihood are barely asymmetric, while the difference between the median and best-fit value is not significant. We found that the standard deviation gives more conservative uncertainties than the CI at 68.3%. We therefore used the standard deviation as an error bar in order to remain conservative, especially for cases with inaccurate derivation of the orbital parameters. For incompletely covered orbits, the MCMC distributions have non-Gaussian tails: a confidence interval equivalent to 3-σ is almost certainly not simply three times broader than that for 1-σ.

For six stars, HD 5470, HD 153376, HD 193554, HD 207992, HD 212735 and BD+212816, SOPHIE and SOPHIE+ data points are not separately sufficiently numerous to derive a meaningful solution, with Δγ = γ_S, + − γ_S being on the order of 50 m s⁻¹ at most. In those cases, we fixed Δγ = 0 to derive the solution.

All the results of the Keplerian fits are summarised in Tables B.5, B.6, and B.8. We find 51 binary systems and 3 triple systems. These can subsequently be divided into two categories of companions, BD candidates and M-dwarfs. We have 11 binary BD candidates in the mass range 15–90 M_J (Tables B.5), and 40 binary companions in the M-dwarf regime (Table B.6). These are presented in more detail in Sect. 4.6.1 below. The results for the three triple systems are presented in Table B.8 and in Sect. 4.6.2. They include one BD candidate and two M-dwarfs, and in both cases, a drift which requires a companion mass above 5 M_J. For the targets with additional public data from the literature, residuals O–C and RV centre-of-mass offset γ of each additional dataset are given in Table B.7. Finally for the objects in the 15–90 M_J mass range, Keplerian solutions and residuals are plotted in Fig. C.1. The Keplerian solutions and residuals for triple systems are plotted on Figs. C.4. And for the objects beyond 90 M_J the solutions are shown individually in Fig. C.3.

In general, the fits are accurate with better than 7% precision on the orbital elements in 90% of the cases, and a median precision of at most 1%. A few cases show however a highly inaccurate derivation of orbital elements that is due mainly to an incomplete covering of the full orbital phase. For HD 85533, although the uncertainty on the period is ~100%, the given value is a lower limit, and the companion should be at least as massive as 450 M_J. On the other hand, the eccentricity is surprisingly accurate with an error of only 20%. This results from a better coverage of an inflexion in the RV curve that leads to a good fit only for eccentricities larger than 0.44. This stands also for HD 13014, as well as HD 40647, HD 60846, and HD 146735, for which the period, RV amplitude, and companion mass, already in the M-dwarf domain, are likely underestimated, while the eccentricity is conversely better constrained.

Finally, the O–C residuals lie below 10 m s⁻¹ except for a few active sources that show much larger dispersion of residuals close to 40 m s⁻¹. We discuss the distribution of the residuals in more detail below in Sect. 4.4, especially comparing the SOPHIE and SOPHIE+ datasets, and describe comparisons with activity indices for the observed sources in Sect. 4.5.

Fig. 2 FWHM (orange) and bissector span (blue) variations of the cross-correlation function of HD 77712 spectra. The error bars are calculated from the RV uncertainties, multiplied by the factors given in Eq. (2).

4.4 Comparing SOPHIE and SOPHIE+

Analysing the O–C residuals of the Keplerian fits allows us to verify the accuracy of SOPHIE data and in particular to compare the quality of the measurements before and after the instrument upgrade in June 2011. The standard deviation of the residuals can give the actual precision of the RV measurements, because the targets in the sample are of similar spectral type, with similar CCF shape and FWHM below ~10 km s⁻¹, as shown in Table B.4. Moreover, with 23 targets observed with both SOPHIE and SOPHIE+ instruments, we can characterise the typical RV offsets between the two datasets. Summarising the data gathered in Tables B.5 and B.6, we show the distribution of O–C values in SOPHIE and SOPHIE+ configurations separately, and the RV offset distribution in Fig. 3.

The distribution of RV offsets between SOPHIE and SOPHIE+ is centred on (3)

This is compatible with the results found by Bouchy et al. (2013) that bound the RV shift due to the upgrade to 0–50 m s⁻¹. To obtain this distribution, we assumed that the offset should not exceed 100 m s⁻¹, in which case the offset should be better explained by a slow drift, due to a third companion. This led us to consider a few systems as being multiple rather than binary, as shown in Sect. 4.6.2 below.

The core of the distribution of the standard deviation of O–C residuals leads to the following estimation of RV accuracy for SOPHIE and SOPHIE+ configurations: (4)

These values are in line with the results obtained by Hébrard et al. (2016) where a median RV accuracy of about 7 m s⁻¹ is derived for observations before June 2011, and of about 3.5 m s⁻¹ for observations taken after the instrument update. The values derived here are higher, which might be explained by the activity index greater than −4.65 for about half of our sample, or the absence of instrumental drift (~m s⁻¹) corrections, as derived in Courcol et al. (2015), in the present study. Moreover, our observed sample also includes spectra with S/Ns down to 30, while the Hébrard et al. (2016) sample only includes spectra with S∕N > 50.

Comparing the O–C individually for every target confirms that the general tendency is of a reduction in the RV dispersion after the upgrade of the instrument towards a value close to 5 m s⁻¹. If this is not the case, it should be explained in terms of supplementary signal in the RV due to either activity jitter or planetary signal. Among our sample, only HD 23965, HD 40647, and HD 161479 have large O–C dispersion >14 m s⁻¹ in both datasets. This is however most likely explained by their significant activity index (Table B.3). We can therefore exclude that more planetary signals are hidden in any residuals beyond an amplitude of about 8 m s⁻¹.

Fig. 3 Upper panel: O–C distribution of SOPHIE (blue) and SOPHIE+ (orange) datasets. Lower panel: RV offset distribution between SOPHIE and SOPHIE+ datasets.

4.5 Residual dispersion and magnetic activity

The residuals of the Keplerian fit, σ_v, can be compared to the activity index to verify whether or not activity can explain the amplitude of the residuals in general and in specific cases where this amplitude is exceptionally large. Saar et al. (1998) and Santos et al. (2000) previously studied the correlation between RV dispersion and magnetic activity. Here, we followed the same procedure as these latter authors in order to make a point comparison. We first exclude datasets with less than 7 pts; we then quadratically subtract the mean internal RV error of all SOPHIE exposures ⟨σ_i ⟩ from the dispersion of the RV residuals for every target, . This should leave only variations from the instrument itself and magnetic activity. Sources for which σ_v is smaller than ⟨σ_i ⟩ were excluded from this analysis. Figure 4 plots and as derived from our sample and compares these to the relation obtained by Santos et al. (2000).

We observe that the dispersion of the residuals correlates well with the magnetic activity, with only a few outlying points. However, we see a discrepancy between our values and the relation derived for G-type stars in the CORALIE sample by Santos et al. (2000), where they find that . In our case, the slope is stronger, with a linear fit of the log–log relation instead leading to (5)

The uncertainty of the fit is σ_fit = 0.3 dex. The slope is closer to the relation obtained by Saar et al. (1998), . After the exclusion of the F and K type stars of our sample, keeping 0.6 < B–V < 0.8, there remain 18 G-type stars. This leads to a similar relation but a larger fit uncertainty of 0.4 dex.

The most significant outlier in Fig. 4 at ~4.75 and ~ 30 m s⁻¹ is HD 207992. We collected 11 RV points in the SOPHIE configuration, but only 2 with SOPHIE+ for this source. The RV curve in Fig. C.3 indeed shows variability in the residuals, which could be due to a supplementary signal for this relatively low-activity star. In Table B.4 we do not see any significant BIS or FWHM variations. We conclude that this signal could be tentative evidence of a third object in the system of HD 207992.

One other case is HD 161479 with σ_v,S+ = 36 m s⁻¹ and σ_v,S− = 42 m s⁻¹. This residual dispersion is large, but might be compatible with magnetic activity since = −4.42 for this K0 star. Moreover, according to Table B.4 the bissector and FWHM variations are relatively significant. We measure a p-value of 0.001 for the no-var model of the FWHM in the SOPHIE+ dataset, and a strong correlation of −0.95 for the bissector in the SOPHIE dataset, although based on only four points. We conclude that the supplementary RV variability of HD 161479 is most likely due to magnetic activity.

Fig. 4 (in m s−1) vs. for SOPHIE (in blue) and SOPHIE+ (in orange) datasets. The symbol size is proportional to B− V with values between 0.5 and 1.2. The black lines represent the relation derived by Santos et al. (2000) for G-type stars with a fit uncertainty of 0.18 dex. The red lines represent the relation derived from the datasets of this work, with a fit uncertainty of 0.3 dex.

Fig. 5 Period vs. M sin i for the systems studied in this paper. The crosses show the new results of this paper. The 13 and 90 MJ limits are drawn as dotted lines. The plain white symbols represent stars with two companions. HD 71827 b and c are represented as green triangles, HD 212735 b and c as red stars, and BD+212816 b and c as blue squares. Most error bars are smaller than the symbols.

4.6 Results of the Keplerian fit

In total, we characterised 54 massive companions in 54 different systems. We report the Keplerian orbit and M sini measurements of 12 BD candidates in the extended range 15–90 M_J. One among the 12 is part of a triple system, HD 71827, whose discovery is reported here. We also characterised the orbit of 42 stellar companions with masses in the M-dwarf regime 90 M_J–0.52 M_⊙. Two BD candidates lie in the grey zone between the classical upper limit of BDs and the lower limit of M-dwarfs, 80–90 M_J.

Moreover, we recall that the constraint on the mass obtained from velocimetry is only a lower limit because of the uncertainty on the inclination of the systems implying an unknown value of sin i. We see in Sects. 5 and 6 that thanks to HIPPARCOS and Gaia astrometry we are able to add constraints on the inclination and thus the true mass for 46 of the candidates presented here.

The M sin i-period diagram summarising the results is shown in Fig. 5. The periods of the derived orbits are large in general, with only nine companions having a period of below 100 days. Among these nine, one is member of a triple system and has an M sini within the BD regime. Eccentricities are also large, with only seven orbits with e < 0.1. The eccentricities are dispersed around 0.42 ± 0.27. Figure 6 shows the period–eccentricity distribution of our results, and compares it to that of the massive planets collected in the Exoplanet.eu database with M > 4M_J. We selectedsystems exclusively compatible with the constraints of our survey (δ > 0°, +0.35 < B–V < +1, d < 60 pc, ±2 mag from MS). The period–eccentricity distribution of the BDs reported in this work agrees with that of giant exoplanets. The eccentricities of giant exoplanets are fully compatible on average with the eccentricities of BDs with e_GP ~ 0.42 ± 0.22. This is in line with the conclusions of Sozzetti & Desidera (2010), who find strong similarities in terms of eccentricity distributions between massive planets and BDs.

One candidate stands apart with a short period and large eccentricity, BD+362641, for which the orbit is however not well constrained due to the small number of points (N_RV = 9). However, the large RV variation observed of ~40 km s⁻¹ places it in the M-dwarf regime with a mass most likely larger than 200 M_J.

Fig. 6 Period vs. eccentricity for the systems studied in this paper. The blue circles show the new results of this paper, while the green circles represent the Exoplanet.eu database with primaries verifying that δ > 0°, +0.35 < B–V <+1, d < 60 pc. The symbol size is proportional to the logarithm of M sin i.

4.6.1 Binary companions in the BD and M-dwarf regime

Among the 12 detected BD candidates, 11 are components of a binary system. They have orbital periods shorter than 30 yr, or semi-major axis smaller than 10 au. Eight of these BD candidates are new discoveries, among which we report six with M sini strictly below 80 M_J. This is a significant increase of the number of known BD candidates. We notice that the orbital period of all these companions is larger than 100 days, even though massive companions with a minimum mass close to but larger than 90 M_J with an orbital period as low as 40 days are also reported. Interestingly, adding BD detections around solar-like stars in the solar neighbourhood that are reported in previous papers tend to confirm this distribution. We discuss all the consequent improvements these new detections bring on the statistics of objects in the BD regime in Sect. 8.

Four objects, HD 28635, HD 210631, HD 211681, and HD 217850, were already published as BD candidates. Improvements on their orbital parameters and M sin i are summarised below.

HD 28635

Also known as vB 88, HD 28635 was reported to host a BD companion with an approximate spectroscopic mass of 70 M_J using Keck/HIRES data (Paulson et al. 2004). Adding 13 SOPHIE data points and 3 Elodie data points, we find that the RVs are compatible with a BD-mass companion at a period of 2636.8 ± 2.2 days with M₂ sini ~ 77.1 ± 2.7M_J and a₂ ~ 4.014 ± 0.068 au.

HD 210631

Latham et al. (2002) reported a 82 ± 6 M_J companion in this system. Adding SOPHIE data, we confirm this result, finding compatible minimum mass of 83.4 ± 6.9 M_J, with a period of 4030 ± 40 days at a separation of 4.976 ± 0.085 au.

HD 211681

The companion of this subgiant G5 was reported as a low-mass star with a minimum mass in the range 72–100 M_J by Patel et al. (2007) using Keck/HIRES data. Adding 30 SOPHIE and 12 ELODIE measurements, we are able to constrain the M sini range of the companion to 77.8 ± 2.6 M_J with a period of 7612 ± 131 days at a semi-major axis of 8.28 ± 0.16 au.

HD 217850

The radial velocity variations of this G8-type star were reported to be compatible with an 11 M_J companion by Butler et al. (2017) who used an incomplete coverage of the orbit with Keck/HIRES data. Adding 41 SOPHIE data points we find the lowest-mass BD of our sample with an orbital period of 3508.2 ± 2.6 days, M₂ sini ~ 22.27 ± 0.77 M_J, and a₂ ~ 4.672 ± 0.079 au. This is the candidate BD with the lowest M sin i in our sample.

Finally, among the 42 massive companions in the M-dwarf regime, 40 form a binary system with their host star. For 24 of them, to our knowledge, this is the first publication of an RV orbital solution. For the 6 systems for which an RV orbit was already published, the last 2 columns of Table B.2 summarises the improvement on the M sini for these stellar companions.

4.6.2 Triple systems

We found evidence for a secondary drift signal in the RV data of three stars, BD+212816, HD 71827, and HD 212735. The result of fitting a single Keplerian and a drift for each system are summarised in Table B.8 and Figs. C.4. In order to derive a minimum estimation of the mass of the second companion, we fitted the drift signal with a Keplerian with the shortest period possible compatible with a drift. For the three sources, the RV offset between SOPHIE or ELODIE, and SOPHIE+ datasets is significantly larger than 100 m s⁻¹. It should be on the order of 10 ± 30 m s⁻¹ between SOPHIE and SOPHIE+ measurements (see Sect. 4.4 below) and on the order of 50–100 m s⁻¹ between ELODIE and SOPHIE+ (Bouchy et al. 2013). This is the sign of a real drift due to a third companion in the system. We had to fix the RV offset to γ_S,+ − γ_S = 0 km s⁻¹ in order to derive a Keplerian solution with a supplementary linear drift. Since the Keplerian of the drift signals cannot be constrained, we only report them here and do not include them in any other analysis in the rest of the paper.

HD 71827

This is a triple system composed of an F8-primary surrounded by one BD and a low-mass star. The 26 M_J-BD stands at a short period of 15 days. This is the only BD in our present sample with a period shorter than 100 days, and it is interesting to note that it is also part of a triple system with possible dynamical interaction. There is clear evidence in SOPHIE+ data of a second signal with a large period and we can confirm the presence of a cubic drift. The shortest-period orbit found compatible with the drift leads to a minimum mass ~163 ± 7 M_J for a companion on a an orbit of 20 yr.

HD 212735

Apart from an obvious 38-day period signal, the RVs of this system display a significant linear drift during the 10 yr of data. Fitting a second long-period Keplerian to the drift signal leads to a minimum estimation of the period and the mass beyond 20 yr and 47 M_J, respectively, for the tertiary. The outer companion is therefore possibly a BD, but is most likely an M-dwarf with a much larger period.

BD+212816

The secondary companion of this K0-type star is an M-dwarf, but a supplementary long-period signal might be present as a drift. However, this linear drift is compatible at 2σ with a constant. It should be considered as a possible yet unconfirmed triple system. The mass of the outer companion could be as low as 5 M_J, but is likely much higher.

For all these triple systems, future Gaia data releases or direct imaging could help to probe for the third companion. In every case, the semi-major axis of the outer orbit is larger than 7 au, with a parallax on the order of 20 mas, and could therefore be seen with adaptive optics that can probe down to about 100 mas in the neighbourhood of stars.

5 HIPPARCOS astrometry

Complementary to the RV orbital derivation, the HIPPARCOS astrometry can allow us to constrain the inclination of the systems, as was performed by for example Sahlmann et al. (2011); Díaz et al. (2013), and Wilson et al. (2016).

For all 54 systems of our sample, the new HIPPARCOS reduction catalogue (van Leeuwen 2007) provides informations on the type of fitting solution (“5” for standard, “X” for stochastic, and “G” for accelerated solutions; see e.g. Perryman et al. 1997 or Lindegren et al. 1997) for number of field-of-view (FoV) transit, measurement time-span, and abscissa measurement errors. A summary of this information is presented in Table B.9.

After a preliminary analysis of all systems, we found 16 of them to show indications of significant orbital motions in the Intermediate Astrometric Data (IAD), plus 1 system, HD 193554, already solved in the HIPPARCOS double star catalogue (ESA 1997), and 18 systems for which it could be possible to derive an upper limit on the astrometric motion due to the massive companion. Outliers in the IAD had to be removed because they can substantially alter the outcome of the astrometric analysis. The result of the Keplerian fit of the HD 193554 astrometric motion analysis done by the HIPPARCOS team is given in Table 1; it compares well with our RV derivation. The true mass estimation for the companion is beyond the 90 M_J limit.

For the 16 systems with a significant orbital solution, we fitted the astrometric measurements with a seven-parameter model, in which the free parameters are the inclination i, the longitude of the ascending node Ω, the parallax ϖ, and offsets to the coordinates (Δα^⋆, Δδ) and proper motions (, Δμ_δ). The other orbital parameters are fixed according to the RV results given in Tables B.5 and B.6. A two-dimensional grid in i and Ω was searched for its global χ²-minimum. The statistical significance of the derived astrometric orbit was determined with a permutation test employing 1000 pseudo-orbits (Sahlmann et al. 2011).

For all 16 sources except two, we detect the astrometric orbit with a significance > 2σ. Those are listed in Table B.11 with their orbital solution. Table B.11 lists updated parallax, proper motion, coordinate offset, inclination, and ascending node of the orbits. The updated parallaxes are compared to the DR2 parallaxes given in Table B.1. Moreover, the updated proper motions are compared to the Tycho-Gaia Astrometric Solution (TGAS) proper motion, which should be closer to the actual proper motion of systems since it is based on a 24-yr baseline of astrometric data. Finally, Figs. C.15 and C.16 show the significant orbits.

In general, the updated HIPPARCOS-2 parallaxes are not compatible with the Gaia DR2 parallax at the 1-σ level. For two systems, HD 133621 and HD 155228, the discrepancy is larger than 3-σ. This shows that accounting for the orbital motion can lead to strong corrections of the published parallax on the order of ~10%. In addition, the comparison of the HIPPARCOS-2 proper motion corrections, after fitting the orbital motion, with the Gaia DR1 proper motion shows good agreement in general, validating the solutions and corrections proposed in Tables B.10 and B.11. Indeed, we expect the proper motion derived in the TGAS sample of the DR1 to be closer to the true linear proper motion of the system, since for those sources the astrometric solution takes into account HIPPARCOS-2 and Gaia measurements along a 24-yr baseline.

In only one case, HD 87899, there is no strong agreement of the proper motion corrections. With a long period of 4.2 yr, the phase coverage of the HIPPARCOS-2 measurements is only partial, and the proper motion corrections are quite uncertain with a given error of 2.8 mas yr⁻¹. Thus, the derived orbital parameters for HD 87899 should be considered only conservatively within their 3-σ error bars. If the Gaia DR1 proper motion is correct, the semi-major axis as derived with the HIPPARCOS-2 data should rather be around 10–15 mas and the mass closer to 0.2 M_⊙, which is more in line with what is derived using Gaia data only in Sect. 6.

For the two stars HD 110376 and HD 155228, the F-test of the orbital model and the permutation test yield significantly discrepant results. The F-test indicates orbit detection whereas the permutation test is inconclusive. Usually, this is caused by strong fit-parameter correlations that skew the average semi-major axis estimation and therefore the result of the permutationtest. For HD 110376 and HD 155228, however, this is not the case and the exact reason for the failure of the permutationtest is unclear. Because the orbit sizes are relatively large, the F-test null probabilities are very small (2.2 × 10⁻¹² and 3.2 × 10⁻⁷, respectively), and the acceptable i-Ω parameter space is well constrained, we present orbit solutions for these two sources as well. As can be seen in Table B.11 the significance is lower than 1σ for these two systems. We also note that the parallax change caused by fitting the orbit model is large (almost 3 mas) for HD 110376.

For HD 225239 and HD 62923, the derived secondary mass is larger than the primary mass, which could be caused by light contribution by the secondary, shifting the position of the photocentre out of the primary star centre. Indeed, our model assumes that the companion is dark (Sahlmann et al. 2011). We are developing supplementary methods to treat these cases and will report results in an upcoming publication. Here, we note that the orbit detection in both cases is significant but that the semi-major axis a refers to the photocentric orbit and that the derived secondary masses are incorrect. Other possibilities could be that the companions of these stars are actually massive dead stars such as white dwarfs, neutron stars, or black holes; they could also be couples of low-mass stars. This could be consistent with the log(g) of the two primaries (Table B.3) that have a lower surface gravity (~4.1–4.2) than that expected for G2-3 dwarf stars (~4.4–4.5). Therefore, these two primaries might rather be more evolved subgiants. Precise astrometry with Gaia and imaging can allow the determination of the exact mass and nature of these companions.

Two BD candidates, BD+210055 and HD 210631, have their mass re-evaluated above 90 M_J. For BD+210055, as guaranteed by the good coverage of the orbit (N_orb = 0.9), the fit by the astrometric model is excellent, with a significance close to 100%. This leads to a real mass that is significantly larger than the M sini ~ 85M_J derived thanks to RV only, with M₂ between 140 and 290 M_J, well within the M-dwarf regime. For HD 210631, the orbital coverage of 0.3 is not ideal, owing to the long (11 yr) orbital period. Still, the fit of the astrometric motion could catch some significant acceleration in HIPPARCOS data points and lead to a 2σ-detection. This shows that the M sini ~ 83M_J of the companion of HD 210631 derived by RV was strongly underestimated compared to its real mass, here constrained to lie between 140 M_J and 1.5 M_⊙. We emphasise that the upper mass range (about >0.6 M_⊙) neglects the fact that in such a domain the secondary might contribute light and even produce a secondary peak in the CCF that we do not detect in our observations.

Finally, we derived upper limits on the astrometric semi-major axis of the primary and the mass of the companion for 18 sources. Provided that at least about 80% of the orbit is covered, the upper-limit of an undetected semi-major axis can be deduced from the value of the median measurement precision σ_Λ. The formula is the one used in Sahlmann et al. (2011), but moreover assumes the most unfavourable case of an edge-on orbit whose projection on the plane of the sky only presents its minor-axis: (6)

The values of the upper limits on semi-major axis of the primary and the corresponding companion mass are added to Table B.11. For the triple systems HD 71827, HD 212735, and BD+212816, only the inner companion b orbit was considered, since the outer companion is not constrained by RVs. Unfortunately, because of the loose constraints on the mass of the orbiting companions, we cannot exclude that all these systems are stellar binaries in the M-dwarf regime. We see in Sect. 6 that using Gaia published astrometric data allows us to tighten the constraints on the mass for several of these systems, including systems for which N_orb < 0.8.

6 Gaia astrometry

To overcome the large uncertainties on the true masses obtained using HIPPARCOS, we also cross match our sample with the Gaia catalogue. We found 51 of our 54 targets in the Gaia DR1 catalogue (Gaia collaboration 2016). Among these 51 systems, we were able to measure the companion mass for 33 of them, and derive upper limits for a further 6. The companion mass of the 12 that remain out of the 51 systems could not be constrained further with Gaia data. Their masses were nonetheless already bounded from below thanks to RV, and this is well within the M-dwarf domain.

The DR1 of Gaia does not provide the individual positional measurements for the whole of our sample of stars. However, two binarity indicators are published in the released catalogues: the astrometric excess noise ϵ and the TGAS discrepancy factor ΔQ (Lindegren et al. 2012; Michalik et al. 2014; Rey et al. 2017).

The astrometric excess noise ϵ is a measure of scatter around the five-parameter astrometric solution as resolved by the Gaia reduction software. The RV orbital parameters and the dates of Gaia data collection for the DR1 can be used to derive an estimation of the inclination of the system, and thus of the true mass of the companion. To do so, we applied an MCMC method, presented in the following Sect. 6.1, that is able to output possible inclinations for a given astrometric excess noise and fixed orbital parameters.

As published in the DR1, the dimensionless quantity ΔQ calculates the difference between the proper motion published in the HIPPARCOS-2 catalogue and the proper motion derived in the TGAS sample by Gaia (Lindegren et al. 2016). We note that this differs from the original ΔQ definition as given in Michalik et al. (2014). The proper motion derived in the TGAS is based on a 24-yr baseline of astrometric measurements, the 4-yr monitoring of HIPPARCOS-2, and the 14-month monitoring of Gaia for the DR1. Measuring a significant long-term astrometric displacement, it can be used as a binarity diagnostic (Michalik et al. 2014; Lindegren et al. 2016). Comparing the value of ΔQ for every source in our sample with the typical value obtained for any source in the DR1 allows us to determine if a system is a likely astrometric binary. This analysis is performed in Sect. 6.2

The DR1 archive provides both quantities, while only ϵ can be found in the DR2 (Gaia collaboration 2018). Moreover, the excess noise values from the DR1, although based on a shortertimeline of astrometric measurements (25 July 2014 – 16 September 2015, or 416 days), are more reliable than in DR2 because of the so-called “DOF-bug” that directly affected the measurement of the dispersion of the final astrometric solution (Lindegren et al. 2018). For these reasons, we only used the DR1 results for ϵ and ΔQ, as extracted from Gaia DR1 archives². They are presented in Table B.12.

6.1 GASTON: Gaia Astrometric Noise Simulation to derive Orbit incliNation

As such, it isnot possible to directly interpret ϵ as a measure of the semi-major axis of an astrometric orbit, because it highly depends on the inclination of the orbit of the system that is seen projected on the plane of the sky. Since the fit of the RVs leads to precise orbital parameters however, the inclination is also the only remaining free parameter that could have an impact on the value of ϵ.

We introduce here the new GASTON method based on Gaia data simulation to derive the inclination of the system from the measurement of astrometric excess noise. The photocentre semi-major axis and secondary masses derived using this method are given in Table B.13.

6.1.1 Basic principle

The principle of this method is to simulate Gaia photocentre measurements along the derived RV orbits presented in Sect. 8. Measurement epochs and Gaia along-scan (AL) axis orientations are randomised along the RV orbit, bounded by the DR1 data collection epochs. Real measurement epochs and AL axis orientations available online³ compare well with random values. Considering random epochs and AL-axis orientation is therefore sufficient for applying the method we present here, which makes use of the excess noise, a quantity that cannot be considered as accurate.

Different inclinations can be tested, each leading to a simulated astrometric excess noise ϵ_s. We then constrain the different possible inclinations by comparing the whole set of ϵ_s with its actual measurement in the DR1, ϵ_DR1. As a result of applying this method on our target sample, we found that the astrometric excess noise follows a one-to-one correspondence with inclination, owing to the increase of the photocentre semi-major axis with decreasing inclination.

A few effects introduce scatter into this relation. First of all, the DR1 excess noise may incorporate bad spacecraft attitude modelling, which means that the value of ϵ does not account only for binary motion (Lindegren et al. 2012). The amplitude of the bad attitude modelling within ϵ could be estimated from its median value in the full sample of objects observed with Gaia (Lindegren et al. 2016), ϵ_med = 0.5 mas. Any value of excess noise below that value cannot be trusted to be genuinely astrophysical, although it could be considered an upper-limit. Conversely any value of ϵ above that level likely contains true binary astrometric motion. To take this effect into account, we added a bad-attitude-modelling noise of 0.5 mas to Gaia measurements in the simulation.

Secondly, the astrometric motion of sources whose orbit has a period close to 1 yr could be modelled by an excess parallax if the orientation of the system coincides with that of parallax motion. Moreover, slow orbital motion with a period on the order of the Gaia–HIPPARCOS baseline (~24 yr) can be absorbed into an erroneous proper motion. These effects cannot be properly taken into account in our simulations, since we have no prior knowledge of the orientation Ω for any of our targets. Thus, the simulated ϵ_s could be overestimated compared toϵ_DR1 for the sources with periods closeto 1 yr or larger than ~20 yr. This tends to lead to underestimation of the inclination and overestimation of the exact mass of the companion.

We first describe the method to calculate the orbital model and the simulated along-scan Gaia measurements, and then how to derive a simulated excess noise for a given inclination. Once this relation is established, we can derive an interval of inclinations compatible with a given value of ϵ_DR1.

6.1.2 Modelling of the along-scan data

In a fixed non-accelerating reference frame, any source has a position vector u_⋆ = (x_⋆, y_⋆) in the plane of the sky. We assume that proper motion, annual parallax, and attitude of the spacecraft have been properly modelled and subtracted with only the orbital motion of the star remaining. Using the RV orbital model derived in this paper, it is fairly easy to obtain the projection of the star’s orbit with a given inclination on the plane of the sky, and to derive u_⋆ with respect to Keplerian parameters and inclination I_c: (7)

where u_x is an arbitrary direction in the plane of the sky, the direction u_z is orthogonal to the plane of the sky and oriented toward the observer, and the remaining direction u_y is directly oriented with respectto u_x and u_z composing the triad (u_x, u_y, u_z). The position vector k(t|P, a, e, T_p) at epoch t is that of a Keplerian orbit whose periastron is oriented along the u_x direction before applying the rotation matrices R_z(ω) and R_x(I_c).

A number of locations are randomly selected along the above orbit by drawing random epochs between the bounds of Gaia DR1 data collection (t_min = 2 456 863.0, t_max = 2 457 282.0). The number of these locations is given by the number of “matched observations” in the DR1; it is the total number of field-of-view N_FoV CCD transits of a given star captured by Gaia. This value is given in Table B.12. At each epoch, there are between 1 and 9 measurements of the AL angle η per FoV transit (Lindegren et al. 2016). For simplicity, we assume a uniform number of CCD transits per epoch, given by N_rec = round(N_tot∕N_obs). Therefore, for each of these locations, we simulate N_rec measurements of η along the AL direction.

The AL axis u_AL(θ) is defined independently at each location with a random orientation θ. The measurements are picked randomly along this axis accounting for the uncertainty on the spacecraft direction at any epoch (0.5 mas; see the preceding section), and the uncertainty on AL measurements during the transit of the target on the CCD (~0.4 mas; see Lindegren et al. 2018). Thus, for any FoV transit observation indexed i, at an epoch t_i, and N_i CCD measurements indexed j, the simulated read-out AL angles are given by (8)

where u_⋆ is projected on the AL direction, and ξ_inst,i and ξ_AL,j are the instrumental and AL errors introduced above.

When dealing with an unresolved binary star, Gaia actually measures the position of the photocentre on the plane of the sky. The photocentre motion has the same orbital parameters as the primary, except for the semi-major axis. At a fixed system inclination, we use the RV orbital solution and a mass-luminosity empirical model for both binary components in order to calculate the photocentre semi-major axis as described in the Appendix A. Since no secondary peaks were seen in any of the CCFs for all targets, we always assumed that M_V,2 − M_V,1 had to be greater than 2.5, and thus the luminosity fraction in the optical range is <10%.

Using the photocentre semi-major axis in Eq. (8), along with all other orbital parameters from the RV Keplerian solution, leads to the final simulated Gaia measurements. Examples are given in Figs. 7 and 8. From these simulated data, we are now able to calculate an astrometric excess noise.

Fig. 7 Left: inclination vs. simulated excess noise for BD+192536. The red lines mark the value of for this star ± 10%. Top-right: example of a simulation of Gaia measurements (blue points) for the peculiar case of BD+192536 that fits the value of ϵDR1 = 1.02 mas. The true proper motion and true parallax are assumed to be subtracted. The purple line is the residual proper motion (moving centroid) fitted to the simulated measurements and the red cross marks the true centre of gravity of the system. Bottom-right: residuals with respect to the moving centroid.

6.1.3 Simulated excess noise

The astrometric excess noise is obtained by estimating the χ² of its η_AL-residuals around the five-parameter solution derived by Gaia reduction software (Lindegren et al. 2012). In the simulations, we did not account for the true proper motion or the parallax, assuming them to have been modelled out already. This results in only two remaining parameters to model out of our simulations, namely the (x, y)-position of the photocentre on the plane of the sky.

The “average” target position published in the DR1 is given by the centroid of η_AL measurements u_c = (x_c, y_c). Assuming that most systematic positional errors have been accounted for, with only the uncertainties ξ_inst,i and ξ_AL,j introduced above remaining, the centroid position is found by minimising the squared sum of residuals R_ℓ with ℓ a given observation

Fig. 8 Same as Fig. 7 but for HD 71827. Here the value of ϵDR1 is larger than any of the simulations provided that the secondary companion is dark (i.e. not visible as a secondary peak in the CCF). The assumption that the inner companion is responsible for the astrometric motion of HD 71827 is likely to be incorrect (see text for explanation).

(9)

This leads to a simple system of four linear equations, which can be inverted, solving for (x_c, y_c). Once u_c is derived, this expression also leads to the χ² of the residuals. In Lindegren et al. (2012) the expression of the χ² and of the excess noise with respect to residuals and the AL uncertainty is given by (10)

Here we assume that the down-weighting factors w_ℓ = 1 since we are only interested in the good AL measurements (with w_ℓ ~ 1). The χ² should follow a χ² distribution with a mean value equal to the number of degrees of freedom, which is the total number of points minus the number of parameters of the astrometric model derived by Gaia, thus N_DOF = N_tot − 5. Therefore, at a given inclination I_c, and assuming a uniform value of σ_AL along all observations, we should solve (11)

The above equation can be solved for I_c by performing the simulations at various inclinations and comparing the right-hand side of Eq. (11) to the value of that is measured by Gaia in the DR1. We sampled the inclination on a grid of 10 000 values uniformly distributed between 0 and π∕2. Each time, the full set of inclinations compatible with ϵ (±10%) leads to a range of possible values of the semi-major axis of the photocentre and the companion mass.

The bounds γ_± of a given parameter γ compatible with ϵ are obtained by solving the following Bayesian equation for diverse posterior probabilities p: (12)

with for example p = 0.68 leading to the 1-σ bound and p = 0.5 the median. To solve this equation, we assumed that ϵ_DR1 is conservatively known at ±10% . The prior P(γ > γ_±) is calculated by assuming that the unknown inclination is uniformly distributed between 90° and the inclination at which the secondary is too massive to be observed in the spectra, that is, verifying M_V,2 − M_V,1 = 2.5. In this case (13)

The likelihood P(ϵ|γ > γ_±) sums all simulations compatible with ϵ(±10%) divided by the total number of simulations such that γ > γ_±. Finally, the marginal probability P(ϵ) is the sum of all simulations compatible with ϵ(±10%) divided by the total number of simulations.

6.1.4 A few caveats

– To begin with, as was mentioned in the preceding section, we did not account for the true proper motion and the parallax in the simulations. In reality, with an additional accelerated motion unaccounted for in the five-parameter model, the Gaia reduction software could have derived an excess of proper motion and an excess of parallax. The excess parallax modelled is maximal when the orbital motion is aligned with the parallax direction and the orbital period is close to 365 days. Unfortunately, the orientation of the orbits of our targets compared to the parallax direction is generally unknown, so this cannot be properly taken into account. This could lead to underestimation of the photocentre semi-major axis and thus the mass of companions with P ~ 365 days. This concerns five systems with orbital periods within 25% of 365 days. On the other hand, the issue of excess proper motion is only relevant for those sources that are members of the secondary dataset of the DR1, i.e. not members of the TGAS sample. Indeed, for those, the time baseline of the astrometric measurement is not 24 yr but rather <416 days, which is less than the duration of the DR1 campaign. In these cases however, the proper motion can be fitted out easily since it is purely linear. This is done by slightly modifying the system of Eqs. (9) with a moving centroid u_c (t) = (x_c + μ_xt, y_c + μ_yt) and inverting the system solving for four parameters rather than two. We incorporated this correction for nine sources in the secondary dataset: BD+132550, BD+210055, BD+680971, HD 147847, HD 155228, HD 207992, HD 225239, HD 24505, and HD 62923.

– For triple systems, the Keplerian solution of the outer long-period companion being unknown, we could not simulate its effect on the motion of the photocentre. We could only simulate the astrometric excess noise derived by Gaia by assuming that the reflex motion of the primary star is mainly explained by the innermost and better-constrained companion. This could be an incorrect assumption (see below), but still allows us to derive a strict upper limit on the mass of the inner companion.

– Finally, we must warn that the excess noise is intrinsically sensitive to outliers of which there are probably quite a few in the DR1 state of the Gaia processing. Added to the issues of attitude modelling errors, the astrometric excess noise measured by Gaia might be in a few cases overestimated. On the other hand, the level of poor calibration in DR1 could also be underestimated, for the targets with the fewest effective epochs and in cases in which the companion is stellar in nature and contributes light. Nevertheless, we found in general that the value given for the excess noise in the DR1 is relevant and agrees with the χ² published in the DR2 accounting for a longer baseline of astrometric monitoring (see following section).

6.1.5 A comparison with Gaia DR2

We added in Table B.12 the DR2 normalised unit weight errors, with ν the number of DOF, as defined in Lindegren et al. (2018) and accounting for an average renormalization factor 1/ due to the DOF-bug on bright (G < 11) targets. Likeϵ, the RUWE is a measurement of the astrometric scatter around the five-parameter solution. As the DR2 is based on 670 days of astrometric measurements, we would expect the astrometric noise to be of better quality than in the DR1, or at least generally agree with the excess noise measurements in the DR1. Unfortunately, we cannot use the values of RUWE for individual sources to directly apply the GASTON method, since neither the individual renormalization factor nor the typical excess attitude noise are known.

A global comparison of the DR1 excess noise to the DR2 RUWE for the sources of our sample reveals a positive correlation, as shown in Fig. 9. This strengthens the reliability of ϵ_DR1 as a measurement of the astrometric scatter.

Here, we focused only on short-period systems with P < 670 days and on sources of the DR1 that are members of the TGAS sample. For long-period orbits, the astrometric displacement as seen by Gaia in the DR2, with a 670-day baseline, could have been modelled by “instantaneous” proper motion, thus biasing the measurement of astrometric scatter. This also stands for sources of the DR1 that are not members of the TGAS sample, since for those only Gaia measurements were used and the time baseline is 416 days. On the contrary, for sources in the TGAS sample, the time baseline is 24 yr, and for short-period systems the DR2 could not have confused orbital motion with proper motion.

In this figure, the outlier HD 71827 is a triple system, with an inner companion at a period of 15 days and an outer companion with a period of greater than 20 yr. The large discrepancy between the DR2 RUWE and the DR1 excess noise shows that Gaia caught the motion of the star due to the outer companion rather than the inner one. In the DR2, this long-term motion could have been confused with proper motion and thus modelled out, while in the DR1, since HD 71827 is a member of the TGAS sample, the whole motion contributes to the scatter in the excess noise. Despite this issue, we continue to assume that the excess noise of triple systems is due to the inner companion in order to derive an upper limit on its mass.

Fig. 9 Comparison of Gaia DR1 excess noise with Gaia DR2 renormalized unit weight error (see text for explanation). Only systems with a short orbital period (<670 days) and which are part of the TGAS sample are compared. For longer period systems or for those of which HIPPARCOS positioning was not taken into account, the orbital motion could have been modelled out of the DR2 or the DR1 data by fitting the proper motion. The outlying HD 71827 case is discussed in the text.

6.1.6 Results

In Table B.13, we give the 1-σ bounds of semi-major axis of the photocentre and companion mass, as defined in the preceding section. We only consider systems for which a well-constrained Keplerian was derived from the RV data, that is for which the uncertainty on the orbital parameters does not exceed 10%. This led to the rejection of five more stars among the 51 considered in this analysis, HD 13014, HD 146735, HD 40647, HD 60846 and HD 85533, all of which have a companion well within the M-dwarf regime. We therefore only consider now the remaining 46 targets that have a well-constrained RV orbit. Moreover, for 13 out of them, ϵ is less than the typical “normal” ϵ value for a single star in Gaia DR1, which is 0.5 mas. In these cases, the derived masses should only be considered as upper limits.

For 43 out of the 46 systems, we found that the marginal probability of ϵ_DR1 ± 10% being produced by simulations with any inclinations is larger than 0.001. This is a positive sign that the GASTON method produces sensible values of the astrometric excess noise compatible with real Gaia measurements, in more than 90% of the cases. Conversely, for three systems, HD 62923, HD 71827 and HD 156728, the value of ϵ_DR1 was difficult to produce, and required strong fine tuning of the simulated data. The excess noise is either incorrectly estimated by Gaia, one of our assumptions is incorrect, or the astrometry is polluted by one of the companions.

Indeed, HD 71827 is a triple system for which the Gaia’s DR1-DR2 comparison in Sect. 6.1.4 suggested that the astrometry recorded by Gaia was rather due to the outer companion. Our assumption that only the inner companion contributes to the excess noise was therefore certainly incorrect. HD 62923 was shown to be massive in Sect. 5 with a companion that could not be assumed as dark. Finally, the phase curve of HD 156728 is not fully covered by radial velocity measurements with a derived orbital period of ~4100 days. Our solution is possibly inexact, or the value of ϵ_DR1 = 0.48 mas is underestimated.

The case of the triple system HD 71827 is worthy of further comment. For this source, we could not produce many ϵ_s as large as what is measured by Gaia for this system: ϵ_DR1 = 1.56 mas. The largest excess noise simulated is obtained for an inclination of 3.3° leading to M₂ = 0.6 M_⊙. Beyond this limit, the magnitude difference between the companion and the primary must be less than 2.5, implying a secondary peak present in the CCF of HD 71827 spectra, which is not the case. This confirms that the excess noise measured in the DR1 is more likely due to the outer companion rather than the inner one. The inner companion mass M₂ = 0.6 M_⊙ should therefore be interpreted as an upper limit only, and the inclination of the system is likely (much) larger than 3.3°.

Figure 10 summarises the results, including a comparison between the semi-major axis obtained with Gaia and those obtained using the HIPPARCOS data (Sect. 5). We also found ground-based speckle interferometry for HD 106888 that led to a separation of 32 ± 3 mas and a magnitude difference of ΔM_V = 1 between thetwo components of this system (Tokovinin 2014). If we assume that the two components were at apastron, this is equivalent to a semi-major axis of the photocentre of about 1.6 ± 0.15 mas.

The comparison with HIPPARCOS and interferometry is quite satisfying. In all the cases where a significant non-edge-on inclination is measured and a corresponding HIPPARCOS solution is derived, the revised semi-major axis of the photocentre a_ph tends to always be much closer to the HIPPARCOS result than the semi-major axis derived with RV results only. This is emphasised in Fig. 11. Most importantly, the GASTON method always leads to a value of a_ph that is larger than the HIPPARCOS and interferometric measurements. Therefore it looks relatively safe to consider the results of this method as an improved measurement of the inclination and of the true mass of the companion compared to RV fitting only.

We tooka closer look at the BD candidate BD+210055 b, which was well-constrained using HIPPARCOS astrometry to be an M-dwarf, with M₂ = 140–290 M_J at 3σ. We found here that the large value of ϵ_DR1 = 1.3 mas measured for this system led with GASTON to derive a companion mass of 96–110 M_J at 1σ. Although not exactly compatible with the HIPPARCOS result, it is remarkable that we reach the same conclusion concerning the real stellar nature of this object.

We found 12 systems that were not already constrained with HIPPARCOS and for which GASTON led to an inclination significantly different from 90° at 3-σ. These are BD+362641, HD 23965, HD 48679, HD 73636, HD 77712, HD 103913, HD 106888, HD 130396, HD 144286, HD 153376, HD 156111, and HD 217850. Five of these latter are particularly interesting to us since their companion was determined to be in the BD mass regime thanks to radial velocities: HD 23965, HD 48679, HD 77712, HD 130396 and HD 217850. For all five, the simulations were able to produce many ϵ_s values compatible with the DR1 excess noise with P(ϵ_DR1 ± 10%) > 1%. The measurements of their inclination thanks to the GASTON method is thus relatively robust. We can safely state that the companions of HD 77712, HD 130396, and HD 217850 are M-dwarfs with masses well above 80 M_J. We also recall that the M sin i of HD 77712 b was underestimated in Sect. 4.6.1 due to the deformation of the CCFs by a hidden component. HD 77712 b is thus well within the M-dwarf domain. On the other hand, for HD 48679 and HD 23965, while the inclination is significantly different from 90°, the companion mass does not exceed 90 M_J. These are likely to be BDs.

Finally, we measured that the mass of seven companions among the 46 considered here are compatible with the BD regime at 1σ: BD+291539 b, HD 23965 b, HD 28635 b, HD 48679 b, HD 71827 b, HD 82460 b, and HD 211681 b. We discuss these in more detail in Sect. 7.

Although it is not free of possible systematic error, we conclude that the GASTON method is able to derive reliable estimations of system inclination without the use of the definitive Gaia intermediate data. It proves to be a useful method for the characterisation of binary mass and to discard massive companions with short periods in exoplanet RV surveys. We are now applying this method to other catalogues of RV-detected binary stars and exoplanets in order to remove the inclination degeneracy on their M sini measurements, and thus constrain their true masses. This should show that some bodies now considered as exoplanets are actually face-on binaries.

For the largest-period orbits, if virtually nothing can be said using this method, the discrepancy between HIPPARCOS and Gaia DR1 proper motions, theΔQ factor, will be more relevant to these cases, with a time baseline larger than 24 yr. We explore this option in the following section.

Fig. 10 Semi-major axis of the photocentre as measured with Gaia using the GASTON method (blue squares). The circles represent the semi-major axis derived if an inclination of 90° (edge-on system) is assumed. Green points are HIPPARCOS measurements or upper limits, as explained in Sect. 5. The red squares correspond to a measured Gaia excess noise with a marginal probability smaller than 0.001. The cyan point is the result found in the HIPPARCOS double and multiple systems catalogue. The magenta square is the speckle measurement of Tokovinin (2014) for HD 106888.

Fig. 11 Direct comparison of Gaia and HIPPARCOS derivation of the semi-major axis of the photocentre when both are available. The colour code is the same as in Fig. 10. Grey points indicate upper limits derived with HIPPARCOS. The red dotted line represents the equality aHIP = aGaia. The outlyingpoint with aHIP = 32 mas is HD 225239. This case is discussed in Sect. 5.

6.2 The TGAS discrepancy factor ΔQ

While it waspointed out that ΔQ, as produced in DR1, does not take into account the perspective acceleration (Michalik et al. 2014), the stars in our sample are too distant and the proper motions too small for perspective acceleration to be significant. In principle, a value of ΔQ, typically larger than 90% of the Gaia primary sample, that is ΔQ > 10 (Lindegren et al. 2016), could be considered as significant, and we should conclude that a non-zero acceleration is detected, supporting the binarity of the system.

As shown in Table B.12, 19 targets have a value of ΔQ larger than 10, while 7 haveΔQ > 100 and 1 has ΔQ > 1000. The value ofΔQ must be related to the amplitude of the orbital motion, and thus should present a correlation with the semi-major axis of the primary star a₁. Indeed, ΔQ should be more sensitive to large orbital periods (P > 4 yr; the HIPPARCOS baseline) that lead to larger differences between the proper motions measured on a 24-yr baseline and those measured on a 4-yr baseline; while larger companion masses also increase the astrometric acceleration. Because of the degeneracy on the inclination of the systems in RV solutions, only the minimum estimation a₁ sini is known. Figure 12displays therelation between ΔQ and a₁ sini for the present sample, only considering binaries and excluding triples. We find that values of ΔQ larger than 100 are exclusively found for primaries with a semi-major axis of at least ~0.7 au. Moreover, values of a₁ sini greater than 1 au systematically lead to ΔQ > 100. On the other hand, below 1 au the values of ΔQ are scattered uniformly between 0 and 100. We conclude that only values of ΔQ > 100 should be trusted as a positive detection of binarity, with a₁ > 1 au.

The system with ΔQ > 1000 is HD 156728. Its period is larger than 10 yr, and the companion mass stands in the stellar domain above 126 M_J. The primary semi-major axis derived from the RV solution is greater than 0.73 au. The detection of a large value of ΔQ suggests an underestimation of a₁ and of the mass of the companion in this system that is most likely seen nearly face-on. The non-detection in ϵ_DR1, lower than 0.5 mas, could be compatible with this result since it allows the mass to be as large as 250 M_J, as shown in Table B.13.

The other six systems with ΔQ > 100 are HD 108436, HD 13014, HD 153376, HD 60846, HD 71827, and HD 85533. Their long-period companions stand beyond 100 M_J and all have a₁ sini > 1 au. No HIPPARCOS astrometric solution could be derived for any of these sources owing to the short span of the HIPPARCOS measurements, namely smaller than 0.4 orbital periods for all of them. Interestingly, the astrometric excess noise of these six systems is significantly greater than 0.5 mas. However, for four them, the RV orbit is not well constrained, having a large period unknown at more than 10% uncertainty. Moreover, the orbital phase is in most cases not fully spanned by the RV measurements. Only HD 108436 and HD 153376 are well fitted and the application of the GASTON method leads to a re-evaluation of the mass at a larger value, and to a semi-major axis of the photocentre larger than 20 mas. This is approximately the displacement ΔQ able to be detected, since the precision of HIPPARCOS is ~10 mas.

The case of the triple system HD 71827 is interesting, since a clear motion is detected by Gaia, with ϵ_DR1 = 1.6 mas, while ΔQ = 358. This suggests that the motion of the star under the influence of its companions is detected by both indicators. While it remains possible that they do not detect the motion due to the same companion, we found in the preceding section that the DR2 χ² and DR1 excess noise agree if ϵ_DR1 measures the astrometric motion of the long-period outer companion. Undoubtedly, the same motion was measured by ΔQ. With a photometric semi-major axis on the order of 1 mas, the motion due to the inner companion is clearly out of the detection zone of ΔQ.

We conclude that ΔQ is a useful binarity indicator, provided that ΔQ > 100, leading to the detection of a primary star motion with a semi-major axis greater than 1 au. It has however of limited use since it does not allow for the exact mass of the companion to be derived.

Fig. 12 TGAS discrepancy factor ΔQ vs. the a1 sini (in AU), the minimum semi-major axis of the primary star orbit derived from RV. The solid red line indicates the ΔQ = 100 limit, and the dotted black line represents the 1 au limit. Some error bars are smaller than the size of the symbols.

7 Detailsof the seven BD companions

Among the initial sample of 12 BD candidates derived by RV in Sect. 4, we excluded 5 of them using astrometric data of HIPPARCOS and Gaia. HD 210631 b and BD+210055 b were found to be M-dwarfs using the Sahlmann et al. (2011) method on HIPPARCOS data in Sect. 5. Moreover, the masses of HD 130396 b, HD 217850 b, and HD 77712 b could be constrained beyond 90 M_J thanks to the GASTON method applied on Gaia DR1 astrometric excess noise in Sect. 6.

Most importantly, we derived in Sect. 6 that the mass of the seven remaining companions could be constrained below 90 M_J. All seven of these companions are thus likely BDs. We list some details on their detection below.

BD+291539

For this G-type star we observed RV variations compatible with a 60-M_J BD companion on a 176-day orbit at a semi-major axis of 0.6 au. The fit of the 17 RV points is of good quality with a residual dispersion of 4.7 m s⁻¹. Probably owing to the short period, the values of ϵ and ΔQ are too small to indicate any significant astrometric motion in Gaia data. Unsurprisingly it was not detected either by HIPPARCOS. This companion is likely a BD with a maximum mass about 69 M_J.

HD 211681

From the orbital parameters and minimum mass of the secondary companion in this system, M₂ sini = 77.8 ± 2.6 M_J with P = 7612 ± 131 days and a₂ = 8.28 ± 0.16 au, we deduce an inferior limit of 7 mas for the semi-major axis of the astrometric orbit of the primary. Neither HIPPARCOS nor Gaia detect any significant motion. Moreover, the comparison between HIPPARCOS and Gaia astrometry is barely significant with ΔQ = 42, which is not surprising considering the HIPPARCOS precision of about 10 mas. We conclude that given the metallicity of HD 211681 (Fe/H ~ 0.36), its companion is likely an object probing a mass regime between star and BD, around 80 M_J.

HD 23965

Using 84 RV measurements obtained with SOPHIE, we derived for this active F-type star a Keplerian compatible with a 40-M_J BD candidate on an 11-yr orbit at 5 au from the star. The large dispersion of the residuals of ~30 m s⁻¹ is compatible with the strong activity that is measured for this source, with = −4.47. The RV jitter tends to magnify the uncertainties of the derived parameters, but they remain known with a precision better than 10%. However, since the full orbital phase has not been covered yet, the period is still not constrained above 3974 days. The Gaia DR1 astrometry, measuring ϵ = 0.6 mas and an insignificant ΔQ, is suggestive of a system close to edge-on. Applying GASTON on this system, assuming the 3974-day period, leads to an inclination of 76 ± 5° and a companion mass of 42 ± 1 M_J. HD 23965 b is thus a strong BD candidate.

HD 28635

Paulson et al. (2004) argued that the mass of the companion could be significantly higher than that found with RVs (~77 M_J). These latter authors proposed 0.86 ± 0.31 M_⊙. Evidence for a small inclination was concluded from the v sin i ~ 1 km s⁻¹ of the primary compared to the estimation of its true rotation velocity. However, the astrometric data presented here, with ϵ = 0.51 mas and ΔQ = 20, do not tend to confirm this result. The value of ϵ, given the RV solution derived, rather lead to an inclination of 66–80° and a companion mass of 82–88 M_J at 1σ. With a Fe/H ~ 0.16, this companion is located slightly above the BD-M dwarf limit. Nevertheless, the phase coverage with Gaia is only partial along the seven-year orbit. Excess proper motion could tend to lower the value of ϵ that was measured for the DR1.

HD 48679

According to the 26 RV measurements obtained with SOPHIE, the companion of this G0 star is a BD candidate with an M₂ sini ~ 36M_J, an orbital period of 1111 days, and a semi-major axis of 2.1 au. The Keplerian fit is of good quality with a small residual dispersion of 4.6 m s⁻¹. For this star, the astrometric excess noise measured by Gaia of 0.8 mas leads to an inclination of between 41 and 65°, a true companion mass of 41–57 M_J, and a_ph ~ 1.3 mas. Moreover, the small extent of the astrometric motion of the photocentre is compatible with a non-detection from comparing Gaia and HIPPARCOS, with ΔQ = 3. Thus, HD 48679 b is a likely BD companion.

HD 71827

The inclination of this triple system and the true mass of the inner companion (M sini = 26 M_J and P = 15 days) could not be constrained by astrometry. Indeed, it was shown in Sect. 6.1.5 that the astrometric scatter derived in Gaia DR2 and DR1 strongly disagree. Moreover, it was also difficult to model with GASTON an excess noise as large as that measured by Gaia in the DR1 for this system up to inclinations that could not be compatible with a dark companion. This shows that the astrometric scatter measured by both the Gaia DR1 and DR2 for this system cannot be explained by the inner companion, but rather by a long-period object such as the outer M-dwarf companion of HD 71827 (P > 20 yr). It follows that the real mass of HD 71827-b is not constrained and could still be compatible with a mass in the BD domain, although it could also still be more massive.

HD 82460

Using 17 RV measurements obtained with SOPHIE, for this early G-type star we report the detection of a BD candidate with M sini ~ 73-M_J on a 590-day orbit at 1.4 au from the primary. The Keplerian fit is of relatively good quality with medium residual dispersion of 7 m s⁻¹ and orbital elements known with errors smaller than 5%. HIPPARCOS intermediate data do not lead to astrometric motion detection, and only allow an upper-limit on the mass of about 270 M_J to be derived. On the other hand, Gaia measures ϵ = 0.51 mas, which is also barely significant. The Gaia-HIPPARCOS discrepancy factor is not particularly informative, with ΔQ < 10. This is not surprising since the orbital period is short compared to the Gaia-HIPPARCOS baseline of 25 yr. Applying the GASTON method on the value of ϵ leads to a mass of the companion next to the classical hydrogen-burning limit at ~80 M_J. This is mostlikely an upper limit on the mass. Therefore, HD 82460 b is likely a BD.

8 Discussion

The present results allow us to complete the statistics of BD companion candidates around solar-like stars. In the following, we consider the M sin i of BD candidates rather than the exact mass, because considering only companions for which the true mass is derived would introduce bias, favouring the inclusion of transiting edge-on systems and systems closer to the Sun for which astrometric motion is easier to measure. Moreover, many companions detected as exoplanets in RV surveys might actually have a true mass in the BD regime and the derived statistics would miss them. Already published and new companions in the northern sky were compiled in Wilson et al. (2016), Sozzetti & Desidera (2010), and the SB9 catalogue (Pourbaix et al. 2004). We included all companions of these publications, even truly stellar, for which the M sini is within 13.5–90 M_J.

Our initial sample consists in a selection of 2350 targets which have δ > 0°, d < 60 pc, +0.35 < B–V < +1, and are located at less than ± 2 mag from the main sequence (Dalal et al., in prep.). These criterions were also used on the additional previously published data. Table B.15 gives a summary of all additional systems used here. Adding those to the new detections in this work, we obtain the M sin i-period diagram plotted in Fig. 13. This diagram shows a clear lack of detection of companions with periods below ~80 days (0.4 au semi-major axis), for the whole 15–90 M_J mass regime. On the contrary, the detected companions are more uniformly distributed in mass above this limit.

Fig. 13 M sin i–P diagram of all companions reported in this study (blue circles) compared to RV-detected BDs of Table A.1 in Wilson et al. (2016) (purple squares) and giant exoplanets (M sin i > 1MJ) from the exoplanet.eu database (green circles). We selected only objects with δ > 0° and those satisfying the constraint of +0.35 < B–V < +1, and d < 60 pc.

8.1 Brown dwarf frequency and a desert below the 80-day period limit

The work presented here adds a significant number of new BD companion candidates with orbital periods shorter than 10 000 days and M sini within 13.5–90 M_J. The full sample of main sequence FGK stars within 60 pc in the northern sky contains ~2950 stars in the new HIPPARCOS catalogue (van Leeuwen 2007). Apart from the 12 new BD candidates reported in this paper, we collected 32 other BD candidates in the literature that are companions to main sequence FGK stars at less than 60 pc from the Sun in the northern sky. This leads to a minimum of 44 BD candidates among the 2950 systems identified by HIPPARCOS.

The monitoring program for the search of Giant planets with Sophie already gathered more than 3 RV points per star for 2050 of them. About 300 sources still have less than 3 RV points and are still uncharacterised, but this number decreases yearly. With more observing time being devoted to sources with more than three RV points that present interesting variations, less observing time is available to complete the monitoring of a random set of stars.

Inspecting the RV variations of the 2050 systems with at least three RV points, apart from those published in this paper, we found that as many as ~30 more BD candidates with a period up to 10 000 days remain to be characterised. The time span of RV measurements in this sample of 2050 systems ranges from 2 to 4200 days for 99% of them, with a median at 850 days. Among these 30 yet-unconstrained companions, about 15 have a long unconstrained period orbit, for which RV measurements probe a drift-like variation on a baseline of 300 to 5000 days. Their periods and M sini could be considerably larger than these time spans. We therefore estimate that between 0 and at least 30 BD companions within the 2950 main sequence and nearby FGK systems of the northern hemisphere remain to be discovered in addition to the 44 BDs gathered here.

This corresponds to a lower limit on the detection frequency of BD candidates within M sini = 13.5–90 M_J and with orbital periods of less than 10 000 days of 2.0 ± 0.5%. This value remains compatible with the upper limit of this frequency obtained by Guenther et al. (2005) with f_BD < 2% for BD companions in the Hyades cluster with a semi-major axis (sma) <8 AU, but larger than the estimation of Sahlmann et al. (2011), f_BD = 1.3% for candidates with sma < 10 AU.

Being obtained from RV and M sin i only, this frequency is overestimated due to the uncertainty on the inclination of the system. Sahlmann et al. (2011) proposed a correction to this number by only considering companions that were not found to be real M-dwarf using astrometry, which decreases this determination to f_BD,corr < 0.6%. This compared well to the rate estimation <0.5% of Marcy & Butler (2000) and the one derived by Santerne et al. (2016) of 0.29 ± 0.17% for transiting BDs with orbital periods smaller than 400 days observed with Kepler (Borucki et al. 2010). Applying the same procedure as in Sahlmann et al. (2011), we find that 35 companions of the 44 considered here are compatible with the BD domain. This leads to a revised lower-limit on the BD frequency of f_BD,corr > 1.7 ± 0.5%, which remains higher than expected. Still, this number is most likely overestimated since 25 of the 32 additional systems in Table B.15 have not yet been constrained by astrometry. We are now applying the GASTON method systematically on all systems with BD candidates to constrain their mass.

Below an 80-day period, we find that the detection frequency drops significantly, with only six BD candidates found in this region. Within the 2050 systems mentioned above, we found only two possibly missing BDs with periods less than 80 days. This leads to a lower-limit f_{BD, low} ~ 0.24 ± 0.04%, a factor of eight lower than above 80 days. This is in line with the findings of Ma & Ge (2014) that BD companions tend to avoid a mass–period region bounded by the limits 30–60 M_J and P < 100 days. Guillot et al. (2014) showed that this scarcity of detections below the 100-day period limit, especially for G-dwarfs, can be explained by the dynamical interactions between the stars and close-in companions, including tidal interactions, stellar evolution, magnetic braking of stars, and tidal dissipation by gravity waves.

8.2 The companion mass distribution beyond the 80-day period limit

Studying the M sin i-histogram of detections beyond the 80-day period limit in more detail, Fig. 14 shows that the M sini distribution is possibly not exactly uniform, with a decrease in detection rate below ~50 M_J. This could be the sign of a lower bound in the true companion mass distribution that was suggested by the work of Halbwachs et al. (2000). We therefore tried to reproduce this M sini distribution out of simulated companions whose masses are uniformly distributed from some lower-bound up until 0.52 M_⊙. We introduced several different lower bounds on the mass from 15 to 100 M_J, with a uniform density function above that limit and zero below. A random inclination is assigned to all these simulated objects, from which we can deduce M sini. We compare these to the 38 detections between 15 and 90 M_J and with P > 80 days. We randomly selected 38 simulated values of M sin i and performed a two-sample Kolmogorov-Smirnov test between the simulated sample and the observed one. This simulation was performed 1000 times. We then counted the number of simulated samples that were incompatible with the observed data, assuming the null-hypothesis rejection for a p-value < 0.05. Two samples drawn from the same distribution should reject the null-hypothesis 5% of the time on average. If a mass distribution leads to a good modelling of the detection statistics, then about 95% of the simulated samples out of this mass distribution should accept the null-hypothesis. We obtain that any uniform mass-distribution with a cut-off between 15 and 100 M_J leads to compatibility with observations for more than 95% of the simulations. The M sini histogram and cumulative plot of detections are plotted in Figs. 14 and 15 and compared to the best of the simulations with the 15-M_J lower bound model.

Conversely to the intuition based on the M sin i distribution, the actual mass distribution of BDs beyond 0.4 au might thus be uniform all the way down to 15 M_J, as is found at wider separation (Reid et al. 2002). We do not find evidence of a lower mass limit in the BD companion population, but cannot exclude the existence of such a bound in the mass distribution (Halbwachs et al. 2000; Sahlmann et al. 2011).

Moreover, a second distribution coming from lower masses, that is the massive planets, does not appear necessary to explain the actual M sin i distribution beyond the 80-day period limit. This suggests that the mass distribution of massive planets does not spread significantly within the BD domain. The tail of the distribution probably stops around 20 M_J but not significantly beyond this value. Independent direct imaging surveys of long-period BDs (Brandt et al. 2014) reached similar conclusions, with low-mass BDs, even those below the deuterium-burning limit, more likely arising from gravitational collapse in discs or fragmenting cloud, as in more massive objects.

Fig. 14 Histogram of the detections with M sin i between 15 and 80 MJ and P > 80 days. The bins are logarithmically spaced. This is compared, in red, to the best-fitting distribution obtained by simulating the effect of the inclination on a random population of objects uniformly distributed in mass from 15 MJ to 0.52 M⊙.

Fig. 15 Cumulative distribution of detections with M sin i between 15 and 80 MJ and P > 80 days. The red curve corresponds to the best-fitting distribution plotted in Fig. 14.

9 Summary and conclusions

We report here the detection of 54 companions to FGK stars in the neighbourhood of the Sun using RV observations of the SOPHIE spectrograph. Among them, 12 were detected as BD candidates according to their projected mass M sini, and 42 as M-dwarfs.

Using HIPPARCOS and Gaia, we were able to study the values of the mass derived for several of the companions, as summarised in Table B.14. We introduced a new method, GASTON, to derive inclination by combining the astrometric excess noise published in the Gaia DR1 and the Keplerian orbit derived from RV variations. This allowed us to reconsider the mass of the 12 BD candidates derived thanks to SOPHIE radial velocities. We found that five of them actually stand in the M-dwarf regime, and confirmed that the seven remaining companions are likely BDs. The stars BD+291539 b, HD 23965 b, HD 48679 b, and HD 82460 b are strongly constrained below 90 M_J. While HD 28635 b, HD 211681 b, and HD 71827 b are possible BDs, they remain as candidates. Moreover, we obtained a stellar mass for the companion of HD 210631 that was previously published as a BD candidate by Latham et al. (2002).

Our 12 BDs detected with RV, added to those reviewed and discovered by Wilson et al. (2016), those in the SB9 (Pourbaix et al. 2004), and those found by Sozzetti & Desidera (2010), place a strong lower limit on the orbital period of BDs at 80 days. The short-period region below 80 days appears to be four times less populated than at a larger period. On the other hand, above 80 days the M sini detection density is well reproduced by a flat distribution of mass that goes down to as low as 15 M_J. Moreover, a long extension of the massive planet distribution beyond 20 M_J does not appear necessary to reproduce the detected M sin i of all companions in the BD regime.

These conclusions should however not be understood as definitive. The statistics are not yet completely free of bias, since the monitoring of the volume-limited sample is not complete. However, we estimate that only a few short-period radial-velocity BDs could have been missed in our sample, and a conservative number of up to 30 companions could still be missing in the 15–90 M_J regime at orbital periods shorter than 10 000 days.

As was demonstrated in Sahlmann et al. (2011), accounting for the exact mass rather than minimum mass can strongly change the picture. Moreover, we showed that only being based on M sini cannot lead to a firm conclusion on the presence or absence of a lower bound on mass for BDs. Only a systematic search for constraints on both mass and inclination for every system with a candidate BD will lead to the ultimate unbiased mass distribution of BDs.

In this work, we showed that such a goal will be achieved thanks to the combination of Gaia and high-resolution spectroscopy. By 2022, it is likely that acceleration solutions and actual orbital solutions for many of the detected companions presented in this work (and many new BD and stellar companions) will become available in the Gaia DR3. In the meantime, the GASTON method presented here allows the already published data release of Gaia to be used to constrain the inclination and mass of many systems and companions, including planets, BDs, and binary stars.

Acknowledgements

We are thankful to the anonymous referee for their constructive comments that helped to significantly improve the quality of this paper. We thank all the staff of Haute-Provence Observatory for their support at the 1.93 m telescope and on SOPHIE. This work is based on observations made with SOPHIE in the context of the programme “Recherche de Planétes Extrasolaire” (PI: I. Boisse, F. Bouchy) and the Programme National de Physique Stellaire (PNPS), in particular the programme "Exploring the Brown Dwarf Desert around FGK stars in the Solar neighborhood" (PI: F. Kiefer; IDs: 2017B_PNPS008, 2018A_PNPS008, 2018B_PNPS013). F.K. acknowledges support by a fellowship grant from the Centre National d’Etude Spatiale (CNES). This work was also supported by the PNPS of CNRS/INSU co-funded by CEA and CNES. A.V-M. and F.K. thank the CNES for financial support. P.A.W. acknowledge support from the European Research Council under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 694 513. V.B. and J.H. acknowledge support by the Swiss National Science Foundation (SNSF) in the frame of the National Centre for Competence in Research PlanetS. V.B. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Four Aces; grant agreement No 724427). This work was supported by FCT – Fundaço para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalizaço by these grants: UID/FIS/04434/2013, UID/FIS/04434/2019 & POCI-01-0145-FEDER-007672; PTDC/FIS-AST /28953/2017 & POCI-01-0145-FEDER-028953 and PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113. N.A-D. acknowledges the support of FONDECYT project 3 180 063. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

Appendix A Photocentre semi-major axis

To fix the value of the photocentre semi-major axis, we use the following formula (van de Kamp 1975): (A.1)

The period P is expressed in years, the masses are given in M_⊙, the parallax ϖ unit is mas. Also, β is the luminosity fraction β = L₂∕(L₁ + L₂) and B the mass fraction B = q∕(1 + q). The mass ratio q can be found by solving the equation of the mass function (see e.g. Halbwachs et al. 2014) (A.2)

with M_RV = . The solutionof this equation, fixing q₀ = M_RV∕sin(i) M_⋆, is (A.3)

Knowing the primary mass M₁ from Table B.3 and deducing M₂ from q, the luminosity fraction at optical wavelength is derived by (A.4)

The visual luminosity fraction can then be calculated thanks to the empirical relation existing between absolute visual magnitude and star mass (Kroupa et al. 1993) which can be approached by (A.5)

Since no secondary peaks were seen in any of the CCFs for all targets, we assumed in the simulations that M_V,2 − M_V,1 had to be greater than 2.5, and thus β < 0.1. This led us to discard certain values of I_c implying overly large values of the luminosity fraction.

Appendix B Additional tables

Appendix C Additional figures

Fig. C.1 Orbital solutions of the RV variations of the 11 binary systems with a BD candidateas secondary. Radial velocity vs. time are presented in the left panel, with O–C residuals below, and RV vs. phase in the right panel. Data point colours are fixed with respect to the order given on the upper-right corner of each figure, with alternatively red, blue, yellow, and pink. Diminutives “SP” stands for SOPHIE-plus, “Soph” for SOPHIE and “Elo” for Elodie. “Keck” and “SB9” are self-explanatory and references of the corresponding data can be found in Table B.2.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.2 Drift+Keplerian orbital solutions of the RV variations of the three triple systems. Radial velocity vs. time is presented in the left panel, with O–C residuals below, and RV vs. phase for the Keplerian solution in the right panel. “SPEloFix” means that the SOPHIE+ and Elodie datasets offset is fixed to 0 m s−1. “SPSophFix” means that the SOPHIE+ and SOPHIE datasets offset is fixed to 0 m s−1. Details are given in Sect. 4.6.2.

Fig. C.3 Orbital solutions of the RV variations of the 39 stars with an M-dwarf companion. Radial velocity vs. time are presented on the left panel, with O–C residuals below, and RV vs. phase in the right panel. The colour code is explained in the caption of Fig. C.1. “SPSophFix” means that the SOPHIE+ and SOPHIE datasets offset is fixed to 0 m s−1.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.3 continued.

Fig. C.4 HIPPARCOS astrometry 3-σ detections in nine systems. Top panels: modelled astrometric orbits projected on the sky. North is up and east is left. The solid red line shows the model orbit and open circles mark the individual HIPPARCOS measurements. Bottom panels: O–C residuals for the normal points of the orbital solution (filled blue circles) and of the five-parameter model without companion (open squares).

Fig. C.5 HIPPARCOS astrometry 2-σ detections in five systems. Top panels: modelled astrometric orbits projected on the sky, cf. Fig. C.4.

References

All Tables

All Figures

	Fig. 1 P-value of the χ2 -test vs. the Pearson correlation coefficient, R, for both FWHM (orange) and BIS (blue) indicators. The radii of the symbols are linearly scaled with the number of points considered. The dashed lines show the limits of significance for R (±0.5) and for the p-value at 3σ (<0.01).
In the text

	Fig. 2 FWHM (orange) and bissector span (blue) variations of the cross-correlation function of HD 77712 spectra. The error bars are calculated from the RV uncertainties, multiplied by the factors given in Eq. (2).
In the text

	Fig. 3 Upper panel: O–C distribution of SOPHIE (blue) and SOPHIE+ (orange) datasets. Lower panel: RV offset distribution between SOPHIE and SOPHIE+ datasets.
In the text

	Fig. 4 (in m s−1) vs. for SOPHIE (in blue) and SOPHIE+ (in orange) datasets. The symbol size is proportional to B− V with values between 0.5 and 1.2. The black lines represent the relation derived by Santos et al. (2000) for G-type stars with a fit uncertainty of 0.18 dex. The red lines represent the relation derived from the datasets of this work, with a fit uncertainty of 0.3 dex.
In the text

	Fig. 5 Period vs. M sin i for the systems studied in this paper. The crosses show the new results of this paper. The 13 and 90 MJ limits are drawn as dotted lines. The plain white symbols represent stars with two companions. HD 71827 b and c are represented as green triangles, HD 212735 b and c as red stars, and BD+212816 b and c as blue squares. Most error bars are smaller than the symbols.
In the text

	Fig. 6 Period vs. eccentricity for the systems studied in this paper. The blue circles show the new results of this paper, while the green circles represent the Exoplanet.eu database with primaries verifying that δ > 0°, +0.35 < B–V <+1, d < 60 pc. The symbol size is proportional to the logarithm of M sin i.
In the text

Fig. 7 Left: inclination vs. simulated excess noise for BD+192536. The red lines mark the value of for this star ± 10%. Top-right: example of a simulation of Gaia measurements (blue points) for the peculiar case of BD+192536 that fits the value of ϵDR1 = 1.02 mas. The true proper motion and true parallax are assumed to be subtracted. The purple line is the residual proper motion (moving centroid) fitted to the simulated measurements and the red cross marks the true centre of gravity of the system. Bottom-right: residuals with respect to the moving centroid.

In the text

	Fig. 8 Same as Fig. 7 but for HD 71827. Here the value of ϵDR1 is larger than any of the simulations provided that the secondary companion is dark (i.e. not visible as a secondary peak in the CCF). The assumption that the inner companion is responsible for the astrometric motion of HD 71827 is likely to be incorrect (see text for explanation).
In the text

Fig. 9 Comparison of Gaia DR1 excess noise with Gaia DR2 renormalized unit weight error (see text for explanation). Only systems with a short orbital period (<670 days) and which are part of the TGAS sample are compared. For longer period systems or for those of which HIPPARCOS positioning was not taken into account, the orbital motion could have been modelled out of the DR2 or the DR1 data by fitting the proper motion. The outlying HD 71827 case is discussed in the text.

In the text

Fig. 10 Semi-major axis of the photocentre as measured with Gaia using the GASTON method (blue squares). The circles represent the semi-major axis derived if an inclination of 90° (edge-on system) is assumed. Green points are HIPPARCOS measurements or upper limits, as explained in Sect. 5. The red squares correspond to a measured Gaia excess noise with a marginal probability smaller than 0.001. The cyan point is the result found in the HIPPARCOS double and multiple systems catalogue. The magenta square is the speckle measurement of Tokovinin (2014) for HD 106888.

In the text

	Fig. 11 Direct comparison of Gaia and HIPPARCOS derivation of the semi-major axis of the photocentre when both are available. The colour code is the same as in Fig. 10. Grey points indicate upper limits derived with HIPPARCOS. The red dotted line represents the equality aHIP = aGaia. The outlyingpoint with aHIP = 32 mas is HD 225239. This case is discussed in Sect. 5.
In the text

	Fig. 12 TGAS discrepancy factor ΔQ vs. the a1 sini (in AU), the minimum semi-major axis of the primary star orbit derived from RV. The solid red line indicates the ΔQ = 100 limit, and the dotted black line represents the 1 au limit. Some error bars are smaller than the size of the symbols.
In the text

	Fig. 13 M sin i–P diagram of all companions reported in this study (blue circles) compared to RV-detected BDs of Table A.1 in Wilson et al. (2016) (purple squares) and giant exoplanets (M sin i > 1MJ) from the exoplanet.eu database (green circles). We selected only objects with δ > 0° and those satisfying the constraint of +0.35 < B–V < +1, and d < 60 pc.
In the text

	Fig. 14 Histogram of the detections with M sin i between 15 and 80 MJ and P > 80 days. The bins are logarithmically spaced. This is compared, in red, to the best-fitting distribution obtained by simulating the effect of the inclination on a random population of objects uniformly distributed in mass from 15 MJ to 0.52 M⊙.
In the text

	Fig. 15 Cumulative distribution of detections with M sin i between 15 and 80 MJ and P > 80 days. The red curve corresponds to the best-fitting distribution plotted in Fig. 14.
In the text

Fig. C.1 Orbital solutions of the RV variations of the 11 binary systems with a BD candidateas secondary. Radial velocity vs. time are presented in the left panel, with O–C residuals below, and RV vs. phase in the right panel. Data point colours are fixed with respect to the order given on the upper-right corner of each figure, with alternatively red, blue, yellow, and pink. Diminutives “SP” stands for SOPHIE-plus, “Soph” for SOPHIE and “Elo” for Elodie. “Keck” and “SB9” are self-explanatory and references of the corresponding data can be found in Table B.2.

In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

Fig. C.2 Drift+Keplerian orbital solutions of the RV variations of the three triple systems. Radial velocity vs. time is presented in the left panel, with O–C residuals below, and RV vs. phase for the Keplerian solution in the right panel. “SPEloFix” means that the SOPHIE+ and Elodie datasets offset is fixed to 0 m s−1. “SPSophFix” means that the SOPHIE+ and SOPHIE datasets offset is fixed to 0 m s−1. Details are given in Sect. 4.6.2.

In the text

	Fig. C.3 Orbital solutions of the RV variations of the 39 stars with an M-dwarf companion. Radial velocity vs. time are presented on the left panel, with O–C residuals below, and RV vs. phase in the right panel. The colour code is explained in the caption of Fig. C.1. “SPSophFix” means that the SOPHIE+ and SOPHIE datasets offset is fixed to 0 m s−1.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

	Fig. C.3 continued.
In the text

Fig. C.4 HIPPARCOS astrometry 3-σ detections in nine systems. Top panels: modelled astrometric orbits projected on the sky. North is up and east is left. The solid red line shows the model orbit and open circles mark the individual HIPPARCOS measurements. Bottom panels: O–C residuals for the normal points of the orbital solution (filled blue circles) and of the five-parameter model without companion (open squares).

In the text

	Fig. C.5 HIPPARCOS astrometry 2-σ detections in five systems. Top panels: modelled astrometric orbits projected on the sky, cf. Fig. C.4.
In the text