The CORALIE survey for southern extrasolar planets XIX. Brown dwarfs and stellar companions unveiled by radial velocity and astrometry

Context. A historical search for exoplanets among a sample of 1647 nearby southern main sequence stars with the CORALIE spectro-graph at La Silla Observatory has been underway since 1998, with a backup subprogram dedicated to the monitoring of binary stars. Aims. We reviewed 25 years of CORALIE measurements and search for Doppler signals consistent with stellar or brown dwarf companions to produce an updated catalog of both known and previously unpublished binary stars in the planet-search sample. We assessed the binarity fraction of the stellar population and survey the prospects for more precise searches for planets in the binary sample. Methods. We performed a new analysis on the CORALIE planet-search sample’s radial velocity measurements, searching for stellar companions and obtaining orbital solutions for both known and new binary systems. We performed simultaneous radial velocity and proper motion anomaly fits on the subset of these systems for which H IPPARCOS and Gaia astrometry measurements are available, obtaining accurate estimates of true mass for the companions. Results. We found 218 stars in the CORALIE sample to have at least one stellar companion, 130 of which are not yet published in the literature and for which we present orbital solutions. The use of the proper motion anomaly allowed us to derive true masses for the stellar companions in 132 systems, which we additionally used to estimate stability regions for possible planetary companions on circumprimary or circumbinary orbits. Finally, we produced detection-limit maps for each star in the sample and obtained occurrence rates of 0 . 43 + 0 . 23 − 0 . 11 % and 12 . 69 + 0 . 87 − 0 . 77 % for brown dwarf and stellar companions, respectively, in the CORALIE sample.


Introduction
Since June 1998, the historical CORALIE exoplanet-search survey has been continuously monitoring a southern hemisphere volume-limited sample composed of 1647 main sequence stars located within 50 pc from the Sun and having spectral types ranging from F8 to K0 (Queloz et al. 2000;Udry et al. 2000).As The radial velocity measurements and additional data products discussed in this paper are available on the DACE web platform at https: //dace.unige.ch/radialVelocities.A copy of the data is also available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/.Appendices A and B, containing orbital solution and detection limit plots, are available in the online version of this paper.
Based on observations collected with the CORALIE spectrograph mounted on the 1.2 m Swiss telescope at La Silla Observatory of the time of writing, the survey has collected more than 60000 radial velocity measurements using the CORALIE Echelle spectrograph mounted on the Euler Telescope at La Silla Observatory, with average measurement precision of ∼5 ms −1 and an average timespan of ∼7600 d.This uniquely long and continuous survey is especially suited for the detection of giant planets with semimajor axes as large as 10 au (Tamuz et al. 2008;Ségransan et al. 2010;Marmier et al. 2013;Rickman et al. 2019) and brown dwarfs (Udry et al. 2002;Santos et al. 2002;Rickman et al. 2019), as well as contributing to statistical studies of the frequency of planetary companions and its dependence on stellar properties (see e.g.Santos et al. 2001;Udry & Santos 2007;Mayor et al. 2011), making the almost 25-year long CORALIE survey an invaluable asset to the field of exoplanetology.Finally, it is worth remark-ing that the continuous monitoring of the less-active stars in the volume-limited sample also makes the CORALIE survey a fertile ground for the search of low-mass exoplanetary candidates suitable for follow-up study with higher-precision instruments such as HARPS (Pepe et al. 2000;Mayor et al. 2003) and ESPRESSO (Pepe et al. 2014), expanding its contributions to the search of exoplanetary bodies toward the realm of terrestrial companions.
During the sample selection process (see Udry et al. 2000), known large amplitude binary stars were collected in a lowpriority subprogram within the planet-search survey, due to the disruptive influence that a close-in stellar companion would have on the stability of the inner region of a planetary system (Holman & Wiegert 1999;Musielak et al. 2005;Turrini et al. 2005;Marzari & Gallina 2016) and in order to limit both the blending effect produced by double-lined spectroscopic binaries (SB2) and the observational effort necessary to disentangle the planetary signal from the higher-amplitude stellar contribution.On the other hand, longer-period binary companions producing only linear trends in radial velocity were instead considered still to be good candidates for the planet search survey, both due to the weak gravitational effect that the distant stellar companion would produce on the inner regions of planetary systems and the fact that such linear trends can easily be corrected for (such as Gl86 and HD41004AB, see Queloz et al. 2000;Santos et al. 2002).This selection strategy also reflected the initial bias against low-separation binaries in favour of stellar environments similar to that of the Solar System (Eggenberger & Udry 2010;Quarles et al. 2020).Still, the large number of measurements collected during almost 25 years of observations with CORALIE have unveiled the binary nature of a non-negligible portion of the stars selected over a wide range of orbital periods, and an updated assessment of the binary population in the sample is a necessary step for advancing the analysis of the CORALIE exoplanet-search survey.
The long-term search for stellar companions within the CORALIE survey represents a key contribution to the current endeavours to further our understanding of stellar formation.Both observational and theoretical studies have now shown that most stars form with at least one stellar companions, and more specifically that about half of FGK stars are th part of a binary system (see e.g.Moe & Di Stefano 2017;Halbwachs et al. 2018;Offner et al. 2022) in which the main-sequence companions appear to follow a lognormal separation distribution peaking around 40 au and roughly uniform mass-ratio distributions (Duquennoy & Mayor 1991;Melo 2003;Raghavan et al. 2010;Tokovinin 2014).It has also been shown that tighter solar-type binaries seem to favor larger mass ratios (Lucy & Ricco 1979;Tokovinin 2000;Moe & Di Stefano 2017), suggesting a common formation and evolution history in a shared circumbinary disc; the fact that wider (a > 200 au) binaries also feature a small but significant fraction of high mass ratio systems (El-Badry et al. 2019) is similarly an indication that at least some wide stellar companion form at intermediate separation and undergo outward migration in later stages of their dynamical evolution.Many different theoretical models have been proposed to explain the formation and observed characteristics of binary systems, such as the fragmentation of filaments and cores in star-forming regions (Könyves et al. 2015;Pineda et al. 2015;Guszejnov & Hopkins 2015;Guszejnov et al. 2017) and of massive accretion discs around individual forming stars (Bonnell 1994;Gammie 2001; Kratter et al. 2010;Harsono et al. 2011), and continuous study of the statistics of binary systems is essential in deepening our understanding of stellar formation.
Considering instead the brown dwarf companion population around solar-type stars, robust study of its demographics is hindered by the currently low number of detection, as less than 100 brown dwarf companions are currently known to orbit such stars (see e.g.Ma & Ge 2014;Grieves et al. 2017), but recent work suggest that only ∼ 4% of solar-type stars have brown dwarf companions (Offner et al. 2022).However, a notable characteristic is a clear paucity of brown dwarf companions around solar-type primary stars on close-in orbits in what is commonly referred to as the brown dwarf desert and that could be explained by post-formation migration processes (see e.g.Grether & Lineweaver 2006;Sahlmann et al. 2011b;Shahaf & Mazeh 2019;Kiefer et al. 2019).
The importance of the characterization of binary stellar systems within the scope of a radial velocity survey aimed at the detection of planetary companions is clear by virtue of the effect that the presence or absence of an additional stellar companion has on the formation and stability of planetary bodies is a fundamental theme in exoplanetary science.While works such as Roell et al. (2012) found the binarity rate among planethosting stars to be about four times smaller than for single solartype stars, the high sample heterogeneity and observational bias are still impediments toward a full understanding of exoplanet demographics in the binary environment (see e.g.Thebault & Haghighipour 2015;Quarles et al. 2020, and references therein for reviews on the subject).More recently, Ngo et al. (2017) found no evidence that host binarity alters the distribution of planet properties in systems characterized by radial velocity observations, while Su et al. (2021) reports a positive correlation between planetary multiplicity and stellar orbital separations in circumprimary planetary systems.
In this paper we present the results of a new analysis of the CORALIE measurements of the 1647 stars in the sample, specifically aimed at the search of radial velocity signals comparable with stellar or brown dwarf companions of the target stars.More precisely, in this study we focus on a specific region of the binary companion parameter space, namely a region limited both in orbital separation as a result of the 25 yr duration of the CORALIE survey and in mass regimes, as we focus on companion having minimum mass higher then 40 M Jup .Companions populating the rest of the parameter space will be the main focus of future papers in this series.We find a total of 218 stars in the sample to have at least one such companion, among which 130 are previously unpublished ones and 88 are instead already known and for which we present updated orbital solutions.Additionally, we present further refined orbital solution for a subset of 132 binary stars in the sample using astrometry constraints provided by Hipparcos (Perryman et al. 1997) and Gaia Early Data Release 3 (Gaia EDR3, Gaia Collaboration et al. 2021) proper motion measurements.
Our paper is organised as follows: in Sect. 2 we describe the physical characteristics of the stars in host stars in our sample.In Sect. 3 we present an overview of the CORALIE observational campaign and of the search for radial velocity signals compatible with brown dwarf and stellar companions, while in Sect. 4 we obtain estimates of dynamical masses for a subset of the presented companions using Hipparcos and Gaia proper motion measurements.In Sect. 5 and 6 we respectively discuss a few systems especially worthy of note and the prospects for follow-up search for exoplanets in the systems comprising our sample, while in Sect.7 we derive occurrence rate values for brown dwarfs and stellar companion in the CORALIE exoplanetary search sample, before concluding and discussing the results of this work in Sect.8.

