3D Selection of 167 Substellar Companions to Nearby Stars

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 Hayashi & Nakano (1963), Kumar (1963), Burrows et al. (2001), Spiegel et al. (2011), Baraffe et al. (2015), and Marley et al. (2021), we use the hydrogen-burning and deuterium-burning mass limits of 13 M_Jup 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, BDs can form through both channels. Hence, BDs can be considered as a unique population, bridging planets and stars, deserving intensive scientific investigations.

BDs can have four different spectral types: M, L, T, and Y. GJ 229 B (Nakajima et al. 1995) and Teide 1 (Rebolo et al. 1995) were the first two BDs to be 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 BD 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 them being classified as individual "free-floating" objects (e.g., objects that are gravitationally bounded only to the galaxy's central potential, rather than as a companion in a binary). Luhman (2007) found that stars and BDs are mixed homogeneously, based on their spatial kinematics being indistinguishable in Chamaeleon I, a young star-forming region, which further supports the understanding that both stars and BDs have the same formation mechanism.

For directly imaged BDs, whether free-floating or having a stellar companion, the cooling model is typically used to provide an indirect estimate of their mass (e.g., Baraffe et al. 2015; Marley et al. 2021). Due to the lack of hydrogen fusion in the cores 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 the BD mass, metallicity, cloud coverage, etc. However, cooling models are diverse, and sometimes estimate a mass inconsistent with the dynamical mass constrained by the direct imaging of BDs. Hence, a large sample of BD companions to stars with known masses are essential for testing various cooling models. The BD host star provides a natural reference point, where age might potentially be accurately determined, for us to better quantify the formation and evolution of BDs.

While we have few direct mass determinations of the thousands of free-floating BDs, we have a much smaller population of BDs orbiting around stars, whose masses are well constrained from orbital fits. Hereafter, we will refer to these BDs orbiting stars as "circumstellar BDs." Radial velocity (RV) surveys find that the occurrence rate of circumstellar BDs is measured as being 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 RV surveys compiled results (e.g., Campbell et al. 1988; Marcy & Benitz 1989; Marcy & Butler 1995, 2000), although the significant observational biases of the RV and subsequent transit surveys were considered as plausible causes. Through an adaptive optics imaging survey, Metchev et al. (2009) found that there were more BDs at wide orbits than BDs in the BD desert. Recently, more and more BDs have been found to reside in this desert (e.g., Carmichael et al. 2019, 2020; Persson et al. 2019; Acton et al. 2021). As suggested by Shahaf & Mazeh (2019), the characterization of the shape of the BD desert in a period–mass diagram by means of a large sample of circumstellar BDs can improve our understanding of its origin.

The correlation between the inner companion planets and the wide-orbit companions may be the key to improving 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 (CJs) were formed by core accretion or gravitational instability (Chabrier et al. 2014). In this manner, our work aims to bolster the sample of detected systems hosting CJs in order to provide improved constraints for theoretical formation models. According to Zhu et al. (2012), when BDs form due to gravitational instability, multiple CJs and BDs are expected to form simultaneously. Therefore, it is reasonable to consider that CJs formed inside the orbit of BDs formed by gravitational instability are both formed through the same process. We expect to be sensitive to such systems in this project.

To constrain the BD cooling models, to test the BD desert hypothesis, to find the correlation between BDs, CJs, and other types of planets, and to provide a large sample of candidates for direct imaging, we need a larger circumstellar substellar sample covering the BD mass range. This requires the improvement of the BD detection sensitivity through a synergy of all the 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, 2021) for nearby stars. Various groups have used this approach to estimate the dynamical mass of directly imaged BDs (Snellen & Brown 2018; Brandt et al. 2019; Kervella et al. 2019; Xuan et al. 2020; Kiefer et al. 2021; Kervella et al. 2022).

We follow the approach developed by Feng et al. (2019a), which uses 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 the Gaia–Hipparcos time difference, it is more sensitive to companion-induced acceleration of the primary star than the proper-motion difference, which is a linear function of time. In other words,

$$\begin{eqnarray}&&{\rm{\Delta }}r=\displaystyle \frac{1}{2}g{\rm{\Delta }}{t}^{2},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}\mu =g{\rm{\Delta }}t,\end{eqnarray} \tag{ 2 }$$

where Δr is the amplitude of the positional change, Δμ is the proper-motion change, Δt is the difference between the reference times of Gaia Early Data Release 3 (EDR3) and the Hipparcos catalog, and g is the acceleration of the primary star induced by a companion. The combined analysis of both proper-motion and positional differences has been 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 could to some extent remove systematics 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 overfitting or underfitting problems (e.g., Foreman-Mackey et al. 2015; Feng et al. 2016).

This paper is structured as follows. The RV and astrometry data are introduced in Section 2. The combined modeling of the 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 our conclusion in Section 8.

In this work, we use the RV data of the University College London Echelle Spectrograph (Diego et al. 1990) mounted on the Anglo-Australian Telescope (AAT), the Automated Planet Finder (APF; Vogt et al. 2014), the Levy Spectrometer at the Lick Observatory, the CORALIE spectrometer (Udry et al. 2000) installed at the Swiss 1.2 m 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) on the ESO La Silla 3.6 m 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 the Very Large Telescope (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.93 m telescope of Observatoire Haute-Provence, the ESO UV-visual echelle spectrograph (UVES) on Unit Telescope 2 of the VLT array, and the High Resolution Spectrograph (HRS; Tull 1998) mounted on the Hobby-Eberly Telescope (Ramsey et al. 1998).

The HARPS data are reduced by Trifonov et al. (2020), using the SERVAL pipeline (Zechmeister et al. 2018). There is a known offset in the RV zeropoint for the post-2015 data set (Lo Curto et al. 2015). Hence, we label the pre-2015 data set "HARPSpre" and the post-2015 data set 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) are published by Bowler et al. (2018). This data set is denoted by "APF1," while the other archived APF data are labeled "APF2." For β Pic (or HD 39060), the RV data reduced by Lagrange et al. (2019) are used, and labeled "AL19." We use the published RV data from the Lick Hamilton spectrograph and label the versions resulting from various updates as Lick6, Lick8, and Lick13 (Fischer et al. 2014). Since the first operation of CORALIE at 1998, it has had major upgrades in 2007 (Ségransan et al. 2010) and in 2014. Hence, we use COR98, COR07, and COR14 to denote the three versions of the 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 ¹⁷ 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, as 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, 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. For HD139357, we use the RVs measured by the coudé échelle spectrograph mounted on the 2 m Alfred Jensch Telescope (AJT) of the Thueringer Landessternwarte Tautenburg (Döllinger et al. 2009). For HD 106515A, we use the RV data measured by the high-resolution spectrograph SARG at TNG (Desidera et al. 2012). The new RV data for all targets are shown in the online figure set in the Appendix; the new data results are also available as data behind the figure.

For a given target with both RV and revised Hipparcos catalog data (van Leeuwen 2007), we use the gaiadr2.tmass_best_neighbour crossmatching catalog in the Gaia data archive to find the Gaia Data Release 2 (DR2) source identity and use the gaiaedr3.dr2_neighbourhood crossmatching catalog to find the EDR3 data. ¹⁸ For a target without Gaia counterparts in the crossmatching catalog, we select the Gaia sources within 0.1° of their 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 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 taken from the TESS input catalog (Stassun et al. 2019), unless it can be derived from the combined analyses of the RV, astrometric, and imaging data.

To compare the significance of companion-induced proper-motion (Δμ) and positional (Δr) differences between Gaia and Hipparcos, we calculate the signal-to-noise ratios (η) for Δμ and Δr as follows:

$$\begin{eqnarray}&&{\eta }_{\mu }={\rm{\Delta }}\mu /{\sigma }_{\mu },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\eta }_{r}={\rm{\Delta }}r/{\sigma }_{r}.\end{eqnarray} \tag{ 4 }$$

In the above equations, the proper-motion (Δμ) and positional (Δr) differences and their uncertainties (σ_μ and σ_r ), as well as the error-weighted proper motion (μ_α and μ_δ ), are

