3-D selection of 167 sub-stellar companions to nearby stars

We analyze 5108 AFGKM stars with at least five high precision radial velocity points as well as Gaia and Hipparcos astrometric data utilizing a novel pipeline developed in previous work. We find 914 radial velocity signals with periods longer than 1000\,d. Around these signals, 167 cold giants and 68 other types of companions are identified by combined analyses of radial velocity, astrometry, and imaging data. Without correcting for detection bias, we estimate the minimum occurrence rate of the wide-orbit brown dwarfs to be 1.3\%, and find a significant brown dwarf valley around 40 $M_{\rm Jup}$. We also find a power-law distribution in the host binary fraction beyond 3 au similar to that found for single stars, indicating no preference of multiplicity for brown dwarfs. Our work also reveals nine sub-stellar systems (GJ 234 B, GJ 494 B, HD 13724 b, HD 182488 b, HD 39060 b and c, HD 4113 C, HD 42581 d, HD 7449 B, and HD 984 b) that have previously been directly imaged, and many others that are observable at existing facilities. Depending on their ages we estimate that an additional 10-57 sub-stellar objects within our sample can be detected with current imaging facilities, extending the imaged cold (or old) giants by an order of magnitude.


INTRODUCTION
Brown dwarfs (BDs) are failed stars, which were unable to initiate nuclear fusion of hydrogen and helium. They are not considered planets as they can induce the fusion of other light elements such as Deuterium and, for very massive BDs, the fusion of Lithium (e.g., Burgasser 2008). Following the studies of Kumar (1963); Hayashi & Nakano (1963); Burrows et al. (2001); Spiegel et al. (2011); Baraffe et al. (2015); Marley et al. (2021), we use the hydrogen-burning and deuterium-burning mass limits of 13 and 75 M Jup to define the range of BD masses, although we are aware of the issue of using such an observationally ambiguous distinction (Boss 1996). While stars are typically formed through core collapse and planets form in circumstellar disks, brown dwarfs could form through both channels. Hence BDs can be considered as an unique population bridging planets and stars, deserving intensive scientific investigations.
BDs can have four different spectral types M, L, T, and Y. GJ 229 B (Nakajima et al. 1995) and Teide 1 (Rebolo et al. 1995) are the first two BDs unambiguously discovered through direct imaging in 1995. The former is a T-type companion BD to an M dwarf, while the latter is a free-floating M-type brown dwarf located in the Pleiades open star cluster. Since 1995, astronomers have discovered more than 2500 BDs (e.g., Kirkpatrick et al. 2021) with the vast majority of those classified as individual "free-floating" objects (e.g., objects which are gravitationally bounded only to the galaxy central potential rather than as a companion in a binary). Luhman (2007) found stars and brown dwarfs are mixed homogeneously based on their spatial kinematics being indistinguishable in Chamaeleon I, a young star forming region, which further supports that both stars and brown dwarfs have the same formation mechanism.
For directly imaged BDs, whether free-floating or having a stellar companion, the cooling model is typically used to provide an indirect estimate of their mass (e.g., Baraffe et al. 2015;Marley et al. 2021). Due to the lack of hydrogen fusion in the core of BDs, they cannot sustain their high temperature and brightness and thus cool down over time. The cooling process as a function of time depends on BD mass, metallicity, cloud coverage, etc. However, cooling models are diverse and sometimes estimate a mass inconsistent with the dynamical mass constrained by direct imaging of BDs. Hence a large sample of BD companions to stars with known masses are essential to test various cooling models. The brown dwarf host star provides a natural reference point where age might potentially be accurately determined for us to better quantify the formation and evolution of BDs.
While we have few direct mass determinations of the thousands of free-floating BDs, we have a much smaller population of BDs orbiting around stars, whose masses are well constrained from orbital fits. Hereafter, we will refer to these BDs orbiting stars as "circumstellar BDs". Radial velocity surveys find the occurrence rate of circumstellar BDs is measured to be about 0.5 − 2% from samples of thousand stars (e.g., Vogt et al. 2002;Patel et al. 2007;Sahlmann et al. 2011;Grieves et al. 2017;Kiefer et al. 2019). The BD desert hypothesis was proposed to explain the low detection rate of circumstellar BDs. This hypothesis was formulated in the late 1980s, when the first precise radial velocity surveys compiled results (e.g., Campbell et al. 1988;Marcy & Benitz 1989;Marcy & Butler 1995, although the significant observational biases of radial velocity (RV) and latter transit surveys were considered as plausible causes. Through an adaptive optics imaging survey, Metchev et al. (2009) found that there are more BDs at wide orbits than BDs in the brown dwarf desert. Recently, more and more BDs have been found to reside in this desert (e.g., Persson et al. 2019;Carmichael et al. 2019Carmichael et al. , 2020Acton et al. 2021). As suggested by Shahaf & Mazeh (2019), characterization of the shape of the brown dwarf desert in the period-mass diagram by a large sample of circumstellar BDs can improve our understanding of its origin.
The correlation between inner companion planets and the wide-orbit companions may be the key to improve our understanding of planet formation. For example, one of the puzzles is whether the existence of a companion in a wide orbit affects the formation of inner companions (Fontanive et al. 2019;Ziegler et al. 2021). In this work, we can detect both inner companions and BDs with wide orbits. Therefore, our sample allows the study of the correlation between the wide-orbit companion and the inner companion. Another long-lasting puzzle is whether cold Jupiters (hereafter CJs) were formed by core accretion or gravitational instability (Chabrier et al. 2014). In this manner, our work aims to bolster the sample of detected systems hosting CJs in order to provide improved constraints for theoretical formation models. According to Zhu et al. (2012), when BDs form due to gravitational instability, multiple CJs and BDs are expected to form simultaneously. Therefore, it is reasonable to consider that CJs formed inside the orbit of BDs formed by the gravitational instability are both formed through the same process. We expect to be sensitive to such systems in this project.
To constrain the BD cooling models, to test the BD desert hypothesis, to find the correlation between BDs, CJs and other types of planets, and to provide a large sample of candidates for direct imaging, we need a larger circumstellar sub-stellar sample covering the BD mass range. This requires the improvement of BD detection sensitivity through the synergy of all available detection techniques. A good approach is to combine the RV and the astrometric difference between Hipparcos (Perryman et al. 1997;van Leeuwen 2007) and Gaia (Gaia Collaboration et al. 2016, 2018 for nearby stars. Various groups have used this approach to estimate the dynamical mass of directly imaged brown dwarfs (Snellen & Brown 2018;Brandt et al. 2019;Kervella et al. 2019;Xuan et al. 2020;Kiefer et al. 2021;Kervella et al. 2022).
We follow the approach developed by Feng et al. (2019a) to use both the position and proper motion differences between Gaia and Hipparcos to constrain the orbits of long period companions. Because the positional difference between Gaia and Hipparcos after subtracting the linear stellar motion is proportional to the square of Gaia-Hipparcos time difference, it is more sensitive to companion-induced acceleration of primary star than the proper motion difference, which is a linear function of time. In other words, where ∆r is the amplitude of the positional change, ∆µ is the proper motion change, ∆t is the difference between the reference times of Gaia EDR3 and Hipparcos catalog, g is the acceleration of the primary star induced by a companion. The combined analysis of both proper motion and positional differences is found to be optimal (Feng et al. 2021) through a comparison of different approaches (e.g., Brandt et al. 2021b) to the Gaia data when considering a sample of low mass companions similar to those considered here. Although a global calibration can remove systematics to some extent a priori (e.g., Cantat-Gaudin & Brandt 2021), we prefer using astrometric jitters and offsets to model the known and unknown systematics a posteriori in order to avoid over-fitting or under-fitting problems (e.g., Foreman-Mackey et al. 2015 andFeng et al. 2016). This paper is structured as follows. The RV and astrometry data are introduced in section 2. The combined modeling of RV and astrometry is described in section 3. The BDs that are discovered and confirmed by our work are listed in section 4. The following section 5 explains the statistics of this sample. The detectability of this sample by the current imaging facilities is discussed in section 6. The dynamical stability of the systems are investigated in section 7. Finally, we present conclusion in section 8.

DATA
In this work, we use the RV data of the University College London Echelle Spectrograph (UCLES; Diego et al. 1990) mounted on the Anglo-Australian Telescope (AAT), the Automated Planet Finder (APF; Vogt et al. 2014) and Levy Spectrometer at the Lick Observatory, the CORALIE spectrometer  installed at the Swiss 1.2-metre Leonhard Euler Telescope at ESO's La Silla Observatory, the ELODIE spectrograph (Baranne et al. 1996) of Observatoire de Haute-Provence, the High Accuracy Radial velocity Planet Searcher (HARPS; Pepe et al. 2000) at the ESO La Silla 3.6m telescope, the HARPS for the Northern hemisphere (HARPS-N or HARPN; Cosentino et al. 2012) installed at the Italian Telescopio Nazionale Galileo (TNG), the HIRES spectrometer (Vogt et al. 1994) at the Keck observatory, the Lick Observatory Hamilton echelle spectrometer (Vogt 1987), the Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2010) installed on VLT, the Magellan Inamori Kyocera Echelle (MIKE) spectrograph (Bernstein et al. 2003) and the Carnegie Planet Finder Spectrograph (PFS; Crane et al. 2010) on the Magellan Clay Telescope, the SOPHIE spectrograph (Perruchot et al. 2008) at the 1.93m telescope of Haute-Provence Observatory, the ESO UV-visual echelle spectrograph (UVES) on the Unit Telescope 2 of the VLT array, and the high resolution spectrograph (HRS; Tull 1998) mounted on the Hobby-Eberly Telescope (HET; Ramsey et al. 1998).
The HARPS data is reduced by Trifonov et al. (2020), using the SERVAL pipeline (Zechmeister et al. 2018). There is a known offset in the RV zero point for the post-2015 dataset (Lo Curto et al. 2015). Hence we label the pre-2015 data set as "HARPSpre", and the post-2015 data as "HARPSpost". The AAT, APF, MIKE, PFS, and UVES data are reduced using the pipeline developed by Butler et al. (1996). The APF data for HD 182488 (or GJ 758) is published by Bowler et al. (2018). This data set is denoted by "APF1" while the other archived APF data is labeled "APF2". For β Pic (or HD 39060), the RV data reduced by Lagrange et al. (2019) is used and labeled "AL19". We use the published RV data from the Lick Hamilton spectrograph and label the versions due to various updates by Lick6, Lick8, and Lick13 (Fischer et al. 2014). Since the first operation of CORALIE at 1998, it had major upgrades in 2007 (Ségransan et al. 2010) and in 2014. Hence we use COR98, COR07, COR14 to denote the three versions of data sets. The ELODIE data for HIP 63762 and the SOPHIE data for HIP 94931, HIP 14729, and HIP 22203 are downloaded from the SOPHIE/ELODIE archives 1 and reduced using the SERVAL pipeline (Zechmeister et al. 2018). We also use the RVs for HD 10697, HD 136118, HD 190228, HD 23596, HD 28185, HD 38529, HD 72659, and HD 95128 measured by the 2.7 m Harlan J. Smith Telescope (HJS) and/or HRS at the McDonald Observatory (Wittenmyer et al. 2009). For HD 14067, we use the data published by Wang et al. (2014), including: the data from the High Dispersion Spectrograph (HDS; Noguchi et al. 2002) installed on the Subaru telescope, and the RV data measured by the High Dispersion Echelle Spectrograph (HIDES) at the Okayama Astrophysical Observatory (OAO), and the data from the High Resolution Spectrograph attached at the Cassegrain focus of the 2.16 m telescope at Xinglong Observatory (XINGLONG). For HD139357, we use the RVs measured by the coudééchelle spectrograph mounted on the 2 m Alfred Jensch Telescope (AJT) of the Thueringer Landessternwarte Tautenburg (Döllinger et al. 2009). For HD 106515A, we use the RV data measured by the high-resolution spectrograph SARG at TNG (Desidera et al. 2012). The new RV data for all targets are shown in the figures in the appendix.
For a given target with both RV and revised Hipparcos catalog data (van Leeuwen 2007), we use the ga-iadr2.tmass best neighbour cross-matching catalog in the Gaia data archive to find the Gaia DR2 source identity and use the gaiaedr3.dr2 neighbourhood cross-matching catalog to find the EDR3 data 2 . For a target without Gaia counterparts in the cross-matching catalog, we select the Gaia sources within 0.1 degree from its Hipparcos ICRS coordinates and with a parallax differing from the Hipparcos one by less than 10%. For stars with both DR2 and EDR3 data, we use the difference between the revised Hipparcos catalog and the Gaia EDR3 to constrain the orbits of companions. For stars with DR2 but without EDR3 data, we use the Hipparcos-DR2 difference. The Hipparcos and Gaia data used in this work are shown in Table 2. The stellar mass for each star is from the TESS input catalog (Stassun et al. 2019) unless it can be derived from the combined analyses of RV, astrometric and imaging data.
To compare the significance of companion-induced proper motion (∆µ) and positional (∆r) differences between Gaia and Hipparcos, we calculate the signal-to-noise ratios (SNRs) for ∆µ and ∆r as follows: In the above equations, the proper motion (∆µ) and positional (∆r) differences and their uncertainties (σ µ and σ r ) as well as the error-weighted proper motion (µ α and µ δ ) are where {α G , δ G , µ G α , µ G δ } and {α H , δ H , µ H α , µ H δ } are respectively the Gaia and Hipparcos astrometry data, including right ascension (RA), declination (DEC), and proper motion in RA and DEC. The uncertainty of the astrometry data is denoted by σ with corresponding subscripts. The correlation between various astrometric parameters are not considered here although they are considered in our full modeling of astrometric data. The SNR µ and SNR r for the stellar sample in this work will be presented in section 4. In the calculation of positional difference, the linear motion due to proper motion is subtracted. Although the perspective acceleration is not considered in the calculation of SNR, it is taken into account in our rigorous modeling of astrometry that will be introduced in the following section.

METHOD
The past and current RV surveys provide a great legacy for the discovery of CJs, BDs, and low mass stellar companions. While RV alone is unable to give the true mass of a companion, astrometry data can break the degeneracy between mass and inclination and fully constrain the mass and orbit. As a successor of the Hipparcos astrometric survey, the Gaia Early Data Release (EDR3) is a third epoch release of astrometric data for more than 1 billion stars. The comparison between the Hipparcos and Gaia catalogs for common stars provides an additional constraint on potential accelerations induced by companions on primary stars.
Rather than just relying on Gaia and Hipparcos proper motions, the frequently used approach is to compare Gaia and Hipparcos positions to derive a third proper motion to calibrate the Gaia and Hipparcos catalogs, as done by Michalik et al. (2015) and Brandt (2018). Nevertheless, such a third proper motion is still biased for targets hosting massive companions with orbital periods comparable or longer than 24 years. Instead of calibrating Gaia and Hipparcos catalogs a priori, we fit the calibration parameters and signal-model parameters simultaneously to the raw catalog data to avoid overfitting problems caused by conducting calibration before signal search (Foreman-Mackey et al. 2015;Feng et al. 2017). This method is optimal compared with the previous methods in terms of constraining the masses and orbits of small planets according to Feng et al. (2021).
Considering that the models and numerical methods are introduced by Feng et al. (2019a) and Feng et al. (2021), we only briefly describe the methodology in this section. The radial velocity model consists of Keplerian components and red noise components that are modeled using the moving average model (MA; Tuomi et al. 2013 andFeng et al. 2016). We select the optimal order of the MA model in the Bayesian framework. Specifically, we calculate the maximum likelihood, L max , for the q th order MA model (or MA(q)) to derive the Bayesian information criterion (BIC) of a model, i.e. (q + 1) ln n − 2 ln L max , where n is the number of data points. We select the optimal order q opt if BIC of MA(q opt − 1) relative to MA(q opt ) is higher than 10 and the BIC of MA(q opt ) relative to MA(q opt + 1) is lower than 10. A detailed description of noise model comparison is given by Feng et al. (2017).
Reflex motion describes the orbital motion of a primary star around the barycenter of the system. The full astrometric model for a star consists of three components: stellar reflex motion, proper motion and parallax of the barycenter of the target system. Though Gaia synthetic data can be generated by GOST 3 , we cannot access the real intermediate data for a reliable detection of short period companions. Because our aim is to detect long period companions, we treat the position (α, δ) and proper motion (µ α , µ δ ) at the reference epochs of Gaia and Hipparcos as the instantaneous astrometry data to be modeled by the combination of the motion of the target system barycenter and the reflex motion. Considering a typical systematic RV of 10 au yr −1 for the barycenter of a system, the change of parallax is about 0.001/d , where d is the heliocentric distance of the target system in pc. For a star 10 pc away from the Sun, the parallax change over 24 years is about 0.01 mas, which is far below the current precision of Gaia parallaxes. Thus we ignore the change of parallax.
For systems with both long-period (P>1000 d) and short-period planets (P<1000 d), we only use the Gaia-Hipparcos astrometry to constrain the long period signals while leaving the inclination I and longitude of ascending node Ω of the short-period companion unconstrained. As is described in Feng et al. (2019a), we model the instantaneous proper motion and position of the target star at the Gaia epoch by combining the barycenter proper motion and the stellar reflex motion at the Gaia reference epoch. We then propagate the proper motion and position of the barycenter to the Hipparcos epoch, and add the stellar reflex motion at the Hipparcos reference epoch to model the proper motion and position of the target star at the Hipparcos epoch. In the calculation of likelihood, we use offsets and jitters to account for unknown systematics in the Gaia and Hipparcos catalog data.
To obtain posterior samples, we use the adaptive and parallel Markov Chain Monte Carlo (MCMC) developed by Haario et al. (2001) and Feng et al. (2019b). All time-related parameters such as orbital periods and correlation time scale follow a log-uniform prior distribution. The inclination I is sampled from a uniform prior distribution over sin I. Note-The column of "LPS" shows the number of long period signals. The number of companions is given for two difference cases: the conservative case in brackets defined by <20% relative mass error (σm c /mc < 20%) and the optimal case defined by <100% relative mass error (σm c /mc < 100%).
The other parameters have a uniform prior distribution. We use BIC> 10 or its equivalent ln(Bayes Factor)> 5 (Kass & Raftery 1995;Feng et al. 2016) to determine whether additional companions are necessary to explain the data. To find signals efficiently, we first conduct an RV-only analysis following the method developed by Feng et al. (2020a). Then we apply the full modeling of RV and astrometry to the systems with P > 1000 d signals to fully constrain the dynamics of the system. For companions with shorter orbital periods (P < 1000 d), the Hipparcos and Gaia EDR3 astrometry significantly deviate from the instantaneous astrometry at reference epochs. Because these short-period companions do not induce significant astrometric signals, they are only constrained by the RV data. For a system with both short and long period companions, the short period companions are constrained only by the RV data while the long period ones are constrained both by RV and by astrometric data. Therefore, the combined analysis of both the RV and the astrometric data for a multi-companion system would give reliable orbital solutions for all companions.

RESULTS OF COMBINED ANALYSES
Among all of the available RV targets from various RV surveys, we select 5108 stars with each star having more than five high precision RV data points. Following the methodology above (with a ln(Bayes Factor)> 5; Kass & Raftery 1995;Feng et al. 2016) and using the RV data alone we find 869 of these stars show long period signals (LPSs; P > 1000 d). By our combined analyses of RV and astrometry, 167 of them are confirmed as companions with masses from 5 to 120 M Jup and with relative mass uncertainty of less than 100%.
The stellar mass and astrometry data used for the sample of 161 stars that host the 167 companions are shown in Table 2. For the companions which are directly imaged, the mass of their hosts are inferred together with other parameters a posteriori. The Hertzsprung-Russell diagram for the sample of stars based on the Gaia BP-RP color is shown in Fig. 1. Most of the hosts of companions are main-sequence AFGKM stars while few hosts are sub-giants.
We also show the distribution of stellar mass, G magnitude, number of RV points, effective temperature, SNR r , SNR µ , and heliocentric distance in Fig. 2. In the figure, we divide the total sample of 161 companions into the known ones found in literature and the new ones found in this work. The major difference between these two population is that the new companions are typically identified from fewer RV data points than the known companions. This is expected because the the astrometric data help to constrain the long period signals that the RV data alone cannot confirm. As shown in Fig. 2, the majority of our sample are brigher than G=8 mag and less than 100 pc away from the Sun. It is apparent that the SNR of positional difference is one order of magnitude higher than the SNR of proper motion difference.
Around the 161 hosts of cold sub-stellar companions with a mass from 5 to 120 M Jup and an orbital period longer than 1000 d, we also find 63 other types of companions, including 60 planets, 1 sub-stellar companion and 2 stellar companions. To count the number of different types of multi-companion systems, we define the mass ranges for planets, BDs, and stars to be < 13 M Jup , [13, 75] M Jup , and > 75 M Jup , respectively. The number of stars and different types of companions are shown in Table 1. Our sample of cold giants is 10 times larger than the current sample of 17 companions with parameters estimated to a similar precision, as shown in Fig. 3.
The fitting results for targets with direct imaging data are shown in Fig. 4 and 5 while example fittings for targets without imaging data are shown in Fig. 6 and 7. For a target with imaging data and shown in Fig. 4 and 5, the first panel from left to right shows the optimal fit to the RV data. The second and third panels respectively shows the fit of the reflex motion induced by wide-orbit companions (P > 1000 d) to the proper motions and positions at the Hipparcos and Gaia epoch after subtracting the barycentric motion. The fourth panel shows the companion's binary orbit to its relative astrometry derived from imaging data. For targets without imaging data and shown in Fig. 6 and 7, the fourth panel is not shown. The fit to the Gaia-Hipparcos difference for a target is based on a prediction of the reflex motion with the optimal parameters over a time-span ranging from the Hipparcos epoch to the Gaia epoch. For a system consist of multiple wide-orbit companions, the reflex motion could be complex (e.g. HD 7449 in Fig. 5, GJ 676 A and HD 105811 in Fig. 7). We present Hipparcos and Gaia EDR3 astrometry for the stars in Table 2 and the orbital parameters for stellar planetary companions based on combined RV and astrometric analyses in Table 3. Because this catalog records dynamical mass and orbital parameters through combined analyses of RV and astrometry data, it is not directly comparable with proper motion anomaly catalogs that are more comprehensive but do not break the degeneracy between companion mass and orbital parameters (e.g., Kervella et al. 2022). We proceed to discuss the individual targets that are found to host companions in previous studies.    Figure 3. The confirmed companions and detected companions found in this work, successively compared with previous companions with relative mass uncertainty less than 100%. The left panel shows the distribution over mass and semi-major axis while the right panel shows the distribution over companion-host mass ratio and semi-major axis. The companions are selected from this work and from the catalog from exoplanet.eu (denoted by black squares and labelled "Literature") if their relative mass uncertainties are less than 100%, inclinations and eccentricities are estimated, host-star masses are higher than 0.2 M , orbital periods are longer than 1000 days. The mass ratio for some previous companions are missing in the right panel plot because their host-star masses are not given.  (2) -2.02(2) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · Note- Table 2 is published in its entirety in the electronic edition of the Astrophysical Journal Supplement Series. A portion is shown here for guidance regarding its form and content.
Note-The masses of stars are from the TESS input catalog (Stassun et al. 2019). The superscripts "gaia" and "hip" are respectively used to denote the Gaia EDR3 and Hipparcos right ascension (α), declination (δ), parallax ( ), proper motion in right ascension (µα), and the proper motion in declination (µ δ ). The hosts with † superscript have been directly imaged and the stellar mass is infered from the relative astrometry a posteriori. For a star without EDR3 but with DR2 data, a star symbol is added behind its name. · · · · · · · · · · · · · · · · · · · · · · · · · · · Table 3 continued on next page  Table 3 is published in its entirety in the electronic edition of the Astrophysical Journal Supplement Series.
A portion is shown here for guidance regarding its form and content. The systems are sorted by the host names listed in the first column while the companions in each system are sorted by their orbital periods. For companions with orbital period less than 1000 days, the inclination is not given and the mass reported in this table is mc sin I. For companions with mass higher than 75 MJup , we use capital letters to label them. Among the parameters, directly inferred parameters are P (orbital period), K (RV semi-amplitude), e (eccentricity), ω (argument of periastron), M0 (mean anomaly at the minimum epoch of RV data), I (inclination), Ω (longitude of ascending node). The derived parameters are Tp (periastron epoch), mc (companion mass), and a (semi-major axis). The dynamical stability of Earth-like planets in the habitable zones of their host stars are determined numerically, as described in section 7. The median and the 1-σ quantiles (i.e. 16% and 84% quantiles) are used to measure the uncertainty of each parameter. The companions with † superscript have been imaged and their relative astrometry data are analyzed in combination with other types of data in the solutions. The targets with superscript have wide binary companions that are identified by El-Badry et al.    For stars with multiple companions, only the astrometry fits for companions with periods longer than 1000 days are shown because the contribution of short-period companions to astrometry is minor. The right panels show the orbits of all companions in a system as well as the relative astrometry derived from imaging data.
For previously published planets, we discuss our new results on an object-by-object basis as follows.
•       • GJ 494 (or Ross 458) is an M-dwarf binary (Heintz 1994;Beuzit et al. 2004) hosting a young BD known as GJ 494 c which has a mass of 11.3±4.5 M Jup (Dupuy & Kraus 2013). Because GJ 494 c is on an extremely wide orbit, it does not change the motion of the inner binary much and is thus not considered in our solution. Based on combined analyses of the RV and Gaia-Hipparcos data as well as the relative astrometry data from Mann et al.  (Heintz 1994;Beuzit et al. 2004).
• GJ 676 is a binary and the primary hosts four planets (Forveille et al. 2011;Anglada-Escudé & Butler 2012 • HD 100939 (HIP 56640) is a K-type giant hosting a companion with a minimum mass of 3.67 ± 0.14 M Jup ). Our combined analyses estimate a dynamical mass of 5.30 +3.309 −0.00 M Jup , putting it into the super-Jupiter category.
• HD 106515 A (GJ 9398) is solar type star hosting a planet with a minimum mass of 9.33 ± 0.16 M Jup (Mayor et al. 2011;Desidera et al. 2012). Our combined analyses give a dynamical mass of 9.52 +6.39 −0.13 M Jup .
• HD 10697 (109 Psc • HD 11506 is a G0V-type star hosting a warm Saturn and a CJ (Fischer et al. 2007). Thanks to the long baseline between Hipparcos and Gaia, we detect an extremely cold super-Jupiter, HD 11506 d, with a period of 40.31 +7.7 −7.5 yr and a mass of 7.38 +2.02 −1.09 M Jup . This makes the system unique in terms of a solar analog hosting a warm Saturn and two CJs.
• HD 120084 is a G star hosting a companion with a minimum mass of 4.5 M Jup (Sato et al. 2013 • HD 125612 is a G star hosting three companions with minimum masses of 3.0, 0.058, and 7.2 M Jup and with orbital periods of 502, 4.1547, and 3008 d, respectively (Fischer et al. 2007;Lo Curto et al. 2010;Ment et al. 2018). Through combined analyses of RV and astrometric data, we are able to constrain the mass of the Jupiter-like planet, HD 125612 d, to be 7.18 +0.93 −0.45 M Jup .
• HD 126614 is a binary consisting of a G-type primary star and a M dwarf companion with the G star hosting a planet with a minimum mass of 0.38±0.04 M Jup (Howard et al. 2010;Stassun et al. 2017 (2010) while a minimum mass of 36 M Jup is given by Halbwachs et al. (2000). We estimate a mass of 47.41 +3.70 −3.53 M Jup , much more precise than previous estimations. This system also hosts a Neptune-sized planet with a minimum mass of about 27 M ⊕ . Its inclination is consistent with the inclination of HD 127506 B, supporting a coplanar architecture.
• HD 129191 (BD-04 37333 or HIP 71803) is a star hosting a companion candidate with a minimum mass of 6.8 M Jup (Hinkel et al. 2019). The use of Gaia-Hipparcos data allow us to break the degeneracy between mass and inclination and confirm this candidate as a BD with a mass of 76.06 +16.22 −14.60 M Jup .
• HD 13724 (HIP 10278) is a G-type star hosting HD 13724 b or (HD 13724 B), a brown dwarf with a mass of 50.5 +3.3 −3.5 M Jup (Rickman et al. 2020). We estimate a mass of 36.32 +1.48 −1.60 M Jup after considering the 10% uncertainty of the mass of the primary star. Without using direct imaging data, we constrain the BD mass to a precision similar to that given by (Rickman et al. 2020) (labeled "ER20") based on combined RV and direct imaging data analysis. We also find a warm Saturn with a minimum mass of about 0.2 M Jup or 72 M ⊕ . Without a well constrained inclination for HD 13724 b, we cannot determine the mutual inclination between HD 13724 b and B.
• HD 139357 is a K-type star hosting a companion with a minimum mass of 9.76±2.15 M Jup (Döllinger et al. 2009).
With combined analyses, we confirm this companion to be a BD with a mass of 16.38 +7.88 −0.00 M Jup .
• HD 14067 is a G9III star hosting a sub-stellar companion with a minimum mass of 9.0 M Jup (Wang et al. 2014). Our combined analyses constrain the mass of the companion to be 9.49 +13.30 −0.00 M Jup .
• HD 142 is a wide binary consist of a F-type star and an M dwarf (Tokovinin & Kiyaeva 2016). The primary hosts at least two planets Wittenmyer et al. 2012 −3.2 M Jup and I c = 101 +31 −33 deg, indicating a strong misalignment. Such a discrepancy may be due to the parameter degeneracy in their constraint of four orbital parameters (Ω b , Ω c , I b , and I c ) by using only the proper motion difference between Gaia and Hipparcos (equivalent to 2 data points).
• HD 145934 hosts a Jupiter analog with a minimum mass of 2.28±0.6 M Jup (Feng et al. 2015). We confirm the outer companion identified by Feng et al. (2015) as a (sub-)stellar companion, HD 145934 B, which is on a wide orbit and has a mass of 87.87 +70.23 −13.93 M Jup . Though well separated from the primary star, this companion cannot be resolved by Gaia due to the small parallax of the system (about 4.36 mas).
• HD 167665 (HIP 89620) is an F-type star hosting a BD (Patel et al. 2007 • HD 214823 (HIP 111928) is a G-type star hosting a companion (Díaz et al. 2016a;Ment et al. 2018) with a minimum mass of 20.56±0.32 M Jup and an orbital period of 1853.9±1.6 d (Luhn et al. 2019). With Gaia and Hipparcos data, we are able to break the degeneracy between mass and inclination, and to constrain the absolute mass to 18.61 +4.14 −1.07 M Jup and the period to be 5.078 +0.004 −0.004 yr.
• HD 217786 (HIP 113834) is a F star hosting a companion with a minimum mass of 13±0.8 M Jup (Moutou et al. 2011) and a stellar companion on an extremely wide orbit El-Badry et al. 2021). Because the RV and astrometric variation of the primary star caused by the wide stellar companion is insignificant, we only model the reflex motion due to the substellar companion and estimate a mass of 13.85 +1.27 −1.31 M Jup . We also find a hot super-Earth with an orbital period of 2.5 d. • HD 27894 (HIP 20277) is a K star hosting three companions with masses of 5.42, 0.16, and 0.67 M Jup (Moutou et al. 2005;Anglada-Escudé et al. 2010;Kürster et al. 2015;Trifonov et al. 2017). However, we only find strong evidence for the biggest two companions. By constraining the inclination of the biggest companion using astrometry, we find a dynamical mass of 6.49 +0.99 −0.35 M Jup .
• HD 28185 (HIP 20723) is a G-type star hosting a companion (Santos et al. 2001;Wittenmyer et al. 2009) with a minimum mass of 6.7 M Jup and an orbital period of 379±2 d (Minniti et al. 2009). With both RV and astrometry data, we constrain the mass of HD 28185 b to be 7.07 +1.29 −0.79 M Jup . Moreover, we detect a BD companion with a mass of 19.64 +2.27 −2.14 M Jup .
• HD 28192 (HIP 20752) is a G star hosting a stellar companion on an extremely wide orbit and with a mass of about 0.4 M (Tokovinin 2014b;El-Badry et al. 2021). This companion at most induces a reflex motion of 0.1 m s −1 yr −1 and is thus not considered in our combined analyses. We find another stellar companion with a mass of 94.78 +8.51 −7.52 M Jup on a 10 au-wide orbit and a planet with a minimum mass of 0.31 +0.02 −0.03 M Jup and an orbital period of about 14 d.
• HD 29461 (HIP 21654) is a G star hosting a companion with a minimum mass of about 0.08 M (Griffin 2012;Bouchy et al. 2016). Our combined analyses estimate a mass of 92.96 +12.61 −4.93 M Jup on a 5 au-wide orbit. The small separation between this companion and the primary star explains null detection of it by previous imaging surveys.
• HD 30177 (HIP 21850) is a G star hosting two companions with minimum masses of 3±0.3 and 8.07±0.12 M Jup Wittenmyer et al. 2017;Barbato et al. 2018). The null detection of these companions in the direct imaging survey conducted by Zurlo et al. (2018) leads to an upper limit of 28-30 M Jup and a minimum inclination of 15 • . Our combined analyses estimate masses of 8.40 +1.24 −0.49 M Jup and 6.15 +1.31 −0.34 M Jup , consistent with the constraints given by previous imaging and RV data analyses. This system is one of few systems hosting two super-Jupiters in our sample.
• HD 39060 (β Pic) is a well-studied system, hosting two giant planets (Lagrange et al. 2009(Lagrange et al. , 2019. We present the parameters of β Pic b and c given by the most recent studies in Table 4. Thanks to the valuable data collected by previous studies, we analyze the relative RV data for b (Snellen et al. 2014), the updated relative astrometry data from Lacour et al. (2021), the recently released Gaia EDR3, and the RV data used by Lagrange et al. (2020) and Vandal et al. (2020) with different reprocessing to remove stellar activity noise. We model the relative astrometry for multiple companions by (1) calculating the reflex motion of the host due to the innermost companion; (2) calculating the reflex motion of the barycenter of the host and the innermost companion due to the outer companion; (3) repeating the above steps until the reflex motion due to all companions are modeled; (4) calculating the position of an outer companion relative to the host star by converting the barycenter of the host and inner companions to the host position. This procedure is introduced by Lacour et al. (2021) in detail and is also used by Brandt et al. (2021a).
We find that our solution is quite sensitive to the RV data. Although both the RV datasets used are from the same HARPS observations different corrections are made to stellar activity particularly from pulsations. We use the Lagrange et al. (2020) dataset (dubbed "AL20") to find β Pic b to be 7.56 +1.35 −1.69 M Jup and that of c to be 8.94 +0.75 −0.78 M Jup while using Vandal et al. (2020) dataset (dubbed "TV20") we find a mass of 11.75 +2.34 −2.15 M Jup and 10.15 +1.20 −1.07 . In Table 4 we present our solutions along with the range of solutions from the literature and note the considerable scatter and discrepancy from the masses of 3.2 M Jup and 5.6±1.5 M Jup measured by Lagrange et al. (2020) and Nowak et al. (2020) respectively for β Pic b using the AL20 RVs when they adopt an uninformative prior. For Table 3 we adopt our TV20 solution since it agrees better with the independent astrometric solution for β Pic c based on interferometric data by Lacour et al. (2021) and provides a solution more consistent with the higher mass predicted by various cooling models (e.g. Baraffe et al. 2003, Spiegel & Burrows 2012. We can anticipate that further evolution in the processing of the RV data for stellar activity as well as the incorporation of Gaia intermediate data into analysis solutions will be useful to resolve the discrepancies that we find in mass measurement from different RV data sets. • HD 39091 (π Men) is a G star hosting at least two companions with masses of 0.015 and 13 M Jup (Jones et al. 2002;Huang et al. 2018;Damasso et al. 2020b) and a possible third companion found recently by Hatzes et al. (2022). The smaller companion is a transiting planet while the bigger companion is found to be significantly misaligned with the inner one (Damasso et al. 2020b;De Rosa et al. 2020;Xuan et al. 2020;Kunovac Hodžić et al. 2021). With combined analyses of the RV data from AAT, CORALIE, ESPRESSO, HARPS, and PFS as well as the Gaia-Hipparcos data, we are able to constrain the mass of this companion to be 12.33 +1.19 −1.38 M Jup and  The solution shown in this table is the so-called "coplanar fit" by EN20 who adopt a Gaussian prior centered on zero and with a standard deviation of 1 • . Without such assumption, the inclination of c has an error of 13 • , leading to significant mass uncertainty for c. c In AL20, the inclination of β Pic c is assumed to be equal to β Pic b because the relative astrometry of β Pic c is not used to constrain it. The mass of b reported here is estimated by AL20 using a Gaussian prior of 14 ± 1 MJup . Without such informative prior, the mass is 3.2 MJup . d In TV20, the stellar mass is given a priori by Wang et al. (2016).
e The mass of b reported here is estimated by MN20 using a Gaussian prior of 15 ± 3 MJup . Without such informative prior, the mass is 5.6±1.5 MJup . f HDR2 represents the data of proper motion difference beween Hipparcos and Gaia DR2.
g HEDR3 represents the data of both proper motion and positional difference between Hipparcos and Gaia EDR3. the inclination to be 54.44 +5.94 −3.72 deg. This inclination of π Men b differs from the 90 deg inclination of the π Men c by 6σ, suggesting significant misalignment as proposed by previous studies. However, we fail to confirm the third companion found by Hatzes et al. (2022).
• HD 39213 (HIP 27491) is a K star hosting a BD with a minimum mass of 0.07±0.01 M (Jenkins et al. 2015). Our combined analyses determine a dynamical mass of 70.77 +10.21 −5.43 M Jup , putting it around the boundary between BD and stellar object. Follow-up direct imaging of this object is needed to further characterize the companion.
• HD 4113 A (HIP 3391) is a Sun-like star hosting a BD (HD 4113 A b or HD 4113 b) and a sub-stellar companion (HD 4113 C) (Tamuz et al. 2008;Cheetham et al. 2018) as well as HD 4113 B, an M dwarf on an extremely wide orbit (Mugrauer et al. 2014). HD 4113 A b has a minimum mass of 1.602 +0.076 −0.075 M Jup based on RV analysis (Cheetham et al. 2018). HD 4113 C has a dynamical mass of 65.8 +5.0 −4.4 M Jup based on analyses of direct imaging data (labeled by "AC18") and an isochronal mass of 36±5 M Jup based on cooling models (Cheetham et al. 2018). Our combined analysis of the RV, astrometry, and imaging data constrains the mass of HD 4113 C to be 51.91 +0.60 −0.46 M Jup , relaxing the previous tension between dynamical and isochronal masses without invoking binarity in HD 4113 C.
• HD 42581 (GJ 229 A) hosts the first imaged BD, GJ 229 B (Nakajima et al. 1995) as well as two planets (Feng et al. 2020a). Recently, Brandt et al. (2021b) estimated a mass of 71.4 ± 0.6 M Jup for GJ 229 B based on combined analyses of RV, imaging data (labeled "MB21"), Gaia EDR3 and Hipparcos data. This mass is in tension with the mass predicted by cooling models. However, with nearly the same data but with additional constraint from Hipparcos-Gaia positional difference, our combined analyses estimate a mass of 60.42 +2.34 −2.38 M Jup , consistent with the 64.8 ± 0.1 M Jup predicted by cooling models (Brandt et al. 2021b). This suggests that the use of positional difference between Hipparcos and Gaia might be important to avoid potential bias by using proper motion difference alone.
• HD 43197 is a Sun-like star hosting a warm Jupiter (Naef et al. 2010). We detect a cold super-Jupiter, HD 43197 c, with a mass of 7.9 ± 1.7 M Jup on a wide orbit with a period of 27±9 yr and a nearly face-on inclination (11.42 +5.39 −3.07 deg). However, the inclination of the inner companion is not well constrained. Assuming a coplanar configuration, HD 43197 b would have a mass of about 4 M Jup .
• HD 65430 (HIP 39064) is a spectroscopic binary with an orbital period of 3138 d (Allen et al. 2012). Our combined analyses estimate a mass of 105.40 +8.37 −8.95 M Jup , confirming its stellar origin.
• HD 66428 b and c were discovered by Butler et al. (2006) and Rosenthal et al. (2021), respectively. In addition to the Keck data used by Rosenthal et al. (2021), we analyze the HARPS data reduced by Trifonov et al. (2020) and • HD 72659 (HIP 42030) is a G-type star hosting a companion Wittenmyer et al. 2009 • HD 72892 is a Sun-like star hosting a super-Jupiter (Jenkins et al. 2017). In addition to this companion, we also detect another companion HD 72892 B with a mass of 77.12 +41.76 −35.48 M Jup on an edge-on and eccentric orbit. The high eccentricity of HD 72892 b (e = 0.419 ± 0.003) might be caused by the strong perturbations from HD 72892 B, which is on an orbit with an eccentricity of 0.38±0.06.
• HD 73267 is a solar type star hosting a Jupiter-like planet, HD 73267 b (Moutou et al. 2009). Through combined RV and astrometry data analyses, we identify an additional companion named "HD 73267 c". It is a super-Jupiter with a mass of 5.13 +0.91 −0.28 M Jup and with an orbital period of 46.74 +2.15 −2.98 yr. The orbits of the two companions are probably misaligned though with large uncertainty.
• HD 74014 (HIP 42634) is a star hosting a BD companion (Patel et al. 2007  • HD 81040 (HIP 46076) is a G star hosting a Jupiter-like companion with a mass of 7.24 +1.0 −0.37 M Jup (Stassun et al. 2017;Sozzetti et al. 2006;Li et al. 2021). Our combined analyses give a mass of 6.77 +1.10 −0.87 M Jup . Though we use the same RV data sets as Li et al. (2021), we model the correlated RV noise using the MA(1) model. This makes our estimation of the dynamical mass more uncertain but more conservative than the values given by Li et al. (2021).
• HD 81817 is a K-type star hosting a substellar companion (HD 81817 b) with a minimum mass of 27.1 M Jup (Bang et al. 2020). With both RV and astrometry data, we are able to constrain its mass to 24.13 +9.83 −0.71 M Jup . We also find another BD in this system (HD 81817 c) although it was diagnosed as an activity signal by Bang et al. (2020) due to a dubious overlap with powers in the periodograms of Hα. However, we confirm HD 81817 c as a BD because this signal shows a unique power (see Fig. 8) in the BFP and is strictly periodic and quite circular based on MCMC posterior samplings. The evidence strongly support a Keplerian origin instead of an activity origin though its inclination is not well constrained due to its short orbital period.
• HD 86264 (HIP 48680) is a K star hosting a companion with a minimum mass of 7 ± 1.6 M Jup (Fischer et al. 2009). According to the solution based on our combined analyses, the dynamical mass of the companion is 9.81 +11.71 −1.95 M Jup .
• HD 8673 (HIP 6702) is a double star system including a F star and an early M-dwarf with a mass of 0.33-0.45 M (Roberts et al. 2015b). The primary F star hosts a sub-stellar companion with a minimum mass of 14.2±1.6 M Jup . Our combined analyses constrain the mass of the sub-stellar companion to be 13.25 +1.70 −1.42 M Jup .
• HD 87883 (HIP 49699) is a K star hosting a companion with a minimum mass of 6.31 +0.31 −0.32 M Jup (Fischer et al. 2009;Stassun et al. 2017;Li et al. 2021). Our solution estimates a mass of 5.3±0.7 M Jup . Compared with the bi-modal posterior distribution of inclination given by Li et al. (2021) based on RV and Gaia-Hipparcos proper motion difference, we use both proper motion and position differences between Gaia and Hipparcos so that we are able to break the degeneracy between I and I + π, and to constrain the inclination to be 25.45 +1.61 −1.05 deg.
• HD 95544 (HIP 54203) is a G star hosting a companion with a minimum mass of 6.84±0.31 (Demangeon et al. 2021). According to our combined analyses of RV and astrometric data, the dynamical mass is 6.02 +1.62 −0.26 M Jup and the inclination is 86.50 +31.19 −24.60 deg, indicating an edge-on configuration.
• HD 984 (or HIP 1134) is an F star hosting a brown dwarf with a mass of 61.0 ± 4.0 M Jup (Johnson-Groh et al. 2017b;Franson et al. 2022). Through combined analyses of the RV, Gaia-Hipparcos, and the imaging data collected by Franson et al. (2022) and is labeled "KF22", we find a dynamical mass of 40.37 +24.33 −18.27 M Jup .
• HD 98649 (HIP 55409) is a G type star hosting a companion with a minimum mass of 6.79 +0.5 −0.3 M Jup (Marmier et al. 2013;Rickman et al. 2019). By combined analyses of RV and Gaia-Hipparcos astrometry, Li et al. (2021) estimate a mass of 9.7 +2.3 −1.9 M Jup . Using both proper motion and position differences between Gaia and Hipparcos, we estimate a mass of 6.76 +3.61 −0.00 M Jup , consistent with and more precise than the mass given by Li et al. (2021).
• HIP 22203 (HD 30246) is a G star hosting a brown dwarf with a minimum mass of 55.1 +20.3 Although the inclination is not well constrained, the companion orbit is unlikely face-on as Kiefer et al. (2021) conclude. Considering that the astrometric signal of this companion is insignificant, we cannot be sure about whether this companion is sub-stellar or stellar.
• HIP 67537 (HD 120457) is a red giant branch star hosting a planet with a minimum masses of 11.1 +0.4 −1.1 M Jup . Our combined analyses constrain the companion's mass to be 10.88 +7.78 −0.00 M Jup .
• HIP 67851 (HD 121056) is a K-type giant star hosting two companions with minimum masses of 5.98±0.76 M Jup and 1.38±0.15 M Jup (Jones et al. 2015b,c;Wittenmyer et al. 2015). Our combined analyses constrain the mass of the bigger companion to be 6.94 +2.06 −0.52 M Jup .
• HIP 78395 (WDS 16003-0148) is a K star hosting two companions (Mason et al. 2001). By combined analyses of the RV, Gaia-Hipparcos astrometry, and the relative astrometric data provided by Mason et al. (2001), we constrain the mass to be 68.09 +8.65 −8.06 M Jup , putting the companion around the boundary between sub-stellar and stellar categories.
• HIP 97233 (HD 186641) is a K star hosting a companion with a minimum mass of 20±0.4 M Jup (Jones et al. 2015c). Based on our analyses, the mass of this companion is 19.19 +3.67 −0.32 M Jup and is on a nearly edge-on orbit.

Mass distribution and occurrence rate
Through our analysis of the 5108 stars with each star having more than 5 high precision RV data point, we find 869 stars with 914 long period signals (>1000 d). Of these 167 of them are confirmed as companions with masses from 5 to 120 M Jup by our combined analyses of RV and astrometry. The relative mass uncertainty of this sample is less than 100%. The masses of 113 companions are constrained to a precision of better than 20%. Without correcting for detection bias, the occurrence rate of the wide-orbit BDs is about 1.3%, consistent with previous estimation (e.g. Grieves et al. 2017;Kiefer et al. 2019).
We define the sample with relative mass uncertainty less than 100% as the "optimistic sample" and the sample with relative mass error less than 20% as the "conservative sample". We show the distribution of the sample over mass and mass ratio in Fig. 9. There are at least three features seen in the mass distribution of the optimistic sample: (1) there is a lack of BDs with a mass around 40 M Jup , consistent with the so-called low-mass and high-mass BD boundary identified by Ma & Ge (2014); (2) there is also a 2-σ valley around the 75 M Jup boundary between stars and BDs; (3) a sharp decrease of companions around the 13 M Jup planet-BD boundary is followed by a shallow decrease from 13 M Jup to 40 M Jup . While the first two features remain in the conservative sample, the third feature becomes insignificant in the conservative sample. In the distribution of mass ratio (right panels of Fig. 9), we see a valley around 0.3-0.4 but fail to find any significant features around the star-BD boundary (0.07 if assuming the host mass to be unit solar mass). Because the detection bias is only significant for cold super-Jupiters, some of the features seen above are not likely to disappear after considering detection bias. In particular, the valley around 40 M Jup is robust to the choice of sample size and the normalization of companion mass. By investigating the distribution over mass (or mass ratio) and semi-major axis (Fig. 3), we observe that the 40 M Jup valley gradually disappear beyond 10 au. Nevertheless, a bias-corrected distribution of sub-stellar companions over mass and semi-major axis is necessary to confirm the above patterns in the sample. Such an investigation will be left to a subsequent study of this sample while this paper is focused on companion detection.

Multiplicity
In the sample of 161 wide-companion hosts, 61 hosts have multiple planet or stellar companions. Among them, 29 have stellar companions from the EDR3 wide-binary catalog given by El-Badry et al. (2021), 6 of them have both companions identified in this work and companions from the wide-binary catalog. Without referring to the EDR3 wide-binary catalog, there are 38 multi-companion systems. Among the 61 multi-companion systems, there are 3 systems that contain planets, BDs, and stellar companions, 12 contain planets and BDs, 21 contain planets and stellar companions, 8 contain BDs and stellar companions, 12 contain multiple planets, 1 contains multiple BDs, and 4 contain multiple stellar companions.
All multi-companion systems are shown in Fig. 10. The apparent impression is that the widest companion in a system that is wider than the Neptune's orbit tend to be stellar. This is due to the incompleteness of sub-stellar companions on extremely wide orbits (e.g. >100 au). On the other hand, the architecture of the inner system seems to be insensitive to the separation between the outer companion and the primary star. This is either due to the incompleteness of the inner companions or due to the insignificant impact of extremely wide companion on inner system.
There is a controversy over whether hosts of hot Jupiters or BDs tend to have high rates of widely separated companions (Fontanive et al. 2019;Ziegler et al. 2021;Moe & Kratter 2021). While close binaries definitely suppress S-type planets, it is unclear whether wide companions could influence the formation of inner planets significantly. Because our sample is mainly massive companions on wide orbits, we would like to assess the influence of widecompanions on inner planets. By cross-matching the EDR3-based wide binary sample given by El-Badry et al. (2021) and our sample of massive companions, we find 29 out of the 161 companion hosts have wide stellar companions, indicating a stellar multiplicity rate of at least 18±3%. Considering that only two host stars in our sample have distances larger than 200 pc, the incompleteness of the identified wide companions is mainly caused by decreasing Gaia completeness below an angular resolution of 2 (El-Badry et al. 2021).
The detection rate of a wide binary around a companion host is P(DWB,s|CH), where s is binary separation. It is derived from the occurrence rate of the wide binary around companion host P(WB,s|CH) and the incompleteness of the EDR3 wide binary sample, P(DWB,s|WB,s), according to P(DWB,s|CH)=P(DWB,s|WB,s)P(WB,s|CH). To calculate the detection rate, P(DWB,s|WB,s), we sample s from 20 to 10,000 au. For each s, we repetitively draw 100,000 samples from the parallaxesω of all companion hosts and select the ones that could be resolved by Gaia, i.e. sω > 2 . The sample becomes significantly incomplete for s < 150 au or ln(s/au) < 5. After correcting for this incompleteness, the occurrence rate of wide binaries for companion hosts or P(WB|CH) is 32±6%, consistent with the binary fraction of 36±2% for binary separation from 20 to 10,000 au (Fontanive et al. 2019). Hence we do not find any preference of multiplicity for massive companions on wide orbits. Based on the calculation of P(DWB,s|CH) and P(DWB,s|WB,s), we derive the wide-binary occurrence rate as a function of separation by using P(WB,s|CH)=P(DWB,s|CH)/P(DWB,s|WB,s). We find that the identified wide binaries approximately follow a log-normal distribution centered around ln (s/au) = 5.8 or s = 330 au. A similar peak around 250 au is also found by Fontanive et al. (2019) for wide binaries hosting hot-Jupiters. After considering the detection bias, we find a power-law distribution of P (WB, s|CH) ∝ s −1.4 to be optimal to model the distribution of the occurrence rate over binary separation. This power-law distribution over separation is consistent with the monotonic decreasing of binary fraction with separation beyond 3 au for the whole EDR3 wide binary sample (El-Badry et al. 2021). Considering that the power-law distribution is found after considering detection bias, it is probably intrinsic to the wide binary sample and thus the peak around 250 au in the separation distribution found by Fontanive et al. Figure 10. The 61 multi-companion systems identified in this work. The yellow dots represent stellar companions with mass higher than 75 MJup , the brown dots represent BDs with mass range from 13 to 75 MJup , and the blue dots represent planet companions with mass lower than 13 MJup . The semi-major axis is used as a proxy for the separation from host star for the companions identified in this work while the binary separation is used for the EDR3 wide-binary catalog given by El-Badry et al. (2021). The Solar System planets are put on the black horizontal line for reference. The masses of the Solar System planets are amplified by 10 times for better visualization. The size of dots represents companion mass and the size for stellar companions is truncated to the largest size.  Figure 11. Orbits of 16 companions with separation larger than 500 mas from their host stars. For a multi-giant system, the black line represents the giant companion on the widest orbit while the red line represents the other companions with orbital period longer than 1000 d in the system. The positions of the companions at some epochs are also denoted by the corresponding years. The black cross in each panel represents the host star.
(2019) is probably due to detection bias. A detailed analysis of the whole EDR3 wide binary sample is needed to understand the intrinsic distribution of wide binaries and is beyond the scope of this paper.

CANDIDATES FOR DIRECT IMAGING
Because the companions identified in this work are massive, nearby, and on wide orbits, they are appropriate targets for direct imaging by the currently available facilities and the next-generation instruments. We find 30 super-Jupiters and BDs with average separation from their hosts larger than 0.5 . The orbits of 16 of them are shown in Fig. 11. Considering that the ages of the host stars are typically not well constrained, we assume ages of 0.1, 1, and 10 Gyr for each star to predict the J, H, K magnitudes and the total luminosity (L) of wide super-Jupiters and BDs (less than 75 M Jup ) based on the cooling models introduced by Phillips et al. (2020). The contrast ratio between the companion and their host stars are shown in Fig. 12. For the modern coronagraphs such as SCExAO/CHARIS installed on the Subaru telescope (Jovanovic et al. 2015), the typical inner working angle is around 0.2 and the contrast limit is about 10 −6 (Currie et al. 2020). Assuming such a detection sensitivity and an age of 1 Gyr for all host stars, 41%, 35%, and 33% companions are detectable in the J, H and K bands. The proportions of detectable companions increase to 62%, 61%, and 61% for all bands if the stars have an age of 0.1 Gyr, and decrease to 16% for J band, 12% for H band, and 11% for K band if the stars have an age of 10 Gyr. Considering 0.1 and 10 Gyr as the lower and upper limits of the ages of host stars respectively, there are 10-57 sub-stellar objects detectable by the current imaging facilities.  Table 5 is published in its entirety in the electronic edition of the Astrophysical Journal Supplement Series. A portion is shown here for guidance regarding its form and content. The contrasts in different bands are denoted by "BandAge" where "Band" is H, J, or K and "Age" is 0.1, 1, or 10 Gyrs. The contrast or flux ratio is in a base-10 log scale. This table only lists the mean values of the separation ρ and the contrast in different bands for each companion. The separation of the companion to the host vary over time and the exact values for different epochs are shown in Fig. 11. The luminosity is either given by the Gaia DR2 or by the mass-luminosity relation provided by Eker et al. (2015). The J, H, and K magntiudes of each star are obtained from the Simbad database (Wenger et al. 2000). The J, H, and K magnitudes as well as the luminosity for each companion are derived by the ATMO 2020 cooling models (Phillips et al. 2020). We performed a large number of N-body simulations to study the dynamical stability of both the planets themselves, as well as the systems' habitable zones (HZs). The majority of our computations utilize the M ercury6 hybrid integrator (Chambers 1999), however we also employed a more direct, Bulirsch-Stoer methodology to accurately simulate systems with high eccentricity planets that make excessively low perihelia passages around their host stars. In general, our simulations are designed to perform a broad, first-order analysis of the long-term dynamical evolution of each planets' orbit given the uncertainties reported in Table 1. Therefore, our work should be viewed as a reasonable measure of the stability of our sample of systems (and a validation of the orbital determinations described earlier in this manuscript), rather than a comprehensive and detailed interrogation of all possible trajectories. While we consider variations within the determined values of the most dynamically significant properties of each planets' orbit (namely their eccentricities and semi-major axes), we do not investigate the possibility of perturbations from other, undetected massive bodies in each system that might perturb the orbital paths of the detected bodies. Thus, while our simulations cannot definitively prove that our sample of systems are stable, we can argue with a high degree of confidence that planets exhibiting regular behavior regardless of the orbital parameters being varied are stable.

Multi-planet system stability
We ran a series of ≥225, 1 Myr dynamical simulations to gauge the orbital stability of each multi-planet system in Table 1. In all cases we consider a grid of five eccentricities and three semi-major axes that span the range of uncertainties reported in Table 1 for each planet (angular orbital elements not determined through our orbit fitting are determined by sampling from uniform distributions of angles). Each simulation leverages a time-step of ∼5% the orbital period of the innermost planet. Through this analysis, we determine that each planet reported in this manuscript are stable in at least 90% of our numerical integrations (the unstable cases typically occur at larger eccentricities). For some multi-companion systems we performed extended (10 Myr) simulations. As an example, the two detected bodies in the HD 205158 system are plotted in Fig. 13. While the brown dwarf companion HD 205158B drives large secular oscillations in the inner Jupiter analog's eccentricity, the pair evolve on stable orbits for the duration of the simulation in spite of the system's uncharacteristically large mutual inclination.

Stability of an Earth analog in the HZ
For each system reported in Table 1, we investigated the stability of an Earth-mass planet on a nearly circular, co-planar orbit situated in the center of the HZ (as determined via relations from Kopparapu et al. 2013). Each of these models utilized the nominal orbital parameters for all planets, the Mercury package's hybrid symplectic integrator, and a total simulation time of 1.0 Gyr. Fig. 14 and Table 3 summarize the results from this batch of lengthy integrations. Unsurprisingly, systems with stable HZs tend to possess planets on longer-period orbits with lower eccentricities. Additionally, the majority of the detected companion bodies in our sample tend to destabilize the  Table 2. The color of each point identifies whether the Earth analog survived the simulation (green points), or was lost via ejection or collision (blue points). The size of each point corresponds to the mass of each planet in our various simulations.
HZ of their host stars. While these findings are by no means novel, we present these simulations to demonstrate the reasonable feasibility of habitability in the majority of our new systems.

CONCLUSION
Based on analyses of the proper motion and positional difference between Gaia and Hipparcos data as well as the RV data, we find 167 circumstellar giants in 161 systems. The occurrence rate of wide-orbit BDs is at least 1.3%, consistent with previous studies. There are 61 stars hosting multiple companions, with 12 hosting both planets and BDs, 21 hosting both planets and stellar companions, 8 hosting both BDs and stellar companions, and 3 hosting all types of companions.
Without correction of detection bias, we observe a monotonic decrease of occurrence rate with mass from 5 to 40 M Jup , a valley around 40 M Jup and 75 M Jup . The investigation of these features by reliable correction of observation bias will be left to a subsequent study of this sample. We do not find significant preference of multiplicity for wide-orbit companions with mass from 5 to 120 M Jup . For a system hosting a 5-120 M Jup companion, the probability of finding a stellar companion on a wide orbit (>3 au) is about 30±5% after accounting for detection bias.
Because our sample of cold giants are nearby and well separated from their hosts, they are good targets for direct imaging by the current facilities. By adopting the "ATMO 2020" model to predict the BD temperature and assuming different ages for the systems, we find that 10-57 super-Jupiters and BDs could be directly imaged by the current facilities. According to the exoplanet.eu archive, 161 exoplanets and BDs are imaged, 16 of them are older than 1 Gyr, 20 of them have dynamical mass, and only 4 of them are older than 1 Gyr and with dynamical mass. Because our sample of stars are selected from various RV surveys and typically show stable RV variation, most of the detectable super-Jupiters and BDs are likely to be old. If imaged by follow-up observations, the substellar companions found in this work will extend the imaged cold (or old) giants by an order of magnitude.
We conduct numerical N-body simulations to investigate the dynamical stability of the systems. Our simulations show a stability probability of at least 90% for single-companion systems over 1 Myr and for multi-companion systems over 10 Myr. Our simulations also show the majority of the massive companions identified in this work tend to destabilize the HZ of their host stars. Hence massive companions on wide and eccentric orbits are not friendly to life.
In addition to the above findings, we also expect the following application of this catalog of sub-stellar companions.
• Sub-stellar mass function: By accounting for the RV and astrometric detection bias, the occurrence rate of objects of different mass can be corrected to derive the mass function as well as the occurrence rate as a function of mass and other orbital parameters. This will quantify the boundary of BD desert in the parameter space and test various scenarios for the formation of cold giants.
• Misalignment in multi-companion system: Unlike previous RV surveys, our synergistic survey of nearby sub-stellar companions fully constrain the orbits and masses of companions. By investigating the orbital mis-alignment between different companions in a system and the correlation between misalignment and other orbital parameters, we can discover the causes of this misalignment and constrain evolutionary models of sub-stellar objects.
• Constraint of cooling models for BDs and super-Jupiters: Many of the companions in our sample could be imaged by the current facilities. Imaging data will be used to derive the luminosity, color, and effective temperature of the companions. Because the host stars of these companions are nearby and bright, their ages are typically well constrained. Assuming a co-natal and coeval formation for companions and their hosts, the age, effective temperature, and mass could be used to constrain the parameters of various cooling models such as cloud coverage.
• Correlation between cold giants and other types of planets: Because our sample significantly extends the well-constrained cold-giant sample, it makes the statistical study of the correlation between these cold giants and inner planets more robust. Unlike the previous work based on RV samples of exoplanets with only a minimum mass available, the study of planet correlations based on our sample will avoid the mass ambiguity caused by the limitation of the RV-only method.
In summary, we detect and confirm 167 sub-stellar companions on wide orbits with well constrained mass, extending the current sample of similar objects by more than one order of magnitudes. This catalog is used to study the correlation between substellar companion and wide binaries, to provide dozens of candidates for direct imaging by the current facility, and to investigate the influence of cold giants on the HZs of their hosts. Future works based on our sample would test the brown dwarf desert hypothesis and shed light on the formation and evolution of planets and sub-stellar companions.
In future, we will detect short period companions through combined analyses of Hipparcos intermediate data and the synthetic Gaia data generated by GOST. By pushing the detection limit of our combined approach towards the low mass regime, we will detect and characterize Jupiter-mass and Saturn-mass planets on wide orbits for the study of population synthesis, dynamical origin of orbital misalignment, and provide a unique sample of CJs and Saturns for direct imaging by James-Webb Space Telescope (JWST; Danielski et al. 2018) and Chinese Space Station Telescope (CSST; Gong et al. 2019).