We analyze 5108 AFGKM stars with at least five highprecision 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 days. Around these signals 167 cold giants and 68 other types of companions are identified through combined analyses of radial velocity astrometry and imaging data. Without correcting for detection bias we estimate the minimum occurrence rate of the wideorbit brown dwarfs to be 1.3 and find a significant browndwarf valley around 40 M SUBJupSUB. We also find a powerlaw 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 substellar 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 1057 substellar objects within our sample can be detected with current imaging facilities extending the imaged cold or old giants by an order of magnitude.
['Astrometric exoplanet detection', 'Radial velocity', 'Exoplanet detection methods', 'Exoplanet astronomy', 'Brown dwarfs', 'Exoplanet catalogs', '2130', '1332', '489', '486', '185', '488', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
3D Selection of 167 Substellar Companions to Nearby Stars
{'3-D selection of 167 sub-stellar companions to nearby stars': "Fabo Feng, 1, 2 R. Paul Butler, 3 Steven S. Vogt, 4 Matthew S. Clement, 3 C.G. Tinney, 5 Kaiming Cui, 1 Masataka Aizawa, 1 Hugh R. A. Jones, 6 J. Bailey, 5 Jennifer Burt, 7 B.D. Carter, 8 Jeffrey D. Crane, 9 Francesco Flammini Dotti, 1, 10, 11, 12 Bradford Holden, 4 Bo Ma, 13 Masahiro Ogihara, 1, 2 Rebecca Oppenheimer, 14 S.J. O'Toole, 15 Stephen A. Shectman, 9 Robert A. Wittenmyer, 8 Sharon X. Wang, 16 D.J. Wright, 8 and Yifan Xuan 1  \n1 Tsung-Dao Lee Institute, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, People's Republic of China 2 School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, People's Republic of China 3 Earth and Planets Laboratory, Carnegie Institution for Science, Washington, DC 20015, USA  \n4  \nUCO/Lick Observatory, University of California, Santa Cruz, CA 95064,USA  \nSchool of Physics and Australian Centre for Astrobiology, University of New South Wales, Sydney 2052, Australia  \n6  \nCentre for Astrophysics Research, University of Hertfordshire, College Lane, AL10 9AB, Hatfield, UK  \nJet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove drive, Pasadena CA 91109  \n8  \nUniversity of Southern Queensland, Centre for Astrophysics, USQ Toowoomba, QLD 4350 Australia  \n- 9 Observatories of the Carnegie Institution for Science, 813 Santa Barbara St., Pasadena, CA 91101 10 Department of Physics, Xi'an Jiaotong-Liverpool University, 111 Ren'ai Rd.,  \nSuzhou Dushu Lake Science and Education Innovation District, Suzhou Industrial Park, Suzhou 215123, P.R. China 11 Department of Mathematical Sciences, University of Liverpool, Liverpool L69 3BX, UK  \n- 12 Astronomisches Rechen-Institut, Zentrum fur Astronomie, University of Heidelberg, Monchhofstrasse 12-14, 69120, Heidelberg, Germany\n- 13  \nDepartment of Astronomy, Sun Yat-sen University, 74 Zhongshan Er Road, Guangzhou 510085, People's Republic of China 14 Astrophysics Department, American Museum of Natural History, Central Park West at 79th Street, New York, NY 10024, USA 15 Australian Astronomical Optics, Faculty of Science and Engineering, Macquarie University, 105 Delhi Rd, North Ryde, NSW 2113,  \nAustralia  \n16 Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China  \n(Dated: September 14, 2022; Received xxx; Revised xxx; Accepted xxx)  \nSubmitted to xx", 'ABSTRACT': 'We analyse 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 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.  \n5  \n7  \nKeywords: Exoplanet astronomy (486), Radial velocity (1332), Exoplanet detection methods (489), Brown dwarfs(185), Astrometric exoplanet detection(2130), Exoplanet catalogs(488)', '1. 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.  \nBDs 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.  \nFor 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.  \nWhile 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, 2000), 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. 2019, 2020; Acton 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.  \nThe 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  \nby the gravitational instability are both formed through the same process. We expect to be sensitive to such systems in this project.  \nTo 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, 2020) 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).  \nWe 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,  \n∆ r = 1 g ∆ t 2 ,  \n∆ µ = g ∆ t ,  \n2 (1) (2)  \nwhere ∆ 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 and Feng et al. 2016).  \nThis 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.", '2. 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 (Udry et al. 2000) 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).  \nThe 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  \nby 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'egransan 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'e 'echelle spectrograph mounted on the 2 m Alfred Jensch Telescope (AJT) of the Thueringer Landessternwarte Tautenburg (Dollinger 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.  \nFor a given target with both RV and revised Hipparcos catalog data (van Leeuwen 2007), we use the gaiadr2.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.  \nTo 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:  \nSNR µ = ∆ µ/σ µ , (3)  \nSNR r = ∆ r/σ r . (4)  \nIn 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  \n∆ µ = √ ( µ G α -µ H α ) 2 +( µ G δ -µ H δ ) 2 , (5)  \n∆ r = √ { ( α G -α H ) cos [( δ G + δ H ) / 2] -µ α ∆ t } 2 +( δ G -δ H -µ δ ∆ t ) 2 , (6)  \nσ µ = √ σ 2 µ H α + σ 2 µ H δ + σ 2 µ G α + σ 2 µ G δ , (7)  \nσ r = √ σ 2 α H + σ 2 δ H + σ 2 α G + σ 2 δ G , (8)  \nµ α = µ H α /σ 2 µ H α + µ G α /σ 2 µ G α 1 /σ 2 µ H α +1 /σ 2 µ G α , (9)  \nµ δ = µ H δ /σ 2 µ δ,H + µ G δ /σ 2 µ G δ 1 /σ 2 µ H δ +1 /σ 2 µ G δ , (10)  \nwhere { α 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  \ndata 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.", '3. 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.  \nRather 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).  \nConsidering 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 and Feng 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).  \nReflex 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 /dglyph[pi1] , 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.  \nFor 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.  \nTo 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 .  \nTable 1.  \nFeng et al. Number of companions identified in this work.  \nNote -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 /m c < 20%) and the optimal case defined by < 100% relative mass error ( σ m c /m c < 100%).  \nThe 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.", '4. 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%.  \nThe 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.  \nWe 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.  \nAround 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.  \nThe 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  \nFigure 1. The Hertzsprung-Russell diagram for the stellar sample in this work. The x axis is the BP-RP color while the y axis is the absolute G magnitude derived from the Gaia EDR3 catalog. The 869 stars showing long period RV signals are shown in blue color. Of these 161 of them are found to host companions with mass from 5 to 120 M Jup and are shown in red color.  \n<!-- image -->  \nHipparcos 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).  \nWe 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.  \nFigure 2. Corner plot for the distribution of the 161 companion hosts over various parameters. The diagonal panels show the histograms of corresponding parameters. Each of the lower panels shows the distribution of the sample over a parameter pair. The red and light green points respectively represent the hosts of known companions and the new companions identified in this work. From left to right on the top, the parameters are stellar mass ( M ), G magnitude, number of RV points ( N rv ), effective temperature ( T eff ), log(SNR r ), log(SNR µ ), and heliocentric distance ( d ). The nonlinear positional difference ∆ r is derived by subtracting the positional difference caused by proper motion from the total positional difference. The logarithmic scales are defined using base 10.  \n<!-- image -->  \nFigure 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 glyph[star] , 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.  \n<!-- image -->  \nTable 2. Hipparcos and Gaia EDR3 catalog astrometry for the sample of stars.  \nNote -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.  \nNote -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 ( glyph[pi1] ), 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.", 'T able 3 . Orbital parameters for stellar and planetary companions detected in this w ork.': '<!-- image -->  \nT able 3 c ont inue d on next p age  \nT able 3 (c ontinue d)', 'HZ Stab ilit y': "a  \nau  \nm c M Jup  \ne  \nI  \nP  \nHost  \nCompanion  \ndeg  \nyr  \nName  \nName  \nNote -T able 3 is published in its en tire t y in the electronic edition of the Astr ophysic al Journal Supplement Series . A p ortion is sho wn here for guidance regarding its form and con ten t. The systems are s orted b y the ho st names listed in the first column while the companions in eac h system are s orted b y their orbital p erio ds. F or companions with orbital p erio d less than 1000 da ys, the inclination is not giv en and the mass rep orted in this table is m c sin I . F or companions with mass higher than 75 M Jup , w e use capital letters to lab el them. Among the parameters, directly inferred parameters are P (orbital p erio d), K (R V semi-amplitude), e (eccen tricit y), ω glyph[star] (argumen t of p eriastron), M 0 (mean anomaly at the min im um ep o c h of R V data), I (inclination), Ω (longitude of ascending no de). The deriv ed parameters are T p (p eriastron ep o c h), m c (companion mass), and a (semi-ma jor axis). The dynamical stabilit y of Earth-lik e planets in the habitable zones of their host stars are determined n umerically , as describ ed in sectio n 7 . The median and the 1σ quan tiles (i.e. 16% and 84% quan tiles) are used to measure the uncertain t y of eac h parameter. The companions with † sup erscript ha v e b een imaged and their relativ e astrometry data are analyzed in co m bination with other t yp es of data in the solutions. The targets with glyph[star] sup erscript ha v e wide binary companio ns that are iden tified b y El-Badry et al. ( 2021 ) or other studies.  \nReferences -(1) Mann et al. ( 2019 ); (2) Agati et al. ( 2015 ); (3) Bona vita et al. ( 2016 ) ; (4) Burgasser et al. ( 2010 ); (5) Hein tz ( 1994 ); (6) Beuzit et al. ( 2004 ); (7) Goldm an et al. ( 2010 ); (8) Dupuy & Kraus ( 2 013 ); (9) F orv eille et al. ( 2 011 ); (10) Stassun et al. ( 2017 ); (11) Anglada-Escud ' e & T uomi ( 2012 ); (12) Marmier et al. ( 2013 ); (13) Desidera et al. ( 2012 ); (14) Li et al. ( 2021 ); (15) Zuc k er & Mazeh ( 200 0 ); (16) Simpson et al. ( 2010 ); (17) Wit tenm y er et al. ( 2009 ); (18) V ogt et al. ( 2000 ); (19) Hink el et al. ( 2019 ); (20) Ma y or et al. ( 2004 ); (21) Borgniet et al. ( 2014 ); (22) Xuan et al. ( 2020 ); (23) Borgnie t et al. ( 2019 ); (24) Jones et al. ( 2016 ); (25) Fisc her et al. ( 2007 ); (26) Men t et al. ( 2018 ); (27) G iguere et al. ( 2015 ); (28) Sato et al. ( 2013 ); (29) Wilson et al. ( 2016 ); (30) Luhn et al. ( 2019 ); (31) Lo Curto et al. ( 2010 ); (32) Ho w ard et al. ( 2010 ); (33) Martioli et al. ( 2010 ); (34) Fisc her et al. ( 2002 ); (35) Ric kman et al. ( 2019 ); (36) Ric kman e t al. ( 2020 ); (37) Do llinger et al. ( 2009 ) ; (38) W ang et al. ( 2014 ); (39) Tinney et al. ( 2002 ); (40) Wittenm y er et al. ( 2012 ); (41) Wittenm y er et al. ( 2020 ); (42) Butler et al. ( 2003 ); (43) Wittenm y er et al. ( 2007 ); (44) F eng et al. ( 2015 ); (45) D ' ıaz et al. ( 2012 ); (46) P atel et al. ( 2007 ); (47) Marcy et al. ( 1999 ); (48) Pily a vsky et al. ( 2011 ); (49) Naef et al. ( 2001 ); (50) Jones et al. ( 2021 ); (51) Arriagada et al. ( 2 010 ); (52) Thalmann et al. ( 2009 ); (53) Vigan et al. ( 2016 ); (54) Brandt et al. ( 2 021b ); (55) Bo wler et al. ( 2018 ); (56) Marcy et al. ( 20 05 ); (57) F uhrmann ( 2004 ); (58) T ok o vinin ( 2014a ); (59) Rob erts et al. ( 2015a ); (60) P errier et al. ( 2003 ); (61) Sahlmann et al. ( 2011 ); (62) Han et al. ( 2001 ); (63) Reffert & Quirren bac h ( 2011 ); (64) D ' ıaz et al. ( 2016a ); (65) S ' egransan et al. ( 2010 ); (66) Rob ertson et al. ( 2012 ); (67) Moutou et al. ( 2017 ); (68) Ginski et al. ( 2016 ); (69) Moutou et al. ( 2011 ); (70) Kane et al. ( 2019 ); (71) Me lo et al. ( 2007 ); (72) V enner et al. ( 2021 ); (73) Jenkins et al. ( 2017 ); (74) Rosen thal et al. ( 2021 ); (75) Moutou et al. ( 2005 ); (76) T rifono v et al. ( 2017 ) ; (77) Anglada-E scud ' e et al. ( 2010 ); (78) K urster et al. ( 2015 ); (79) San tos et al. ( 2001 ); (80) T ok o vinin ( 2014b ); (81) Griffin ( 2012 ); (82) Bouc h y et al. ( 2016 ); (83) Wittenm y er et al. ( 2017 ); (84) Zurlo et al. ( 2018 ); (85) Fisc her et al. ( 2001 ); (8 6) Lagrange et al. ( 2009 ); (87) Nielsen et al. ( 2 014 ); (88) Lagrange et al. ( 201 9 ); (89) Lacour et al. ( 2021 ); (90) Jones et al. ( 2002 ); (91) Huang et al. ( 2018 ); (9 2) Damasso et al. ( 2020a ); (93) Jenkins et al. ( 2015 ); (94) T am uz et al. ( 2008 ); (95) Che etham et al. ( 2018 ); (96) Mugrauer et al. ( 2014 ) ; (97) F eng et al. ( 2020b ); (98) Nak a jima et al. ( 1995 ); (99) Naef et al. ( 2010 ); (100) Al len et al. ( 2012 ); (101) Butler e t al. ( 2006 ); (102) Barbato et al. ( 2018 ); (103) Moutou et al. ( 2009 ); (104) Naef et al. ( 2004 ); (105) Dum usque et al. ( 2 011 ); (106) Wittenm y er et al. ( 2019 ); (107) Ro digas e t al. ( 2016 ); (108) Demangeon et al. ( 2021 ); (109) Sozzetti et al. ( 2006 ); (110) Bang et al. ( 2020 ); (111) Fisc her et al. ( 2009 ); (11 2) Hartmann et al. ( 2010 ); (113) Johnson-Groh et al. ( 2017a ); (114) W agner et al. ( 2019 ); (115) F ranson e t a l. ( 2022 ); (116) Jones et al. ( 2015c ); (117) Jones et al. ( 2017 ); (118) Jones et al. ( 2015a ); (119) Mason et al. ( 2001 ) .  \nFigure 4. Combined fit to RV, Gaia-Hipparcos astrometry, and direct imaging data for all brown dwarf targets. The left panels show the RV fits, the second and third columns respectively show the fit to the proper motion and positional differences between Gaia and Hipparcos, the right panels show the fit to the direct imaging data collected from literature. For each target, the linear motion is subtracted from the proper motion and position to show planet-induced nonlinear motion. The RV data sets as well as Gaia and Hipparcos astrometric data are encoded by the same colors as shown by the left-hand legends. The covariances of the Gaia and Hipparcos proper motions and positions are denoted by error ellipses. The straight line connects the data point to the best fitting model, which represents where the companion is expected to be. The MAP values of parameters of the companion with the longest orbital period are shown on the right side of the figure. 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.  \n<!-- image -->  \nBarycenter  \nFor previously published planets, we discuss our new results on an object-by-object basis as follows.  \n- · GJ 234 A (or Ross 614) is an M-dwarf hosting a companion with a mass of about 100 M Jup (Agati et al. 2015; Bonavita et al. 2016; Mann et al. 2019). Based on combined analyses of the RV, Gaia-Hipparcos astrometry and the imaging data collected by Mann et al. (2019), we estimate a mass of 94 . 84 +0 . 88 -1 . 37 M Jup .  \nHD39060  \nFigure 5. Similar to Fig. 4 but for other targets.  \n<!-- image -->  \nFigure 6. Combined fit to the RV and Gaia-Hipparcos data for all brown dwarf targets. The left panels show the RV fits, the middle and right ones respectively show the fit to the proper motion and positional differences between Gaia and Hipparcos. For each target, the linear motion is subtracted from the proper motion and position to show planet-induced nonlinear motion. The RV data sets as well as Gaia and Hipparcos astrometric data are encoded by the same colors as shown by the left-hand legends. The covariances of the Gaia and Hipparcos proper motions and positions are denoted by error ellipses. The straight line connects the data point to the best fitting model, which represents where the companion is expected to be. The MAP values of parameters of the companion with the longest orbital period are shown on the right side of the figure. 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.  \n<!-- image -->  \nFigure 7. Similar to Fig. 6 but for other targets.  \n<!-- image -->  \n- · 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. (2019) (labeled 'AM19'), we estimate a mass of 0.511 ± 0.004 M glyph[circledot] for GJ 494 A and a mass of 88 . 92 +1 . 84 -2 . 84 M Jup for GJ 494 B. According to the mass-luminosity relation derived by Mann et al. (2019), GJ 494 B is an M7.4-type dwarf while GJ 494 A is an M1.5-type dwarf, differing from previous classifications based on less accurate mass estimation (Heintz 1994; Beuzit et al. 2004).\n- · GJ 676 is a binary and the primary hosts four planets (Forveille et al. 2011; Anglada-Escud'e & Butler 2012). Our solution constrains the inclinations of GJ 676 A b and c to be 48 . 92 +3 . 31 -2 . 78 and 33 . 69 +1 . 36 -1 . 32 · . Our solutions determines their masses as 5.79 +0 . 47 -0 . 48 M Jup and 13 . 49 +1 . 05 -1 . 13 1 M Jup , respectively. With a projected separation of 761 au from GJ 676 A (El-Badry et al. 2021), GJ 676 B may play a role to warp the planetary disk such that the orbits of GJ 676 A d and e are slightly misaligned.\n- · HD 100939 (HIP 56640) is a K-type giant hosting a companion with a minimum mass of 3 . 67 ± 0 . 14 M Jup (Jones et al. 2021). Our combined analyses estimate a dynamical mass of 5 . 30 +3 . 309 -0 . 00 M Jup , putting it into the superJupiter category.\n- · 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 .\n- · HD 10697 (109 Psc) is a G star hosting a companion with a minimum mass of 6.38 ± 0.08 M Jup (Vogt et al. 2000; Butler et al. 2006; Wittenmyer et al. 2009; Luhn et al. 2019). Our analyses estimate a mass of 5 . 74 +1 . 01 -0 . 29 M Jup .\n- · HD 111031 (GJ 3746) is a G-type star hosting a planet candidate with a minimum mass of 3.7 M Jup (Hinkel et al. 2019). Based on combined analyses of the RV and astrometry data, we estimate a mass of 54 . 17 +5 . 32 -6 . 15 M Jup .\n- · HD 111232 (HIP 62534) is a G star hosting a companion with a minimum mass of 6.80 M Jup (Mayor et al. 2004; Stassun et al. 2017). Our analyses estimate a mass of 7 . 97 +1 . 13 -0 . 48 M Jup . We also find a BD companion with a mass of 18 . 06 +4 . 21 -1 . 61 M Jup on an 17 au-wide orbit. The HD 111232 system is one of few systems hosting multiple sub-stellar companions in our sample.\n- · HD 113337 hosts a disc with an inclination of about 13 +10 -9 deg (Xuan et al. 2019). Xuan et al. (2019) estimate an inclination of 31 +5 -4 deg based on analyses of Gaia DR2 and Hipparcos data. However, our solution determines an inclination of 57 . 48 +4 . 55 -3 . 65 deg. Such discrepancy is probably due to our use of the EDR3 catalog.\n- · HD 11343 hosts a super-Jupiter (Jones et al. 2016) and a wide-orbit M-dwarf companion (El-Badry et al. 2021). We find a new Jupiter analog, HD 11343 c, with a minimum mass of about 0.5 M Jup . The astrometric signal induced by this HD 11343 c is weak and thus the astrometric data does not constrain its inclination well. The mass of HD 11343 b is 7 . 71 +0 . 73 -1 . 19 M Jup .\n- · 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.\n- · HD 120084 is a G star hosting a companion with a minimum mass of 4.5 M Jup (Sato et al. 2013). Our combined analyses estimate a mass of 5 . 76 +4 . 62 -0 . 31 M Jup .\n- · HD 122562 is a G star hosting a companion with a minimum mass of 24 ± 2 M Jup (Wilson et al. 2016). Our solution gives a dynamical mass of 38 ± 4 M Jup . In addition to this companion, we find a warm super-Neptune, HD 122562 c, with an orbital period of 84 d. On the basis that there seem to be other planets in the system we remain consistent with previous naming of the more massive companion and label it HD 122562 b. Its mass is 37 . 20 +3 . 67 -2 . 60 M Jup .  \n- · HD 125390 hosts a companion with a minimum mass of 22.16 ± 0.96 M Jup (Luhn et al. 2019). Our combined analyses estimate a mass of 27 . 20 +4 . 13 -0 . 31 M Jup .\n- · 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 .\n- · 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). With combined analyses, we estimate a mass of 0 . 34 +0 . 20 -0 . 02 M Jup for the planetary companion and 81 . 13 +7 . 78 -7 . 92 M Jup for the M-dwarf companion.\n- · HD 127506 (GJ 554) is a K-type star hosting a brown dwarf HD 127506 B (Halbwachs et al. 2000; Zucker & Mazeh 2001). The BD is further constrained by Sozzetti & Desidera (2010) and Reffert & Quirrenbach (2011) using Hipparcos intermediate astrometric data. A mass of 45 ± 21 M Jup is estimated by Sozzetti & Desidera (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.\n- · 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 .\n- · HD 136118 (HIP 74948) is an F-type star hosting a BD with a mass of 42 +11 -18 mj (Fischer et al. 2002; Wittenmyer et al. 2009). Our solution estimates a mass of 13 . 10 +1 . 35 -1 . 27 M Jup .\n- · 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.\n- · HD 139357 is a K-type star hosting a companion with a minimum mass of 9.76 ± 2.15 M Jup (Dollinger 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 .\n- · 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 .\n- · 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 (Tinney et al. 2002; Wittenmyer et al. 2012). Our analyses confirm HD 142 A d, a warm Saturn identified by Wittenmyer et al. (2020). The astrometric data put a strong constraint on the inclination of HD 142 A c, leading to a mass estimation of 10 . 90 +1 . 28 -0 . 94 M Jup .\n- · HD 145675 (14 Her) hosts two companions (Butler et al. 2003; Wittenmyer et al. 2007; Bardalez Gagliuffi et al. 2021). Our combined analyses estimate m b = 8 . 05 +1 . 63 -0 . 88 M Jup , m c = 5 . 03 +0 . 87 -1 . 07 M Jup , I b = 144 . 65 +6 . 28 -3 . 24 deg and I c = 120 . 10 +6 . 26 -29 . 05 deg for HD 145675 b and c, consistent with an aligned orbital configuration. However, Bardalez Gagliuffi et al. (2021) give a solution of I b = 32 . 7 +5 . 3 -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).  \n- · 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).\n- · HD 167665 (HIP 89620) is an F-type star hosting a BD (Patel et al. 2007) with a minimum mass of 50.6 ± 1.7 M Jup (Wilson et al. 2016). Based on our combined RV and astrometry analysis, HD 167665 B has a mass of 52 . 71 +5 . 11 -4 . 40 M Jup after considering stellar mass uncertainty. Its inclination is 95 . 20 +7 . 47 -9 . 86 deg, suggesting a non-zero transit probability.\n- · HD 170707 (HIP 90988) is a K giant hosting a planet with a minimum mass of 1.96 ± 0.11 M Jup (Jones et al. 2021). Our combined analyses estimate a minimum mass of 2 . 10 +0 . 08 -0 . 09 M Jup for this companion. We also find a stellar companion with a mass of 101 . 79 +5 . 91 -6 . 11 M Jup .\n- · HD 181234 (HIP 95015) is a G-type star hosting a companion with a minimum mass of 8.37 +0 . 34 -0 . 36 M Jup (Rickman et al. 2019). Based on combined RV and astrometry analysis, we constrain the mass of this companion to be 8 . 24 +2 . 15 -0 . 00 M Jup and the inclination to be 98 . 73 +24 . 47 -41 . 61 deg.\n- · HD 182488 hosts a BD companion named 'HD 182488 B' or 'GJ 758 B' Thalmann et al. (2009) and imaged by Vigan et al. (2016). It is further studied by multiple groups (e.g., Bowler et al. 2018 and Brandt et al. 2021b) through combined analyses of astrometry, imaging data (labeled 'BB18') and RV data. In particular, Brandt et al. (2021b) estimate a mass of 38 . 0 ± 0 . 8 M Jup and an eccentricity of 0.24 ± 0.11. By analyzing the same data set as that used by Brandt et al. (2021b) as well as the Gaia-Hipparcos positional difference, we constrain the mass to be 36 . 39 +1 . 21 -1 . 08 M Jup and the eccentricity to be 0 . 37 +0 . 08 -0 . 08 , consistent with the solution given by Brandt et al. (2021b). With astrometry for both components of the binary, we also estimate a mass of 0.93 ± 0.03 M glyph[circledot] for the primary a posteriori .\n- · HD 185414 (HIP 96395; WDS 19359+5659) is a spectroscopic binary hosting a wide-orbit companion detected by 2MASS and Gaia (Mason et al. 2001; Fuhrmann 2004; Tokovinin 2014b; Roberts et al. 2015b; El-Badry et al. 2021). The spectroscopic binary is unresolved. Our solution estimates a mass of 117 . 59 +9 . 03 -9 . 55 M Jup for this unresolved companion.\n- · HD 190228 is a G star hosting a Jupiter-like companion, HD 190228 b (Perrier et al. 2003; Wittenmyer et al. 2009). The inclination of HD 190228 b is 4.5 ± 2.1 deg based on the analyses of the Hipparcos intermediate astrometric data (Reffert & Quirrenbach 2011). Our solution does not put a strong constraint on the inclination because the orbital period is comparable with the Hipparcos or Gaia EDR3 baselines. A detailed analyses of intermediate data in a future Gaia data release will strongly constrain the mass.\n- · HD 191806 is a K star hosting a companion with a minimum mass of 8.52 ± 0.63 M Jup (D'ıaz et al. 2016b). Our combined analyses constrain the mass to be 9 . 33 +0 . 92 -0 . 85 M Jup .\n- · HD 203473 B was discovered by Ment et al. (2018) based on analysis of KECK RV data. With both KECK and HARPS data as well as the Hipparcos-Gaia astrometry, we estimate a mass of about 96 M Jup , much higher than the minimum mass of 7.84 ± 1.15 M Jup estimated by Ment et al. (2018). Our solution estimates an orbital period of 8 . 10 +0 . 01 -0 . 02 yr, much longer than the period of 4.25 yr estimated by Ment et al. (2018), who may have identified a half-period harmonic of the true signal due to limited RV baseline.\n- · HD204313 is claimed to host three planets based on the RV data collected by the 2.7 m Harlan J. Smith Telescope and CORALIE spectrometer (S'egransan et al. 2010; Robertson et al. 2012). However, D'ıaz et al. (2016a) do not confirm HD 204313 d based on their analyses of previous and new HARPS data. Our solution agrees with D'ıaz et al. (2016a) that there is no significant signal around a period of 2800 d as claimed by Robertson et al. (2012). However, we confirm the trend identified by D'ıaz et al. (2016a) as a new long period signal with a period of 20 . 06 +1 . 10 -1 . 01 yr based on combined RV and astrometry analysis. This new companion is dubbed 'HD 204313 e' to be distinguished from the dubious candidate HD 204313 d. It corresponds to a BD with a mass of 15 . 32 +4 . 89 -5 . 18 M Jup .  \n- · HD 211847 is a G star hosting a sub-stellar companion with a minimum mass of 19.2 ± 1.2 M Jup (Sahlmann et al. 2011; Moutou et al. 2017). Our analyses constrain the mass to be 55 . 32 +1 . 34 -18 . 49 M Jup .\n- · 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.\n- · 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 (Ginski et al. 2016; 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.\n- · HD 217958 (HIP 113948) is a G star hosting a possible stellar companion (Kane et al. 2019). Our combined analyses of the RV, Gaia-Hipparcos astrometry and the imaging data given by (Kane et al. 2019) do not show strong evidence for this candidate companion on a wide orbit. However, we discover a companion with a mass of 81 . 82 +34 . 17 -38 . 20 M Jup and a planet with a mass of 0 . 52 +0 . 35 -0 . 06 M Jup on a Jupiter-like orbit.\n- · HD 219077 (HIP 114699) is a G star hosting a companion with a minimum mass of 13.4 ± 0.78 M Jup (Marmier et al. 2013; Kane et al. 2019). Our combined analyses gives a dynamical mass of 9 . 62 +1 . 00 -0 . 73 M Jup .\n- · HD219828 (HIP 115100) is a G star hosting two companions with minimum masses of 15.1 and 0.0661 M Jup (Melo et al. 2007; Ment et al. 2018). Our combined analyses estimate a mass of 16 . 05 +3 . 52 -1 . 27 M Jup for the bigger companion on a Jupiter-like orbit.\n- · HD 221420 (GJ 4340) is a G-type star hosting a BD (Kane et al. 2019) with a wide-orbit M-dwarf companion (Venner et al. 2021; El-Badry et al. 2021). Based on combined RV and the HGCA catalog of calibrated Gaia DR2 and Hipparcos proper motions (Brandt 2018), Venner et al. (2021) estimated a mass of 20 . 35 +1 . 99 -3 . 43 M Jup , consistent with 20 ± 3 M Jup estimated in this work.\n- · HD 224538 (HIP 118228) is an F star hosting a companion with a minimum mass of 5.97 ± 0.42 M Jup (Jenkins et al. 2017). Our combined analyses estimate a dynamical mass of 6 . 53 +1 . 84 -0 . 35 M Jup .\n- · HD 23596 (HIP 17747) is an F star hosting a sub-stellar companion with a mass of 8.1 M Jup (Perrier et al. 2003; Wittenmyer et al. 2009; Stassun et al. 2017) and a wide-binary companion (El-Badry et al. 2021). Our combined analyses show a dynamical mass of 11 . 91 +0 . 99 -1 . 77 M Jup .\n- · HD 25015 (HIP 18527) is a K star hosting a companion with a minimum mass of 4.5 ± 0.3 M Jup (Rickman et al. 2019). Our combined analyses show a mass of 9 . 08 +1 . 16 -1 . 81 M Jup .\n- · HD 26161 (HIP 19428) is a G-type star hosting a companion with a minimum mass of 13.5 +8 . 5 -3 . 7 M Jup (Rosenthal et al. 2021). Based on the combined RV and astrometry analysis, the dynamical mass of this companion is 28 . 46 +20 . 05 -0 . 22 M Jup and the orbital period is about 39 years. HD 26161 also hosts an M dwarf companion with a projected separation of 561 au (El-Badry et al. 2021).\n- · 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'e et al. 2010; Kurster 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 .\n- · 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 .  \n- · 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 glyph[circledot] (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.\n- · HD 29461 (HIP 21654) is a G star hosting a companion with a minimum mass of about 0.08 M glyph[circledot] (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.\n- · 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 (Tinney et al. 2002; 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.\n- · HD 38529 (HIP 27253) is a G star hosting two companions with masses of 18 ± 3 and 0.17 M Jup (Fischer et al. 2003; Wittenmyer et al. 2009; Barbato et al. 2018; Benedict et al. 2010; Xuan et al. 2020). Our combined analyses estimate a mass of 10 . 38 +1 . 03 -0 . 88 M Jup and an inclination of 104 . 56 +6 . 39 -8 . 72 deg for HD 38529 c, consistent with the the solution given by Xuan et al. (2020) but with higher precision.\n- · HD 39060 ( β Pic) is a well-studied system, hosting two giant planets (Lagrange et al. 2009, 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).  \nWe 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.  \n- · 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ˇzi'c 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  \nTable 4. Parameters for the β Pic system based on analyses of various data sets in literature and in this work  \nthe 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).  \n- · HD 39213 (HIP 27491) is a K star hosting a BD with a minimum mass of 0.07 ± 0.01 M glyph[circledot] (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.\n- · 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.\n- · 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  \nin 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.  \n- · 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 .\n- · 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.\n- · 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 estimate a mass of 10.6 ± 2.8 M Jup for HD 66428 b and 3.1 ± 1.7 M Jup for HD 66428 c. Compared with an orbital period of 107 +153 -49 yr and a minimum mass of 0 . 085 +0 . 069 -0 . 053 M Jup estimated for HD 66428 c by Rosenthal et al. (2021), our estimation of a period of 28 . 69 +9 . 21 -5 . 35 yr and 1 . 76 +3 . 40 -0 . 04 M Jup is much more precise. However, further constraint on the inclination of HD 66428 c is needed to determine whether the orbits of the two companions are misaligned.\n- · HD 72659 (HIP 42030) is a G-type star hosting a companion (Butler et al. 2003; Wittenmyer et al. 2009) with a minimum mass of 3.15 ± 0.14 M Jup (Moutou et al. 2011). Based on the combined analysis of RV and astrometry, HD 72659 b is found to have a mass of 2 . 99 +2 . 59 -0 . 10 M Jup . We also find a new companion HD 72659 c with a mass of 18 . 81 +4 . 44 -4 . 80 M Jup . Assuming a coplanar configuration between b and c, HD 72659 b would have a mass of about 15 M Jup .\n- · 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.\n- · 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.\n- · HD 74014 (HIP 42634) is a star hosting a BD companion (Patel et al. 2007) with a minimum mass of 49.0 ± 1.7 M Jup (Sahlmann et al. 2011). Our combined RV and astrometry analyses constrain the mass to 61 . 54 +5 . 61 -5 . 77 M Jup .\n- · HD 74156 (HIP 42723) is a G star hosting two sub-stellar companions with masses of 1.78 and 8.00 M Jup (Naef et al. 2004; Wittenmyer et al. 2009). Our combined analyses lead to estimated masses of 1 . 71 +0 . 96 -0 . 06 and 8 . 67 +1 . 39 -0 . 47 M Jup .\n- · HD 7449 A b and HD 7449 A c (or HD 7449 B) are detected by Dumusque et al. (2011) and Rodigas et al. (2016), respectively. Through combined analyses of RV, Gaia-Hipparcos astrometry, and the relative astrometry derived from the imaging data collected by Rodigas et al. (2016) (labeled 'TR16'), we estimate inclinations of 171 . 63 +2 . 61 -3 . 74 deg and 68 . 40 +4 . 10 -3 . 89 deg for Ab and B, indicating significant misalignment between the two companions. The mass of HD 7449 A b is 8 . 17 +3 . 06 -2 . 70 M Jup . The mass of HD 7449 B is 178 . 15 +16 . 61 -13 . 66 M Jup , consistent with the value of 0.23 +0 . 22 -0 . 05 M glyph[circledot] estimated by Rodigas et al. (2016) based on photometry.\n- · HD 80869 (HIP 46022) is a G star hosting a planetary companion with a minimum mass of 4.86 +0 . 65 -0 . 29 M Jup (Demangeon et al. 2021). The use of Gaia-Hipparcos astrometry allows us to break the inclinationmass degeneracy and estimate a mass of 5 . 07 +2 . 54 -0 . 56 M Jup .  \nFigure 8. BFPs and phase curves for the 600 and 1000d signals detected in the RV data for HD 81817. Upper panels: the left panel shows the BFP for the raw data while the right panel shows the BFP for the residual with the 600-day signal subtracted. The red lines in the upper panels indicate the periods of the two signals. The dashed and dotted lines show the thresholds of lnBF=5 and 0, respectively. Lower panels: the red curves in the lower panels show the best fit corresponding to the parameter values at the maximum a posteriori. The black error bars represent 10 binned data points.  \n<!-- image -->  \n- · 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).\n- · HD81817 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.  \n- · 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 .\n- · 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 glyph[circledot] (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 .\n- · 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.\n- · 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.\n- · 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 .\n- · 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).\n- · HIP 22203 (HD 30246) is a G star hosting a brown dwarf with a minimum mass of 55 . 1 +20 . 3 -8 . 2 M Jup (D'ıaz et al. 2012). Our combined analyses constrain the mass to be 51 ± 5 M Jup and the inclination to be 50 . 76 +4 . 73 -4 . 54 deg.\n- · HIP 65891 (HD 117253) is a red giant branch K star hosting a planet with a minimum mass of about 6 M Jup (Jones et al. 2015c). Using both RV and Gaia astrometric excess noise, Kiefer et al. (2021) estimate a mass ranging from 168.4 to 713.7 M Jup for HD 117253 b. However, based on Gaia-Hipparcos data, we find a mass of 5 . 89 +0 . 76 -0 . 28 M Jup and an inclination of 88 . 97 +16 . 21 -17 . 33 deg. 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.\n- · 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 (Jones et al. 2017). Our combined analyses constrain the companion's mass to be 10 . 88 +7 . 78 -0 . 00 M Jup .\n- · 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 .\n- · 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.\n- · 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.", '5. STATISTICS OF THE COMPANION SAMPLE': "5.1. Mass distribution and occurrence rate  \nThrough 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  \n5 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).  \nWe 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.", '5.2. 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.  \nAll 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.  \nThere 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).  \nThe 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  \nFigure 9. Mass distributions of the sub-stellar companions detected in this work. The samples are divided into 10 bins. The upper and lower panels respectively show the distributions of the companions with relative mass uncertainty of less than 100% and 20% . The left and right panels respectively show the distributions over companion mass and companion-host mass ratio. The error bars are determined by assuming that the number of detected companions follows a Poisson distribution.  \n<!-- image -->  \nbinary 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.  \nBased 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 = 330au. 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.  \nFigure 10. The 61 multi-companion systems identified in this work. The yellow dots represent stellar companions with mass higher than 75 M Jup , the brown dots represent BDs with mass range from 13 to 75 M Jup , and the blue dots represent planet companions with mass lower than 13 M Jup . 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.  \n<!-- image -->  \nFigure 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.  \n<!-- image -->  \n(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.", '6. 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  \nFigure 12. Companion-host contrast in J, H, and K bands by assuming ages of 0.1, 1, and 10 Gyr for the systems.  \n<!-- image -->  \nand 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.  \nTable 5 . Mean separation and base-10 logarithmic flux ratio of companions to host stars.  \nNote -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).", '7. DYNAMICAL STABILITY': "Figure 13. An example evolutionary scheme from one of our dynamical simulations studying the low mass companion and Jupiter analog in HD 205158. The pericenter and apocenter for each body is plotted with red (HD 205158 B) and blue (HD 205158 b) lines.  \n<!-- image -->  \nWe 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 Mercury 6 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.", '7.1. 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.", '7.2. 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  \n2  \nFigure 14. Summary of our simulations studying the stability of an Earth-mass planet in the HZ of our various systems identified in 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.  \n<!-- image -->  \nHZ 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.", '8. 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.  \nWithout 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.  \nBecause 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.  \nWe 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.  \n- · 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.\n- · 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-  \nignment 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.  \n- · 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.\n- · 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.  \nIn 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.  \nIn 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).", 'ACKNOWLEDGEMENTS': "This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has also made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration. This research has also made use of the services of the ESO Science Archive Facility, NASA's Astrophysics Data System Bibliographic Service, and the SIMBAD database, operated at CDS, Strasbourg, France. Support for this work was provided by Chinese Academy of Sciences President's International Fellowship Initiative Grant No. 2020VMA0033. The authors acknowledge the years of technical support from LCO staff in the successful operation of PFS, enabling the collection of the data presented in this paper. We would also like to acknowledge the many years of technical support from the UCO/Lick staff for the commissioning and operation of the APF facility atop Mt. Hamilton. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Specifically, it used the Bridges2 system, which is supported by NSF award number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC: Nystrom et al. 2015). F.F.D. was supported by the Research Development Fund (grant RDF-16-01-16) of Xi'an Jiaotong-Liverpool University (XJTLU) and was supported by the SPP 1992 exoplanet diversity program and Berlin Technical University. B.M. acknowledges financial support from NSFC grant 12073092 and the science research grants from the China Manned Space Project (No. CMS-CSST-2021-B09). This paper is partly based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes: 0100.C-0097,0100.C-0414,0101.C-0232,0101.C-0379,072.C-0488,074.C0364,075.C-0332,075.D-0800,076.C-0155,077.C-0101,077.C-0364,079.C-0657,079.C-0927,081.C-0148,081.C-0802,082.C-  \n0212,082.C-0427,084.C-0228,085.C-0019,085.C-0063,086.C-0284,087.C-0368,087.C-0831,088.C-0662,089.C-0497,089.C0732,090.C-0395,090.C-0421,091.C-0034,091.C-0866,091.C-0936,091.D-0469,092.C-0721,093.C-0409,094.C-0901,095.C0551,095.D-0026,096.C-0460,097.C-0090,098.C-0366,098.C-0739,099.C-0458,180.C-0886,183.C-0437,183.C-0972,183.D0729,188.C-0265,190.C-0027,191.C-0873,192.C-0224,192.C-0852,196.C-1006,198.C-0836,60.A-9036, and 60.A-9700.  \nSoftware: Astropy (Astropy Collaboration et al. 2013, 2018), Astroquery (Ginsburg et al. 2019), Numpy (van der Walt et al. 2011; Harris et al. 2020), Mercury (Chambers 1999), magicaxis (Robotham 2016), fields (Nychka et al. 2018), MASS (Ripley et al. 2013), minpack.lm (Elzhov et al. 2016)", 'APPENDIX': 'Fig. Set 1. Combined fit for the 144 systems for online access.'}