Host stars characteristics
Out of the 1647 stars composing the CORALIE exoplanetsearch sample, we first exclude known SB2 identified as such either by archive query or by identifying double peaks in the cross-correlation function (CCF) of the radial velocity spectra.
From this we further identify a subset of 218 stars for which we detect a robust radial velocity signal hinting at the presence of a massive companion having minimum mass Msin i 40 M Jup ; the radial velocity analysis that led to this selection is fully detailed in Sect.3.For clarity and simplicity, we'll refer to this subsam-ple of 218 stars hosting stellar companions that is the main focus of this work simply as "the binary sample" throughout this paper, while the larger CORALIE exoplanet-search sample will be referred to as "the CORALIE sample".
In order to have an updated and homogeneous characterization of the physical properties of every star in the sample, we fit the stellar Spectral Energy Distribution (SED) of each star, using the MESA Isochrones and Stellar Tracks (MIST) (Dotter 2016;Choi et al. 2016) via the IDL suite EXOFASTv2 (Eastman et al. 2019).With this method, the stellar parameters are simultaneously constrained by the SED and the MIST isochrones, since the SED primarily constrains the stellar radius R and effective temperature T eff , while a penalty for straying from the MIST evolutionary tracks ensures that the resulting star is physical in nature (see Eastman et al. 2019, for more details on the method).
For each star, we fitted all available archival magnitudes from Tycho B T and V T bands (Høg et al. 2000), Johnson's B, V and 2MASS J, H, K bands from the UCAC4 catalog (Zacharias et al. 2012), WISE bands (Cutri et al. 2021), and Gaia G, G BP and G RP bands (Gaia Collaboration et al. 2016), imposing Gaussian priors on each star's effective temperature T eff and metallicity [Fe/H] based on their respective values in the Anders et al. (2019) catalog, as well as on the stellar parallax based on Gaia EDR3 astrometric measurement (Gaia Collaboration et al. 2021).
The stellar parameters derived from the SED fitting for the binary sample are listed in Table 1, along with each star's archival spectral type, while the distribution of a few selected parameters for both the binary and CORALIE sample are plotted in Fig. 1.The median value of host star mass in our sample is of 0.94 M , and the average relative error on this parameter is 10%; median values and average errors for other stellar parameters of interest are 1.02 R and 4% for stellar radii, 5985 K and 3.72% for effective temperature, 4.41 and 1.41% for surface gravity log g, -0.27 dex and 95% for metallicity.In order to compare the distributions of the binary sample with those of the larger CORALIE sample we perform a Kolmogorov-Smirnov test for each stellar parameter derived from the SED fits, finding p-values < 0.05 for M (p = 0.008) R (p = 0.009) and log g (p = 0.014), suggesting that the underlying population of the binary sample is not the same as the overall CORALIE sample.Indeed, we find the median values of stellar mass, radius and surface gravity in the CORALIE sample to be 0.91 M , 0.95 R and 4.45, suggesting therefore that the underlying population of our binary sample consists of slightly more massive and larger stars than the underlying population of the overall exoplanetary search sample.