$$\begin{eqnarray}&&{\rm{\Delta }}\mu =\sqrt{{({\mu }_{\alpha }^{G}-{\mu }_{\alpha }^{H})}^{2}+{({\mu }_{\delta }^{G}-{\mu }_{\delta }^{H})}^{2}}\,,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}r\\ =\sqrt{\{{({\alpha }^{G}-{\alpha }^{H})\cos [({\delta }^{G}+{\delta }^{H})/2]-{\mu }_{\alpha }{\rm{\Delta }}t\}{}^{2}+({\delta }^{G}-{\delta }^{H}-{\mu }_{\delta }{\rm{\Delta }}t)}^{2}}\,,\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\sigma }_{\mu }=\sqrt{{\sigma }_{{\mu }_{\alpha }^{H}}^{2}+{\sigma }_{{\mu }_{\delta }^{H}}^{2}+{\sigma }_{{\mu }_{\alpha }^{G}}^{2}+{\sigma }_{{\mu }_{\delta }^{G}}^{2}}\,,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\sigma }_{r}=\sqrt{{\sigma }_{{\alpha }^{H}}^{2}+{\sigma }_{{\delta }^{H}}^{2}+{\sigma }_{{\alpha }^{G}}^{2}+{\sigma }_{{\delta }^{G}}^{2}}\,,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\mu }_{\alpha }=\displaystyle \frac{{\mu }_{\alpha }^{H}/{\sigma }_{{\mu }_{\alpha }^{H}}^{2}+{\mu }_{\alpha }^{G}/{\sigma }_{{\mu }_{\alpha }^{G}}^{2}}{1/{\sigma }_{{\mu }_{\alpha }^{H}}^{2}+1/{\sigma }_{{\mu }_{\alpha }^{G}}^{2}}\,,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\mu }_{\delta }=\displaystyle \frac{{\mu }_{\delta }^{H}/{\sigma }_{{\mu }_{\delta ,H}}^{2}+{\mu }_{\delta }^{G}/{\sigma }_{{\mu }_{\delta }^{G}}^{2}}{1/{\sigma }_{{\mu }_{\delta }^{H}}^{2}+1/{\sigma }_{{\mu }_{\delta }^{G}}^{2}}\,,\end{eqnarray} \tag{ 10 }$$

where $\{{\alpha }^{G},{\delta }^{G},{\mu }_{\alpha }^{G},{\mu }_{\delta }^{G}\}$ and $\{{\alpha }^{H},{\delta }^{H},{\mu }_{\alpha }^{H},{\mu }_{\delta }^{H}\}$ are the Gaia and Hipparcos astrometry data, respectively, including R.A., decl., and proper motion in R.A. and decl. The uncertainty of the astrometry data is denoted by σ, with corresponding subscripts. The correlations between the various astrometric parameters are not considered here, although they are considered in our full modeling of the astrometric data. The η_μ and η_r for the stellar sample in this work will be presented in Section 4. In the calculation of the positional difference, the linear motion due to proper motion is subtracted. Although the perspective acceleration is not considered in the calculation of η, it is taken into account in our rigorous modeling of the astrometry, which will be introduced in the following section.

Past and current RV surveys have provided 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, Gaia EDR3 is the third-epoch release of astrometric data for more than one billion stars. Comparison between the Hipparcos and Gaia catalogs for common stars provides an additional constraint on the potential accelerations induced by companions on primary stars.

Rather than just relying on Gaia and Hipparcos proper motions, the frequently used approach is to compare the 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 to or longer than 24 yr. Instead of calibrating Gaia and Hipparcos catalogs a priori, we fit the calibration parameters and signal model parameters to the raw catalog data simultaneously, to avoid overfitting problems caused by conducting the calibration before the signal search (Foreman-Mackey et al. 2015; Feng et al. 2017). This method is optimal, compared with the previous methods, in terms of constraining the masses and orbits of small planets, according to Feng et al. (2021).

Considering that the models and numerical methods are introduced by Feng et al. (2019a) and Feng et al. (2021), we only briefly describe the methodology in this section. The RV 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 qth-order MA model (or MA(q)), to derive the Bayesian information criterion (BIC) of a model, i.e., $(q+1)\mathrm{ln}n-2\mathrm{ln}{{ \mathcal L }}_{\max }$, where n is the number of data points. We select the optimal order q_opt if the 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 the noise model comparison is given by Feng et al. (2017).

Reflex motion describes the orbital motion of a primary star around the barycenter of the system. The full astrometric model for a star consists of three components: stellar reflex motion, proper motion, and the parallax of the barycenter of the target system. Though Gaia synthetic data can be generated by GOST, ¹⁹ 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's barycenter and the reflex motion. Considering a typical systematic RV of 10 au yr⁻¹ for the barycenter of a system, the change of parallax is about 0.001/d ϖ, where d is the heliocentric distance of the target system in parsecs. For a star 10 pc away from the Sun, the parallax change over 24 yr is about 0.01 mas, which is far below the current precision of Gaia parallaxes. Thus, we ignore the change of parallax.

For systems with both long-period (P > 1000 days) and short-period planets (P < 1000 days), we only use the Gaia–Hipparcos astrometry to constrain the long-period signals (LPSs), while leaving the inclination I and longitude of the ascending node Ω of the short-period companion unconstrained. As described in Feng et al. (2019a), we model the instantaneous proper motion and position of the target star at the Gaia epoch by combining the barycenter proper motion and the stellar reflex motion at the Gaia reference epoch. We then propagate the proper motion and position of the barycenter to the Hipparcos epoch, and add the stellar reflex motion at the Hipparcos reference epoch to model the proper motion and position of the target star at the Hipparcos epoch. In the calculation of likelihood, we use offsets and jitters to account for unknown systematics in the Gaia and Hipparcos catalog data.

To obtain posterior samples, we use the adaptive and parallel Markov Chain Monte Carlo (MCMC) developed by Haario et al. (2001) and Feng et al. (2019b). All time-related parameters, such as orbital periods and correlation timescale, follow a log-uniform prior distribution. The inclination I is sampled from a uniform prior distribution over $\sin I$. The other parameters have a uniform prior distribution. We use BIC > 10 or its equivalent ln (Bayes factor) > 5 (Kass & Raftery 1995; Feng et al. 2016) to determine whether additional companions are necessary to explain the data. To find signals efficiently, we first conduct an RV-only analysis, following the method developed by Feng et al. (2020). Then we apply the full modeling of RV and astrometry to the systems with P > 1000 day signals to fully constrain the dynamics of the system. For companions with shorter orbital periods (P < 1000 days), 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 astrometric data. Therefore, the combined analysis of both the RV and astrometric data for a multicompanion system would give reliable orbital solutions for all companions.

Among all of the available RV targets from the various RV surveys, we select 5108 stars, with each star having more than five high-precision RV data points. Following the methodology above (with ln (Bayes Factor) > 5; Kass & Raftery 1995; Feng et al. 2016), and using the RV data alone, we find that 869 of these stars show LPSs (P > 1000 days). From our combined analyses of RV and astrometry, 167 of them are confirmed as companions with masses from 5 M_Jup to 120 M_Jup, and with relative mass uncertainty of less than 100%.

The stellar mass and astrometry data used for the sample of 161 stars that host the 167 companions are shown in Table 2. For the companions that are directly imaged, the masses 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 Figure 1. Most of the hosts of companions are main-sequence AFGKM stars, while few hosts are subgiants.

Zoom In Zoom Out Reset image size

We also show the distribution of the stellar mass, G magnitude, number of RV points, effective temperature, η_r , η_μ , and heliocentric distance in Figure 2. In the figure, we divide the total sample of 161 companions into the known ones found in the literature and the new ones found in this work. The major difference between these two populations is that the new companions are typically identified from fewer RV data points than the known companions. This is expected, because the astrometric data help to constrain the LPSs, which the RV data alone cannot confirm. As shown in Figure 2, the majority of our sample is brighter than G = 8 mag and less than 100 pc away from the Sun. It is apparent that the η of the positional difference is one order of magnitude higher than the η of the proper-motion difference.

Zoom In Zoom Out Reset image size

Around the 161 hosts of cold substellar companions with masses ranging from 5 M_Jup to 120 M_Jup and orbital periods longer than 1000 days, we also find 63 other types of companions, including 60 planets, one substellar companion, and two stellar companions. To count the number of different types of multicompanion 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 Figure 3.

Zoom In Zoom Out Reset image size

Table 1. Number of Companions Identified in This Work

Relative Uncertainty	Full Sample	LPS (Systems)	Super-Jupiter (Systems)	BDs (Systems)	Stellar Companions (Systems)
			5–13 MJup	13–75 MJup	>75 MJup
${\sigma }_{{m}_{c}}/{m}_{c}\lt 100$%	5108	914 (869)	71 (68)	67 (66)	29 (29)
${\sigma }_{{m}_{c}}/{m}_{c}\lt 20$%	5108	914	40 (37)	49 (48)	24 (24)

Note. The "LPS" column 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 (${\sigma }_{{m}_{c}}/{m}_{c}\lt 20$%) and the optimal case, defined by <100% relative mass error (${\sigma }_{{m}_{c}}/{m}_{c}\lt 100$%).

Download table as:  ASCIITypeset image

The fitting results for targets with direct imaging data are shown in Figures 4 and 5, while example fittings for targets without imaging data are shown in Figures 6 and 7. For a target with imaging data that is shown in Figures 4 and 5, the first panel (from left to right) shows the optimal fit to the RV data. The second and third panels, respectively, show the fit of the reflex motion induced by wide-orbit companions (P > 1000 days) to the proper motions and positions at the Hipparcos and Gaia epochs, after subtracting the barycentric motion. The fourth panel shows the companion's binary orbit to its relative astrometry derived from imaging data. For the targets without imaging data that are shown in Figures 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 consisting of multiple wide-orbit companions, the reflex motion could be complex (e.g., HD 7449 in Figure 5, GJ 676 A and HD 105811 in Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We present Hipparcos and Gaia EDR3 astrometry for the stars in Table 2 and the orbital parameters for the 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, which 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.

