The SOPHIE search for northern extrasolar planets -- XVII. A wealth of new objects: Six cool Jupiters, three brown dwarfs, and 16 low-mass binary stars

Distinguishing classes within substellar objects and understanding their formation and evolution need larger samples of substellar companions such as exoplanets, brown dwarfs, and low-mass stars. In this paper, we look for substellar companions using radial velocity surveys of FGK stars with the SOPHIE spectrograph at the Observatoire de Haute-Provence. We assign here the radial velocity variations of 27 stars to their orbital motion induced by low-mass companions. We also constrained their plane-of-the-sky motion using HIPPARCOS and Gaia Data Release 1 measurements, which constrain the true masses of some of these companions. We report the detection and characterization of six cool Jupiters, three brown dwarf candidates, and 16 low-mass stellar companions. We additionally update the orbital parameters of the low-mass star HD 8291 B, and we conclude that the radial velocity variations of HD 204277 are likely due to stellar activity despite resembling the signal of a giant planet. One of the new giant planets, BD+631405 b, adds to the population of highly eccentric cool Jupiters, and it is presently the most massive member. Two of the cool Jupiter systems also exhibit signatures of an additional outer companion. The orbital periods of the new companions span 30 days to 11.5 years, their masses 0.72 Jupiter mass to 0.61 Solar mass, and their eccentricities 0.04 to 0.88. These discoveries probe the diversity of substellar objects and low-mass stars, which will help constrain the models of their formation and evolution.


Introduction
The first exoplanet orbiting a "solar-type star," 51 Pegasi b (Mayor & Queloz 1995), was discovered using the radial velocity (RV) technique. This technique is currently the second-most successful planet detection method, having as of Jan 21, 2021, detected 913 exoplanets ranging from Earth-mass to more massive than Jupiter (www.exoplanet.eu). It remains much more efficient than the transit method at detecting long-period planets, and thanks to improved RV precision and increasingly long temporal baselines, RV surveys are ideal for discovering analogs of Jupiter and Saturn. The ELODIE and SOPHIE exoplanet survey at Observatoire de Haute Provence (France) is one of the longest operating RV surveys, with a total time baseline of over 25 years which extends from before the discovery of 51 Pegasi b.
This paper presents discoveries of new substellar and lowmass stellar companions of solar-type stars -cool Jupiters (CJs), brown dwarfs (BDs), and M-dwarfs or low-mass stars (SCs)from two ongoing exoplanet detection surveys with the SOPHIE spectrograph (Section 2). These substellar objects are classified according to their masses. CJs are defined as massive planets with masses above 0.3 M J and orbital periods above 100 days (Wittenmyer et al. 2020). BDs are substellar objects in the 13-75 M J mass range and occupy the domain between massive planets and stars. SCs are low-mass stars with masses above approxie-mail: shweta.dalal@iap.fr mately 75 M J (0.072 M ). The mass limit that separates BDs from massive planets is conventionally the minimum mass required to fuse deuterium in the core of the substellar object, that is 13 M J (Boss et al. 2003;Chabrier et al. 2014), and the limit between BDs and SCs comes from the minimum mass for hydrogen fusion in the core, that is 75 M J (Chabrier et al. 2000). Both boundaries depend to a small extent on the metallicity of the object (Chabrier & Baraffe 1997;Spiegel et al. 2011).
Despite the discovery of thousands of giant planets and BDs, the statistics on their occurrence and properties are still very incomplete. Most statistical studies of massive planets are for hot and warm Jupiters (Howard et al. 2010;Fernandes et al. 2019), and CJs, which are suspected to be more abundant than hot Jupiters, are poorly characterized as a population (Wittenmyer et al. 2020). Like our Solar System's Jupiter and Saturn, they are dynamically dominant in their system and influence the formation and evolution of any interior planets, including habitable worlds (Raymond et al. 2006;Morbidelli et al. 2012;Raymond & Izidoro 2017;Bryan et al. 2019). Detecting more CJs will help to provide detailed statistics on massive planets with a range of periods from a few days to a decade.
The two leading planet formation models are core accretion (Pollack et al. 1996;Mordasini et al. 2009;Guilera et al. 2010) and disk instability (Cai et al. 2010;Boss 2011). Previous studies Narang et al. 2018;Schlaufman 2018) suggest that giant planets might divide into two distinct giant planet Article number, page 1 of 27 arXiv:2105.09741v1 [astro-ph.EP] 20 May 2021 A&A proofs: manuscript no. aa populations. The metallicity of the host stars of giant planets with a mass above 4 M J is, on average, lower than that of the host stars of giant planets with mass under 4 M J . This hints towards possibly two distinct planet formation scenarios for these two populations. Planet-formation models can probably form bodies up to 40 M J (Ida & Lin 2004;Alibert et al. 2005;Mordasini et al. 2009), which suggests that low-mass BDs might form like massive planets, through the disk gravitational instability scenario. High mass BDs, by contrast, are likely to form like stellar binary systems, through molecular cloud fragmentation (Ma & Ge 2014).
Brown dwarfs are interesting as the transition between the formation mechanisms of giant planets and stars probably runs through their population. One interesting characteristic of the BD population is that few are detected at short orbital periods, a feature known as the BD desert (e.g., Halbwachs et al. 2000;Grether & Lineweaver 2006;Sahlmann et al. 2011b). New detections have shrunk this desert in recent years (Csizmadia & CoRot Team 2016;Wilson et al. 2016;Kiefer et al. 2019;Šubjak et al. 2020), but there is still a detection deficit for orbital periods under 100 days and masses between 30 and 60 M J (Ma & Ge 2014;Ranc et al. 2015;Kiefer et al. 2019). One obvious path towards a better understanding of substellar mass objects is to detect additional objects in and around the desert to better characterize its shape, which in turn will constrain the giant planet and BD formation and evolution processes.
This paper combines the RV technique with astrometry to detect substellar objects and constrain their mass. The RV technique can, in most cases, only determine the minimum mass of the companion (m sin i , where m is the true mass and i the orbital inclination of the planet). As mass is what distinguishes giant planets from BDs, and from SCs, the inclination ambiguity must be lifted to asserting the true nature of a substellar companion (e.g., Díaz et al. 2012;Curiel et al. 2020;Kiefer et al. 2021). As just one example, HD 5388 b was first announced as a likely gas giant ) but turned out to be a BD companion when Sahlmann et al. (2011a) detected its astrometric signature in the HIPPARCOS measurements of HD 5388. Here as well, we combine astrometry and RV measurements to overcome the sin i ambiguity and constrain the true mass of the companions (McArthur et al. 2010;Tokovinin & Latham 2017).
Section 2 presents the two SOPHIE surveys which produced these new detections, while Sect. 3 explains how the observations were performed and the data reduced. Sections 4, 5, and 6 respectively discuss the spectral analysis, the stellar activity, and the RV analysis of the SOPHIE observations. In Sect. 7 we analyze the HIPPARCOS and Gaia astrometric measurements. In Sect. 8 we review the new detections of CJs, BDs and SCs. Finally, Section 9 discusses and summarizes our results.

