The detection and characterisation of 54 massive companions with the SOPHIE spectrograph -- 7 new brown dwarfs and constraints on the BD desert

Brown-dwarfs are substellar objects with masses intermediate between planets and stars within about 13-80Mjup. While isolated BDs are most likely produced by gravitational collapse in molecular clouds down to masses of a few Mjup, a non-negligible fraction of low-mass companions might be formed through the planet formation channel in protoplanetary disks. The upper mass limit of objects formed within disks is still observationally unknown, the main reason being the strong dearth of BD companions at orbital periods shorter than 10 years, a.k.a. the BD desert. To address this question, we aim at determining the best statistics of secondary companions within the 10-100Mjup range and within 10 au from the primary star, while minimising observational bias. We made an extensive use of the RV surveys of FGK stars below 60pc distance to the Sun and in the northern hemisphere performed with the SOPHIE spectrograph at the Observatoire de Haute-Provence. We derived the Keplerian solutions of the RV variations of 54 sources. Public astrometric data of the Hipparcos and Gaia missions allowed constraining the mass of the companion for most sources. We introduce GASTON, a new method to derive inclination combining RVs Keplerian and astrometric excess noise from Gaia DR1. We report the discovery of 12 new BD candidates. For 5 of them, additional astrometric data led to revise their mass in the M-dwarf regime. Among the 7 remaining objects, 4 are confirmed BD companions, and 3 others are likely in this mass regime. We also report the detection of 42 M-dwarfs within 90Mjup-0.52Msun. The resulting Msin(i)-P distribution of BD candidates shows a clear drop in the detection rate below 80-day orbital period. Above that limit, the BD desert reveals rather wet, with a uniform distribution of the Msin(i). We derive a minimum BD-detection frequency around Solar-like stars of 2.0+/-0.5%.