Table 2. Hipparcos and Gaia EDR3 Catalog Astrometry for the Sample of Stars

Star	Mass	αH	δH	${\mu }_{\alpha }^{H}$	${\mu }_{\delta }^{H}$	αG	δG	${\mu }_{\alpha }^{G}$	${\mu }_{\delta }^{G}$
	(M⊙)	(deg)	(deg)	(mas yr−1)	(mas yr−1)	(deg)	(deg)	(mas yr−1)	(mas yr−1)
GJ 2030	0.96 ± 0.125	50.82374837	−7.79356453	2.2 ± 0.58	−219.27 ± 0.34	50.82376065334	−7.79507336136	1.295 ± 0.179	−219.292 ± 0.124
GJ 234 * †	0.195 ± 0.001	97.34581159	−2.81247675	705.28 ± 2.66	−611.92 ± 2.4	97.35069832373	−2.81701377835	750.14 ± 1.648	−802.947 ± 1.475
GJ 3222	0.89 ± 0.114	50.89679287	−40.07649183	40.76 ± 0.61	43.06 ± 0.62	50.89726944179	−40.07622273995	63.985 ± 0.024	35.175 ± 0.036
GJ 494 †	0.511 ± 0.004	195.19563106	12.37576423	−616.31 ± 1.51	− 13.59 ± 1	195.19112542751	12.37559492615	−628.715 ± 0.184	−33.472 ± 0.133
GJ 676A	0.626 ± 0.063	262.54769928	−51.63652557	−260.02 ± 1.34	−184.29 ± 0.82	262.54483233204	−51.63780340759	−258.759 ± 0.034	−185.119 ± 0.025
GJ 680	0.456 ± 0.046	263.80642761	−48.68197465	82.62 ± 3.84	454.25 ± 2.84	263.80723698262	−48.67877743079	74.069 ± 0.024	470.166 ± 0.015
GJ 864	0.576 ± 0.058	339.04021326	−0.84015108	48.8 ± 2.01	−628.16 ± 1.46	339.04061610661	−0.84440195456	55.904 ± 0.04	−630.315 ± 0.044
GJ 9714	0.67 ± 0.076	315.41221712	−32.5241822	244.41 ± 1.74	−121.19 ± 1.23	315.41421729626	−32.52501976182	246.267 ± 0.019	−122.028 ± 0.018
HD 100939	1.75 ± 0.101	174.20069282	−37.03903769	25.23 ± 0.58	−15.59 ± 0.46	174.2008982268	−37.03913688429	23.86 ± 0.024	−14.544 ± 0.016
HD 105618	1.01 ± 0.129	182.40284617	11.21158266	−103.19 ± 1.35	−2.33 ± 0.64	182.40211396767	11.21156927773	−104.35 ± 0.025	−2.024 ± 0.02

Note. The masses of the stars are taken from the TESS input catalog (Stassun et al. 2019). The reference epochs of the Gaia and Hipparcos data are J1991.25 and J2016.0, respectively. The superscripts "gaia" and "hip" are used to denote, respectively, the Gaia EDR3 and Hipparcos R.A. (α), decl. (δ), parallax (ϖ), proper motion in R.A. (μ_α ), and proper motion in decl. (μ_δ ). The hosts with a "†" superscript have been directly imaged, and the stellar mass is inferred from the relative astrometry a posteriori. For a star without EDR3 data, but with DR2 data, a star symbol is added behind its name.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 3. Orbital Parameters for the Stellar and Planetary Companions Detected in This Work

Companion	Host	P	I	e	mc	a	HZ Stability	Reference
Name	Name	(yr)	(deg)		(MJup)	(au)
GJ 2030 b	HD 21019	${0.007}_{-4\times {10}^{-7}}^{+4\times {10}^{-7}}$	⋯	${0.239}_{-0.058}^{+0.075}$	${0.015}_{-0.002}^{+0.002}$	${0.034}_{-0.002}^{+0.001}$	Stable
GJ 2030 c	HD 21019	${69.906}_{-6.641}^{+6.979}$	${16.999}_{-2.535}^{+2.935}$	${0.041}_{-0.008}^{+0.003}$	${12.803}_{-2.136}^{+2.335}$	${16.761}_{-1.337}^{+1.294}$	Stable
GJ 234 B	Ross 614	${16.586}_{-0.004}^{+0.004}$	${52.918}_{-0.016}^{+0.016}$	${0.382}_{-1\times {10}^{-4}}^{+1\times {10}^{-4}}$	${94.837}_{-1.370}^{+0.880}$	${4.187}_{-0.009}^{+0.008}$	Stable	1, 2, 3
GJ 3222 b	HD 21175⋆	${0.029}_{-6\times {10}^{-6}}^{+8\times {10}^{-6}}$	⋯	${0.929}_{-0.080}^{+0.042}$	${0.036}_{-0.006}^{+0.011}$	${0.091}_{-0.004}^{+0.004}$	Unstable
GJ 3222 c	HD 21175⋆	${3.915}_{-0.015}^{+0.016}$	${2.081}_{-0.160}^{+0.186}$	${0.537}_{-0.038}^{+0.065}$	${52.213}_{-4.946}^{+5.113}$	${2.391}_{-0.109}^{+0.097}$	Unstable
GJ 494 B†	BD+132618	${13.678}_{-0.036}^{+0.032}$	${130.677}_{-0.200}^{+0.203}$	${0.245}_{-0.001}^{+0.001}$	${88.918}_{-2.844}^{+1.836}$	${4.757}_{-0.016}^{+0.014}$	Stable	1, 4, 5, 6, 7, 8
GJ 676 A d	CD-5110924	${0.010}_{-6\times {10}^{-7}}^{+5\times {10}^{-7}}$	⋯	${0.192}_{-0.079}^{+0.073}$	${0.012}_{-1\times {10}^{-3}}^{+0.001}$	${0.039}_{-0.001}^{+0.001}$	Unstable	9, 10, 11
GJ 676 A e	CD-5110924	${0.097}_{-4\times {10}^{-5}}^{+4\times {10}^{-5}}$	⋯	${0.152}_{-0.060}^{+0.086}$	${0.021}_{-0.002}^{+0.002}$	${0.181}_{-0.006}^{+0.006}$	Unstable	9, 10, 11
GJ 676 A b	CD-5110924	${2.879}_{-0.001}^{+0.001}$	${48.919}_{-2.781}^{+3.312}$	${0.319}_{-0.003}^{+0.003}$	${5.792}_{-0.477}^{+0.469}$	${1.735}_{-0.060}^{+0.056}$	Unstable	9, 10, 11
GJ 676 A c	CD-5110924	${38.115}_{-4.157}^{+3.391}$	${33.690}_{-1.324}^{+1.362}$	${0.295}_{-0.049}^{+0.033}$	${13.492}_{-1.127}^{+1.046}$	${9.726}_{-0.793}^{+0.629}$	Unstable	9, 10, 11
GJ 680 b	CD-4811837⋆	${47.315}_{-10.730}^{+13.232}$	${21.984}_{-5.062}^{+3.313}$	${0.381}_{-0.145}^{+0.095}$	${25.100}_{-11.149}^{+6.158}$	${10.138}_{-1.696}^{+1.836}$	Stable
GJ 864 b	HD 214100	${14.236}_{-1.142}^{+1.070}$	${15.788}_{-1.064}^{+0.257}$	${0.524}_{-0.044}^{+0.041}$	${59.962}_{-6.078}^{+8.858}$	${4.931}_{-0.332}^{+0.296}$	Stable
GJ 9714 b	HD 199981⋆	${52.678}_{-12.360}^{+22.871}$	${155.326}_{-5.041}^{+4.804}$	${0.155}_{-0.041}^{+0.073}$	${9.403}_{-1.695}^{+3.055}$	${12.319}_{-2.230}^{+3.623}$	Stable
HD 100939 b	HD 100939	${7.325}_{-0.298}^{+0.502}$	${85.093}_{-43.671}^{+49.211}$	${0.118}_{-0.076}^{+0.062}$	${6.770}_{-1.213}^{+1.859}$	${4.550}_{-0.150}^{+0.221}$	Stable
HD 105618 b	HD 105618	${0.027}_{-3\times {10}^{-6}}^{+4\times {10}^{-6}}$	⋯	${0.455}_{-0.038}^{+0.041}$	${0.080}_{-0.008}^{+0.007}$	${0.091}_{-0.004}^{+0.004}$	Stable
HD 105618 c	HD 105618	${211.127}_{-28.121}^{+63.068}$	${31.968}_{-3.897}^{+3.441}$	${0.126}_{-0.080}^{+0.098}$	${26.491}_{-2.827}^{+4.340}$	${35.731}_{-3.639}^{+6.866}$	Stable
HD 105811 b	HD 105811	${3.218}_{-0.020}^{+0.019}$	${0.474}_{-0.075}^{+0.177}$	${0.766}_{-0.045}^{+0.071}$	${55.534}_{-11.899}^{+13.075}$	${2.699}_{-0.052}^{+0.054}$	Unstable
HD 105811 B	HD 105811	${6.763}_{-0.012}^{+0.004}$	${81.410}_{-8.027}^{+6.232}$	${0.608}_{-0.003}^{+0.006}$	${206.640}_{-7.718}^{+9.641}$	${4.574}_{-0.089}^{+0.082}$	Unstable
HD 106515 A b	HD 106515 A⋆	${9.820}_{-0.093}^{+0.162}$	${88.922}_{-39.902}^{+42.184}$	${0.565}_{-0.025}^{+0.021}$	${9.518}_{-0.128}^{+6.385}$	${4.491}_{-0.187}^{+0.187}$	Stable	12, 13, 14
HD 10697 b	109 Psc	${2.944}_{-0.002}^{+0.002}$	${86.116}_{-20.530}^{+19.957}$	${0.104}_{-0.008}^{+0.009}$	${5.743}_{-0.289}^{+1.011}$	${2.051}_{-0.087}^{+0.079}$	Stable	15, 16, 17, 18
HD 109988 b	HD 109988⋆	${46.905}_{-5.998}^{+10.794}$	${33.992}_{-4.296}^{+1.773}$	${0.056}_{-0.034}^{+0.036}$	${21.821}_{-4.155}^{+2.265}$	${12.518}_{-1.319}^{+1.517}$	Stable
HD 110537 b	HD 110537	${66.994}_{-21.792}^{+27.441}$	${51.891}_{-5.422}^{+9.617}$	${0.401}_{-0.152}^{+0.115}$	${27.099}_{-8.990}^{+12.755}$	${16.714}_{-3.907}^{+4.538}$	Stable
HD 111031 b	HD 111031	${46.740}_{-4.870}^{+3.194}$	${169.682}_{-0.582}^{+0.404}$	${0.031}_{-0.012}^{+0.001}$	${54.167}_{-6.149}^{+5.319}$	${13.101}_{-1.092}^{+0.750}$	Stable	19