Description of the SOPHIE surveys
The SOPHIE volume-limited survey for giant planets and BDs observes a catalog of about 2300 FGK stars in the northern sky (δ > +00 : 00 : 00). These targets are within 60 pc of the Sun and have B−V between 0.35 and 1.0. Around 2000 of those stars have SOPHIE observations at this point Kiefer et al. 2019). Two ongoing programs with the SOPHIE spectrograph contribute new substellar companions to this paper.

Giant planets survey
The goal of this ongoing volume-limited program is to increase the number of detections of giant planets orbiting nearby FGK stars, and to identify candidates for follow-up studies: multiplanetary systems for dynamics and transiting systems for structure characterization. This survey constrains the distributions of exoplanet parameters, which helps understand the diversity of planetary systems. The new CJs presented in this paper are the continuation of work by Boisse et al. (2010), Moutou et al. (2014), Díaz et al. (2016), andHébrard et al. (2016).

Brown dwarfs survey
The giant planet survey stops observing any star with a companion which is clearly outside the planetary mass range (i.e., ≥ 13 M J ), and these stars are transferred to the BD survey which has looser RV precision requirements. The goal of the BD program is to obtain an unbiased inventory of companions within and about the BD mass regime for orbital periods up to 10 yrs. This includes stellar companions with mass (or m sin i ) > 75 M J (0.072 M ), because detection of these stellar companions is inevitable while aiming for completeness for BD and because they probe the connections between massive BDs and low-mass stellar companions. The new BDs and SCs presented in this paper are the continuation of work by Díaz et al. (2012), Bouchy et al. (2016), Wilson et al. (2016), and Kiefer et al. (2019).

Observations
We present new observations of 27 stars with the SOPHIE spectrograph, a cross-dispersed, environmentally stabilized echelle spectrograph at the 1.93 m telescope of Observatoire de Haute Provence (OHP). SOPHIE has been in operation since 2006 and covers the 3872 to 6943 Å wavelength range (Perruchot et al. 2008;Bouchy et al. 2009a). The spectrograph is fed through two optical fibers, one of which is always illuminated by starlight from the telescope. Our observations illuminate the second fiber with light from the background sky to estimate its contribution to the on-star spectrum (obj AB mode), and they were carried out in the high-resolution (R = 75000) mode of the spectrograph. Wavelength calibrations and drift measurements were obtained approximately every two hours during the night, as well as at the beginning and end of each night. In 2011, the circular-section fiber was replaced with octagonal-section fiber in the fiber link to improve the stability of the spectrograph illumination (Perruchot et al. 2011;Bouchy et al. 2013), and the pre and post upgrade data have distinct characteristics. This work, therefore, distinguishes two SOPHIE datasets, labeled SOPHIE and SO-PHIE+, depending on whether the spectra were taken before or after this SOPHIE upgrade.