Introduction
According to the classical convention, brown dwarfs are substellar objects whose mass is too small to maintain hydrostatic equilibrium thanks to hydrogen-based nuclear reactions, while massive enough to ignite Deuterium nuclear reactions in the e-mail: flavien.kiefer@iap.fr core, at least for few million years. Following this definition, the brown-dwarf domain is framed within the mass range 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 many debates (in e.g. Saumon et al., 1996, Chabrier & Baraffe 2000, Luhman et al. 2007, Luhman 2012, Chabrier et al. 2014) out of which it is Article number, page 1 of 51 arXiv:1909.00739v1 [astro-ph.SR] 2 Sep 2019 A&A proofs: manuscript no. article_BD_accepted proposed that the existence of nuclear reactions in its core is not the crucial parameter to define the nature of a substellar body.
The observation of objects with masses as low as 5 M J in young stellar clusters is a strong evidence that molecular cloud fragmentation is not limited in mass and can form objects in the brown dwarf and giant planet mass regime (see e.g. de Marchi, Paresce & Portegies Zwart 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 brown dwarf domain (see e.g. Burgasser et al. 2007 and reference therein). On the other hand, a dearth of detections of brown dwarf companions with orbital periods shorter than 10 yrs, 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 a few tens of Jupiter masses.
Planet formation pathways, such as disk instability or core accretion, could in principle allow to form bodies up to 40 M J within protoplanetary disks (Pollack 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 thus help constraining 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 the recent years have led to abundant detections of brown dwarf companions with diverse instrumental methods (Sozzetti & Desidera 2010, Sahlmann et al. 2011, Diaz et al. 2012, Ranc et al. 2015, Wilson et al. 2016. A difficulty arises due to the dearth of brown dwarf companion detected at short orbital periods, i.e. the brown dwarf desert. For some reason, the presence of substellar companions is forbidden at close distance of a more massive primary star. This implies that the mass-period distribution of brown dwarf companions to sun-like stars are deformed by possibly several perturbing effects, such as tidal interactions, magnetic braking and tidal dissipation (Guillot et al. 2014). This strongly bias the determination of the real mass distribution of giant planets and very low-mass stars.
It is thus necessary to constrain the minimum orbital period above which this effect becomes negligible. Mixing the results from several surveys done with diverse detection methods, Ma & Ge et al. (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 BD 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 tends to confirm the framing of the desert at periods lower than 100 days. But the use of detections arising from several diversely biased or incomplete surveys is perilous. To our knowledge, there exists no fully complete non-biased statistical sample of detected brown-dwarfs companions.
It would be most valuable to achieve a survey of browndwarf companions that is non-biased, or at least for which the selection function of the followed-up sample is well known and allow deriving a meaningful statistics of BD population. Some of the most problematic issues with gathering detections from multiple surveys, apart from instrumental bias, are the diverse observers own interests. Observations are usually stopped as soon as the followed-up target is not anymore of interest regarding the given study. Typically, on one side, sources with a companion that is not within the planetary mass domain, beyond about 20 M J , are not continued and not always published. On the other side, orbits and mass ratio of obvious stellar binaries are easily characterised and published. It follows that brown dwarfs within the BD desert and especially at periods larger than 1 year are under-sampled.
The volume-limited FGK stars survey program for searching 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. in prep). To this date, around 2050 sources were observed at least 3 epochs each, with an aimed signal-to-noise ratio (SNR) per spectrum of at least 50.
The reflex-motion due to brown-dwarfs within the desert leads to RV amplitudes larger than 100 m s −1 . Since the SOPHIE spectrograph is able to detect RV signals as low as a few m s −1 (Courcol et al. 2015) on a time baseline of 13 years, brown dwarfs companions can easily be detected around nearby bright stars. Therefore, we expect to reach eventually almost 100 % completion of RV-detected brown-dwarf candidates with orbital period less than 10,000 days around FGK stars in this volumelimited sample of stars, which 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 brown dwarf desert with SOPHIE, we present here the latest results of this radial velocity survey on 54 solar-like sources with spectral types ranging from K5 to F5. With RV only, we report 12 BD candidates with M sin i within 15-90 M J , among which 8 never published.
We conservatively extend the BD-domain above 80 M J in order to include objects in the grey zone 80-90 M J , separating Mdwarf from brown dwarfs. We think this is justified for essentially three reasons. First, there is always an uncertainty (up to few M J ) on the M sin i derived with RV. Second, the mass limit for Hydrogen-burning is not a strict one, and may vary according to e.g. metallicity from 83 to 75 M J within M/H∼[-1 ; 0] (Chabrier & Baraffe 1997). And third, extending towards lowmass M-dwarfs allows exploring the tail of the BD-mass distribution on the stellar side.
Although velocimetry is an efficient mean to detect companions, either stellar, sub-stellar or planetary, to stars, 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 determine 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, Gaia's intermediate data 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 Section 2 we present the target selection. In Section 3 we review the observations performed and the targets observed. In Section 4, the spectroscopic analysis of the SOPHIE's observations is discussed, including the result of Keplerian fitting to the radial velocity variations. In Sections 5 and 6 we study the astrometric measurements made with Hipparcos and Gaia. In Section 7 we review the 7 discovered brown dwarfs. And finally, in Section 8, we discuss the implication of the presented results on the brown dwarf desert localisation. We conclude in Section 9.

Target selection
The goal of the programme in which this study takes place, is to complete a meaningful unbiased statistic of companions detected within and about the brown-dwarf mass regime, and up to 10 yrs period. Extracting brown dwarf candidates out of a sample of stars which selection function 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 sin i>15 M J , on a time span larger than 10 years. This could allow the determination of the orbit of BD companions with period as large as 10 years. In order to gather the largest possible number of browndwarfs in the BD desert, and to be able to compare the browndwarf population to the low-mass star population, we were especially focused on sources with companion masses in the broad M sin i range of 20-150 M J . This range includes the whole BD regime, from the upper end of the giant planets domain, but also extends up to the late M-dwarfs domain. We thus continued the radial velocity monitoring of the sources that present any clue 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 at a signal-to-noise ratio (SNR) per spectrum of at least 30. Table B.1 summarises the basic informations on the 54 targets covered by the present study. It includes only sources for which we gathered more than 6 RV data and for which a meaningful Keplerian solution of the RV variations, or a lower mass limit beyond 150 M J , could be derived. We excluded SB2 sources from this publication.

Observations
The observations were performed with the SOPHIE spectrograph, fiber-fed from the Cassegrain focus of the 1.93-m telescope at the Haute-Provence Observatory (OHP, France). It 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 highresolution 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 fiber, the instrument is also exposed to the background sky in a second fiber allowing subtraction of scattered light contamination in the science spectrum. The exposure time was varied to reach a signal-to-noise ratio (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 lines removal, sky spectrum removal, CTI correction, CCF computation, and barycentric earth radial velocity correction. In the reduction software, the CCFs are fitted by Gaussians to calculate the sources radial velocities (Baranne et al. 1996;Pepe et al. 2002). Moreover, the bisector spans (BIS) and full-width-athalf-maximum (FWHM) of each CCF are computed following Queloz et al. (2001). We did not correct the seasonal RV zero point variation, from standard stars variation (see e.g. Courcol et al. 2015) that are being tiny (∼m/s) compared to the expected velocity variations amplitude (∼km/s).
In June 2011, the SOPHIE spectrograph hexagonal fibers 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 −1 in the measured velocities was observed on standard stars (Bouchy et al. 2013). We therefore separated the data about June 2011 (JD 2455731.5). Before that date, the data will be referred to as SOPHIE, and after that date, they will be referred to as SOPHIE+. Moreover a systematic noise of 5 m s −1 was quadratically added to the measured RV uncertainty of the SOPHIE data before June 2011 (Hebrard et al. 2016).
Additional non-SOPHIE data were found in the literature, with occasionally, already 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 companions parameters. Some of the public data were found in the SB9 catalogue 1 (Pourbaix et al. 2004).
Data points with less than half of the median S/N and large uncertainty on the radial velocity 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 years, with a median at 8 years.

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, which number is given in 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 onedimensional Kurucz model atmospheres. All derived parameters are given in Table B.3. We estimated the stellar mass and radius of the primary star using the Torres et al. (2010) empirical relation. The log(g) were corrected in order to be calibrated on log(g) derived using asterosismology, following equation (4) in Mortier et al. (2014): The host stars presented in this paper are of type K5 to F8 (from the SIMBAD catalog) on the main sequence with metallicities [Fe/H] ranging from −0.3 to +0.3 dex. HD24505, HD109157 and HD204613 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 HD204613 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 appended to the table the average activity indicator log R HK calculated using all spectra of each target, with the SO-PHIE reduction software (Boisse et al. 2010). 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 log R HK in SOPHIE spectra (Boisse et al. 2010). Our targets show medium stellar activity level in general, with 17 targets having log R HK <−4.75, 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 log R HK >−4.5. Nevertheless, the amplitude of the derived RVs are all larger than a few hundred m s −1 , while activity is expected to influence RV measurements only at the scale of a few tens of m s −1 (Campbell et al. 1991, Saar & Donahue 1997, Saar et al. 1998, Santos et al. 2000, Boisse et al. 2010. All the detections presented in this paper are securely those of true companions but the magnetic jitter of the most active stars will add scatter and imply larger uncertainty to the measurement of orbital parameters for the hosted companions.

Search for activity and binarity indicators in the CCF
We calculated the full width at half maximum (FWHM), the bissector span (BIS) and the signal-to-noise ratio (SNR) variations for all sources. This allows us to check whether the radial velocity 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 SNR 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 the 2 indicators being the variations of the FWHM and of the bissector span (Santerne et al. 2014, Santerne et al. 2015. For each indicator, we performed a χ 2 -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, 8 with less than 4 spectra were not analysed regarding these diagnostics, since any variations would hardly be meaningful.
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) and Santerne et al. (2015). These multiplications factors can be refined here, comparing the scatter of the FWHM or BIS to the median RV uncertainty for every sources. The median factors found lead to Thus for the BIS we confirm the result of Santerne et al. (2015). For the FWHM, we found that the multiplication factor rather stands higher. These corrected factors were used eventually to calculate the χ 2 test and the Pearson correlation coefficient that are summarised in Table B.4 and presented in Figure 1.  Only 10 systems show an FWHM or BIS dispersion that is significant with p-values lower than the 3σ limit. But a single source shows also a strong correlation coefficient of FWHM with RV variations, R(FWHM, RV)=-0.82. This system is HD77712, a K-type star with a medium activity level at log R HK ∼−4.7. On the other hand, it presents no significant variations of the bissector. This is similar to the case of a triple system studied in Section 2.9 of Santerne et al. (2015) with the pollution of the CCF from a weak secondary peak always present at fixed radial velocity. A possible explanation is therefore that HD77712 is a triple system with a long period binary which secondary is polluting the CCF and a shorter period binary with a dark companion. The RV amplitude and M sin i we report for this system are likely underestimated.

Keplerian orbits fitting
The yorbit software (Segransan et al. 2011) was used to calculate the solution by genetic algorithm refining of initial pa-Article number, page 4 of 51 F. Kiefer et al.: 7 new brown dwarfs rameters for a Levenberg-Marquardt optimisation. This leads to priors for an MCMC estimation of errorbars following Díaz et al. (2014Díaz et al. ( ,2016. MCMC was applied on 1,000 iterations. The varied parameters are the period P, the RV amplitude K, the eccentricity e, the angle of periastron ω, the periastron passage time T 0 , and the offsets γ S and γ S,+ for SOPHIE and SOPHIE+ datasets respectively. Specific additional offsets are used for supplementary datasets as indicated in Table B.2. When the RV errorbars 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 parameters 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 errorbars defined by the confidence interval (CI) at 68.3% around the best-value are barely asymmetric, while the difference between median and best-value is not significant. We found that the standard deviation gives more conservative uncertainties than the CI at 68.3%. We thus uses the standard deviation as errorbar in order to keep on the conservative side, especially for cases with inaccurate derivation of the orbital parameters. For incompletely covered orbits the MCMC distributions have non-Gaussian tails: the interval for a confidence interval equivalent to 3-σ is almost certainly not just 3 times broader than that for 1-σ.
For 7 stars, HD5470, HD7747, HD153376, HD193554, HD207992, HD212735 and BD+212816, SOPHIE and SO-PHIE+ data are not sufficiently numerous separately to derive a meaningful solution with ∆γ=γ S,+ − γ S on the order of 50 m s −1 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 divide subsequently in 2 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 details in Section 4.6.1 below. The results for the 3 triple systems are presented in Table B.8 and in Section 4.6.2. They include 1 BD candidate and 2 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 center-of-mass offset γ of each additional dataset are given in Table B In general, the fits are accurate with precision on the orbital elements better than 7% 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 HD85533, 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 HD13014, as well as HD40647, HD60846 and HD146735, 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 lies below 10 m s −1 except for a few active sources that show much larger dispersion of residuals close to 40 m s −1 . We discuss the distribution of the residuals in more details, especially comparing the SOPHIE and SOPHIE+ datasets below in Section 4.4, and confronting to activity indices for the observed sources in Section 4.5.

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 comparing 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 −1 , 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 2 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 Figure 3.
The distribution of RV offsets between SOPHIE and SO-PHIE+ is centred on This is compatible with the results found by Bouchy et al. (2013) that bounds the RV shift due to the upgrade to 0-50 m s −1 . To obtain this distribution, we assumed that the offset should not exceed 100 m s −1 , in which case the offset should be better explained by a slow drift, due to third companion. This led to consider a few systems as rather multiple than binary, as shown in Section 4.6.2 above.
The core of the distribution of O-C residuals standard deviation leads to the following estimation of RV accuracy for SO-PHIE and SOPHIE+ configurations: 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 −1 is derived for observations before June 2011, and of about 3.5 m s −1 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, the absence of instrumental drift (∼m s −1 ) corrections, as derived in Courcol et al. (2015), in the present study. Moreover, our observed sample also includes spectra with SNRs down to 30, while the Hébrard et al. (2016) sample only includes spectra with SNR>50.
Comparing the O-C individually for every targets 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 −1 . If 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 HD23965, HD40647, and HD161479, admit large O-C dispersion >14 m s −1 in both datasets. This is however most likely explained by their significant activity index log R HK >−4.5 (Table B.3). We can thus exclude that more planetary signals are hidden in any residuals beyond an amplitude of about 8 m s −1 .

Residuals dispersion and magnetic activity
The residuals of the Keplerian fit, σ v , can be compared to the activity index log R HK to verify if it can explains the amplitude of the residuals, in general and in specific cases where it is exceptionally large. Previous study of the correlation between RV dispersion and magnetic activity was done in Saar et al. (1998) and in Santos et al. (2000). We followed here the same procedure in order to make a point comparison. We first exclude datasets with less than 7 pts ; then we quadratically subtract the mean internal RV error of all SOPHIE exposures < σ i > from the dispersion of the RV residuals for every targets, σ v = σ 2 v − < σ i > 2 . This should let 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 log R HK and σ v as derived from our sample and compares to the relation obtain in Santos et al. (2000).
We observe that the dispersion of the residuals correlates well with the magnetic activity, with only few outlying points. But 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 σ v,G = 7.8 × (10 5 R HK ) 0.55 . In our Fig. 4. σ v (in m s −1 ) plotted vs 10 5 R HK for SOPHIE (in blue) and SO-PHIE+ (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 σ v = 7.8 (10 5 R HK ) 0.55 with a fit uncertainty of 0.18 dex. The red lines represent the relation derived from the datasets of this work, σ v = 2.6 (10 5 R HK ) 1.0 with a fit uncertainty of 0.3 dex.
case, the slope is stronger, with a linear fit of the log-log relation leading rather to The uncertainty of the fit is σ fit =0.3 dex. The slope is closer to the relation obtained by Saar et al. (1998) HK . After the exclusion of the F and K type stars of our sample, keeping 0.6<B-V<0.8, there remains 18 G-type stars. It leads to a similar relation σ v,G =3.6 × (10 5 R HK ) 0.9 but a larger fit uncertainty of 0.4 dex.
The most significant outlier in Fig. 4 at log R HK ∼4.75 and σ v ∼30 m s −1 is HD207992. We collected 11 RV points in the SO-PHIE configuration, but only 2 with SOPHIE+ for this source. The RV curve in Fig. C.3 shows indeed 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 nor FWHM variations. We conclude that this signal could be a tentative evidence of a third object in the system of HD207992.
One other case is HD161479 with σ v,S + =36 m s −1 and σ v,S − =42 m s −1 . This residuals dispersion is large, but might be compatible with magnetic activity since log R HK =-4.42 for this K0 star. Moreover, according to Table B.4 the bissector and FWHM variations are relatively significant. We measure a pvalue of 0.001 for the no-var model of the FWHM in the SO-PHIE+ dataset, and a strong correlation of -0.95 for the bissector in the SOPHIE dataset, although based on only 4 points. We conclude that the supplementary RV variability of HD161479 is most likely due to magnetic activity.

Results of the Keplerian fit
In total, we characterised 54 massive companions in 54 different systems. We report the Keplerian orbit and M sin i measurements of 12 brown dwarf candidates in the extended range 15-90 M J . One among the 12 is part of a triple system, HD71827, which discovery is reported here. We also characterised the orbit of 42  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 will see in Sections 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 period of the derived orbits are large in general, with only 9 companions below 100-days period. Among the latters, one is member of a triple system and has an M sin i within the BD regime. Eccentricities are large as well, with only 7 orbits with e<0.1. The eccentricities are dispersed around 0.42±0.27. Fig. 6 shows the period-eccentricity distribution of our results, and compare it to the massive planets collected in the Exoplanet.eu database with M>4 M J . We selected systems 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 brown dwarfs reported in this work agrees with that of giant exoplanets. The eccentricities of giant exoplanets are fully compatible in average with the eccentricities of brown dwarfs with e GP ∼0.42±0.22. This is in line with the conclusions of Sozzetti & Desidera (2010) finding strong similarities in terms of eccentricity distributions between massive planets and BD.
One candidate stands apart at small period and large eccentricity, BD+362641, for which the orbit is actually not well constrained because of the small number of points (N RV =9). However, the large radial velocity variation observed of ∼40 km s −1 places it in the M-dwarf regime with a mass most likely larger than 200 M J .

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 brown dwarf candidates are brand new discoveries, among which we report 6 of them with M sin i 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 solarlike 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 Section 8.
Four objects, HD28635, HD210631, HD211681, and HD217850 were already published as brown dwarf candidates. Improvements on their orbital parameters and M sin i are summarised here: HD28635. Also known as "vB 88", it was reported hosting a BD companion with an approximate spectroscopic mass of 70 M J using Keck/HIRES data (Paulson et al. 2004). Adding 13 SOPHIE and 3 Elodie data, we find that the RVs are compatible with a BD-mass companion at a period of 2636.8±2.2 days with M 2 sin i∼77.1±2.7 M J and a 2 ∼4.014±0.068 au. 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.

HD211681.
The companion of this sub-giant 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 SO-PHIE and 12 ELODIE measurements, we are able to narrow down the M sin i 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.
HD217850. The radial velocities variations of this G8-type star were reported to be compatible with an 11 M J compan-Article number, page 7 of 51 A&A proofs: manuscript no. article_BD_accepted ion in Butler et al. (2017) using an incomplete coverage of the orbit with Keck/HIRES data. Adding 41 SOPHIE data we find the lowest mass BD of our sample with an orbital period of 3508.2±2.6 days, an M 2 sin i∼22.27±0.77 M J and a 2 ∼4.672±0.079 au. This is the candidate BD with 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 sin i for these stellar companions.

Triple systems
We found evidence for a secondary drift signal in the RV data of 3 stars, BD+212816, HD71827, and HD212735. The result of fitting a single Keplerian and a drift for each system are summarised in Table B.8. 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 3 sources, the RV offset between SOPHIE or ELODIE, and SOPHIE+ datasets is significantly larger than 100 m s −1 . It should be on the order of 10±30 m s −1 between SO-PHIE and SOPHIE+ measurements (see Section 4.4 below) and on the order of 50-100 m s −1 between ELODIE and SOPHIE+ (Boisse 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 −1 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.

HD71827.
It 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 are clear evidences in SOPHIE+ data of a second signal with a large period and 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 20 yrs orbit.

HD212735.
Apart from an obvious 38-days period signal, the RVs of this system display a significant linear drift during the 10 years 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 yrs and 47 M J for the tertiary. The outer companion is thus possibly a brown dwarf, but 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 probing for the third companion. In every cases, the semi-major axis of the outer orbit is larger than 7 au, with a parallax on the order of 20 mas. Thus it could be seen with adaptive optics that can probe down to about 100 mas in the neighbourhood of stars.

Hipparcos astrometry
In complement to the RV orbital derivation, the Hipparcos astrometry can allow to constrain the inclination of the systems as was performed in e.g. 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 catalog (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), number of field-of-view transit, measurement time span, and abscissa measurement errors. A summary of these informations is presented in Table B.9.
After a preliminary analysis of all systems, we found 16 of them for which there are indications of significant orbital motions in the Intermediate Astrometric Data (IAD), plus 1 system, HD193554, already solved in the Hipparcos double star catalog (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 HD193554 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 radial velocity results given in Tables B.5 and B.6. A two-dimensional grid in i and Ω was searched for its global χ 2 -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.10 with their orbital solution. Table B.11 lists updated parallaxes, proper motion, coordinates offsets, 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 that should be closer to the actual proper motion of systems since based on a 24-years baseline of astrometric data. Finally, Figures C.4, and C.5 show the significant orbits.
In general, the updated Hipparcos-2 parallax are not compatible with the Gaia DR2 parallax at the 1-σ level. Even for 2 systems, HD133621 and HD155228, 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%. Besides, the comparison of the Hipparcos-2 proper motion corrections, after fitting the orbital motion, with the Gaia DR1 proper motion shows a global nice agreement, 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 24years baseline.
In only one case, HD87899, there is no strong agreement of the proper motion corrections. With a long period of 4.2 years, 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 HD87899 should be considered only conservatively within their 3-σ errorbars. If the Gaia DR1 proper motion is correct, the semi-major axis as derived with the Hipparcos-2 data should be rather around 10-15 mas and the mass closer to 0.2 M , which is more in line with what will be derived using Gaia data only in Section 6.
For the two stars HD110376 and HD155228, 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 permutation test. For HD110376 and HD155228, however, this is not the case and the exact reason for the failure of the permutation test is unclear. Because the orbit sizes are relatively large, the F-test null probabilities are very small (2.2 10 −12 and 3.2 10 −7 , 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.10 the significance is lower than 1σ for these 2 systems. Note also that the parallax change caused by fitting the orbit model is large (almost 3 mas) for HD110376.
For HD225239 and HD62923, 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 photocenter out of the primary star center. 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 dwarf, neutron stars or black holes ; they could also be couples of low mass stars. This can be consistent with the log(g) of these two primaries (Table B.3) that are small (∼4.1-4.2) compared to the expected value of surface gravity for G2-3 dwarf stars (∼4.4-4.5). Therefore these 2 primaries might rather be more evolved sub-giants. Precise astrometry with Gaia and imaging can allow determining the exact mass and nature of these companions.
Two brown dwarf candidates, BD+210055 and HD210631, 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%. It leads to a real mass that is significantly larger than the M sin i∼85 M J derived thanks to RV only, with M 2 between 140 and 290 M J , well within the M-dwarf regime. For HD210631, the orbital coverage of 0.3 is not ideal, owing to the long 11-yrs orbital period. Still, the fit of the astrometric motion could catch some significant acceleration in Hipparcos data points and led to a 2σ-detection. It shows that the M sin i∼83 M J of HD210631's companion 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 that in such domain the secondary might contribute light and even produce a secondary peak in the CCF that we actually do not detect in our observations. Finally, we derived upper limits on the astrometric semimajor 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 which projection on the plane of the sky only presents its minor-axis: The value of the upper limits on semi-major axis of the primary and the corresponding companion mass are added to Table B.10. For the triple systems HD71827, HD212735 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-dwarfs regime. We will see in Section 6, that using Gaia's published astrometric data allows tightening up the constraints on the mass for several of these systems, including systems for which N orb <0.8.

Gaia astrometry
To overcome the large uncertainties on the true masses obtained using Hipparcos, we also cross match our sample with the Gaia catalog. We found 51 of our 54 targets in the Gaia DR1 catalog (Gaia collaboration et al. 2016). Among these 51 systems, we could measure the companion mass for 33 of them, and derive upper-limits for 6 others. The companion mass of the 12 that remain out if the 51 systems could not be constrained further with Gaia data. Their mass was nonetheless already bounded from below, thanks to RV, well within the M-dwarf domain.
The DR1 of Gaia does not provide the individual positional measurements for all our sample stars. However, two binarity indicators are published in the released catalogues, the astrometric Article number, page 9 of 51 excess noise and the TGAS discrepancy factor ∆Q (Lindegren et al. 2012, Michalik et al. 2014, Rey et al. 2015.
The astrometric excess noise is a measure of scatter around the 5-parameters astrometric solution as resolved by the Gaia Reduction Software. Given the knowledge of the RV orbital parameters and the dates of Gaia data collection for the DR1, it 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 Section 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 catalog 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-year baseline of astrometric measurements, the 4-years monitoring of Hipparcos-2 and the 14-months monitoring of Gaia for the DR1. Measuring a significant long-term astrometric displacement, it can be used as a binarity diagnosis (Michalik et al. 2014, Lindegren et al. 2016. Comparing the value of ∆Q for every sources 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 Section 6.2 The DR1 archive provides both quantities, while only can be found in the DR2 (Gaia collaboration et al. 2018). Moreover, the excess noise values from the DR1, although based on a shorter timeline 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 have only used the DR1 results for and ∆Q, as extracted from Gaia's DR1 archives 2 . They are presented in Table B.12.

GASTON: Gaia Astrometric Noise Simulation To derive Orbit incliNation
As such, it is not 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. But since the fit of the RVs leads to precise orbital parameters, 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 measure of astrometric excess noise. The photocenter semimajor axis and secondary masses derived using this method are given in Table B.13.

Basic principle
The principle of this method is to simulate Gaia photocenter measurements along the derived RV orbits presented in Section 8. Measurements 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 on-line 3 compare well with random values. Considering random epochs and AL-axis orientation is thus sufficient for applying the method we present here that 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 targets sample, we found that the astrometric excess noise follows a one-to-one correspondence with inclination, owing to the increase of the photocenter 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-attitudemodelling noise of 0.5 mas to Gaia's measurements in the simulation.
Second of all, 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 period on the order of Gaia-Hipparcos baseline (∼25 yrs) 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 period close to 1 yr or larger than ∼20 years. This will tend to underestimate the inclination and overestimate the exact mass of the companion.
We will 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 established, we will be able to derive an interval of inclinations compatible with a given value of DR1 .

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 will 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 : 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 ; the remaining direction u y is directly oriented with respect to 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 which 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 Article number, page 10 of 51 DR1 data collection (t min =2456863.0, t max =2457282.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 along-scan (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: 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 photocenter on the plane of the sky. The photocenter motion has the same orbital parameters than 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 photocenter semi-major axis, as described in the Appendix A. Since no secondary peaks was seen in any CCF 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 photocenter semi-major axis, along with all other orbital parameters from the RV Keplerian solution, in Equation 8 leads to the final simulated Gaia's measurements. Examples are given in Figs. 7 and 8. Out of these simulated data, we are now able to calculate an astrometric excess noise.

Simulated excess noise
The astrometric excess noise is obtained by estimating the χ 2 of its η AL -residuals around the 5-parameters solution derived by Gaia's reduction software (Lindegren et al. 2012). In the simulations, we did not account for the true proper motion and the parallax, assuming they have been already modelled out. It results in only 2 remaining parameters to model out of our simulations, i.e. the (x, y)-position of the photocenter 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 remaining the uncertainties ξ inst,i and ξ AL, j introduced above, the centroid position is found by minimising the squared sum of residuals R with a given observation 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 leads also to the χ 2 of the residuals. In Lindegren et al. (2012) the expression of the χ 2 and of the excess noise with respect to residuals and the AL uncertainty is given by Here we will assume that the down-weighting factors w =1 since we are only interested in the good AL measurements (with w ∼1). The χ 2 should follow a χ 2 distribution with a mean value equal to the number of degrees of freedom, i.e. 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 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 2 +σ 2 AL 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 photocenter 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 with e.g. 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% ∼ DR1 / √ N tot . 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 for not being observed in the spectra, i.e. verifying M V,2 − M V,1 =2.5. In this case 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.

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 5-parameters model, the Gaia's 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 Article number, page 11 of 51 A&A proofs: manuscript no. article_BD_accepted the orbits of our targets compared to the parallax direction are generally unknown, so this cannot be taken into account properly. This could lead to underestimate the photocenter semi-major axis and thus the mass of companions with P∼365 days. This concerns 5 systems with orbital period within 25% of 365 days. On the other hand, the issue of excess proper motion is only relevant for those sources that are member 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 years but rather <416 days, i.e. less than the duration of the DR1 campaign. But in these cases, the proper motion can be fitted out easily since it is purely linear. This is done by slightly modifying the system of equations (9) with a moving centroid u c (t) = (x c + µ x t, y c + µ y t) and inverting the system solving for 4 parameters rather than 2. We incorporated this correction for 8 sources in the secondary dataset: BD+132550, BD+210055, BD+680971, HD147847, HD155228, HD207992, HD225239, HD24505, and HD62923.
-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 photocenter. 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 better constrained companion. This could be a wrong assumption (see below), but still allows deriving 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 bad 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 χ 2 published in the DR2 accounting for a longer baseline of astrometric monitoring (see following section).

A comparison with Gaia DR2
We added in Table B.12 the DR2 normalised unit weight errors, RUWE= χ 2 /2, as defined in Lindegren et al. (2018) and accounting for an average renormalization factor 1/ √ 2 due to the DOF-bug on bright (G<11) targets. Like , the RUWE is a measurement of the astrometric scatter around the 5-parameters so-  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). Most likely the assumption that the inner companion is responsible for the astrometric motion of HD71827 is wrong (see text for explanation).
lution. 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 nice positive correlation, as shown in Figure 9. This strengthen 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-days long baseline, could have been modelled by "instantaneous" proper motion, thus biasing the measurement of astrometric scatter. This stands also 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 years, and for short period systems the DR2 could not have confused orbital motion with proper motion.
In this figure, the outlier HD71827 is a triple system, with an inner companion at a period of 15 days and an outer companion with a period larger than 20 years. The large discrepancy between the DR2 RUWE and the DR1 excess noise shows that Gaia rather 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 HD71827 is a member of the TGAS sample, the whole motion participates to the scatter accounted in the excess noise. Despite this issue, we will keep on assuming that the excess noise of triple systems is due to the inner companion in order to derive an upper-limit on its mass.

Results
In Table B.13, we give the 1-σ bounds of semi-major axis of the photocenter 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, i.e. for which the uncertainty on the orbital parameters does not exceed 10%. This led to reject 5 more stars among the 51 considered in this analysis, HD13014, HD146735, HD40647, HD60846 and HD85533, all of which admit a companion well within the M-dwarf regime. We will thus only consider now the remaining 46 targets with well-constrained RV orbit. Moreover, for 13 out of them, is less than the typical "normal" value for a single star in Gaia's DR1, i.e. 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 produce sensible values of the astrometric excess noise compatible with real Gaia measurements, in more Article number, page 13 of 51 A&A proofs: manuscript no. article_BD_accepted Fig. 9. A comparison of Gaia DR1 excess noise with Gaia DR2 renormalized unit weight error (see text for explanation). Only systems with 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 HD71827 case in discussed in the text. than 90% of the cases. Conversely for 3 systems, HD62923, HD71827 and HD156728, the value of DR1 was difficult to produce, and needed strong fine tuning of the simulated data. The excess noise is either wrongly estimated by Gaia, either one of our assumptions is incorrect, either the astrometry is polluted by one of the companion.
Indeed, HD71827 is a triple system for which the Gaia's DR1-DR2 comparison in Section 6.1.4 suggested that the astrometry recorded by Gaia was rather due to the outer companion. Our assumption that only the inner companion participates to the excess noise was therefore certainly wrong. HD62923 was shown to be massive in Section 5 with a companion that could not be assumed as dark. Finally, the HD156728's phase curve 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 HD71827 is worth commenting further. For this source, we could not produce many s as large as what 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 2 =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 HD71827 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 2 =0.6 M has thus to 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. It includes a comparison between the semi-major axis obtained with Gaia to those obtained by using the Hipparcos Data (Section 5). We also found ground-based speckle interferometry for HD106888 that led to a separation of 32±3 mas and a magnitude difference of ∆M V =1 between the 2 components of this system (Tokovinin et al. 2014). If we assume that the 2 components were at apoastron, this is equivalent to a semi-major axis of the photocenter 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 in- clination is measured and a corresponding Hipparcos solution is derived, the revised semi-major axis of the photocenter a ph tends to always be much closer to the Hipparcos result than the semimajor 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 can take a particular look at the BD candidate BD+210055 b that was well-constrained using Hipparcos astrometry to be an actual M-dwarf, with M 2 =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 to 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 than 90 • at 3-σ. These are BD+362641, HD23965, HD48679, HD73636, HD77712, HD103913, HD106888, HD130396, HD144286, HD153376, HD156111, and HD217850. Five of them are particularly interesting to us since their companion was determined to be  in the BD mass regime thanks to radial velocities: HD23965, HD48679, HD77712, HD130396 and HD217850.
For all five, the simulations could 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 pretty robust. We can safely state that the companions of HD77712, HD130396 and HD217850 are M-dwarfs with masses well above 80 M J . We also recall that the M sin i of HD77712 b was underestimated in Section 4.6.1 due to the deformation of the CCFs by a hidden component. HD77712 b is thus well within the M-dwarf domain. On the other hand, for HD48679 and HD23965, while the inclination is significantly different than 90 • , the companion mass does not exceed 90 M J . These are likely to be brown dwarfs.
Finally, we measured that the mass of 7 companions, among the 46 considered here, are compatible with the BD regime at 1σ. They are BD+291539 b, HD23965 b, HD28635 b, HD48679 b, HD71827 b, HD82460 b, and HD211681 b. We discuss them in more details in Section 7.
Although it is not free of possible systematics, we conclude that the GASTON method is able to derive reliable estimations of systems inclination without the use of the definitive Gaia intermediate data. It proves to be a useful method allowing the characterisation of binaries mass and discarding massive companions with short periods in exoplanet RV survey. We are now applying that method to other catalogues of RV-detected binary stars and exoplanets in order to remove the inclination degeneracy on their M sin i measurements, and thus constrain their true masses. This should show that some bodies now considered as exoplanets actually are face-on binaries.
For the largest period orbits, if virtually nothing can be said using this method, the discrepancy between Hipparcos and Gaia's DR1 proper motions, the ∆Q factor, will be more rele-vant to these cases, with a time baseline larger than 25 yrs. We explore this option in the following section.

The TGAS discrepancy factor ∆Q
While it was pointed 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 Gaia primary sample, i.e. ∆Q>10 (Lindegren et al. 2016), could be considered as significant, and we should conclude that a nonzero acceleration is being detected, advocating for binarity of the system.
As showed in Table B.12, 19 targets have a value of ∆Q larger than 10, while 7 have ∆Q>100 and 1 of them 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 1 . Indeed, ∆Q should be more sensible to large orbital period (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 mass also increase the astrometric acceleration. Because of the degeneracy on the inclination of the systems in RV solutions, only the minimum estimation a 1 sin i is known. Figure 12 display the relation between ∆Q and a 1 sin i 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 1 sin i 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 >1 au.
The system with ∆Q>1000 is HD156728. Its period is larger than 10 yrs, and the companion mass stands in the stellar domain above 126 M J . The primary semi-major axis derived from the radial velocity solution is greater than 0.73 au. The detection of a large value of ∆Q advocates for an underestimation of a 1 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 6 systems with ∆Q>100 are HD108436, HD13014, HD153376, HD60846, HD71827, and HD85533. Their long-period companions stand beyond 100 M J and have all a 1 sin i>1 au. No Hipparcos astrometric solution could be derived for any of these sources owing to the short span of the Hipparcos measurements, smaller than 0.4 orbital periods for all of them. Interestingly, the astrometric excess noise of these 6 systems is significantly greater than 0.5 mas. However, for 4 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 HD108436 and HD153376 are well fitted and the application of the GASTON method leads to re-evaluate the mass at a larger value, and to semi-major axis of the photocenter larger than 20 mas. This is on the order of magnitude of the displacement ∆Q could be able to detect since the precision of Hipparcos is ∼10 mas.
The case of the triple system HD71827 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 Article number, page 15 of 51 A&A proofs: manuscript no. article_BD_accepted Fig. 12. The TGAS discrepancy factor ∆Q vs the a 1 sin i (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.
it remains possible that they do not detect the motion due to the same companion, we found in the preceding section that the DR2 χ 2 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 semimajor 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 a limited usage since it does not allow for deriving the exact mass of the companion.

Detailing the 7 brown-dwarf companions
Among the initial sample of 12 BD candidates derived by RV in Section 4, we excluded 5 of them by using astrometric data of Hipparcos and Gaia. HD210631 b and BD+210055 b were found to be M-dwarfs using the Sahlmann et al. (2011) method on Hipparcos data in Section 5. Moreover, the mass of HD130396 b, HD217850 b, and HD77712 b could be constrained beyond 90 M J thanks to the GASTON method applied on Gaia's DR1 astrometric excess noise in Section 6.
Most importantly, we derived in Section 6 that the mass of the 7 remaining companions could be constrained below 90 M J . All these 7 companions are thus likely brown-dwarfs. We list some details on their detection below.

BD+291539.
We observed for this G-type star radial velocity variations compatible with a 60-M J BD companion on a 176 days 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 −1 . 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 brown-dwarf with a maximum mass about 69 M J . HD211681. From the orbital parameters and minimum mass of the secondary companion in this system, M 2 sin i=77.8±2.6 M J with P=7612±131 days and a 2 =8.28±0.16 au, we deduce an inferior limit of 7 mas for the semi-major axis of the primary's astrometric orbit. 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 HD211681 (Fe/H∼0.36), its companion is likely an object probing a mass regime between star and brown-dwarf, around 80 M J .
HD23965. 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-yrs orbit at 5 au from the star. The large dispersion of the residuals ∼30 m s −1 is compatible with the strong activity that is measured for this source, with log R HK =-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-days period, lead to an inclination of 76±5 • and a companion mass of 42±1 M J . HD23965 b is thus a strong brown dwarf candidate.
HD28635. Paulson et al. (2004) argued that the mass of the companion could be significantly higher than what found with RVs (∼77 M J ). They proposed 0.86±0.31M . Evidence for small inclination was drawn according to the v sin i∼1 km s −1 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 7-years orbit. Excess proper motion could tend to lower the value of that was measured for the DR1.
HD48679. According to the 26 RV measurements obtained with SOPHIE, the companion of this G0 star is a brown dwarf candidate with an M 2 sin i∼36 M 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 −1 . For this star, the astrometric excess noise measured by Gaia of 0.8 mas leads to derive an inclination 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 photocenter is compatible with a non-detection from comparing Gaia and Hipparcos, with ∆Q=3. Thus, HD48679 b is a likely brown dwarf companion.
HD71827. The inclination of this triple system and the true mass of the inner companion (M sin i=26 M J and P=15 days) could not be constrained by astrometry. Indeed, it was shown in Section 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 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 astro-Article number, page 16 of 51 metric scatters 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 HD71827 (P>20 yrs). It follows that the real mass of HD71827b is not constrained and could still be compatible with a mass in the brown-dwarf domain, although it could also still be more massive.
HD82460. Using 17 RV measurements obtained with SOPHIE, we report for this early G-type star the detection of a BD candidate with M sin i∼73-M J on a 590-days orbit, at 1.4 au from the primary. The Keplerian fit is of fairly good quality with medium residual dispersion of 7 m s −1 and orbital elements known with errors smaller than 5%. Hipparcos intermediate data does not lead to astrometric motion detection, and only allow deriving a upper-limit on the mass about 270 M J . On the other han, Gaia measures =0.51 mas, which is also barely significant. The Gaia-Hipparcos discrepancy factor is not much informative with ∆Q<10. This cannot be surprising since the orbital period is short compared to the Gaia-Hipparcos baseline of 25 yrs. 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 most likely an upper-limit on the mass. Therefore HD82460 b should be a brown-dwarf.

Discussion
The present results allow us to complete the statistics of brown dwarfs companions 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 which astrometric motion is easier to measure. Moreover, many companions detected as exoplanets in RV survey 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 catalog (Pourbaix et al. 2004). We included all companions of these publications, even truly stellar, for which the M sin i 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 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 Figure 13. This diagram shows a clear dearth of detection below ∼80-days period (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.

Brown dwarf frequency and a desert below 80 days period
The work presented here add a consequent number of new BD companions candidates, with orbital period shorter than 10,000 days and M sin i within 13.5-90 M J . The full sample of mainsequence FGK stars within 60 pc in the northern sky gathers ∼2950 stars in the new Hipparcos catalog (van Leeuwen et al. 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. More observing time being devoted to sources with more than 3 RV points presenting 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 3 RV points, apart from those published in this paper, we found that there could remain as much as ∼30 more BD candidates with a period up to 10,000 days to characterise. 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, which RV measurements probe a drift-like variation on a baseline of 300 to 5000 days. Their period and M sin i could be considerably larger than these time spans. We thus 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 brown-dwarfs gathered here. This corresponds to a lower limit on the detection frequency of BD candidates within M sin i=13.5-90 M J and with orbital period 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 system's inclination. 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 brown dwarfs 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 over the 44 considered here are compatible with the BD domain. This leads to a revised lower-limit on the BD frequency 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 were not constrained by astrometry yet. We are now applying the GASTON method systematically on all systems with BD candidates to constrain their mass.
Below 80-days period, we find that the detection frequency drops significantly, with only 6 BD candidates found in this region. Within the 2050 systems mentioned above, we found only two possibly missing BDs with period less than 80 days. This leads to a lower-limit f BD, low ∼0.24±0.04%, a factor of 8 lower than above 80 days. This is in line with the findings of Ma & Ge (2014) that brown dwarfs companions tend to avoid a massperiod region bounded by the limits 30-60 M J and P<100 days. Guillot et al. (2014) showed that this dearth of detections below 100-days period, 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.

The companion mass distribution beyond 80-days period
Studying in more details the M sin i-histogram of detections beyond 80-days period, Fig. 14 shows that the M sin i 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 thus tried to reproduce this M sin i distribution out of simulated companions which 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 sin i. We compare these to the 38 detections between 15 and 90 M J and with P>80 days. We thus randomly selected 38 simulated values of M sin i and perform a two-samples Kolmogorov-Smirnov test between the simulated sample and the observed one.
This simulation is performed 1000 times. We then counted the number of simulated samples that are 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 in average 5% of the time. 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 sin i histogram and cumulative plot of detections are plotted in Fig. 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 brown dwarf 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 brown-dwarf companions population, but cannot exclude the existence of such bound in the mass distribution (Halbwachs et al. 2000, Sahlmann et al. 2011. Moreover, a second distribution coming from lower masses, i.e. the massive planets, does not appear necessary to explain the actual M sin i distribution beyond 80 days. This suggests that the mass distribution of massive planets does not spread much within the BD domain. The tail of the distribution probably stops around 20 M J but not much beyond. Independent direct imaging surveys of long-period brown-dwarfs (Brandt et al. 2014) reached similar conclusions, with low-mass BD, even those below the deuterium burning limit, more likely arising, as more massive objects, from gravitational collapse in disk or fragmenting cloud.

Summary and conclusions
We reported here the detection of 54 companions to FGK stars in the neighbourhood of the Sun using radial velocities observations of the SOPHIE spectrograph. Among them, 12 were detected as brown dwarfs candidates according to their projected mass M sin i, and 42 as M-dwarfs.
Using Hipparcos and Gaia, we could reconsider 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's radial veloci-  ties. We found that 5 of them actually stand in the M-dwarf regime, and confirmed that the 7 remaining companions are likely brown-dwarfs. BD+291539 b, HD23965 b, HD48679 b, and HD82460 b are strongly constrained below 90 M J . While HD28635 b, HD211681 b and HD71827 b are possible brown dwarfs, although still remains candidates. Moreover, we obtained a stellar mass for the companion of HD210631 that was previously published as a BD candidate by Latham et al. (2002).
Our 12 BD detected with RV, added to those reviewed and discovered in Wilson et al. (2016), in the SB9 (Pourbaix et al. 2004) and in Sozzetti & Desidera (2010), place a strong lower limit on the orbital period of brown dwarfs at 80 days. The short period region below 80 days appears 4 times less populated than at a larger period. On the other hand, above 80 days the M sin i 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  These conclusions should however not be understood as definitive. This statistics is not completely free of biases yet, since the monitoring of the volume limited sample is not completed. However, we estimated that only few short-period radialvelocity BD 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 period 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 sin i cannot lead to a firm conclusion on the presence or absence of a lower bound on mass for brown dwarfs. 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 brown-dwarfs.
In this work, we showed that such goal will be fulfilled 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 brown dwarf and stellar companions) will become available in the Gaia DR3. In the meanwhile, the GASTON method presented here will allow using the already published data release of Gaia to constraint the inclination and mass of many systems and companions, including planets, brown dwarfs and binary stars.
Acknowledgements. We are thankful to the anonymous referee for his constructive comments that help 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)  To fix the value of the photocenter semi-major axis, we use the following formula (van de Kamp 1975): The period P is expressed in years ; the masses are given in M ; the parallax unit is mas. β is the luminosity fraction β=L 2 /(L 1 + L 2 ) 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) with M RV =M 2/3 K(P/2πG) 1/3 . The solution of this equation, Knowing the primary mass M 1 from Table B.3 and deducing M 2 from q, the luminosity fraction at optical wavelength is derived by Then, the visual luminosity fraction can be calculated thanks to the empirical relation existing between absolute visual magnitude and star mass (Kroupa et al. 1993)  Since no secondary peaks was seen in any CCF 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 too large values of the luminosity fraction.
Article number, page 21 of 51     SOPHIE measurements before and after the instrument upgrade in June 2011 are referred to respectively as S-and S+. Residual O-C and RV γ of supplementary data from different instruments used to calculate these orbits are given in Table B   Notes. (*) We fixed the γ of the S+ and S-dataset at a common value.

Appendix B: Tables
( †) With an eccentricity of 0, the time of periastron passage is ill-defined. In this case, T p indicates the time of primary transit.    Notes.
(a) Doubtful secondary mass estimation, because the Sahlmann et al. (2011) model assumes the companion to be dark. ( †) a rel is the relative semi-major axis, a rel =a 1 +a 2 . Notes.
(a) ∆α and ∆δ, the corrections on the equatorial coordinates of the star in the tangent plane of the sky with respect to Hipparcos-2 catalog. (b) and ∆ , the new parallax and the corresponding corrections with respect to Hipparcos-2 catalog and Gaia DR2 catalog.   Lindegren et al. (2012) and Lindegren et al. (2016). D measures the significance of the value, and should be larger than 2 (Lindegren et al. 2012). N good is the number of reliable measurements taken into account in Gaia's astrometric solution and N FoV is the number of field-of-view transits on the CCD detector. N orb counts the number of orbits covered by the 416-days long DR1 time span.  Table B.13. Photocenter semi-major axis a ph , inclination I c and companion mass M 2 , as derived from Gaia DR1 astrometric excess noise using the GASTON method introduced in Section 6.1. For periods larger than 200 days, these values could be underestimated, as explained in the text. Error bars give 1-σ confidence intervals, but a value of < 0.5 mas leads to an upper-limit on mass and semi-major axis and a lower limit on inclination (see text). Only the sources for which the uncertainty on the period is better than 10% are shown. For triple systems, we calculate the inclination for companion b orbit, assuming that only the inner short-period companion has a measurable effect on Gaia's astrometry. We marked with an * the BD candidates which true mass is estimated to be above 90 M J at 1σ. The definition of the marginal probability is explained in Section 6.1. The minimum inclination I c,min measures the minimum inclination a system can have while the secondary companion remains dark.