Note. The systems are sorted by the host names listed in the first column, while the companions in each system are sorted by their orbital periods. For companions with orbital periods of less than 1000 days, the inclination is not given and the mass reported in this table is ${m}_{c}\sin I$. For companions with mass higher than 75 M_Jup, we use capital letters to label them. Among the parameters, the directly inferred parameters are P (orbital period), K (RV semi-amplitude), e (eccentricity), ω_⋆ (argument of periastron), I (inclination), and Ω (longitude of ascending node). The derived parameters are T_p (periastron epoch), m_c (companion mass), and a (semimajor axis). The dynamical stabilities of Earth-like planets in the habitable zones (HZs) of their host stars are determined numerically, as described in Section 7. The median and 1σ quantiles (i.e., the 16% and 84% quantiles) are used to measure the uncertainty of each parameter. The companions with a "†" superscript have been imaged and their relative astrometry data are analyzed in combination with other types of data in the solutions. The targets with a "⋆" superscript have wide-binary companions that are identified by El-Badry et al. (2021) or other studies.

References. (1) Mann et al. (2019); (2) Agati et al. (2015); (3) Bonavita et al. (2016); (4) Burgasser et al. (2010); (5) Heintz (1994); (6) Beuzit et al. (2004); (7) Goldman et al. (2010); (8) Dupuy & Kraus (2013); (9) Forveille et al. (2011); (10) Stassun et al. (2017); (11) Anglada-Escudé & Tuomi (2012); (12) Marmier et al. (2013); (13) Desidera et al. (2012); (14) Li et al. (2021); (15) Zucker & Mazeh (2000); (16) Simpson et al. (2010); (17) Wittenmyer et al. (2009); (18) Vogt et al. (2000); (19) Hinkel et al. (2019); (20) Mayor et al. (2004); (21) Borgniet et al. (2014); (22) Xuan et al. (2020); (23) Borgniet et al. (2019); (24) Jones et al. (2016); (25) Fischer et al. (2007); (26) Ment et al. (2018); (27) Giguere 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) Howard et al. (2010); (33) Martioli et al. (2010); (34) Fischer et al. (2002); (35) Rickman et al. (2019); (36) Rickman et al. (2020); (37) Döllinger et al. (2009); (38) Wang et al. (2014); (39) Tinney et al. (2002); (40) Wittenmyer et al. (2012); (41) Wittenmyer et al. (2020); (42) Butler et al. (2003); (43) Wittenmyer et al. (2007); (44) Feng et al. (2015); (45) Díaz et al. (2012); (46) Patel et al. (2007); (47) Marcy et al. (1999); (48) Pilyavsky et al. (2011); (49) Naef et al. (2001); (50) Jones et al. (2021); (51) Arriagada et al. (2010); (52) Thalmann et al. (2009); (53) Vigan et al. (2016); (54) Brandt et al. (2021b); (55) Bowler et al. (2018); (56) Marcy et al. (2005); (57) Fuhrmann (2004); (58) Tokovinin (2014a); (59) Roberts et al. (2015a); (60) Perrier et al. (2003); (61) Sahlmann et al. (2011); (62) Han et al. (2001); (63) Reffert & Quirrenbach (2011); (64) Díaz et al. (2016); (65) Ségransan et al. (2010); (66) Robertson et al. (2012); (67) Moutou et al. (2017); (68) Ginski et al. (2016); (69) Moutou et al. (2011); (70) Kane et al. (2019); (71) Melo et al. (2007); (72) Venner et al. (2021); (73) Jenkins et al. (2017); (74) Rosenthal et al. (2021); (75) Moutou et al. (2005); (76) Trifonov et al. (2017); (77) Anglada-Escudé et al. (2010); (78) Kürster et al. (2015); (79) Santos et al. (2001); (80) Tokovinin (2014b); (81) Griffin (2012); (82) Bouchy et al. (2016); (83) Wittenmyer et al. (2017); (84) Zurlo et al. (2018); (85) Fischer et al. (2001); (86) Lagrange et al. (2009); (87) Nielsen et al. (2014); (88) Lagrange et al. (2019); (89) Lacour et al. (2021); (90) Jones et al. (2002); (91) Huang et al. (2018); (92) Damasso et al. (2020); (93) Jenkins et al. (2015); (94) Tamuz et al. (2008); (95) Cheetham et al. (2018); (96) Mugrauer et al. (2014); (97) Feng et al. (2020); (98) Nakajima et al. (1995); (99) Naef et al. (2010); (100) Allen et al. (2012); (101) Butler et al. (2006); (102) Barbato et al. (2018); (103) Moutou et al. (2009); (104) Naef et al. (2004); (105) Dumusque et al. (2011); (106) Wittenmyer et al. (2019); (107) Rodigas et al. (2016); (108) Demangeon et al. (2021); (109) Sozzetti et al. (2006); (110) Bang et al. (2020); (111) Fischer et al. (2009); (112) Hartmann et al. (2010); (113) Johnson-Groh et al. (2017); (114) Wagner et al. (2019); (115) Franson et al. (2022); (116) Jones et al. (2015c); (117) Jones et al. (2017); (118) Jones et al. (2015); and (119) Mason et al. (2001).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

For previously published planets, we discuss our new results on an object-by-object basis as follows.

• 1.  
  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}_{-1.37}^{+0.88}$ M_Jup.
• 2.  
  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_⊙ for GJ 494 A and a mass of ${88.92}_{-2.84}^{+1.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).
• 3.  
  GJ 676 is a binary and the primary hosts four planets (Forveille et al. 2011; Anglada-Escudé & Butler 2012). Our solution constrains the inclinations of GJ 676 A b and c to be ${48.92}_{-2.78}^{+3.31}$° and ${33.69}_{-1.32}^{+1.36}$°, respectively. Our solutions determine their masses as ${5.79}_{-0.48}^{+0.47}$ M_Jup and ${13.49}_{-1.13}^{+1.05}$ 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 in warping the planetary disk, such that the orbits of GJ 676 A d and e are slightly misaligned.
• 4.  
  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}_{-0.00}^{+3.309}$ M_Jup, putting it into the super-Jupiter category.
• 5.  
  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}_{-0.13}^{+6.39}$ M_Jup.
• 6.  
  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}_{-0.29}^{+1.01}$ M_Jup.
• 7.  
  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}_{-6.15}^{+5.32}$ M_Jup.
• 8.  
  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}_{-0.48}^{+1.13}$ M_Jup. We also find a BD companion with a mass of ${18.06}_{-1.61}^{+4.21}$ M_Jup on a 17 au wide orbit. The HD 111232 system is one of a few systems hosting multiple substellar companions in our sample.