Radial velocity observations and analysis
Since its first observations in June 1998, the CORALIE spectrograph went through two significant upgrades in June 2007 and in November 2014 in order to increase the overall efficiency and accuracy of the instrument.Specifically, the 2007 upgrade consisted in the replacement of CORALIE's fiber link and crossdisperser optics (Ségransan et al. 2010), while the 2014 upgrade consisted in replacing CORALIE's fiber link with octagonal fibers (Chazelas et al. 2012) and adding a Fabry-Pérot calibration unit (Cersullo et al. 2017).Both interventions on the instrument introduced small offsets between the radial velocity measurements collected before and after each upgrade, depending on such parameters as the spectral type and systemic velocity of the observed star; during the course of the timeseries analysis we therefore consider CORALIE as three different instruments, each marked by the different upgrades.Dur- Name [cgs] [cgs] [K] HD7320   2: Distributions of orbital period, semimajor axis, minimum mass, eccentricity and mass ratio q of the Msin i >40 M Jup companions identified in the sample and characterized via radial velocity analysis.The minimum masses, eccentricity and q distributions for inner (a < 5 au) and outer (a ≥ 5 au) companions found in the sample are shown in yellow and green respectively.ing the course of this work we'll refer to the original CORALIE dataset as CORALIE-98 (C98), to the dataset collected after the first upgrade as CORALIE-07 (C07) and to the one collected after the most recent upgrade as CORALIE-14 (C14).For selected stars in our sample, we additionally include in the analysis the measurements collected at lower precision (∼300 ms −1 ) with the CORrelation-RAdial-VELocities (CORAVEL) spectrometer (Baranne et al. 1979(Baranne et al. ) between 1981(Baranne et al. and 1998, especially when the CORAVEL data are numerous enough to help identify long-period signals and constraining the orbital parameters of the companions found. As of the time of writing, over the course of almost 25 years of observations on the CORALIE sample we collected a total of 62600 radial velocity measurements for the 1647 stars in the sample, averaging 38 datapoints per star, with median photonnoise uncertainty and timespan of 5.21 ms −1 and 7698 d.In order to search for Doppler signals consistent with the presence of stellar companions in the CORALIE radial velocity timeseries, we follow an iterative process of investigation of successive dominant peaks in the radial velocity periodogram, such as described in Delisle et al. (2016).As mentioned in Sect.2, we once again note that stars found to be SB2s in the CORALIE sample are excluded from the following analysis.
First of all, we consider in our analysis only the 1497 stars for which a total of at least 10 CORALIE measurements have been collected over the years, to ensure robust identification of significant signals.We model instrumental offsets, noise and stellar jitter for each star in the CORALIE sample following the formalism detailed in Díaz et al. (2016) and Delisle et al. (2018), computing false alarm probabilities (FAPs) on the periodogram of the residuals as described in Baluev (2008).The periodogram's main peak is considered significant if characterized by a FAP lower than 0.1% and is modeled as a Keplerian, and the thusly obtained new radial velocity residuals is again investigated for significant signals, re-adjusting the jitter, noise and offsets at each step of the iterative process.This method is however valid only when enough measurements are available to compute a value of FAP; for those cases in which no robust value of FAP is obtained but a clear variation having a scatter in excess of the observation formal errors is present in the radial velocity measurements, we still model the data with a Keplerian model, assessing its significance using the difference between the Bayesian Information Criterion (∆BIC) of the Keplerian and flat models, computed as: with k as the number of model parameters, n the number of datapoints and log L the maximised log-likelihood of the model evaluated following Delisle et al. (2020).Additionally, the longer-period signals found during this search are also modelled as linear or quadratic trends instead, and are included in the final sample that is the main focus of this work only if ∆BIC > 10 in favour of the Keplerian solution.Overall, the signal search process is similar to that undertaken in parallel in Unger et al., in prep., focused instead on identifying new giant planets and brown dwarf companions in the CORALIE sample.
At the end of this analysis, we identified a total of 218 stars featuring at least one signal compatible with a companion minimum mass Msin i 40 M Jup within 1σ, a threshold between giant planets and brown dwarfs we select following the findings reported in Sahlmann et al. (2011a), composing the binary sample representing the main focus of this work.The CORALIE radial velocity dataset for this sample is comprised of a total of 7226 CORALIE measurements, averaging 33 datapoints per star and featuring a median radial velocity uncertainty of 5.31 ms −1 and an observational timespan of 7581 d; all data products are publicly available at the Data and Analysis Center for Exoplanets (DACE) 1 .We run a Markov chain Monte Carlo (MCMC) analysis for each star in our sample based on the algorithm described in Díaz et al. (2016) and Delisle et al. (2016Delisle et al. ( , 2018) )   order to obtain the posterior distribution of the model parameters, using initial conditions drawn from the orbital solutions we obtained during our preliminary iterative signal search and computing each parameter's confidence intervals for a 68.27% confidence level.A summary of the bestfit radial velocity orbital solutions for all companions found in the sample is listed in the left portion of Table 2, the distributions of selected orbital parameters and mass ratio are plotted in the histograms shown in It can be seen that the population of the companions having Msin i >40 M Jup identified by radial velocity in the sample peaks at ∼5.92 au and ∼290 M Jup (around 0.27 M ) and that the orbital elements cover a large variety, with periods ranging from 4 d (for HD196998B) to 65213 d (HD3795B), semimajor axes from 0.045 au (HD196998B) to 36.40 au (HD3795B), minimum masses from 41.60 M Jup (∼0.04 M , HD30774B) to ∼0.71 M (HD181199B) and eccentricities values from fully circular (HD207450B) to 0.95 (HD137763B).It can be also seen that all the companions in the sample have a minimum mass below the solar mass.
Another point of interest is the distinction between single-lined binaries (SB1) and double-lined spectroscopic binaries (SB2) in the sample.Following Halbwachs et al. (2003) we can use the minimum mass ratio parameter q min between secondary and primary component to identify possible SB2s as those having q > 0.8.By doing so we find no binary systems in the sample with such high value of q as characterized by our radial velocity solutions.
To search for differences in the properties of inner (a <5 au) and outer (a ≥5 au) companions we again perform a Kolmogorov-Smirnov test of the orbital elements of the detected companions, finding p-values only marginally lower than 0.05 for minimum mass (p = 0.03) and eccentricity (p = 0.04) sug-10 1 10 0 10 1 a from RV (au) 10 1 10 0 10 1 a from orvara (au) Fig. 4: Comparison between the best-fit semimajor axes obtained by the simultaneous fitting of radial velocity timeseries and proper motions variations and those obtained by the fitting of radial velocities alone.gesting a possible difference in the respective distributions for inner and outer companions.Indeed, as shown in Fig. 2, more companions are found on low-eccentricity inner orbits than outer ones, likely as a result of orbit circularization effects.
Finally, we find three stars in the sample to host more than one companion.As described by our radial velocity solution, HD94340 is a triple star system in which the 1.28 M primary is orbited by a Msin i ∼0.08 M companion on a 6.84 d orbit and by a Msin i ∼0.07 M body with an orbital period of ∼1123 d, while HD206276 hosts a Msin i ∼0.09 M companion on a 32 d orbit and a Msin i ∼0.16 M at 1374 d, while HD196885 hosts both an inner giant planet with minimum mass of 1.95 M Jup and orbital period of 1330 d and an outer stellar companion with Msin i ∼0.26 M on a 14912 d orbit.Two of these systems are already known in the literature (see Tokovinin et al. 2006, 2012for HD94340 and Correia et al. 2008 for HD196885) but in the present work we provide updated orbital parameters for all components, especially with the use of astrometry constrains detailed in Sect.4, while for HD206276 we provide an updated solution for the outer stellar companion, whose presence was already hinted at by astrometric observations, and present a new inner stellar companion (see Sections 5.2-5.4 for more details on the specific systems).