Data reduction
The SOPHIE pipeline extracts the spectra and cross-correlates them with a numerical mask (Bouchy et al. 2009a). The crosscorrelation functions (CCFs) are produced by considering masks corresponding to their stellar type and incorporating all of the spectral orders. The CCFs were then fitted with Gaussians to derive the radial velocities (RVs) (Baranne et al. 1996;Pepe et al. 2002). The exposure time was adjusted to reach a signal-to-noise ratio (S/N per pixel at 550 nm) of at least 50 for the giant planet Article number, page 2 of 27 survey and 30 for the BD survey (observations obtained while a BD target was still in the giant planet survey targeted S/N=50), under diverse weather conditions. Spectra that are significantly contaminated by the Moon were discarded, as were all spectra with less than half of the median S/N for a given target or a large uncertainty on the RV measurements. Spectra for the stars of the giant planet survey and the BD survey have an average S/N of 53.2 and 43.8, respectively. The RVs were also corrected from the CCD charge transfer inefficiency following Bouchy et al. (2009b). Parameters such as Full Width at Half Maximum (FWHM), contrast, stellar rotational velocity (v sin i , where i is the inclination of the star's rotational axis with respect to the line of sight), and the Bisector Inverse Span (BIS), were also derived from the CCF by the SOPHIE reduction software following the method of Boisse et al. (2010). A 5 m s −1 systematic uncertainty is added in quadrature to the uncertainty of the RV measurements obtained before the June 2011 upgrade to account for the poor scrambling properties of the early exposures. The main characteristics of the SOPHIE and SOPHIE+ data sets are summarized in Table A.1.

Spectral analysis
For each star, we performed spectral analyses of an optimally weighted average of all the SOPHIE spectra unaffected by Moon pollution. We used the ARES+MOOG 1 , following closely what was done in Santos et al. (2004) and Sousa et al. (2008) to derive the effective temperature T eff , the surface gravity log g, and the metallicity [Fe/H]. Using those derived spectroscopic parameters as input, stellar masses were derived from the calibration of Torres et al. (2010) with the correction of Santos et al. (2013). Their uncertainties were computed from 10,000 random draws of the stellar parameters within their error bars and assuming Gaussian distributions. Table A.2 lists the resulting stellar parameters, as well as the logR' HK and v sin i, which were obtained following the approach of Boisse et al. (2010).

Stellar activity analysis
Activity in the atmosphere of the star can alter the shape of stellar lines (Queloz et al. 2001), as can face-on binaries (e.g., Díaz et al. 2012;Wright et al. 2013). This gives rise to apparent variations in RV signatures which can mimic a planetary signal. We use various indicators, such as BIS, FWHM, and logR' HK , to probe whether the observed RV signal stems from spectral-line profile changes related to stellar activity.
We evaluate the expected activity-related RV scatter σ a (Table A.3) from our measurement of the mean logR' HK index (Table A.2) using the Santos et al. (2000) σ a vs logR' HK relation. The most active star in the giant planet survey sample is HD 204277, with logR' HK = -4.50 ± 0.11 and σ a 21.8 m s −1 . The dispersion of its measured RVs is on the order of this σ a value, and Sect. 8.1 presents a detailed discussion of the nature of its signal. We exclude it from the rest of our analysis and the next sections are dedicated to the 26 remaining stars. Their σ a range between 5 and 8 m s −1 for the giant planet survey targets and from 5 to 22 m s −1 for the BD survey targets. The dispersion of the measured RVs of these 26 stars is significantly larger than their estimated σ a (see Table A.1 for dispersion of the measured RVs).
To further investigate whether the observed RV signals can be caused by stellar activity, we looked for correlations between the measured RVs and two probes of the line shape, namely the FWHM and BIS. We calculated the Pearson correlation coefficients and the significance of the correlation (p-value) ( Table  A.3) and find that none of them is significant.

Radial velocity analysis
We use the Yorbit software Bouchy et al. 2016) to fit the Keplerian RV signal induced by a companion. Yorbit uses a genetic algorithm to produce starting values for a Levenberg-Marquardt optimization, which in turn provides the priors for a Markov chain Monte Carlo (MCMC) estimation of the error bars following Díaz et al. (2014).
The first step of the RV data analysis is to identify significant periodic signals in the data. This is done by computing the Generalized Lomb-Scargle (GLS) periodogram algorithm of the RV measurements (Zechmeister & Kürster 2009). In the case of giant planet survey targets, we then estimate the false-alarm probabilities (FAPs) of the tentative signals through a bootstrap permutation of the data. The GLS periodogram is however known to fail for the signals induced by companions in highly eccentricity orbits (Zechmeister & Kürster 2009). For two of the stars (BD+631405 and HD 331093), we, therefore, used the Plan-etPack software to compute Keplerian periodograms (Baluev 2013a,b). We ran its kpow command for a range of orbital periods with a frequency step of 0.01 and an upper limit on eccentricity. The latter is needed to control the computational cost, as the closer emax approaches unity, the longer is the computational time.
When a significant period is identified, the RVs are first fit using a single Keplerian orbital model initialized at that period. The following parameters are varied while fitting the single Keplerian model: P the orbital period, K the RV semi-amplitude, e the orbital eccentricity, ω the orbital argument of periastron, T p the time of passage through pericenter, γ S and γ S + the RV offsets for SOPHIE and SOPHIE+ datasets, respectively. To obtain robust confidence intervals for the free parameters, we use 1000 MCMC iterations.
In most cases, a single companion (planet, BD, or stellar) on a Keplerian orbit is a good description of the RV measurements. For HD 124330 and BD+550362, however, linear and quadratic drifts are considered in addition to one Keplerian. Six targets, HD 8291, HD 25603, HD 76332, HD 187057, HD 211961, and HD 352975, have too few RV measurements with either SOPHIE or SOPHIE+ to constrain both γ S and γ S+ . Fitting for both parameters, therefore, produced unrealistically high values for γ S+γ S : they are on the order of 100 m s −1 when the RV offsets between SOPHIE and SOPHIE+ measurements are known to be bounded by 50 m s −1 (Bouchy et al. 2013;Kiefer et al. 2019). We therefore fix γ S+ -γ S to 0 when deriving Keplerian solutions for those six stars.
Section 8 discusses the results of the RV analysis of the newly detected objects in detail, and these objects are divided into CJs, BDs and SCs based on their minimum mass. Three tables summarize the Keplerian orbital elements as well as the m sin i and a derived parameters for each category. Table 1 reports the orbital parameters of the six CJs, while Fig. 1 plots their Keplerian fits as a function of time along with their RV measurements and residuals.

Astrometry analysis
Astrometry can complement the RV orbital information and measure the inclination of the systems. We, therefore, used the HIPPARCOS and Gaia DR1 data to perform the astrometric analysis presented below. The true mass of the companion (either CJ, BD, or SC) is expressed by M c . 7.1. HIPPARCOS astrometry 7.1.1. Selecting the orbit candidates As the input sample for SOPHIE giant planet survey was selected from the HIPPARCOS catalog, all 26 stars listed in Table A.5 were observed by the HIPPARCOS satellite (Perryman et al. 1997). After a preliminary examination of all 26 systems, we selected 12 systems with indications of significant orbital motion 2 . This subset contains all sources with nonstandard HIP-PARCOS solution types ('1' for stochastic solutions, '5' for standard solutions, and '7' or '9' for accelerated solutions). We list upper mass limits in Table A.5 for those of the 26 systems where one could be derived.