• 9.  
  HD 113337 hosts a disk with an inclination of about ${13}_{-9}^{+10}$° (Xuan et al. 2020). Xuan et al. (2020) estimate an inclination of ${31}_{-4}^{+5}$°, based on analyses of Gaia DR2 and Hipparcos data. However, our solution determines an inclination of ${57.48}_{-3.65}^{+4.55}$°. Such a discrepancy is probably due to our use of the EDR3 catalog.
• 10.  
  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 11343c, with a minimum mass of about 0.5 M_Jup. The astrometric signal induced by HD 11343c is weak, and thus the astrometric data does not constrain its inclination well. The mass of HD 11343 b is ${7.71}_{-1.19}^{+0.73}$ M_Jup.
• 11.  
  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.5}^{+7.7}$ yr and a mass of ${7.38}_{-1.09}^{+2.02}$ M_Jup. This makes the system unique in terms of a solar analog hosting a warm Saturn and two CJs.
• 12.  
  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}_{-0.31}^{+4.62}$ M_Jup.
• 13.  
  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 122562c, with an orbital period of 84 days. On the basis that there seem to be other planets in the system, we remain consistent with the previous naming of the more massive companion, and label it HD 122562 b. Its mass is ${37.20}_{-2.60}^{+3.67}$ M_Jup.
• 14.  
  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}_{-0.31}^{+4.13}$ M_Jup.
• 15.  
  HD 125612 is a G star hosting three companions with minimum masses of 3.0 M_Jup, 0.058 M_Jup, and 7.2 M_Jup and orbital periods of 502, 4.1547, and 3008 days, respectively (Fischer et al. 2007; Lo Curto et al. 2010; Ment et al. 2018). Through combined analyses of the RV and astrometric data, we are able to constrain the mass of the Jupiter-like planet, HD 125612 d, to ${7.18}_{-0.45}^{+0.93}$ M_Jup.
• 16.  
  HD 126614 is a binary consisting of a G-type primary star and an 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.02}^{+0.20}$ M_Jup for the planetary companion and ${81.13}_{-7.92}^{+7.78}$ M_Jup for the M-dwarf companion.
• 17.  
  HD 127506 (GJ 554) is a K-type star hosting a BD, 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.53}^{+3.70}$ 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.
• 18.  
  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 allows us to break the degeneracy between mass and inclination, and confirm this candidate as a BD with a mass of ${76.06}_{-14.60}^{+16.22}$ M_Jup.
• 19.  
  HD 136118 (HIP 74948) is an F-type star hosting a BD with a mass of ${42}_{-18}^{+11}$ M_Jup (Fischer et al. 2002; Wittenmyer et al. 2009). Our solution estimates a mass of ${13.10}_{-1.27}^{+1.35}$ M_Jup.
• 20.  
  HD 13724 (HIP 10278) is a G-type star hosting HD 13724 b or (HD 13724 B), a BD with a mass of ${50.5}_{-3.5}^{+3.3}$ M_Jup (Rickman et al. 2020). We estimate a mass of ${36.32}_{-1.60}^{+1.48}$ 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.
• 21.  
  HD 139357 is a K-type star hosting a companion with a minimum mass of 9.76 ± 2.15 M_Jup (Döllinger et al. 2009). With combined analyses, we confirm this companion to be a BD with a mass of ${16.38}_{-0.00}^{+7.88}$ M_Jup.
• 22.  
  HD 14067 is a G9III star hosting a substellar 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}_{-0.00}^{+13.30}$ M_Jup.
• 23.  
  HD 142 is a wide binary consisting of an 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}_{-0.94}^{+1.28}$ M_Jup.
• 24.  
  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}_{-0.88}^{+1.63}$ M_Jup, ${m}_{c}={5.03}_{-1.07}^{+0.87}$ M_Jup, ${I}_{b}={144.65}_{-3.24}^{+6.28}$° and ${I}_{c}\,={129.10}_{-29.05}^{+6.26}$° 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}_{-3.2}^{+5.3}$° and ${I}_{c}={101}_{-33}^{+31}$°, 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 two data points).
• 25.  
  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}_{-13.93}^{+70.23}$ 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).
• 26.  
  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}_{-4.40}^{+5.11}$ M_Jup, after considering the stellar mass uncertainty. Its inclination is ${95.20}_{-9.86}^{+7.47}$°, suggesting a nonzero transit probability.
• 27.  
  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.09}^{+0.08}$ M_Jup for this companion. We also find a stellar companion with a mass of ${101.79}_{-6.11}^{+5.91}$ M_Jup.
• 28.  
  HD 181234 (HIP 95015) is a G-type star hosting a companion with a minimum mass of ${8.37}_{-0.36}^{+0.34}$ M_Jup (Rickman et al. 2019). Based on combined RV and astrometry analysis, we constrain the mass of this companion to be ${8.24}_{-0.00}^{+2.15}$ M_Jup and the inclination to be ${98.73}_{-41.61}^{+24.47}$°.
• 29.  
  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 has been further studied by multiple groups (e.g., Bowler et al. 2018; 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 used by Brandt et al. (2021b), as well as the Gaia–Hipparcos positional difference, we constrain the mass to be ${36.39}_{-1.08}^{+1.21}$ M_Jup and the eccentricity to be ${0.37}_{-0.08}^{+0.08}$, consistent with the solution given by Brandt et al. (2021b). With the astrometry for both components of the binary, we also estimate a mass of 0.93 ± 0.03 M_⊙ for the primary a posteriori.
• 30.  
  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.55}^{+9.03}$ M_Jup for this unresolved companion.
• 31.  
  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 45 ± 21, 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 analysis of the intermediate data in a future Gaia data release will strongly constrain the mass.
• 32.  
  HD 191806 is a K star hosting a companion with a minimum mass of 8.52 ± 0.63 M_Jup (Díaz et al. 2016). Our combined analyses constrain the mass to be ${9.33}_{-0.85}^{+0.92}$ M_Jup.
• 33.  
  HD 203473 B was discovered by Ment et al. (2018), based on analysis of Keck RV data. With both the 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.02}^{+0.01}$ 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 a limited RV baseline.
• 34.  
  HD 204313 is claimed to host three planets, based on the RV data collected by the 2.7 m HJS and the CORALIE spectrometer (Ségransan et al. 2010; Robertson et al. 2012). However, Díaz et al. (2016) do not confirm HD 204313 d, based on their analyses of previous and new HARPS data. Our solution agrees with Díaz et al. (2016) that there is no significant signal around a period of 2800 days, as claimed by Robertson et al. (2012). However, we confirm the trend identified by Díaz et al. (2016) as a new LPS, with a period of ${20.06}_{-1.01}^{+1.10}$ yr, based on combined RV and astrometry analysis. This new companion is dubbed "HD 204313 e," to distinguish it from the dubious candidate HD 204313 d. It corresponds to a BD with a mass of ${15.32}_{-5.18}^{+4.89}$ M_Jup.
• 35.  
  HD 211847 is a G star hosting a substellar 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}_{-18.49}^{+1.34}$ M_Jup.
• 36.  
  HD 214823 (HIP 111928) is a G-type star hosting a companion (Díaz et al. 2016; Ment et al. 2018) with a minimum mass of 20.56 ± 0.32 M_Jup and an orbital period of 1853.9 ± 1.6 days (Luhn et al. 2019). With Gaia and Hipparcos data, we are able to break the degeneracy between mass and inclination, and constrain the absolute mass to ${18.61}_{-1.07}^{+4.14}$ M_Jup and the period to be ${5.078}_{-0.004}^{+0.004}$ yr.
• 37.  
  HD 217786 (HIP 113834) is an 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.31}^{+1.27}$ M_Jup. We also find a hot super-Earth with an orbital period of 2.5 days.
• 38.  
  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 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}_{-38.20}^{+34.17}$ M_Jup and a planet with a mass of ${0.52}_{-0.06}^{+0.35}$ M_Jup on a Jupiter-like orbit.
• 39.  
  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 give a dynamical mass of ${9.62}_{-0.73}^{+1.00}$ M_Jup.
• 40.  
  HD 219828 (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}_{-1.27}^{+3.52}$ M_Jup for the bigger companion on a Jupiter-like orbit.
• 41.  
  HD 221420 (GJ 4340) is a G-type star hosting a BD (Kane et al. 2019) with a wide-orbit M-dwarf companion (El-Badry et al. 2021; Venner et al. 2021). Based on the combined RV and HGCA catalog of calibrated Gaia DR2 and Hipparcos proper motions (Brandt 2018), Venner et al. (2021) estimated a mass of ${20.35}_{-3.43}^{+1.99}$ M_Jup, consistent with the 20 ± 3 M_Jup estimated in this work.
• 42.  
  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}_{-0.35}^{+1.84}$ M_Jup.
• 43.  
  HD 23596 (HIP 17747) is an F star hosting a substellar 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}_{-1.77}^{+0.99}$ M_Jup.