Astrometric constraints
In order to achieve a more complete characterization of the orbital parameters of the stellar companions of the CORALIE sample identified in our radial velocity analysis we use the variations in proper motion measurements δµ in conjunction with the radial velocity data already analysed in Sect 3 to derive precise dynamical mass and orbits for the massive companions in the sample.Specifically, we use the proper motion anomalies between the Hipparcos (Perryman et al. 1997) andGaia (Gaia Collaboration et al. 2016) epochs, which have been shown to be able to detect accelerations as small as a few µas yr −2 (Sahlmann 2016;Brandt 2018).Such an approach has been successfully used more and more commonly in recent years to characterize both stellar and substellar companions (see e.g.Calissendorff & Janson 2018;Snellen & Brown 2018;Kervella et al. 2019aKervella et al. ,b,c, 2020Kervella et al. , 2022;;Brandt et al. 2019Brandt et al. , 2020Brandt et al. , 2021a;;Damasso et al. 2020a,b;Makarov et al. 2021a,b;Venner et al. 2021;Feng et al. 2021;Llop-Sayson et al. 2021), especially with the advent of more and more precise proper motions measurements provided by the different Gaia Data Releases.
In this work we use the proper motion measurements provided by Gaia EDR3, which are on average about 3-4 times more precise than those provided by the previous Gaia DR2 (Gaia Collaboration et al. 2021;Lindegren et al. 2021;Brandt 2021).More specifically, we use the variations ∆µ α,δ from a purely linear motion as reported in the Hipparcos-Gaia Catalog of Accelerations (HGCA, Brandt 2018, 2021), in which the Hipparcos and Gaia EDR3 catalogues have been cross-calibrated to account for systematics and shift all proper motions in the Gaia EDR3.For this purpose, we use the open source code orvara (Brandt et al. 2021b), an MCMC orbit fitting code that is able to fit Keplerian orbits to any combination of HCGA proper motion variations, absolute astrometry, relative astrometry, and radial velocities to obtain precise dynamical masses and orbital elements.We use a version of orvara that has been especially modified to accept priors on orbital periods and semimajor axes 2 , in addition to the priors on primary mass and stellar jitter already included in the original version of the code, to help in achieving better convergence and more constrained orbital elements especially for intermediate-separation binaries.
We find 40 stars in our binary sample not to be included in the HGCA, leaving us with 178 stars for which the simultaneous use of proper motion anomalies and radial velocity is in principle viable.For each of them, we impose priors on the primary mass equal to the stellar masses obtained with the SED fits (see Sect. 2), and on orbital periods equal to the periods retrieved from the radial velocity analysis (see Sect. 3).
As the use of proper motion anomaly is based on two measurements separated by almost 25 yr, the method is clearly more efficient for long orbital periods (Kervella et al. 2019a); when analysing in such a manner the entirety of our sample, which as discussed in Sect. 3 features a large variety of orbital periods, special cautions must be taken.This is particularly evident when comparing the semimajor axes for the companions in the sample as characterized by radial velocities alone with those resulting from the orvara simultaneous fit of proper motion variations and radial velocities.Such a comparison, shown in Fig. 4, highlights a large discrepancy between the values obtained by the two methods for the companions that radial velocity characterized as having shorter separations, and that orvara typically overestimates the orbital periods of such close-in stellar companions.For the purposes of this work, we identify a RV ∼1 au as a reasonable threshold above which astrometry constraints are indeed helpful in constraining the orbital elements of the stellar companions in our sample, and we decide to focus our following analysis only on the 132 stars in our sample for which our radial velocity analysis identified a companion with semimajor axis greater than 1 au.In the cases of HD206276 and HD94340, being both triple star systems hosting one companion below 1 au, we perform the orvara analysis using as radial velocity timeseries the residuals obtained after subtracting the Keplerian signal of the inner stellar companion at 0.19 and 0.08 au respectively, so that the timeseries used for the simultaneous fit contain in principle only the contribution of the outer companions in the systems.Lastly, we must note that binary systems with mass ratios q > 0.6 could, in principle, be affected by the luminosity of the secondary companion shifting the photocenter orbit from the primary orbit.To account for this, we analysed all the CCFs for the systems having The blue dots and histogram respectively show the parameter space position and distribution of semimajor axes and minimum mass as retrieved by the radial velocity analysis alone, while the respective red plot elements refer instead to the results of the simultaneous astrometry and radial velocity fits, and therefore true dynamical mass instead of minimum mass.The components of multiple systems are connected by gray dash-dotted lines, while the horizontal dashed brown, orange and yellow lines respectively indicate the brown dwarf (40 M Jup ), dwarf star (80 M Jup ) and solar-mass (1047.58M Jup ) thresholds.such mass ratios, failing to detect any secondary peak; this corresponds to a lower limit on magnitude ratio of 2.5-3.0, which would lead at most to a barycentric semimajor axis underestimation of of 5-9%.
A summary of the best-fit solution obtained for these 132 systems using orvara is listed in Table 2, while Fig 5 compares the distributions of the stellar companions as characterized by radial velocity alone and the simultaneous use of radial velocity and proper motion anomaly.Generally speaking, this plot shows again the good agreement in orbital separations for the companions in the sample between the two different solutions, with the semimajor axis distribution from the radial velocity fits alone and the one obtained with the addition of proper motion anomalies peaking both at ∼6.30 au.Of clear interest is also the comparison between the mass distributions obtained by the two solutions, as the determination of orbital inclinations derived from the use of astrometry permits the derivation of the true dynamical mass of the stellar companions, shifting the peak of the mass distribution upward to ∼343 M Jup (∼0.33 M ) from the peak value of ∼262 M Jup (∼0.25 M ) for the Msin i distribution derived from radial velocity alone.Fig 6 shows a comparison between the values of selected orbital parameters of Msin i >40 M Jup companions as obtained by radial velocities and those obtained using also astrometry constraints, while Fig 7 shows the distribution of the same parameters as obtained by the simultaneous radial velocity and proper motion anomaly fits.Once again, we generally find good agreement between the orbital separations found from the two solutions and between the eccentricities as well.As done for the orbital parameters obtained from the radial velocity solutions (see Sect. 3) we again perform Kolmogorov-Smirnov tests to search for significant differences in the distributions of the properties of inner (1≤ a <5 au) and outer (a ≥5au) com- panions, finding this time no significant difference between the respective distributions of true mass (p = 0.39) and eccentricities (p = 0.66).
The true masses derived from the simultaneous fit once again show the majority of the companions found in the sample to lie below the solar mass value within the respective uncertainties, although the precise estimate of orbital inclination provided by the astrometry measurements allow some companions to reveal themselves to have a true mass greater than 1 M , namely the ones found orbiting HD27019 (M= 1.74 ± 0.17 M ), HD223084 (M= 0.93 ± 0.06 M ), HD173872 (M= 1.10 +0.47 −0.29 M ), HD181199A (M= 1.18 ± 0.06 M ) and HD3795 (M= 0.90 +0.15  −0.17 M ).As previously done in Sect.3, we follow Halbwachs et al. (2003) and select a threshold of mass ratio q > 0.8 for SB2s in the sample, using this time the companion true mass values  7: Distributions of orbital period, semimajor axis, minimum mass, eccentricity and mass ratio q of Msin i >40 M Jup companions as obtained by the simultaneous radial velocity and proper motion anomaly fits for the 132 stars considered in Sect. 4. The minimum masses, eccentricity and q distributions for inner (1 ≤ a < 5 au) and outer (a ≥ 5 au) companions found in the sample are shown in yellow and green respectively.obtained by the joint radial velocity and proper motion anomalies analysis.By virtue of the true masses derived, this time we find four systems to have q > 0.8, namely HD3795 (q = 0.86), HD39012 (q = 0.81), HD181199A (q = 0.90) and HD223084 (q = 0.90) which are therefore possible undetected blended binaries.
In addition to the companions detected in this sample we also report some results from Unger et al., in prep., in which a Table 2: Best-fit orbital solutions for the binary systems identified in the sample, left side reporting the results from the radial velocity fits (see Sect. 3 and right side referring to the simultaneous radial velocities and proper motions fits (see Sect. 4).We note that systems with mass ratio q > 0.6 could in principle be affected by photocenter bias leading to underestimating the barycentric semimajor axis at most by 5-9%

Notes.
Full table is available at the CDS.A portion is shown here for guidance regarding its form and content.
similar joint analysis of radial velocity and astrometric measurements is performed on the planetary companions detected in the CORALIE sample.Specifically, the possible brown dwarf companions found to be orbiting stars HD162020 and HD112758, having respective minimum masses of 14.96 ± 0.53 M Jup and 32.7±1.9M Jup are found to have true masses of 0.392±0.005M and 0.245 ± 0.001 M .While we shall not go into the details of the orbital solutions of these two companions in the present work as they are thoroughly analysed in Unger et al., in prep., these additional stellar objects in the CORALIE sample will be here considered in the occurrence rates analysis detailed in Sect.13.

Brown dwarfs in the sample
Following the radial velocity analysis detailed in Sect.3, we find 28 companions in the sample to have a minimum mass between 40 and 80 M Jup and that are therefore describable as possible brown dwarfs.While the low number of such companions found in the sample do not allow for robust statistical analysis, it is possible to note that most of them are found at orbital separations larger than 0.5 au from the primary star, in a further example of the brown dwarfs desert observed around solar-type stars (Grether & Lineweaver 2006;Shahaf & Mazeh 2019;Kiefer et al. 2019).As the simultaneous fit of radial velocities and proper motion anomalies performed and discussed in Sect. 4 allow for the determination of orbital inclination and therefore true mass of the companions, it is of clear interest to discuss how many and which of these possible brown dwarfs are confirmed to be as such by this analysis and which are instead revealed to be stellar companions instead.
It is first of all important to note that the primary stars hosting three of these possible brown dwarf companions (namely those found orbiting HD53680, HD153284, HD184860A) are not included in the HGCA, and therefore no astrometric analysis is possible for these stars, and that seven additional possible brown dwarf companions (HD3277, HD28454, HD30774, HD89707, HD151528, HD164427A and HD219709) are found to be on orbits tighter than the 1 au threshold that we have set for reliable analysis with orvara; this therefore leaves us with 18 such companions for which we are able to determine true masses and confirm or reject their nature as brown dwarfs.
Of these 18 objects, our use of proper motion anomalies confirm 7 companions (those orbiting HD4747, HD17289, HD30501, HD74014, HD112863, HD167665 and HD206505) to remain in the brown dwarf mass range, while the remaining 11 (orbiting HIP22059, HD17155, HD20916, HD43848, HD78746, HD87359, HD94340, HD119559, HD154697, HD195010 and HD217580) are found to have a true mass above 80 M Jup and are therefore revealed to be stellar companions.Of the latter, the companion characterized by the larger difference in minimum and true mass is the one found orbiting star HD119559, which starting from a minimum mass of 76.267 +4.149  −4.274 M Jup is revealed by astrometric constraints to have a true mass of 0.35 ± 0.02 M .
Finally, we again report some results from Unger et al., in prep.like done in Sect.4, in this work the usage of Gaia DR3 astrometric measurements allow for four of the aforementioned eleven possible brown dwarf companions to be revealed as stellar companions having a true mass above 80 M Jup , namely the ones found around HD3277 (0.468 +0.023 −0.005 M ), HD89707 (0.100± 0.001 M ), HD151528 (0.1403 +0.0001 −0.0030 M ), and HD164427A (0.3369 +0.003 −0.002 M Jup ).As mentioned before, we refer to Unger et al., in prep.for further details on the respective orbital solution, Fig. 8: Radial velocity best-fit solution and phase folded model curves for triple stellar system HD206276, with CORAVEL, CORALIE98, CORALIE07 and CORALIE14 measurements shown in orange, blue, green and purple respectively.but they will considered in the present work as stellar companions in the occurrence rates analysis detailed in Sect.6.2.

A new triple star system: HD206276 (HIP107143)
The possible presence of a long-period massive companion around K3 V star HD206276 was first noted in Goldin & Makarov (2007), in which the use of a genetic optimizationbased algorithm allowed for obtain additional orbital solutions for a subsample of Hipparcos stars with previous stochastic solutions.In the cited work, HD206276 was reported as possibly hosting a massive companion on a 1338 +171 −73 d orbit with a 0.14 +0.21  −0.11 eccentricity and an orbital inclination of 40 ± 6 deg.However, no radial velocity measurements have been published ever since or used to confirm thet companion or to provide a better orbital solution for the system.
Within the scope of the CORALIE exoplanetary search, we have observed HD206276 over a total of 6932 d collecting 35 radial velocity measurements (divided as 18 C98, 3 C07 and 14 C14); the additional usage of two CORAVEL measurements brings the total of datapoints available for our analysis to 37 over the course of 10584 d.We identify in the timeseries periodogram a highly significant peak (FAP=4.3• 10 −6 ) at 32 d, with an one-A&A proofs: manuscript no.coralie-binaries Keplerian residual peak present at 1363 d with FAP=3.1 • 10 −4 corresponding to the astrometric signal reported in Goldin & Makarov (2007), with no further significant signals present in the two-Keplerian residuals.We therefore present our two-Keplerian bestfit model (shown in Fig. 8) with which we confirm the presence of the outer companions having an orbital period P C = 1374.27± 0.97 d, semiamplitude K C = 3.50 ± 0.03 kms −1 and eccentricity e C = 0.275 ± 0.008, while we also report the detection of an inner companion with P B = 32.005± 0.0002 d, K B = 6.8 ± 0.03 kms −1 and eccentricity e B = 0.255 ± 0.003.By virtue of the primary star having a mass of 0.88 +0.06 −0.05 M (see Sect. 2), we derive values of minimum masses and semimajor axes of 168.89 +7.59 −7.66 M Jup (corresponding to 0.161 ± 0.007 M ) and 2.45±0.05au for the outer companion and of 93.96 +4.09 −4.19 M Jup (corresponding to 0.089±0.004M ) and 0.195±0.004au for the inner one.
As mentioned in Sect.4, when performing the simultaneous fit of radial velocities and proper motion anomaly with orvara we first subtract the Keplerian signal of the inner companion, since its short period of 32 d is not detectable though the use of proper motion measurements at the two Hipparcos and Gaia epochs, therefore using for the joint analysis a radial velocity timeseries containing in principle only the contribution of the outer companion.However, we find the joint orvara to be unable to converge, failing to contraint the orbital elements of the outer stellar companions, likely due to the 32 d stellar companion still producing a significant astrometric contribution to the proper motion anomalies.The same failure to converge is obtained when trying to model the triple system by jointly fitting our radial velocity with the Hipparcos epoch astrometric time series via the kepmodel python package (see Delisle et al. 2016;Delisle & Ségransan 2022) and the samsam MCMC sampler (e.g.Delisle et al. 2018).We additionally note that this triple system is not listed amongst the Gaia DR3 astrometric orbital solutions validated in Holl et al. (2022), further highlighting the difficulties in disentagling the astrometric signature of the inner stellar companion.For this system we are therefore able to provide only the orbital solutions resulting from the radial velocity fit alone.

An updated triple star system: HD94340 (HIP53217)
First hints on the multiple nature of the G4 V star HD94340 were reported in Makarov & Kaplan (2005), in which the presence of a stellar companion to the primary star is suggested by the detection of a large proper motion acceleration between Hipparcos and Tycho-2 measurements.More specifically, using preliminary results from the CORALIE survey, Tokovinin et al. (2006) reports the presence of a stellar companion with an orbital period of 6.8 d and hints of a possible second massive companion at 1200 d; a joint analysis of Hipparcos proper motion anomaly, adaptive optics and speckle interferometry presented in Tokovinin et al. (2012) confirm the triple nature of the system, win an inner companion of 6.8 d and an outer one with an estimated orbital period of 3.3 yr period and 40mas axis unresolved by speckle interferometry.
Over the course of our survey, we collected a total of 63 radial velocity measurements ( 17 with no further residual significant peak and no correspondence between the identified significant peaks and stellar activity signals among the activity indicators analysed.According to our two-Keplerian bestfit model (shown in the top two rows of Fig. 9), we find the inner companion to have an orbital period P B = 6.84 ± 0.01 d, semiamplitude K B = 8.182 ± 0.003 kms −1 and eccentricity e B = 0.009 ± 0.001, while for the outer com-panion we find P C = 1122.96+0.60 −0.48 d, K C = 1.356 +0.008 −0.010 kms −1 and e C = 0.305 ± 0.004; having obtained from the SED fit of the primary star (see Sect. 2) a stellar mass of 1.28 +0.15  −0.19 M we derive values of minimum masses and semimajor axes of 90.08 +8.71  −9.13 M Jup and 0.08 ± 0.01 au for the inner companion and of 77.85 +7.54  −7.93 M Jup and 2.34 +0.11 −0.12 au for the outer companion.As mentioned in Sect. 4 and already done for HD206276, we subtract from the radial velocity timeseries the Keplerian signal of the 6.84 d inner companion, performing the orvara fit using the thusly obtained residuals.From the results of this simultaneous radial velocity and proper motion anomaly fit we find for the outer companion a true mass of M C = 269 +51 −45 M Jup (corresponding to 0.26 +0.05 −0.04 M ), inclination i C = 18.2 +2.8 −2.3 deg, with values of semimajor axis (a C = 2.44 +0.11  −0.12 au) and eccentricity (e C = 0.332 +0.023 −0.022 ) in good agreement with those derived from fitting the radial velocity measurements alone.
We additionaly note that the outer stellar companion in this system is part of the Gaia DR3 astrometric orbital solutions validated in Holl et al. (2022), in which it is characterized by having an orbital period of 1213.8 ± 22.0 d and eccentricity 0.30 ± 0.01.The inclination corresponding to this solution is i C = 20.9+1.1 −1.3 deg, the companion true mass is M C = 0.37 ± 0.03 M , the relative semi-major axis is a C = 2.63 +0.12 −0.13 au.This astrometric-only solution would correspond to a minimum-mass of M C sin i C = 137 ± 15 M Jup or 0.131 ± 0.015 M , which we note is larger by a factor 1.75 than the radial-velocity solution.The parallax is 22.508 ± 0.036 instead of 19.72 ± 0.56 for the Gaia DR3 single-star solution.
Moreover, we additionally analyzed jointly the radial velocity and Hipparcos epoch astrometric time series using the kepmodel python package (see Delisle et al. 2016;Delisle & Ségransan 2022) and the samsam MCMC sampler (e.g.Delisle et al. 2018).While the inclination of the inner 6.8 d companion remains, as expected, unconstrained, the outer 1100 d companion is detected by Hipparcos and its inclination and true mass are constrained.We find an inclination of i C = 13.2 ± 0.8 deg, corrsponding to a true mass of M C = 0.39±0.04M .The relative semi-major axis is a C = 2.5 ± 0.1 au (53.8 ± 2.7mas), and we find an orbital period for the outer companion of 1122.9±0.4 d, an eccentricity of 0.305 +0.004 −0.003 , ω C = 98.9 +0.7 −0.8 deg, Ω = 8.6±4.0 deg and a revised parallax of 21.4 +0.7 −0.8 mas instead of the 19.58 ± 1.46mas reported for the Hipparcos single-star solution.
We find some discrepancy between the joint Hipparcos and radial velocity solution, the orvara solution, and the Gaia DR3 astrometric orbit solution, in particular between the respective outer companion orbital inclination values.We note that the orbital period is slightly longer than the Gaia DR3 timespan, which makes the Gaia DR3 solution typically less accurate.Moreover, the errorbars of the Gaia DR3 solution are probably underestimated, which is confirmed by the poor matching between the Gaia DR3 period and the well determined RV period (about 4σ).Finally, in this regime of period, the orvara solution is expected to be more sensitive to the Gaia scanning-law, which could lead to inaccurate results.We therefore adopt the joint Hipparcos and radial velocity solution for this system.

A planet-hosting system: HD196885 (HIP 101966)
The presence of a possible stellar companion to the F8 V star HD196885 was first reported in Chauvin et al. (2006Chauvin et al. ( , 2007) ) as a result of a near-infrared adaptive optics survey, while Correia et al. (2008) additionally detected a planetary companion of hav-ing minimum mass of 3 M Jup and orbital period of 1349 d using CORAVEL, CORALIE and ELODIE radial velocity measurements.Follow-up relative astrometry and radial velocity observations (Fischer et al. 2009;Chauvin et al. 2011) confirmed the binary nature of the system as well as the orbital parameters of the planetary body found around the primary star.
In the scope of this work, we analyse 18 CORAVEL and 41 CORALIE radial velocity measurements (divided as 33 C98, 1 C07 and 7 C14), over a total of 14574 d of and we are able to provide yet another update on the HD196885 multiple system.The timeseries periodogram features a highly significant region (FAP=1.08• 10 −8 ) at periodicity higher than 16800 d, indicating either a Keplerian signal longer that our observational timespan or a long-term drift trend, and a residual peak at 1323 d with a FAP level of 1 • 10 −3 after subtracting a one-Keplerian solution, with no further significant residual signals.In our two-Keplerian solution, shown in the top row of Fig. 10, we find the outer stellar companion to have orbital period P B = 14912.14+1.13  −1.16 d, semiamplitude K B = 2.102 +0.001 −0.001 kms −1 and eccentricity e B = 0.322 +0.001 −0.003 , for which we then derive using a primary mass of 1.24 +0.12 −0.15 M (see Sect. 2) values of minimum mass M B sin i B =278.19 +21.98  −22.84 M Jup and semimajor axis a B =13.59 +0.50  −0.54 au.Similarly, we characterize the inner planetary companion as having orbital period P b = 1330.64+0.07 −0.43 d, semiamplitude K c = 36.94+1.24 −1.67 ms −1 and eccentricity e b = 0.521 +0.325  −0.085 , and we derive a minimum mass M b sin i b =1.96 +0.34  −0.53 M Jup and semimajor axis a b =2.54 +0.10 −0.11 au.As both companions of the primary star orbit beyond our 1 au threshold, we perform the orvara fit on the original radial velocity timeseries instead of removing one of the companion.As a result, we find for the outer companion a true mass of M B = 334 +26 −27 M Jup (corresponding to 0.32 ± 0.02 M ) and inclination i B = 101.9+1.6 −1.5 deg, and for the inner one M b = 2.67 +1.4 −0.63 M Jup and inclination i C = 89 +42 −44 deg; the bestfit solutions for the proper motions anomaly curves are shown in the bottom row of Fig. 10.
The simultaneous usage of radial velocity and proper motion measurements therefore allow us to confirm the planetary and stellar nature of the inner and outer companion, respectively.Following Tokovinin (1993), we additionally compute the relative orientation of angular momenta ϕ to try and provide further information on the system's configuration.However, due to the highly unconstrained value of the longitude of the ascending node of the inner planetary companion (Ω b = 169 +142 −121 deg), we obtain a relative orientation of ϕ = 87 deg, a value too close to the 90 deg threshold proposed in Tokovinin (1993) to provide any further robust statement on the configuration of the HD196855 system components.

Prospects for exoplanetary search
As discussed in Sect. 1, the binarity of a stellar system has a deep influence on the formation and evolution of planetary companions, especially when the orbital separations between the stellar components of the systems is below a few hundreds au, in which case most studies in the literature suggest that the formation and survival of more massive and close-in exoplanets than those found around single stars (see e.g.Fontanive et al. 2019;Fontanive & Bardalez Gagliuffi 2021;Cadman et al. 2022).Additionally and more closely related to the present work, different studies have highlighted how planetary formation and evolution is especially affected by binary separations less than 100 au (see e.g.Mayer et al. 2005;Moe & Kratter 2021).Fontanive et al. (2019) derives a binary fraction of 79.0 +13.2 −14.7 % for stars hosting giant planets and brown dwarfs on orbits shorter than 1 au, again supporting the critical influence that the presence of stellar companions have on the formation and evolution of such planetary systems.Similarly, based on a literature search for binary companions to exoplanet-hosting stars within 200 pc, Fontanive & Bardalez Gagliuffi (2021) finds that while exoplanets found on circumprimary orbits in very wide binary systems show similar physical properties than those around single stars, tighter binary systems with separations up to a few hundred astronomical units tend to promote instead the formation and survival of more massive giant planets and brown dwarfs with shorter orbital periods and typically in single-planet configurations.The same work also suggests that the properties of close-in exoplanets in wide binary systems are consistent with them being the result of for-mation via fragmentation in a gravitationally unstable disc.This result is further supported by the simulations detailed in Cadman et al. (2022), which find that intermediate separations between the component of a binary stellar system promoting fragmentation is consistent with those featured in the systems displaying an excess of close-in giant planets and brown dwarf (Wang et al. 2014;Kraus et al. 2016;Ngo et al. 2016) While we have detected no new robust exoplanetary signal in the radial velocity analysis of our binary sample the lowto-intermediate orbital separations we derived for the stellar companions discussed in this work, ranging from ∼0.045 au to ∼36.40 au, makes our sample a significant opportunity for the search of exoplanetary companions in binary systems and a suitable testing field for planetary formation theoretical models, provided a larger number of radial velocity measurements with high enough density and precision are collected to successfully disentangle the stellar companion's radial velocity signal from that of any lower-mass body that can orbit the system.

Planetary stability in the binary sample
In order to provide a first assessment of the regions in which exoplanets could be found on stable orbits in the systems of our sample, we use the analytical stability criteria provided in Ballantyne et al. ( 2021) and based on the numerical simulations performed in Holman & Wiegert (1999).
Considering a planet on a circumprimary (or S-type) orbit in a binary system, Ballantyne et al. ( 2021) defines the critical semimajor axis a cS as the maximum stable distance from the primary star: being a bin and e the semimajor axis and eccentricity of the binary, and where: with m p and m s as the masses of the primary and secondary stellar components of the binary.Similarly, for a circumbinary (or P-type) orbit the critical semimajor axis a cP , being the minimum stable distance from the binary system, is given by: While finally considering a planet on a circumsecondaty orbit, the maximum stable distance from the secondary star a cS,sec is given by using again Eq. 2 with instead: As noted in Ballantyne et al. ( 2021), these stability criteria are to be taken as a first-order indication of the circumprimary or circumbinary stability of a planet, since there is a variety of mechanisms that can further enhance or disrupt the stability of such orbits (see e.g.Pilat-Lohinger & Dvorak 2002;Pilat-Lohinger et al. 2003;Parker & Quanz 2013;Lam & Kipping 2018;Quarles et al. 2018Quarles et al. , 2020;;Kong et al. 2021, and references therein).However, since a full dynamical characterization of the binary systems in our sample is beyond the scope of the present work, we assume the validity of the defined stability criteria for the purpose of producing a first estimate of the stability regions of the considered systems.
Additionally, we consider as the minimum stable distance from the primary star the Roche limit of the host star, defined as: where R and ρ are the radius and density of the primary star and ρ pl is the density of the orbiting planet.We use for each primary component the stellar parameters derived from the SED fits detailed in Sect 2, while for the planetary density used to compute the Roche limits we use the Earth density ρ ⊕ as a lower limit scenario.The values of stability limits d Roche , a cS and a cP computed for the systems in the sample are listed in Table 3 We first focus our attention on the binary systems for which we obtained true mass values from the simultaneous fit of radial velocities and proper motion anomalies performed with orvara and detailed in Sect.4, ignoring for the moment the 40 systems absent from the HGCA and the 44 systems with orbital period too short to be detected by proper motion anomalies.We therefore consider in the following analysisc only the 132 systems for which we have derived values of companion true dynamical mass.Additionally we note that the dynamical stability of a planet in the triple systems HD206276 and HD94340 is likely beyond the validity of the Ballantyne et al. ( 2021) criteria and therefore warrants further study and focused analysis such as numerical simulations, and that the stability assessments presented here for these two systems should then be interpreted as firstorder estimates.Fig. 11 shows the system architectures compared with that of the Solar System, the stability regions represented by green bands.
From our stability estimates, we note that only 4 systems in our sample do not allow for stable S-type orbits, namely the HD58696, HD43848, HD139696 and HD118598 systems, due to the high eccentricities (0.98, 0.97, 0.98 and 0.95 respectively) of their secondary components as derived by the simultaneous radial velocity and proper motion analysis.Additionally, we find both triple systems HD206276 and HD94340 to have a very narrow stable region for S-type orbits, spanning below 0.05 au and therefore unlikely to host circumprimary planetary companions, especially by virtue of the presence of the two stellar companions discussed in Sects.5.2-5.3.The system featuring the widest Fig. 11: Overview of the architectures of the 132 binary systems for which true masses were derived using simultaneous radial velocity and proper motion anomaly fits, compared to that of the Solar System.For each system, the companion periastron and apoastron are connected by a thin black line, while the thick green and thin red lines highlight the stable and unstable regions for additional planetary companions respectively as detailed in Sect.6.1; we note that circumbinary stable regions are typically smaller than the companion marker size due to the logarithmic axis used.The mass of each companion is represented by the marker growing size and different color, namely grey for terrestrial (M<2 M ⊕ ), green for Neptune-mass (10 M ⊕ <M<30 M ⊕ ), orange for giant planets (30 M ⊕ <M<40 M Jup ), brown for brown dwarfs (40 M Jup <M<80 M Jup ) and yellow for stellar companions (M>80 M Jup ).Fig. 12: Same as Fig. 11 but for the 86 systems with only radial velocity solutions available, for which the dynamically stable and unstable regions for additional planetary companions are computed using the Msin i value of the detected companions.
S-type stability region in the sample is HD30517, with said region spanning 5.16 au, while the system for which the secondary companion has the larger impact on planetary stability is instead HD3795, in which no planetary orbit appear to be stable from an orbital separation of 5.07 au to 144.91 au from the primary star.The latter system is also the one in the sample characterized by the wider minimum P-type stable orbit in the sample, while the system with the tighter circumbinary stable orbit is HD120559, having a a cP of 2.58 au.Finally focusing on planetary companions orbiting the secondary component of our binary systems, the largest stable circumsecondary region in our sample is found in the HD3795 system with a a cS,sec of 4.65 au.Considering instead the 86 binary systems in the sample for which we have only radial velocity orbital solution and therefore only minimum mass values for the companions detected in these systems, we apply the same stability criterion using the value of Msin i as the m s in Eqs.2-4; the stable regions thus obtained are then to be considered estimates of the maximum range of semimajor axes in which stable orbits for planetary companions are possible.The architecture of these systems is then shown in Fig. 12, from which it can be noted that 3 systems (HD8129, HD89707 and HD137763) do not allow for stable S-type orbits, again by virtue of the large eccentricities (0.95 for all of them) of the detected companions.

Detection limits
In order to investigate the detection capabilities of the CORAVEL and CORALIE data analysed so far we compute the detection limits for the binary sample object of the present work, especially focusing on the substellar (Msin i <80 M Jup ) regime.
To this end we follow a injection and retrieval scheme similar to the one pursued in Barbato et al. (2018), injecting synthetic companion signals into the radial velocity residuals timeseries of each star in the sample, obtained by subtracting from the original radial velocity data the contribution of the companions detected and characterized as described in Sect.3. The synthetic signals are generated over a 40x40 grid of semimajor axes a inj evenly spaced in logarithm from 0.01 to 100 au and minimum masses M inj similarly spaced from 2 M ⊕ to 80 M Jup .For each of the 1600 (a inj , M in j ) realizations we generate and inject into the residuals 500 synthetic radial velocity curves with randomly drawn values of mean longitude λ 0,inj , eccentricity e inj and argument of periastron ω inj .Lastly, we add to each synthetic data thus generated a random Gaussian noise with an amplitude equal to the standard deviation of the instrumental uncertainties of the original timeseries.Each of the resulting 8 • 10 5 synthetic radial velocity timeseries is then fitted with a flat model and a Keplerian one, and the injected signal is considered as detected only if the ∆BIC between the two models is at least 10 points in favour of the Keplerian model.The resulting detection limit maps obtained for each stars in our binary sample are collected in Appendix B, highlighting both the dynamically unstable regions estimated in Sect.6.1 and the parameter space region of additional companions that we can exclude based on the available radial velocity measurements.
We show in Fig. 13 the averaged completeness map for the whole binary sample; while the CORALIE and CORAVEL measurements collected for the sample do not have the precision and sampling necessary to ensure detection completeness for planetary companion below 10 M ⊕ and we have only partially completeness for giant companions below the mass of Jupiter in the explored range of orbital separation, we are instead complete for companions above 1 M Jup within 1.80 au.

Occurrence rates for stellar and brown dwarfs companions in the sample
From the detection limit map produced in Sect.6.2 and shown in Fig 13 it is possible to notice that the radial velocity measurements collected for the sample analysed in the present paper allow us to reach full detection completeness for both brown dwarf companions (40<Msin i <80 M Jup ) with semimajor axis below ∼62 au and stellar (Msin i >80 M Jup ) companions within 100 au from the primary stars.The thorough analysis of the detection limits map produced instead for the whole CORALIE sample will be the focus of a future paper in the series, but for the purposes of the present work it is possible to report that the detection completeness for the aforementioned brown dwarf and stellar companion parameter space is confirmed for the overall CORALIE sample.We exclude from the following analysis the five companions with q > 0.8 identified in Sect. 3 and Sect.4. We can use this information to provide an assessment of the occurrence rate f for brown dwarfs and stellar companions using the binomial distribution: being N the size of the CORALIE search sample and m the number of detections within the parameter space here taken into account.In order to derive f we follow the approach described in Burgasser et al. (2003), McCarthy & Zuckerman (2004) and previously applied in different occurrence rate sudies (such as Sozzetti et al. 2009;Faria et al. 2016;Barbato et al. 2018;Santos et al. 2011, to name a few), where the inverse binomial function p ( f ; m, N) ∝ p(m; N, f ) is normalized to 1 over a range of f values bound between 0 and 1.This yield the result: and the occurrence rate f is then found as the value corresponding to the mode of p .Lastly, upper and lower uncertainty limits f U , f L corresponding to 1σ confidence limits are computed by finding the range covering 68% of the p distribution by numerically solving the relation: Firstly, considering the 209 stellar companions detected (203 from this work and 6 from Unger et al., in prep.)we find an occurrence rate for stellar companions of f = 12.69 +0.87 −0.77 %.If again we distinguish between the 127 inner and 82 outer stellar companions setting the threshold at 5 au, we find occurrence rate values of 7.71 +0.71  −0.60 % and 4.98 +0.59 −0.48 % respectively.It can be noted that these occurrence rates are much lower than the ∼50% computed for the CORAVEL survey presented in Duquennoy & Mayor (1991), but it is important to underline not only the fact that the sample analysed in the cited work was composed of 164 stars and therefore much smaller in size than the whole CORALIE sample we instead considered in the present study, but also that the CORAVEL sample also included a number of SB2 that we instead excluded from our analysis, making a direct comparison between the two studies non-trivial.
Considering the 13 detections with 40<Msin i <80 M Jup for which the astrometric analysis either confirm the brown dwarf nature or is not possible (see Sect. 5.1 and Unger et al.,in prep.)we obtain an occurrence rate for brown dwarf companions in the CORALIE sample of f BD = 0.79 +0.29  −0.16 %; considering only the 8 such companions found within 5 au from the primary star we obtain an occurrence rate of close-in brown dwarfs of 0.49 +0.24  −0.12 %, while for the 5 wider brown dwarf companions we obtain an occurrence rate of 0.30 +0.21  −0.07 %, values that we note to be compatible within 1σ.While this apparent surplus of brown dwarfs on closer orbits might be interpreted as in opposition with the known brown dwarf desert, it is important to remember that, while the CORALIE survey is certainly able to detect the large amplitude signals of brown dwarf companions on wide orbits as evident from Fig. 13, its 25 yr timespan allows only for the robust identification of Keplerian signals corresponding to an orbital separation up to ∼7 au while wider companions would be instead detected as radial velocity trends which have not been considered for analysis in the present work (see Sect. 3): a number of possible long-period brown dwarf companions in the CORALIE sample are therefore possibly still to be detected, leading to an larger occurrence rate difference between the two population, and the same applies to the stellar companions in the sample.
By instead considering only the 7 brown dwarfs in the sample confirmed as such by the joint radial velocity and proper motion analysis we obtain an occurrence rate of f BD = 0.43 +0.23  −0.11 %, which we therefore propose as a lower limit on the occurrence rate of brown dwarfs in the sample.Selecting again a threshold of 5 au between inner and outer brown dwarfs, we find occurrence rate values of 0.12 +0.17 −0.03 % and 0.30 +0.21 −0.07 % respectively.While we now find a lower occurrence rate of inner brown dwarfs, we must again note that these two values are compatible within 1σ.

Discussion and conclusions
We present in this paper the results of the homogeneous analysis performed in search for brown dwarf and stellar companions to the 1647 stars comprising the long-term CORALIE exoplanetary survey, in order to produce an updated catalog of binary objects in the sample using a combination of radial velocity and astrometry measurements.As a result, we find 218 stars in the CORALIE sample to host at least one stellar or brown dwarf companion, 88 of which are already known in the literature and for which we present updated orbital solutions, and 130 of which are instead not known so far and for which we therefore provide first assessment of the orbital parameters.
Furthermore, by combining radial velocity measurements and astrometric accelerations as computed between the Hipparcos and Gaia EDR3 we are able to derive precise dynamical masses of 132 stellar and brown dwarf companions with an orbital separation down to 1 au.Notably, we are also able to confirm the planetary nature of HD196885 b as well as the brown dwarf nature of 7 companions with 40≤Msin i ≤80 M Jup , while we find 11 companions with minimum masses within this range to be revealed as stellar-mass companions, again stressing the power of joint usage of radial velocity and astrometric measurements in painting a full picture of system characterization.
The detection completeness analysis we perform on the sample also allow us to derive occurrence rates f = 12.69 +0.87  −0.77 % and f BD = 0.43 +0.23  −0.11 % for stellar and brown dwarf companions respectively.While our occurrence rates also show an apparent overabundance of stellar and brown dwarf companions below 5 au compared to those found on wider orbits, it is imperative to stress that in the present work we have considered only those companions that are found to be best characterized by Keplerian models instead of linear or parabolic trends, and therefore a possibly large number of wide massive companions are still to be found and fully characterized by continuous observations by spectrographs and especially direct imaging instruments, and will be the subject of future papers in this series.
The binary sample presented and characterized in this work can not only represent an important element for follow-up studies on binary statistics and comparison with formation and evolution theoretical models for both stellar and brown dwarf companions, but also an unparalleled opportunity for the search of exoplanetary bodies in binary systems, as theoretical models have shown that planetary formation and survival is deeply influenced by stellar companions within ∼100 au such as those detailed in the present study.Both the dynamical stability assessment and detection limit maps we have produced show that there is still significant space for the discovery of exoplanets on circumprimary and circumbinary orbits around the stars here analysed, and continued follow-up observations will allow in the near future to deeply probe the exoplanetary discovery space in the sample, allowing this catalog to be used as testing field for models of planetary formation in binary systems.

Fig. 1 :
Fig. 1: Distributions of stellar mass, radius, effective temperature, surface gravity and metallicity of the 218 stars composing the sample discussed in this work (orange) and for the whole CORALIE exoplanetary search sample (green).
Fig.2: Distributions of orbital period, semimajor axis, minimum mass, eccentricity and mass ratio q of the Msin i >40 M Jup companions identified in the sample and characterized via radial velocity analysis.The minimum masses, eccentricity and q distributions for inner (a < 5 au) and outer (a ≥ 5 au) companions found in the sample are shown in yellow and green respectively.

Fig. 3 :
Fig.3: Distribution in the Msin i-a parameter space of the Msin i >40 M Jup companions identified in the sample and characterized via radial velocity analysis.In the main plot, the kernel density estimation of the population is plotted as contour levels, the components of multiple systems are connected by black dash-dotted lines, while the horizontal dashed brown, orange and yellow lines respectively indicate the brown dwarf (40 M Jup ), dwarf star (80 M Jup ) and solar-mass thresholds.The top-left and right-hand histograms show the distribution and kernel density estimation of the semimajor axes and minimum mass of the companions, respectively.
Fig 2, while the companions distribution on the Msin i-a parameter space is shown in Fig 3, and finally the phase folded radial velocity curves for every companion in the sample are collected in Appendix A.

Fig. 5 :
Fig.5: Comparison between the distributions in the Msin i-a parameter space of the Msin i > 40 M Jup companions orbiting the 132 stars in the sample for which we performed simultaneous radial velocity and proper motion anomaly fits.The blue dots and histogram respectively show the parameter space position and distribution of semimajor axes and minimum mass as retrieved by the radial velocity analysis alone, while the respective red plot elements refer instead to the results of the simultaneous astrometry and radial velocity fits, and therefore true dynamical mass instead of minimum mass.The components of multiple systems are connected by gray dash-dotted lines, while the horizontal dashed brown, orange and yellow lines respectively indicate the brown dwarf (40 M Jup ), dwarf star (80 M Jup ) and solar-mass (1047.58M Jup ) thresholds.

Fig. 6 :
Fig.6: Comparison between the best-fit orbital period, semimajor axis, minimum (or true) mass and eccentricity values of Msin i >40 M Jup companions as obtained by the simultaneous fitting of radial velocity timeseries and proper motions variations and those obtained by the fitting of radial velocities alone for the 132 stars considered in Sect. 4.
Fig.7: Distributions of orbital period, semimajor axis, minimum mass, eccentricity and mass ratio q of Msin i >40 M Jup companions as obtained by the simultaneous radial velocity and proper motion anomaly fits for the 132 stars considered in Sect. 4. The minimum masses, eccentricity and q distributions for inner (1 ≤ a < 5 au) and outer (a ≥ 5 au) companions found in the sample are shown in yellow and green respectively.

Fig. 9 :
Fig. 9: Top two rows: same as Fig. 8 but for HD94340.Bottom panel: best-fit solution astrometric orbit induced by HD94340C, Hipparcos intermediate astrometric data shown in red.

Fig
Fig. 10: Best-fit orbital solutions for HD196885.Top row: Radial velocity solution and phase folded model curves, with CORAVEL, CORALIE98, CORALIE07 and CORALIE14 measurements shown in orange, blue, green and purple respectively.Bottom row: Observed and fitted proper motions in right ascension and declination.The thick black line is the best-fit orbit obtained by the simultaneous fit of radial velocities and proper motion anomalies, with 50 orbits randomly drawn from the posterior distributions and color-coded according to the mass of the outer massive companion, and the proper motion measurements from Hipparcos and Gaia EDR3 are shown in orange.Note that in the radial velocity phase folded plot for HD196885 b the low-precision CORAVEL data are not shown in order to better show the radial velocity curve of the low-amplitude inner companion.

Fig. 13 :
Fig. 13: Completeness map of the binary sample, focused on the substellar (2 M ⊕ <Msin i <80 M Jup ) companion regime.Detection frequency contour levels of 10, 50 and 100% are respectively shown as solid, dashed and dotted curves, while the companions detected in our search are shown as white circles.
Fig.A.1: Phase folded radial velocity curves for the companions detected around each star in the binary sample as described in Sect.3.In each plot, the CORAVEL, CORALIE98, CORALIE07 and CORALIE14 measurements are shown in orange, blue, green and purple respectively over the phase folded model radial velocity curve.Companion orbital period and orbit completeness, indicated as the ratio between observational span and period, are noted in each plot's box.

Table 1 :
Stellar parameters for the sample discussed in this work, derived by the SED fits described in Sect. 2 except spectral types, retrieved from Simbad.

Table 3 :
Boundaries of the circumprimary, circumbinary and circumsecondary stability regions of the systems in the binary sample as described in Sect.6.1.The Note column indicates with value of the detected companion mass is used for the stability estimation.
Notes.Full table is available at the CDS.A portion is shown here for guidance regarding its form and content.