Analysis of the HIPPARCOS astrometry
We analyze the Intermediate Astrometric Data (IAD) of the most recent HIPPARCOS reduction (F. van Leeuwen 2007) for signatures of orbital motion, following Sahlmann et al. (2011b) where a detailed description of the method can be found. Fixing the other orbital elements to their values from the RV orbit (Tables 1, 2, and A.4) we adjust seven free parameters to the IAD of each star: the inclination of the orbit I p , the longitude of its ascending node Ω, the parallax , and offsets to both coordinates (∆α , ∆δ) and both proper motion components (∆µ α , ∆µ δ ). We search a two-dimensional grid in the two nonlinear pa-  rameters, I p and Ω, for its global χ 2 -minimum. The false-alarm probability of the derived astrometric orbit was then evaluated through a permutation test employing 1000 pseudo-orbits. The uncertainties in the solution parameters were derived by Monte Carlo simulations which also propagate the uncertainties in the RV parameters. This method has a good track record in reliably detecting orbital signatures in the HIPPARCOS IAD (Sahlmann et al. 2011b,a;Díaz et al. 2012;Sahlmann & Fekel 2013). Table A.5 lists the target names and the basic parameters of the HIPPARCOS observations relevant for the astrometric analysis. The solution type S n indicates the astrometric model adopted by the (F. van Leeuwen 2007) reduction. The code is '5' for the standard five-parameter solution, whereas it is '7' and '9' when the model included proper-motion derivatives of first and second order, respectively. N orb represents the number of orbital periods covered by the HIPPARCOS observation time span, N Hip is the number of astrometric measurements, and σ Λ is their median precision. Outliers in the IAD can very substantially alter the outcome of the astrometric analysis and therefore need to be clipped.
Even when the astrometric data detected no orbital signal (i.e., the derived significance is low), we can set an upper limit on the companion mass by determining the minimum astrometric signal a min that would be detectable for the individual target. When the data cover at least one complete orbit, Sahlmann et al. (2011b,a) showed that an astrometric signal-to-noise of S/N 6 − 7 is required to obtain a detection at the 3 σ level, where S/N = a N Hip /σ Λ and a is the semi-major axis of the detected orbit. We conservatively use a S/N-limit of 8 to derive an upper limit on the astrometric signal of where the factor 1 − e 2 accounts for the most unfavorable orientation of I p = 90 • and Ω = 90 • , in which the astrometric signal is given by the semi-minor axis of the orbit. The last column in Table 3: Solution parameters determined for the significant detections in HIPPARCOS data.   We detect the astrometric orbit of 7 sources with a significance of at least 2 σ, as determined by the permutation test. Those are listed in Table 3, where we also include BD+031552 even though its permutation significance is low. However, the null probability is small, the decrease in residual amplitude is large, the derived orbit looks visually good. Given that this source also has an accelerated HIP solution type, we indicate this solution as viable despite the failure of the permutation test. The 8 sources with orbit determinations then include all 7 sources with nonstandard HIPPARCOS solution types. Table A.6 lists the updated parallaxes, proper motion, inclination, and ascending node of the orbits of those 8 stars. The 8 orbits detected from HIPPARCOS astrometry include no CJ, and we astrometrically detect the orbit of one of the BD candidates, HD 205521 (Fig. 3). Its low inclination however shifts the companion mass from the BD to the SC domain, as discussed in Sect. 8.4. The remaining astrometric orbits are for SCs are plotted in Fig. B.6.

Gaia astrometry
Recent studies have shown that information from the first Gaia Data Release (DR1) can constrain the mass of RV-detected nontransiting companions Mugrauer & Michel 2020). The GASTON code uses the astrometric excess noise, as published in DR1, to determine the amplitude of the astrometric motion of the host star of a known RV-detected companion (Kiefer 2019;Kiefer et al. 2019Kiefer et al. , 2021. This measures, or at least constrains from below, the inclination of the orbit, and thus resolves inclination ambiguity on the true mass. Starting from the RV-derived parameters of the companion, the parallax of the star and an estimation of its mass as priors, GASTON uses a Markov-Chain Monte-Carlo algorithm to explore the space of astrometric orbits compatible with the measured astrometric excess noise, hereafter ε DR1 , and hence to constrain the possible inclinations. It accounts for the systematics in the Gaia DR1 catalog and distinguishes between stars belonging to the TGAS (Tycho-Gaia Astrometric Solution) and the secondary subset of the DR1 catalog (Lindegren et al. 2016;Kiefer et al. 2021). For the targets in the TGAS dataset, the inputs to the 5-parameter fit include a 24-years older Tycho-2 or HIPPARCOS-2 astrometric point in addition to 14 months of Gaia measurements, leading to an improved proper motion accuracy. Targets in the secondary dataset Table 4: Detected orbits with Gaia astrometric excess noise ε DR1 >0.85 mas. All sources belong to the primary dataset of Gaia DR1. When the 3-σ upper-limit on the orbital inclination is above 89.5 • , the constraint on the companion true mass is handled as an upper-limit.  (2021)). Here we additionally impose that the secondary contributes no more than 10% of the total luminosity since the systems we study are single-lined spectroscopic binaries.
Out of 26 targets studied in the present paper, 25 targets figure in the DR1 catalog, with HD 211961 as the sole exception. Table A.7 summarizes for those 25 stars the DR1 inputs to GASTON. For targets in the TGAS (respectively secondary) dataset of the DR1, the reported value of ε DR1 is considered significant beyond noise, if above 0.85 mas (resp. 1.2 mas) (Kiefer et al. 2021). Targets with ε DR1 above these thresholds can lead to, but does not guarantee, a true mass measurement from an allowed orbit inclination range that excludes 90 • . Targets with ε DR1 below the thresholds cannot be distinguished from noise, so can only provide an upper-limit (lower-limit) on the mass (inclination).
Tables 4, A.8, and A.9 summarize the results of this analysis. The results for duplicate sources should be considered with care, as the Gaia observations for a single source can be mistakenly divided between two Gaia "sources" with different IDs (Gaia Collaboration et al. 2016;Lindegren et al. 2016). The astrometric solution in such cases is thus based on an incomplete astrometric data series.
Similar to what happens with the HIPPARCOS astrometry, the 8 orbits detected with the Gaia astrometry include 7 SCs, a single BD candidate, HD 205521, which Gaia DR1 likewise demonstrates is actually an SC (Section 8.4), and no CJs. The remaining detected orbits are detected for the SCs. As discussed at length in Sect. 8, the RV-orbit of the HD 76332 system is incompatible with its observed ε DR1 of 2.19 mas: with the orbital elements of the RV orbit of for HD 76332 as prior, all GASTON Monte Carlo draws produce an astrometric excess noise significantly below 2.19 mas.
Finally, a word of caution. GASTON does not take into account both companions of multiple systems at once, since this would require a totally unpractical number of MCMC steps -and thus a ridiculous time -to converge. We thus examine multiple companions with GASTON, one-by-one, which means that the derived mass of each is overestimated when both contribute significantly to the astrometric signal. In the HD 114762 system ) the wide orbit companion HD 114762 B with a B ∼130 au has little effect on the astrometry of the host star compared to the closer companion with P∼84 days, but a larger mass for the companion would lead to larger perturbations on the star's orbit.