• 44.  
  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.81}^{+1.16}$ M_Jup.
• 45.  
  HD 26161 (HIP 19428) is a G-type star hosting a companion with a minimum mass of ${13.5}_{-3.7}^{+8.5}$ M_Jup (Rosenthal et al. 2021). Based on the combined RV and astrometry analysis, the dynamical mass of this companion is ${28.46}_{-0.22}^{+20.05}$ 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).
• 46.  
  HD 27894 (HIP 20277) is a K star hosting three companions with masses of 5.42 M_Jup, 0.16 M_Jup, and 0.67 M_Jup (Moutou et al. 2005; Anglada-Escudé et al. 2010; Kürster et al. 2015; Trifonov et al. 2017). However, we only find strong evidence for the biggest two companions. By constraining the inclination of the biggest companion using astrometry, we find a dynamical mass of ${6.49}_{-0.35}^{+0.99}$ M_Jup.
• 47.  
  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 days (Minniti et al. 2009). With both the RV and astrometry data, we constrain the minimum mass of HD 28185 b to be ${5.837}_{-0.510}^{+0.486}$ M_Jup. Moreover, we detect a BD companion with a mass of ${19.64}_{-2.14}^{+2.27}$ M_Jup.
• 48.  
  HD 28192 (HIP 20752) is a G star hosting a stellar companion on an extremely wide orbit and with a mass of about 0.4 M_⊙ (Tokovinin 2014b; El-Badry et al. 2021). This companion induces a reflex motion of 0.1 m s⁻¹ yr⁻¹ at most, and is thus not considered in our combined analyses. We find another stellar companion with a mass of ${94.78}_{-7.52}^{+8.51}$ M_Jup on a 10 au wide orbit, and a planet with a minimum mass of ${0.31}_{-0.03}^{+0.02}$ M_Jup and an orbital period of about 14 days.
• 49.  
  HD 29461 (HIP 21654) is a G star hosting a companion with a minimum mass of about 0.08 M_⊙ (Griffin 2012; Bouchy et al. 2016). Our combined analyses estimate a mass of ${92.96}_{-4.93}^{+12.61}$ M_Jup on a 5 au wide orbit. The small separation between this companion and the primary star explains the null detection of it by previous imaging surveys.
• 50.  
  HD 30177 (HIP 21850) is a G star hosting two companions with minimum masses of 3 ± 0.3 M_Jup and 8.07 ± 0.12 M_Jup (Tinney et al. 2002; Wittenmyer et al. 2017; Barbato et al. 2018). The null detections 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}_{-0.49}^{+1.24}$ M_Jup and ${6.15}_{-0.34}^{+1.31}$ M_Jup, consistent with the constraints given by previous imaging and RV data analyses. This system is one of a few systems hosting two super-Jupiters in our sample.
• 51.  
  HD 38529 (HIP 27253) is a G star hosting two companions with masses of 18 ± 3 M_Jup and 0.17 M_Jup (Fischer et al. 2003; Wittenmyer et al. 2009; Benedict et al. 2010; Barbato et al. 2018; Xuan et al. 2020). Our combined analyses estimate a mass of ${10.38}_{-0.88}^{+1.03}$ M_Jup and an inclination of ${104.56}_{-8.72}^{+6.39}$° for HD 38529c, consistent with the the solution given by Xuan et al. (2020), but with higher precision.
• 52.  
  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 as 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 motions due to all companions are modeled; and (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 has also been used by Brandt et al. (2021a). We find that our solution is quite sensitive to the RV data. Although both the RV data sets used are from the same HARPS observations, different corrections are made to the stellar activity, particularly from pulsations. We use the Lagrange et al. (2020) data set (dubbed "AL20") to find the mass of β Pic b to be ${7.56}_{-1.69}^{+1.35}$ M_Jup and that of c to be ${8.94}_{-0.78}^{+0.75}$ M_Jup, while using the Vandal et al. (2020) data set (dubbed "TV20") we find masses of ${11.73}_{-2.14}^{+2.34}$ M_Jup and ${10.14}_{-1.03}^{+1.18}$ M_Jup. 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 for β Pic b measured by Lagrange et al. (2020) and Nowak et al. (2020), respectively, 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 the interferometric data of 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 for resolving the discrepancies that we find in the mass measurements from different RV data sets.
• 53.  
  HD 39091 (π Men) is a G star hosting at least two companions, with masses of 0.015 M_Jup and 13 M_Jup (Jones et al. 2002; Huang et al. 2018; Damasso et al. 2020), 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. 2020; De Rosa et al. 2020; Xuan et al. 2020; Kunovac Hodžić et al. 2021). With combined analyses of the RV data from AAT, CORALIE, ESPRESSO, HARPS, and PFS, as well as the Gaia–Hipparcos data, we are able to constrain the mass of this companion to be ${12.33}_{-1.38}^{+1.19}$ M_Jup and the inclination to be ${54.44}_{-3.72}^{+5.94}$°. This inclination of π Men b differs from the 90° inclination of π Men c by 6σ, suggesting a significant misalignment, as proposed by previous studies. However, we fail to confirm the third companion found by Hatzes et al. (2022).
• 54.  
  HD 39213 (HIP 27491) is a K star hosting a BD with a minimum mass of 0.07 ± 0.01 M_⊙ (Jenkins et al. 2015). Our combined analyses determine a dynamical mass of ${70.77}_{-5.43}^{+10.21}$ 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.
• 55.  
  HD 4113 A (HIP 3391) is a Sun-like star hosting a BD (HD 4113 A b or HD 4113 b) and a substellar 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.075}^{+0.076}$ M_Jup, based on RV analysis (Cheetham et al. 2018). HD 4113 C has a dynamical mass of ${65.8}_{-4.4}^{+5.0}$ M_Jup, based on analyses of direct imaging data (labeled "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.46}^{+0.60}$ M_Jup, relaxing the previous tension between dynamical and isochronal masses, without invoking binarity in HD 4113 C.
• 56.  
  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. 2020). 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 data, imaging data (labeled "MB21"), and Gaia EDR3 and Hipparcos data. This mass is in tension with the mass predicted by cooling models. However, with nearly the same data, but with an additional constraint from the Hipparcos–Gaia positional difference, our combined analyses estimate a mass of ${60.42}_{-2.38}^{+2.34}$ 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 for avoiding potential bias when using proper-motion difference alone.
• 57.  
  HD 43197 is a Sun-like star hosting a warm Jupiter (Naef et al. 2010). We detect a cold super-Jupiter, HD 43197c, 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}_{-3.07}^{+5.39}$°). 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.
• 58.  
  HD 65430 (HIP 39064) is a spectroscopic binary with an orbital period of 3138 days (Allen et al. 2012). Our combined analyses estimate a mass of ${105.40}_{-8.95}^{+8.37}$ M_Jup, confirming its stellar origin.
• 59.  
  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 66428c. Compared with the orbital period of ${107}_{-49}^{+153}$ yr and minimum mass of ${0.085}_{-0.053}^{+0.069}$ M_Jup estimated for HD 66428c by Rosenthal et al. (2021), our estimation of a period of ${28.69}_{-5.35}^{+9.21}$ yr and ${1.76}_{-0.04}^{+3.40}$ M_Jup is much more precise. However, further constraint on the inclination of HD 66428c is needed to determine whether the orbits of the two companions are misaligned.
• 60.  
  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}_{-0.10}^{+2.59}$ M_Jup. We also find a new companion, HD 72659c, with a mass of ${18.81}_{-4.80}^{+4.44}$ M_Jup. Assuming a coplanar configuration between b and c, HD 72659 b would have a mass of about 15 M_Jup.
• 61.  
  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}_{-35.48}^{+41.76}$ 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.
• 62.  
  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.28}^{+0.91}$ M_Jup and an orbital period of ${46.74}_{-2.98}^{+2.15}$ yr. The orbits of the two companions are probably misaligned, though with large uncertainty.
• 63.  
  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.77}^{+5.61}$ M_Jup.
• 64.  
  HD 74156 (HIP 42723) is a G star hosting two substellar companions with masses of 1.78 M_Jup and 8.00 M_Jup (Naef et al. 2004; Wittenmyer et al. 2009). Our combined analyses lead to estimated masses of ${1.71}_{-0.06}^{+0.96}$ M_Jup and ${8.67}_{-0.47}^{+1.39}$ M_Jup.
• 65.  
  HD 7449 A b and HD 7449 A c (or HD 7449 B) were 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}_{-3.74}^{+2.61}$° and ${68.40}_{-3.89}^{+4.10}$° for Ab and B, indicating significant misalignment between the two companions. The mass of HD 7449 A b is ${8.17}_{-2.70}^{+3.06}$ M_Jup. The mass of HD 7449 B is ${178.15}_{-13.66}^{+16.61}$ M_Jup, consistent with the value of ${0.23}_{-0.05}^{+0.22}$ M_⊙ estimated by Rodigas et al. (2016) based on photometry.