Comparison HIPPARCOS and Gaia astrometry results
We detect astrometric orbits for eight stars in each of the HIP-PARCOS and Gaia astrometric analyses, with five being common to both. The results of the two astrometric analyses agree within 3-σ for the true mass of these five companions, with the exceptions of HD 26596. As discussed in Sect. 8.4, the Gaia analysis finds a significantly lower mass for HD 26596 B.

Results
In this section, we discuss each SOPHIE RV detection in the light of the astrometric data from HIPPARCOS and Gaia, and we classify them as CJs, BDs, or SCs based on their minimum mass (or true mass for some companions). After this RV and astrometric analysis, we have six CJs, three BD candidates, and 17 low-mass stellar companions. We also present a detailed analysis of the RV signal from HD 204277, previously presented as a tentative planetary candidate and which we conclude is due to magnetic activity.
8.1. HD 204277 : activity rather than a planet HD 204277 is a V=6.7 magnitude F8V star located at a distance of 33 parsecs from the Sun and has a stellar mass of 1.14 ± 0.08 M . It is an active star which has logR' HK = -4.50 ± 0.11 and therefore has σ a 21.8 m s HIRES spectrograph, listed this system as an SRC 3 with a period of 30 days. The Generalized Lomb-Scargle periodogram of the SOPHIE RV measurements for HD 204277 instead shows a strong signal around a 250 days period. However, the periodograms of both BIS and FWHM peak at the same period as the RV periodogram, and so does the periodogram of logR' HK (Fig. 4). This suggests starspots, magnetic cycles, and other stellar activities as the likely origin of the RV signal, and there is thus a low probability that it is due to a planet. We conclude that the observed RV variations are likely due to the stellar activity, and will discuss the nature of this signal in a dedicated forthcoming paper.

BD+450564
BD+450564 is a K1 type star with mass 0.81 ± 0.07 M and logR' HK of -4.98 ± 0.11. It was only observed after the SOPHIE upgrade and the final dataset has 14 SOPHIE+ RV measurements. Its GLS periodogram (Fig. B.1) has a clear peak around 300 days, with a FAP below 0.1%. The Keplerian fit of the 14 RVs has a period of 307.8 ± 1.5 days and a semi-amplitude of 47.7 ± 2.8 m s −1 , and it corresponds to a minimum mass of m sin i = 1.36 ± 0.12 M J . The orbit of BD+450564 b is quasicircular, with e = 0.12 ± 0.06. Due to the low (and possibly null) eccentricity, we fit for the time of possible transit instead of the time of periastron. The 3.3 m s −1 dispersion of the RV residuals is consistent with the typical uncertainty of the RVs. No significant astrometric orbit is detected in HIPPARCOS data. The Gaia astrometric excess noise ε = 0.25 leads to a 3-σ upper limit on the true mass of the giant planet companion, that is M c < 31.4 M J .
3 Signal Requiring Confirmation by additional data before rising to classification as planet candidate