• 66.  
  HD 80869 (HIP 46022) is a G star hosting a planetary companion with a minimum mass of ${4.86}_{-0.29}^{+0.65}$ M_Jup (Demangeon et al. 2021). The use of Gaia–Hipparcos astrometry allows us to break the inclination–mass degeneracy and estimate a mass of ${5.07}_{-0.56}^{+2.54}$ M_Jup.
• 67.  
  HD 81040 (HIP 46076) is a G star hosting a Jupiter-like companion with a mass of ${7.24}_{-0.37}^{+1.0}$ M_Jup (Sozzetti et al. 2006; Stassun et al. 2017; Li et al. 2021). Our combined analyses give a mass of ${6.77}_{-0.87}^{+1.10}$ 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).
• 68.  
  HD 81817 is a K-type star hosting a substellar companion (HD 81817 b) with a minimum mass of 27.1 M_Jup (Bang et al. 2020). With both RV and astrometry data, we are able to constrain its mass to ${24.13}_{-0.71}^{+9.83}$ 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 81817c as a BD, because this signal shows a unique power (see Figure 8) in the Bayes factor periodogram and is strictly periodic and quite circular based on MCMC posterior samplings. The evidence strongly supports a Keplerian origin, instead of an activity origin, though its inclination is not well constrained due to its short orbital period.
• 69.  
  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}_{-1.95}^{+11.71}$ M_Jup.
• 70.  
  HD 8673 (HIP 6702) is a double star system, including an F star and an early M dwarf with a mass of 0.33–0.45 M_⊙ (Roberts et al. 2015b). The primary F star hosts a substellar companion with a minimum mass of 14.2 ± 1.6 M_Jup. Our combined analyses constrain the mass of the substellar companion to be ${13.25}_{-1.42}^{+1.70}$ M_Jup.
• 71.  
  HD 87883 (HIP 49699) is a K star hosting a companion with a minimum mass of ${6.31}_{-0.32}^{+0.31}$ M_Jup (Fischer et al. 2009; Stassun et al. 2017; Li et al. 2021). Our solution estimates a mass of ${5.41}_{-0.64}^{+0.47}$ M_Jup. Compared with the bimodal posterior distribution of inclination given by Li et al. (2021), based on the RV and Gaia–Hipparcos proper-motion difference, we use both the proper-motion and positional differences between Gaia and Hipparcos so that we are able to break the degeneracy between I and I + π, constraining the inclination to be ${25.45}_{-1.05}^{+1.61}$°.
• 72.  
  HD 95544 (HIP 54203) is a G star hosting a companion with a minimum mass of 6.84 ± 0.31 M_Jup (Demangeon et al. 2021). According to our combined analyses of RV and astrometric data, the dynamical mass is ${6.02}_{-0.26}^{+1.62}$ M_Jup and the inclination is ${86.50}_{-24.60}^{+31.19}$°, indicating an edge-on configuration.
• 73.  
  HD 984 (or HIP 1134) is an F star hosting a BD with a mass of 61.0 ± 4.0 M_Jup (Johnson-Groh et al. 2017; Franson et al. 2022). Through combined analyses of the RV, Gaia–Hipparcos, and the imaging data collected by Franson et al. (2022), which is labeled "KF22," we find a dynamical mass of ${40.37}_{-18.27}^{+24.33}$ M_Jup.
• 74.  
  HD 98649 (HIP 55409) is a G-type star hosting a companion with a minimum mass of ${6.79}_{-0.3}^{+0.5}$ M_Jup (Marmier et al. 2013; Rickman et al. 2019). Through combined analyses of RV and Gaia–Hipparcos astrometry, Li et al. (2021) estimate a mass of ${9.7}_{-1.9}^{+2.3}$ M_Jup. Using both proper-motion and positional differences between Gaia and Hipparcos, we estimate a mass of ${6.76}_{-0.00}^{+3.61}$ M_Jup, consistent with and more precise than the mass given by Li et al. (2021).
• 75.  
  HIP 22203 (HD 30246) is a G star hosting a BD with a minimum mass of ${55.1}_{-8.2}^{+20.3}$ M_Jup (Díaz et al. 2012). Our combined analyses constrain the mass to be ${50.76}_{-4.54}^{+4.73}$ M_Jup and the inclination to be ${84.36}_{-2.56}^{+8.68}$°.
• 76.  
  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 M_Jup to 713.7 M_Jup for HD 117253 b. However, based on Gaia–Hipparcos data, we find a mass of ${5.89}_{-0.28}^{+0.76}$ M_Jup and an inclination of ${88.97}_{-17.33}^{+16.21}$°. 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 substellar or stellar.
• 77.  
  HIP 67537 (HD 120457) is a red giant branch star hosting a planet with a minimum mass of ${11.1}_{-1.1}^{+0.4}$ M_Jup (Jones et al. 2017). Our combined analyses constrain the companion's mass to be ${10.88}_{-0.00}^{+7.78}$ M_Jup.
• 78.  
  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. 2015, 2015c; Wittenmyer et al. 2015). Our combined analyses constrain the mass of the bigger companion to be ${6.94}_{-0.52}^{+2.06}$ M_Jup.
• 79.  
  HIP 78395 (WDS 16003-0148) is a K star hosting two companions (Mason et al. 2001). Through 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.06}^{+8.65}$ M_Jup, putting the companion around the boundary between the substellar and stellar categories.
• 80.  
  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}_{-0.32}^{+3.67}$ M_Jup and it is on a nearly edge-on orbit.

Zoom In Zoom Out Reset image size

Table 4. Parameters for the β Pic System Based on Analyses of Various Data Sets in the Literature and in This Work

Reference a	EN20 b	AL20 c	TV20 d	MN20 e	TB21	SL21	This Paper (AL20)	This Paper (TV20)
m⋆ [M⊙]	1.76${}_{-0.02}^{+0.03}$	1.77 ± 0.03	1.80${}_{-0.04}^{+0.03}$	1.82 ± 0.03	1.83 ± 0.04	1.75${}_{-0.02}^{+0.03}$	1.76 ± 0.02	1.80 ± 0.03
mb [MJup]	8.03${}_{-2.62}^{+2.61}$	11.1 ± 0.8	11.7 ± 1.4	9.0 ± 1.6	9.3${}_{-2.5}^{+2.6}$	11.90${}_{-3.04}^{+2.93}$	${7.56}_{-1.69}^{+1.35}$	${11.73}_{-2.14}^{+2.34}$
mc [MJup]	${9.18}_{-0.87}^{+0.96}$	7.8 ± 0.4	8.5 ± 0.5	8.2 ± 0.8	8.3 ± 1.0	8.89${}_{-0.75}^{+0.75}$	${8.94}_{-0.78}^{+0.75}$	${10.14}_{-1.03}^{+1.18}$
Ib [deg]	88.82${}_{-0.01}^{+0.01}$	89.01 ± 0.01	${88.88}_{-0.03}^{+0.04}$	88.99 ± 0.01	88.94 ± 0.02	${88.93}_{-0.01}^{+0.00}$	88.93${}_{-0.09}^{+0.09}$	${89.01}_{-0.01}^{+0.01}$
Ic [deg]	88.85${}_{-0.71}^{+0.72}$	89.01 ± 0.01	⋯	89.17 ± 0.50	89.1 ± 0.66	${88.95}_{-0.10}^{+0.09}$	89.00${}_{-0.01}^{+0.01}$	${88.95}_{-0.09}^{+0.08}$
AL20 RV		✓		✓			✓
TV20 RV			✓		✓	✓		✓
b RV					✓		✓	✓
b Astrometry	✓	✓	✓	✓	✓	✓	✓	✓
c Astrometry		✓		✓	✓	✓	✓	✓
HDR2 f	✓				✓
HEDR3 g							✓	✓

Notes.

^a EN20: Nielsen et al. (2020); AL20: Lagrange et al. (2020); TV20: Vandal et al. (2020); MN20: Nowak et al. (2020); TB21: Brandt et al. (2021a); SL21: Lacour et al. (2021). ^b The solution shown in this table is the so-called "coplanar fit" by EN20, who adopt a Gaussian prior centered on zero and with a standard deviation of 1°. Without such an assumption, the inclination of c has an error of 13°, leading to significant mass uncertainty for c. ^c In AL20, the inclination of β Pic c is assumed to be equal to β Pic b, because the relative astrometry of β Pic c is not used to constrain it. The mass of b reported here is estimated by AL20 using a Gaussian prior of 14 ± 1 M_Jup. Without such an informative prior, the mass is 3.2 M_Jup. ^d In TV20, the stellar mass is given a priori by Wang et al. (2016). ^e The mass of b reported here is estimated by MN20 using a Gaussian prior of 15 ± 3 M_Jup. Without such an informative prior, the mass is 5.6 ± 1.5 M_Jup. ^f HDR2 represents the data of the proper-motion difference between Hipparcos and Gaia DR2. ^g HEDR3 represents the data of both the proper-motion and positional differences between Hipparcos and Gaia EDR3.

Download table as:  ASCIITypeset image

5.1. Mass Distribution and Occurrence Rate

Through our analysis of 5108 stars, with each star having more than five high-precision RV data points, we find 869 stars with 914 LPSs (>1000 days). Of these, 167 of them are confirmed as companions, with masses ranging from 5 M_Jup to 120 M_Jup through 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 estimations (e.g., Grieves et al. 2017; Kiefer et al. 2019).

We define the sample with relative mass uncertainty of less than 100% as the "optimistic sample," and the sample with relative mass error of less than 20% as the "conservative sample." We show the distribution of the sample over mass and mass ratio in Figure 9. There are at least three features seen in the mass distribution of the optimistic sample: (1) there is a lack of BDs with masses 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 the mass ratio (the right panels of Figure 9), we see a valley around 0.3–0.4, but fail to find any significant features around the star–BD boundary (0.07, assuming the host mass to be units of 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 the companion mass. By investigating the distribution over mass (or mass ratio) and semimajor axis (Figure 3), we observe that the 40 M_Jup valley gradually disappears beyond 10 au. Nevertheless, a bias-corrected distribution of substellar companions over mass and semimajor 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.

Zoom In Zoom Out Reset image size

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), while six 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 multicompanion systems. Among the 61 multicompanion systems, there are three systems that contain planets, BDs, and stellar companions, 12 that contain planets and BDs, 21 that contain planets and stellar companions, eight that contain BDs and stellar companions, 12 that contain multiple planets, one that contains multiple BDs, and four that contain multiple stellar companions.

All multicompanion systems are shown in Figure 10. The apparent impression is that the widest companion in a system that is wider than Neptune's orbit tends to be stellar. This is due to the incompleteness of substellar 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 an extremely wide companion on the inner system.

Zoom In Zoom Out Reset image size

There is a controversy over whether hosts of hot Jupiters or BDs tend to have high rates of widely separated companions (Fontanive et al. 2019; Moe & Kratter 2021; Ziegler et al. 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 wide companions on inner planets. By crossmatching the EDR3-based wide-binary sample given by El-Badry et al. (2021) with our sample of massive companions, we find that 29 out of the 161 companion hosts have wide stellar companions, indicating a stellar multiplicity rate of at least 18% ± 3%. Considering that only two host stars in our sample have distances larger than 200 pc, the incompleteness of the identified wide companions is mainly caused by decreasing Gaia completeness below an angular resolution of 2'' (El-Badry et al. 2021).

The detection rate of a wide binary around a companion host is P(DWB, s∣ CH), where s is the binary separation. It is derived from the occurrence rate of the wide binary around the 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 $\tilde{\omega }$ of all the companion hosts and select the ones that could be resolved by Gaia, i.e., $s\tilde{\omega }\gt 2^{\prime\prime}$. The sample becomes significantly incomplete for s < 150 au or $\mathrm{ln}(s/\mathrm{au})\lt 5$. After correcting for this incompleteness, the occurrence rate of wide binaries for companion hosts or P(WB∣CH) is 32% ± 6%, consistent with the binary fraction of 36% ± 2% for binary separation from 20 to 10,000 au (Fontanive et al. 2019). Hence, we do not find any preference of multiplicity for massive companions on wide orbits.

Based on the calculations 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 $\mathrm{ln}(s/\mathrm{au})=5.8$ or s = 330 au. A similar peak around 250 au is also found by Fontanive et al. (2019) for wide binaries hosting hot Jupiters. After considering the detection bias, we find a power-law distribution of P(WB, s∣CH) ∝ s^−1.4 to be optimal for modeling 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). Given 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. (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.

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 separations from their hosts larger than 0.5''. The orbits of 16 of them are shown in Figure 11.

Zoom In Zoom Out Reset image size

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, and 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 ratios between the companions and their host stars are shown in Figure 12 and Table 5. For the modern coronagraphs, such as SCExAO/CHARIS, installed on the Subaru telescope (Jovanovic et al. 2015), the typical inner working angle is around 02 and the contrast limit is about 10⁻⁶ (Currie et al. 2020). Assuming such a detection sensitivity and an age of 1 Gyr for all host stars, 41%, 35%, and 33% of the companions are detectable in the J, H, and K bands. The proportions of the 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 the host stars, respectively, there are 10–57 substellar objects that are detectable by the current imaging facilities.

Zoom In Zoom Out Reset image size

We performed a large number of N-body simulations to study the dynamical stability of the planets themselves, as well as the systems' habitable zones (HZs). The majority of our computations utilize the Mercury 6 hybrid integrator (Chambers 1999), but 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 planet's 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 article), 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 planet's orbit (namely their eccentricities and semimajor axes), we do not investigate the possibility of perturbations from other, undetected massive bodies in each system, which might perturb the orbital paths of the detected bodies. Thus, while our simulations cannot definitively prove that the systems in our sample 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. Multiplanet System Stability

We ran a series of ≥225 1 Myr dynamical simulations to gauge the orbital stability of each multiplanet system in Table 1. In all cases, we consider a grid of five eccentricities and three semimajor axes, which span the range of uncertainties reported in Table 1 for each planet (the 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% of the orbital period of the innermost planet. Through this analysis, we determine that each planet reported in this article is stable in at least 90% of our numerical integrations (the unstable cases typically occur at larger eccentricities). For some multicompanion systems, we performed extended (10 Myr) simulations. As an example, the two detected bodies in the HD 205158 system are plotted in Figure 13. While the BD 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.

Zoom In Zoom Out Reset image size

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, coplanar orbit, situated in the center of the HZ (as determined via the 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. Figure 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 HZs 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.

Zoom In Zoom Out Reset image size

Based on analyses of the proper-motion and positional differences 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, eight hosting both BDs and stellar companions, and three hosting all types of companions.

Without the correction of detection bias, we observe a monotonic decrease of the occurrence rate with mass, from 5 M_Jup to 40 M_Jup, and a valley between 40 M_Jup and 75 M_Jup. The investigation of these features with 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 M_Jup to 120 M_Jup. For a system hosting a 5–120 M_Jup companion, the probability of finding a stellar companion on a wide orbit (>3 au) is about 30% ± 5%, after accounting for detection bias.

Because the cold giants in our sample are nearby and well separated from their hosts, they are good targets for direct imaging with 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 archive, ²⁰ 161 exoplanets and BDs are imaged, 16 of them are older than 1 Gyr, 20 of them have dynamical mass, and only four of them are older than 1 Gyr and have dynamical mass. Because our sample of stars has been selected from various RV surveys, and typically shows stable RV variation, most of the detectable super-Jupiters and BDs are likely to be old. If imaged by follow-up observations, the substellar companions found in this work will extend the imaged cold (or old) giants by an order of magnitude.

We conducted 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 multicompanion systems over 10 Myr. Our simulations also show that the majority of the massive companions identified in this work tend to destabilize the HZs 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 applications of this catalog of substellar companions.

• 1.  
  Substellar mass function. By accounting for the RV and astrometric detection bias, the occurrence rates 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 the BD desert in the parameter space and test various scenarios for the formation of cold giants.
• 2.  
  Misalignment in multicompanion systems. Unlike previous RV surveys, our synergistic survey of nearby substellar companions fully constrains the orbits and masses of the companions. By investigating the orbital misalignment between different companions in a system, and the correlations between misalignment and other orbital parameters, we can discover the causes of this misalignment and constrain evolutionary models of substellar objects.
• 3.  
  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 conatal and coeval formation for the 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.
• 4.  
  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 previous work based on RV samples of exoplanets, with only the minimum mass available, the study of planet correlations based on our sample will avoid the mass ambiguity caused by the limitation of the RV-only method.

In summary, we detect and confirm 167 substellar companions on wide orbits with well-constrained mass, extending the current sample of similar objects by more than one order of magnitude. This catalog is used to study the correlation between substellar companions and wide binaries, to provide dozens of candidates for direct imaging by current facilities, and to investigate the influence of cold giants on the HZs of their hosts. Future works based on our sample would test the BD desert hypothesis and shed light on the formation and evolution of planets and substellar companions.

In the 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 toward the low-mass regime, we will detect and characterize Jupiter-mass and Saturn-mass planets on wide orbits for the study of population synthesis, the dynamical origin of orbital misalignment, and to provide a unique sample of CJs and Saturns for direct imaging by the James Webb Space Telescope (Danielski et al. 2018) and the Chinese Space Station Telescope (Gong et al. 2019).

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 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 the 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.C-0364, 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-0212, 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.C-0732, 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.C-0551, 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.D-0729, 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.

Software: 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 (https://cran.r-project.org/web/packages/fields/index.html), MASS (https://cran.r-project.org/web/packages/MASS/index.html), minpack.lm (https://cran.r-project.org/web/packages/minpack.lm/index.html).

The complete set of fits for all 161 systems is available in Figure A1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (161 panels)

Tarball of figure set images