BD+550362
BD+550362 is a K3 type star with mass of 0.91 ± 0.10 M and logR' HK = -5.11 ± 0.12. Its 22 RV measurements were all obtained after the SOPHIE upgrade. The GLS periodogram (Fig. B.1) shows two peaks around 260 and 2000 days with FAPs below 0.1% and 1%, respectively. We tested three models for this target: k1 (one-planet Keplerian), k1d2 (one-planet Keplerian and a quadratic drift), and k2 (two-planet Keplerian). An Ftest comparison of the k1 and k1d2 models reveals that the k1d2 model gives a better fit than the k1 model with confidence of 96% (F-value = 24.3). The F-test comparison of the k1d2 and k2 models finds no significant improvement for the k2 model (F-value=3.8, Probability= 84.5%). The k2 model fits the shorter period planet well, but its orbital parameters for the longer period planet are very uncertain because our data do not sample its period well. We, therefore, adopt the k1d2 model for our final fit. The k1d2 Keplerian model fit of the RVs gives a 25.1 ± 1.7 m s −1 semi-amplitude, and a minimum mass of m sin i = 0.72 ± 0.08 M J . The orbit of BD+550362 b has a 265.6 ± 1.0 days period and a significant eccentricity of e = 0.27 ± 0.06. The 3.8 m s −1 , dispersion of its residuals is consistent with the noise of the RV measurements. The RV drift is well fitted by a quadratic trend and we obtained quadratic (d2) and linear (d1) drift coefficients of 1.4365 ± 0.0005 m s −1 yr −2 and 1.76 ± 0.04 m s −1 yr −1 . If one forces a circular orbit for this additional companion, the k2 fit converges to a CJ with an orbital period of at least 3600 days and a mass of at least 2.1 M J . Other orbital solutions, with longer periods, higher eccentricity, and higher masses, are however equally consistent with the measurements.
Additional RV data will therefore be needed to reveal the nature of this outer companion. The 3-σ upper limit on the mass of BD+550362 b from the GASTON astrometric analysis of its host star, M c < 72.45 M J , is not sufficiently tight to nail down its planetary nature.

BD+631405
BD+631405 is a K0 type star with a mass of 0.82 ± 0.08 M and has low magnetic activity (logR' HK = -4.93 ± 0.13). Its 16 RVs were all obtained after the SOPHIE spectrograph upgrade. Figure B.2 shows both a GLS periodogram and a Keplerian periodogram of those velocities. We produced the Keplerian periodogram using the kpow command of the PlanetPack package for orbital periods between 3-1500 days with a frequency step of 0.01 days −1 and an emax=0.91 upper limit on the eccentricity. The GLS periodogram has no obvious peak, but its Keplerian counterpart shows a clear one around 1200 days.
A one-planet Keplerian fit of the RVs finds a highly eccentric orbit (e = 0.88 ± 0.02) of period 1198.5 ± 60.8 days. The semiamplitude, 186.0 ± 14.9 m s −1 of the Keplerian fit corresponds to a planet with minimum mass, m sin i = 3.96 ± 0.31 M J . The poor sampling of the fast RV changes at periastron and the incomplete coverage of the orbital period are the main reasons for the large uncertainty on the orbital period and T p (the time of possible transit is reported in Table 1 rather than the time of periastron). The 4.6 m s −1 dispersion of the residuals is consistent with the noise of the RV measurements. From the ε = 0.39 in the Gaia DR1 astrometry, we set a 3-σ upper limit on the true mass of the giant planet of M c < 40.23 M J .

HD 124330
HD 124330 is a G4IV type star which has a mass of 1.15 ± 0.08 M and logR' HK = -5.27 ± 0.12. Its 58 RV measurements were all secured after the SOPHIE upgrade. The GLS periodogram shows a clear peak around 270 days (Fig. B.1). The k1d1 (oneplanet Keplerian and a linear drift) model fit of the radial velocities gives a significant eccentricity, e = 0.34 ± 0.05 and an orbital period of 270.7 ± 1.2 days. This Keplerian has a semiamplitude of 22.8 ± 1.2 m s −1 and corresponds to a minimum planet mass of m sin i = 0.75 ± 0.06 M J . The 5.4 m s −1 dispersion of the residuals is slightly larger than the noise of the RV measurements.
The 2.95 ± 0.07 m s −1 yr −1 linear drift of the RVs indicates that this system contains an additional outer companion. For a circular orbit and the 3900 days minimum period defined by the extent of the measurements, its mass is at least 0.85 M J , but longer orbital periods and large eccentricities are obviously also compatible with the RV measurements. Neither HIPPAR-COS nor Gaia detects any significant astrometric orbit and we set a 3-σ upper limit on the true mass of the giant planet of M c < 53.76 M J .

HD 155193
HD 155193 is a magnetically quiet F8IV type star with logR' HK = -5.14 ± 0.25 and a 1.22 ± 0.08 M mass. Its 73 RV measurements were all acquired after the SOPHIE upgrade. Their GLS periodogram (Fig. B.1) peaks around 350 days. The one-planet Keplerian fit to the radial velocities has a 352.6 ± 2.6 days orbital period, a 19.5 ± 1.4 m s −1 semiamplitude, a modestly significant eccentricity of e = 0.21 ± 0.08, and it corresponds to a m sin i = 0.75 ± 0.06 M J minimum mass. Although the Keplerian fit appears robust, the nearly one year orbital period induces a large phase gap in our coverage of the planetary orbit, which limits the precision on orbital parameters such as ω and T p . Table 1 reports the time of possible transit instead of the time of periastron. The 7.05 m s −1 dispersion of the residuals is slightly above the estimated accuracy of the measurements. The ε = 0.55 Gaia astrometric excess noise translates into a M c < 63.21 M J 3-σ GASTON upper limit for the true mass of the companion .

HD 331093
HD 331093 is a magnetically quiet K0 type star with logR' HK = -5.10 ± 0.13 and a 1.03 ±0.08 M mass. It was observed before and after the SOPHIE upgrade for a total of 20 RV measurements. Figure B.2 shows both their GLS periodogram and a Keplerian periodogram, which we produced using the kpow command of the PlanetPack package for orbital periods between 3 and 1400 days, with a 0.01 days −1 frequency step, and an emax=0.7 upper limit on the eccentricity. Similar to the also eccentric BD+631405 system, the GLS periodogram is featureless but the Keplerian periodogram shows a clear peak, here around 600 days.
The Keplerian fit of the radial velocities gives a 621.6 ± 16.1 days orbital period, a high eccentricity of e = 0.59 ± 0.03, a 43.6 ± 2.2 m s −1 semi-amplitude, and it corresponds to a m sin i = 1.5 ± 0.1 M J minimum mass. The dispersion of the residuals, 3.3 and 3.2 m s −1 for respectively SOPHIE+ and SOPHIE, is consistent with the typical uncertainty on the RV measurements. The Gaia astrometric excess noise only provides a very loose 3-σ upper limit on the true mass of the companion, M c < 270.5 M J , and therefore contains no additional insight on the nature of HD 331093 b.

Brown dwarfs
The BDs presented in this section are "BD candidates", as the outcome of their astrometric analysis is compatible with a substellar mass but does not demonstrate one.

BD-004475
BD-004475 is a G0 star with a 0.81 ± 0.10 M estimated mass. It was observed both before and after the SOPHIE upgrade, for a total of 13 radial velocities. The Keplerian fit finds a 723 days eccentric orbit (e=0.39) and corresponds to a m sin i = 25-M J BD candidate companion with a 1.48 AU semi-major axis. The dispersion of the residuals of the Keplerian fit is 6.19 m/s. We detect no significant astrometric motion in the HIPPARCOS data, and the 0.72 mas astrometric excess noise measured by Gaia translates into a 125 M J (∼ 0.12 M ) upper limit on the true mass of the companion. BD-004475 b can therefore still be either a BD or an M-dwarf star.

HD 184601
The 15 RVs of HD 184601 were all obtained after the SO-PHIE upgrade. It is a G0 type star with a 0.95 ± 0.07 M mass. HD184601 b has an orbital period of 849 days and an eccentricity of 0.49, and its minimum mass is 60.27 ± 2.15 M J . The dispersion of the residual of the keplerian fit to the RVs is 7.69 m s −1 , and compatible with their measurement noise. Neither HIPPARCOS nor Gaia detect any significant astrometric motion, and the upper limit on the true mass is very loose, 276 M J (∼ 0.2 M ).

HD 5433
HD 5433 is a G5 star with a mass of 0.98 ± 0.07 M . Its 20 RV measurements identify a companion with an orbital period of 576 days, the highest eccentricity (e = 0.81) among the BD candidates presented in this paper, and a minimum mass of 49 ± 3.4 M J . The astrometric analysis of the Gaia DR1 data only sets a loose upper limit on its true mass of 236 M J (∼ 0.23 M ), which leaves the true nature of HD 5433 b undetermined.  Figure B.5 present the Keplerian orbits of the 17 stellar companions with masses (or m sin i ) above 75 M J (or equivalently 0.072 M ). All of these except for HD 8291 are detected for the first time. Their orbital periods range from 30 days to 4198 days and their eccentricities from 0.03 to 0.65, and we discuss a few of the more interesting ones below.

Low-mass stars
HD 205521 is a G5 type star and has a mass of 1.10 ± 0.082 M . The Keplerian orbit has an orbital period of 2032.32 days, an eccentricity of 0.17, and a semi-amplitude of 406.91 ± 5.82 m s −1 . The 2.57 m/s dispersion of its residuals is compatible with the measurement uncertainties of the RVs. The minimum mass of the companion is m sin i = 26.62 ± 1.64 M J and firmly into BD candidate territory, but our analyses of the HIPPAR-COS and Gaia astrometric data both find that the orbit is close to face-on and that the companion is actually a star. The astrometric orbit is firmly detected in the HIPPARCOS time series, with a ∼24 mas semi-major axis of the photocenter orbit around the center of mass. The 2.6 mas Gaia DR1 astrometric excess noise points towards a smaller value of ∼10 +3 −2 mas (1-σ confidence interval) but is compatible with the HIPPARCOS estimate at the 3-σ level. The orbit of the companion is thus within 3 to 12 • of face-on, and its true mass is within the 0.13 to 0.7 M range. While initially classified as a BD candidate from its m sin i , it is a low-mass star.
BD+031552 is a K5 type star with an uncertain stellar mass. Its substellar companion has a 879 days orbital period, an eccentricity of 0.47, and a minimum mass of 133.2 ± 49.0 M J . The analysis of the HIPPARCOS time series of BD+031552 detects its astrometric orbit with low significance but sets a 3-σ upper limit of 0.22 M on the mass of the companion. The upper limit from the Gaia DR1 astrometric excess noise is looser and therefore not informative.
HD 162735 has a stellar companion with an orbital period of 4197 days, an eccentricity of 0.65, and a minimum mass of m sin i = 227.5 ± 13.1 M J . Neither HIPPARCOS nor Gaia detect any significant astrometric motion, and the 3-σ upper limit on the true mass of the companion from the Gaia DR1 excess noise is a loose 0.73 M . HD 162375 B has the highest eccentricity and the longest orbital period of the 17 low-mass stellar companions presented in this work.

Updated parameters for HD 8291b
HD 8291 is a well-known wide binary system with a low-mass (0.073 M ) companion at a 2222 AU projected distance (Baron et al. 2015, and references therein). dos Santos et al. (2017) first reported the detection of a closer-in m sin i = 124.6 ± 2.1 M J companion on an eccentric orbit (e = 0.680 ± 0.009) of 1852.3 ± 3.2 days orbital period using HARPS and SOPHIE data ( Table 3 in dos Santos et al. (2017)). Adding 15 new SO-PHIE and SOPHIE+ RV data to the previous data, we refine this detection to an orbital period of 1862.53 ± 2.9 days, an eccentricity of 0.632 ± 0.001, and m sin i = 99.48 ± 5.86 M J . The poorly covered periastron in dos Santos et al. (2017) led to an overestimated eccentricity and therefore a generally less accurate Keplerian orbit.
This 5-year period M-dwarf companion induces a motion of the star with semi-major axis >4.6 mas, which is indeed detected in Gaia DR1 as a 1.98 mas astrometric excess noise. To 1-σ confidence GASTON finds i c =49±12 • and a companion mass of 142 +37 −21 M J (∼ 0.14 M ), but the 3-σ confidence region includes edge-on orbits.

Incompatible orbits
The companion of HD 76332 has an orbital period of 2489 days, an eccentricity of 0.14, and m sin i of 216.01 ± 12.46 M J . The semi-major axis of the reflex orbit is at least 13 mas, and thus expected to be well detected by both HIPPARCOS and Gaia. They indeed both detect significant nonlinear motion even though their measurements do not span a full orbit, but cannot agree on a common true mass for this companion. The HIPPARCOS time series detect the reflex orbital motion with high significance, and find a a=25±5 mas semi-major axis and a 3-σ range of 0.24-0.84 M for the companion mass. Masses above ∼0.6 M would give rise to an unobserved second peak in the CCFs of this system and can therefore be excluded. The Gaia DR1 astrometric excess noise of 2.2 mas is well below the minimum expected semimajor axis. This is qualitatively expected from Gaia DR1 only covering ∼1/6 of the orbit (HD 76332 is not part of the TGAS dataset), but none of the GASTON simulations, whatever the mass of their companion, could produce a Gaia DR1 astrometric excess noise above ε=1.55 mas. We note that HD 76332 has a duplicate source and that duplication may be the sign that something went wrong during the reduction for this specific source. For non-TGAS sources, the GASTON simulations fit the proper motion through the orbital motion observed with Gaia and we suspect that the best explanation for this discrepancy is an inaccurate proper motion fit in the DR1. If not fitting for the proper motion, GASTON leads indeed to I c within 40-70 • . We adopt the HIPPARCOS result as the best estimation of HD 76332 B's true mass.
HD 26596 has a stellar companion in an eccentric orbit (e=0.4) at 890 days and a minimum mass of 121 ± 7 M J . To 3σ confidence, GASTON derives an upper limit of < 0.19 M for the true mass and the HIPPARCOS analysis derives a range of 0.23-0.43 M . This source again has a duplicate in the DR1 database, and we suspect that its Gaia astrometric excess noise is underestimated. We therefore again prefer the HIPPARCOS confidence range for the mass of this companion to its likely underestimated GASTON upper limit.

Discussion and conclusions
This paper reports the discovery of 6 CJs, 3 BD candidates, and 16 SCs with the SOPHIE spectrograph at OHP. We also present updated orbital parameters for the low-mass star HD 8291 B. We analyzed the HIPPARCOS and Gaia astrometry to constrain the inclinations of the orbits and hence, the true masses of the companions. Figure 5 shows the period-mass distribution of these 26 companions.  (Tamuz et al. 2008), and is today the most massive of these 5 highly eccentric giant planets. Figure 6 shows the periodeccentricity distribution of all known giant planets ( Table 1 in Ma & Ge (2014) and Table A ant planet hosts in this paper are metal-rich, and their average metallicity is 0.07. We compare the metallicities of two samples of giant planet hosts (including the new detections), divided between planetary masses above and below 4 M J . Figure 7 shows the cumulative metallicity distribution for host stars of these two samples, which a Kolmogorov-Smirnov (K-S) test 5 shows has a very low p-value of 9.53 ×10 −8 to be drawn from a single population. This reinforces the previous observations that the average metallicity of the host star decreases as the mass of the planet increases Narang et al. 2018;Swastik et al. 2021). A possible interpretation is that the population of giant planets is bimodal, with giant planets (m p ≤ 4 M J ) forming via core accretion and more massive giant planets (m p > 4 M J ) forming via the disk gravitational instability. The three new BD candidates have minimum masses between 25 and 60 M J and orbital periods between 576 and 850 days. Figure 8 shows the period-eccentricity distribution of all known BDs. All three newly detected BD candidates have significant eccentricities. The lightest of the three, BD-004475, with a minimum mass of 25 M J , has an eccentricity below 0.4, and the two BDs with masses 49 and 60 M J have eccentricities of 0.5 and 0.8. These new detections are consistent with the observation by Ma & Ge (2014) that a significant number of the BDs with (minimum) masses below 42.5 M J and with orbital periods between 300 and 3000 days have eccentricities below 0.4. The heaviest CJ and BD companions in our small sample have higher eccentricities.
We also present RV and astrometric orbit solutions for 16 low-mass stellar companions heavier than 75 M J or 0.072 M . These low-mass stars have orbital periods between 30 days and 11.5 years and have eccentricities between 0.04 and 0.65. These objects were followed with SOPHIE as potential BDs and turned out to be SCs. Information from the Gaia and HIPPARCOS astrometric space missions allowed us to constrain the orbital inclination and hence, the true mass of some of these objects.
The SOPHIE giant planet and BD surveys play an important role in identifying promising targets for future direct imaging detections. Furthermore, the SOPHIE sample is volume-limited, which will help to obtain unbiased statistics for these planetary systems, which is crucial for understanding the formation and evolution of these objects. N is the Number of RV measurements. σ a is the expected activity-related RV scatter which is estimated from the mean logR' HK value for each star. BIS mean is the mean of bisector span and BIS std is the variation in bisector span. FWHM mean is the mean of full width at half maximum and FWHM std is the variation in FWHM. RV-BIS and RV-FWHM is the Pearson correlation coefficient between RV-BIS and RV-FWHM respectively, and the significance of correlation is computed by p-value (For N ≥ 3).

Notes:
The uncertainty in the host star mass is taken into account while obtaining the uncertainty on m sin i and a.   Table A.9: Companion HD 76332 B, whose RV orbit is incompatible with Gaia astrometric excess noise. The minimum and maximum ε that we were able to simulate for this secondary dataset source is given as ε simu,min and ε simu,max .