© E. L. Rickman et al. 2022

1 Introduction

Dynamical masses of companions to main-sequence stars provide a powerful tool for testing mass-luminosity relations and evolutionary models of stellar and substellar objects. They can act as benchmark objects for models that are used to infer masses and radii of isolated objects using their luminosities alone. While there are a number of dynamical masses in the literature for solar mass and higher stars, the known dynamical masses of very low-mass stars (<0.5 M_⊙) with individual mass components are still limited (e.g., Biller et al. 2022). In recent years updated mass-luminosity relationships have been constructed (Chen et al. 2014; Benedict et al. 2016; Mann et al. 2019) that build upon earlier work (e.g., Delfosse et al. 2000). We note the importance of continuing to place observational constraints on these relations in order to constrain them.

A powerful, yet underused, method to determine precise dynamical masses of individual components in stellar binary systems is by combining different detection and measurement techniques that are typical of exoplanet and brown dwarf searches. Radial velocities provide constraints on the minimum mass (m sin i), with an unknown orbital inclination (i), of a companion orbiting a host star, as well as the orbital period. Direct imaging provides the astrometry of a companion relative to a host star. Proper motions from the combination of Gaia (Gaia Collaboration 2016) and HIPPARCOS (ESA 1997) place additional absolute astrometric constraints on the orbit of the system.

Using one of these techniques alone is not sufficient to reveal the dynamical masses of the individual components of a stellar system. However, long-period radial-velocity data, like those from the CORALIE spectrograph (Udry et al. 2000) on the 1.2 m Swiss/EULER telescope and from HARPS (Mayor et al. 2003) on the ESO 3.6 m telescope, both of which are situated at La Silla Observatory, Chile, can be combined with relative astrometry from direct imaging with instruments like VLT/SPHERE (Beuzit et al. 2019) and astrometric accelerations (e.g., Brandt 2021 ; Kervella et al. 2022). This removes the degeneracy of orbital parameters, such as the orbital inclination of the observed system, and reveals model-independent masses. In this way we are able to detect and characterize a number of low-mass companions to bright primary stars, which can serve as benchmark objects to test mass-luminosity relations, as well as as stellar and substellar evolutionary models. It is only recently that the power of combining direct and indirect techniques has been utilized to derive the precise dynamical masses of stellar and substellar objects, as well as exoplanets (e.g., Sahlmann et al. 2011b; Crepp et al. 2012; Dupuy & Liu 2017; Bowler et al. 2018; Cheetham et al. 2018; Brandt et al. 2019, 2021a; Maire et al. 2020; Rickman et al. 2020).

We use radial-velocity data taken from the CORALIE survey as part of the ongoing survey for extra-solar planets (Udry et al. 2000) in the southern hemisphere since 1999. This survey contains a volume-limited sample of solar-type stars within 50 pc. With more than 20 yr of radial-velocity data the CORALIE survey provides a wealth of data for understanding the occurrence rate of planetary, brown dwarf, and stellar companions in the solar neighborhood. From the CORALIE survey, we selected targets that present long-term linear or quadratic drifts in their radial-velocity data, that indicate massive companion candidates amenable to direct imaging observations, and that could potentially be masquerading as exoplanets or brown dwarfs.

We directly imaged selected promising candidates using VLT/SPHERE to reveal their true masses and therefore companion nature. In order to further refine derived orbital parameters and precise dynamical masses, we combined astrometric information from long baseline proper motion anomalies from HIPPARCOS (ESA 1997; van Leeuwen 2007; van Leeuwen & Michalik 2021) and Gaia eDR3 (Gaia Collaboration 2021) from the HIPPARCOS-Gaia catalog of accelerations (HGCA; Brandt 2021) to provide additional constraints on fitted orbital parameters. In order to determine the precise orbital parameters we use the orbital fitting code orvara (Brandt et al. 2021c), which is an efficient code written specifically for the purpose of combining radial-velocity measurements, absolute astrometry, and relative astrometry.

In this paper we report the discovery of new very low-mass stellar companions orbiting HD 142234, HD 143616, and HIP 22059. We also present the first direct detection of HD 92987 B, which was previously detected via radial velocities (Rickman et al. 2019), and update the orbital parameters and dynamical masses of two already known very low-mass stellar companions: HD 157338 and HD 195010, both previously imaged with VLT/NACO by Montagnier (2008). For each of these systems we present precise dynamical mass determinations combining radial velocities, relative astrometry from direct imaging, and HIPPARCOS-Gaia eDR3 accelerations.

The paper is organized as follows. The properties of the host stars are summarized in Sect. 2. In Sect. 3 we present our RV and direct imaging observations and data reduction. In Sect. 4 we discuss the process of determining the astrometry and photometry from the high-contrast images. In Sect. 5 we present the orbit fitting and solution of each of the companions. In Sect. 6 we make a comparison of our derived dynamical masses with low-mass stellar M_H-mass relations. The results are discussed in Sect. 7, and some concluding remarks are presented.

2 Stellar characteristics of the primary stars

For each of the systems, the primary star spectral types are taken from the HIPPARCOS catalog (Perryman et al. 1997), while the V_T band magnitudes and color indices are taken from the Tycho-2 catalog (Høg et al. 2000). The astrometric parallaxes (π) are taken from the Gaia early data release 3 (eDR3; Gaia Collaboration 2021), and the luminosities (L) and effective temperatures (T_eff) are taken from the second Gaia data release (Gaia Collaboration 2018).

The v sin(i) for each primary star is computed using the calibration of the width of CORALIE’s cross-correlation function (CCF; Santos et al. 2001; Marmier 2014). The values for the stellar surface gravities (log g) and metallicities ([Fe/H]) are gathered from literature values across several catalogs (Casagrande et al. 2011; Soto & Jenkins 2018; Sousa et al. 2021).

Isochronal masses, ages, and radii of the primary stars, as well as their uncertainties, are derived using the Geneva stellarevolution models (Ekström et al. 2012; Georgy et al. 2013). The interpolation in the model grid was made through a Bayesian formalism using observational Gaussian priors on T_eff, M_V, and [Fe/H] following the procedure outlined in Marmier (2014). The stellar parameters for each primary star are summarized in Table 1.

3 Observations and data reduction

3.1 Absolute astrometry

Our absolute astrometry comes from the HIPPARCOS-Gaia Catalog of Accelerations (HGCA, Brandt 2018, 2021). The HGCA is a cross-calibration of HIPPARCOS (ESA 1997; van Leeuwen 2007) and Gaia EDR3 (Gaia Collaboration 2016, 2021; Lindegren et al. 2021) that places both on a common reference frame with calibrated uncertainties. Each star in the HGCA has three proper motions: a Hipparcos proper motion near 1991.25, a Gaia proper motion near 2016.0, and the position difference between HIPPARCOS and Gaia scaled by the time baseline between the missions.

Table 2 lists each of our target star’s absolute astrometry from the HGCA. The listed Gaia proper motions match the EDR3 catalog values, but the epochs listed are those that minimize the uncertainty in each star’s position. The epochs differ from 2016.0 because of covariance between position and proper motion at the catalog epoch adopted by Gaia. The Hipparcos proper motions are a linear combination of the ESA (ESA 1997) and HIPPARCOS-2 (van Leeuwen 2007) reductions, corrected for a locally variable frame rotation as discussed in Brandt (2018). All uncertainties are inflated from the catalog values to achieve statistical agreement between the three proper motions for a sample of stars that show no radial-velocity acceleration (Brandt 2021).

All stars listed in Table 2 show highly significant discrepancies between the proper motion measurements, indicative of a tug by an unseen companion. The formal significance of this astrometric acceleration ranges from almost 50σ for HD 142234 to nearly 300σ for HD 92987.

3.2 Radial velocities

The radial-velocity data are taken from the CORALIE Survey for Extra-Solar Planets (Queloz et al. 2000), a planet search survey in progress for more than 20 yr in the southern hemisphere, whose observations began in June 1998. The survey makes use of the CORALIE spectrograph at La Silla Observatory, and monitors a volume-limited sample of 1647 main-sequence stars spanning from spectral types F8 to K0 within 50 pc of the Sun.

In order to increase the efficiency and accuracy of the CORALIE spectrograph over the years, it underwent two major upgrades in June 2007 (Ségransan et al. 2010) and in November 2014. While these changes improved the overall performance of the CORALIE spectrograph, they also introduced small offsets in the measured radial velocities; therefore, we treat the radial-velocity data from the CORALIE spectrograph as three separate instruments corresponding to these upgrades. We refer to the original CORALIE as CORALIE-98 (C98), the upgrade in 2007 as CORALIE-07 (C07), and the latest upgrade in 2014 as CORALIE-14 (C14).

We also include some additional data for a few of the targets from other spectrographs, namely HARPS (Mayor et al. 2003) and HIRES (Vogt et al. 1994). All of the data products presented in this paper are available at the Data and Analysis Center for Exoplanets (DACE)¹.

The RV data were reduced using the CORALIE automated pipeline, which also measures the cross-correlation function (CCF), full width at half maximum (FWHM), the bisector, and the H_α chromospheric activity indicator. These indicators allow us to pinpoint the origin of observed periodic signals, to ensure that we are not observing, for example, significant stellar activity.

As part of the ongoing CORALIE survey, we have been monitoring for long-term trends in the data that could be indicative of widely separated companions that are amenable to high-contrast imaging observations. To monitor these candidates, we assessed our sample of RV data for long-period linear or quadratic trends over the past >20-yr baseline.

Fig. 1 Contours showing the astrometric predictions of each companion relative to the host star in ascension (Δα) and declination (Δδ) using the orvara orbit fitting package (Brandt et al. 2021c), with the measured companion positions shown by the overplotted gray stars. Relative astrometric positions indicate whether a companion can be directly detected with imaging given its projected separation between the host star and the companion itself. The 1, 2, and 3σ contour levels of the predicted positions of each companion are plotted from a full orbit fit using just the radial velocities and the HGCA astrometry, with the astrometric predictions being shown for the year in which each companion was then subsequently directly imaged with VLT/SPHERE in order to visually compare the prior predictions to the direct detections. The measured position of each companion from the VLT/SPHERE imaging are shown by the overplotted gray stars. In the case of HD 142234, the left star corresponds to HD 142234 B and the right star corresponds to HD 142234 C. The measured positions of each companion from direct imaging observations agree well with the predicted positions from the orbit fitting using just the radial velocities and astrometry from the HGCA.

3.3 High-contrast imaging

The radial-velocity data taken from the CORALIE spectrograph (Queloz et al. 2000), as well as HARPS (Mayor et al. 2003) and HIRES (Vogt et al. 1994), are used as precursor observations to assess the detectability of any companion candidates through direct imaging. We combined the radial-velocity data with the astrometric information from the HGCA to determine the predicted relative astrometry in ascension and declination (Δα, Δδ) of each companion on the sky, using the orbit-fitting code orvara (Brandt et al. 2021c). The results of the predicted positions, which show full orbital fits of each system using the radial velocities and the astrometric accelerations from the HGCA, are shown in Fig. 1. The detectability of these objects was assessed by translating these values into a predicted contrast ratio and a projected angular separation for each system. These values were then compared against the on-sky measured performance 5σ contrast curves from VLT/SPHERE.

We selected to observe the most promising targets using the VLT/SPHERE dual-band imaging mode (Vigan et al. 2010) using the InfraRed Dual-Band Imager and Spectrograph (IRDIS; Dohlen et al. 2008) in the H2 and H3 bands (λ_H2 = 1.593 µm, λ_H3 = 1.667 µm). Observations were taken with SPHERE (Beuzit et al. 2019) via the extreme adaptive optics system at the VLT, under programs 0102.C-0236(A) (PI: Rickman), 0104.C-0724(A) (PI: Rickman), and 105.20SZ.001 (PI: Rickman).

The observations were taken using an apodized Lyot coronagraph (Soummer et al. 2003), where the length of exposure time was calculated to be optimal for the expected contrast ratio of the detections, with observations ranging from 0.5 to 3 h. In order to measure the position of the centre of the star behind the coronagraph, several exposures at the beginning and at the end of the observing sequence were taken with two orthogonal sinusoidal modulations applied to the deformable mirror, which generates satellite spots around the star.

The observation sequence also consisted of short exposure images of the primary star moved from behind the coronagraph to obtain the flux and the shape of the point spread function (PSF), making use of a neutral density (ND) filter². Longexposure sky frames were also taken to estimate the background flux and correct for any bad pixels on the detector. The observations were carried out in pupil tracking mode, which allows the observed field of view to rotate during the course of the coronagraphic observations, in order to obtain high-contrast images using the angular differential imaging (ADI) approach.

The VLT/SPHERE data were reduced using the Geneva Reduction and Analysis Pipeline for High-contrast Imaging of planetary Companions (GRAPHIC, Hagelberg et al. 2016). The GRAPHIC pipeline carries out sky subtraction, flat-fielding, bad pixel correction, frame selection, and a correction for the distortion as outlined in Maire et al. (2016). We then performed a classical ADI (cADI) PSF subtraction algorithm before derotating and median combining to produce a final PSF-subtracted image. The resulting images in the H2 band are shown in Fig. 2

Fig. 2 High-contrast images taken with VLT/SPHERE using IRDIS coronagraphic imaging in the H2 band. In each image the star is masked behind the central white circle. Some of the coronagraphic images shown here show some saturation, which is why the flux frames were used to obtain the photometry of each companion. The image for HD 142234 is shown on a wider angular scale than the other images.

4 Astrometry and photometry

Due to saturation in some of the coronagraphic images caused by the brightness of the companions, we opted to use the flux frames of each data set to determine the astrometry and photometry of the companions. These frames use a neutral density filter with no coronagraph and a shorter exposure time than a norma science sequence. Our use of unsaturated flux frames enabled us to perform PSF photometry and astrometry, using the primary star itself as the PSF template.

We first bias-correct and flat-field the images, interpolating over bad pixels. We then perform reference differential imaging (RDI), subtracting a scaled IRDIS stellar PSF taken under the same observing conditions to remove the primary star’s light. For this purpose we use HD 142234 (which has no companions within 1″) and two additional comparison stars without companions from the ESO archive: HIP 82154 and HIP 102979. HD 142234 itself has companions sufficiently far away that RD] is unnecessary. For all other stars, we have six image pairs: three PSF template stars, and two sets of SPHERE/IRDIS coadds for each epoch.

We optimize the position and relative flux of each PSF reference star within one of our flux frames and use it to remove the light of the primary star. We then perform PSF photometry on the residual, which contains the light from the companion, together with noise and low-level residuals from the primary. We use the unsubtracted PSF of the primary as the template for the secondary. In this way, we fit for an offset between the two stars without needing to centroid the image since the stellar PSF and companion PSF are both present in the same image. We also obtain a contrast without having to assume anything about the form of the PSF.

Our approach resulted in six measurements of the contrast and position (three PSF reference images times two sets of coadds) for each epoch. We report the mean and standard deviation of these six epochs as our measurements and uncertainties; we bias-corrected the standard deviations by a factor of . We conservatively did not divide the uncertainties by as we did not have a large enough data set to verify the statistical independence of our measurements.

The separation and position angle detector positions were converted into on-sky separation and position angles by accounting for an anamorphic distortion of 1.0062 ± 0.0002, plate scale of 12.255 ± 0.009 mas pix⁻¹ for the H2 band observations, and 12.250 ± 0.009 mas pix⁻¹ for the H3 band observations (Maire et al. 2016). A true north offset of −1.75 ± 0.08 deg (Maire et al. 2016), and a pupil offset of 135.99 ± 0.11 deg were also applied³.

The resulting astrometry and photometry of the very low-mass stellar companions are given in Table 3. We note that the error on the relative astrometry and photometry of HD 92987 is larger than the other systems. This is due to the HD 92987 system having a higher contrast ratio of the companion to the primary star than the other observed systems, making it more difficult to obtain the same level of precision for the relative astrometry and photometry. Despite this, we still measure values with reasonable uncertainty that place orbital constraints on the HD 92987 system.

In Fig. 1 we compare the predicted and measured relative astrometry of each companion by overplotting the measured astrometry outlined in Table 3 on the predicted positions. From this we see that the positions of each observed companion are in good agreement with each other, where each companion is within at least 2σ of the prediction.

The values from Table 3 are shown in Fig. 3 as a comparison against the measured on-sky performance 5σ contrast curves from H2 VLT/SPHERE observations to assess their direct imaging detectability. The observations were obtained through injection and recovery of PSFs of companions at increasing separation for a number of targeted stars from the CORALIE survey.

Fig. 3 Measured contrast of each companion at the measured separations with stacked 5σ contrast curves corresponding to measured on-sky performance from VLT/SPHERE data for IRDIS H2 bands. The detection probability scale is shown on the right. The darker color corresponds to a lower likelihood of detectability with imaging with VLT/SPHERE IRDIS H2 band, and the lighter colors correspond to a higher likelihood of detectability at a given separation and contrast ratio. The astrometry and photometry for HD 157338 and HD 195010 are each shown as an average of the two 2021 SPHERE H2 band data points.

5 Orbital solutions

In order to obtain orbital fits of our observed systems, we made use of orvara (Brandt et al. 2021c), which is an orbit-fitting code specifically designed to combine radial-velocity data, relative astrometry from direct imaging, and absolute astrometric data from the HIPPARCOS and Gaia. The orvara code (Brandt et al. 2021c) utilizes a comprehensive Markov chain Monte Carlo (MCMC) approach, incorporating astrometric accelerations from HIPPARCOS-Gaia eDR3 using the HIPPARCOS-Gaia catalog of accelerations (HGCA; Brandt 2021) and the intermediate astrometry fitting tool htof (Brandt et al. 2021b).

Using orvara, we were able to obtain robust orbital fits by combining the radial-velocity data with relative astrometry obtained from our VLT/SPHERE observations, as described in Sect. 4, as well as proper motion anomalies from the long baseline between the HIPPARCOS and Gaia missions. In this section we present the orbital solutions for each of the systems. The full list of posteriors of the fitted orbital parameters is given in Table 4, and each of the corner plots showing the posterior distributions is shown in Appendix B.

5.1 HD 92987 (HIP 52472)

HD 92987 is a G2/3V star that has been observed with CORALIE at La Silla Observatory since January 1999. Fifty-three measurements were taken with CORALIE-98, 18 additional RV measurements were obtained with CORALIE-07, and 57 additional RV measurements were obtained with CORALIE-14. The detection of HD 92987 B was originally presented in Rickman et al. (2019) with a minimum mass derived from the radial velocities of M_Jup, interpreted as being either a massive planetary companion, or a very low-mass brown dwarf. The radial-velocity detection of HD 92987 B was also reported independently by Kane et al. (2019), who measured a minimum mass of M_Jup, in agreement with the interpretation presented in Rickman et al. (2019).

HD 92987 was then observed with the VLT/SPHERE program 0102.C-0236(A) (PI: Rickman) on 2019 January 17, with a total observation time of three hours, with a 64-second integration time on the coronagraphic frames using IRDIS, as well as several sky, center, and flux frames taken which was necessary for the post-processing of the data. The direct detection of HD 92987 B is shown in Fig. 2. Since Rickman et al. (2019), we increased the radial-velocity baseline by about three years, continuing to monitor how the radial-velocity trend changes over time, which enabled us to further constrain its orbital parameters.

Combining the radial-velocity measurements with the relative astrometry in Table 3 and with the HIPPARCOS-Gaia eDR3 accelerations, we ran the MCMC orbital fit using orvara (Brandt et al. 2021c) with a log-flat prior on the primary mass and with 500 000 steps in each chain. The resulting plots are shown in Fig. 4 with a burn-in phase of 700 multiplied by every 250th step that is saved on the chain with 500 randomly drawn orbits.

We derive a companion mass of M_Jup, with an orbital period of yr. For the primary mass of HD 92987 A we derive a mass of M_⊙, which is in agreement with the mass we determined using the Geneva stellar-evolution models (Ekström et al. 2012; Georgy et al. 2013), as described in Sect. 2. Our dynamical mass value of the primary star, HD 92987 A, is also in agreement with a recent value from Sousa et al. (2021), who derived a spectroscopic mass of 1.155 ± 0.005 M_⊙.

Kane et al. (2019) presented imaging observations of HD 92987 using the Differential Speckle Survey Instrument (DSSI) operating on the Gemini South telescope, but did not directly detect HD 92987 B. They did, however, use these observations to constrain the upper mass limit of HD 92987 B as a function of the angular separation, and found the upper mass limit to be in the 200–500 M_Jup mass range at approximately 200 mas separation. This is in agreement with the mass derived using our SPHERE observations and absolute astrometry from HIPPARCOS and Gaia, with a mass of M_Jup at 200.81 ± 5.99 mas in the H2 band at 1.593 µm.

More recently, the mass of HD 92987 B was determined by Venner et al. (2021) through combining published radial velocities with astrometry from HIPPARCOS and Gaia eDR3, who derived a mass of 268.4 ± 4.7 M_Jup. For their orbital fit they placed a Gaussian prior on the primary star based on their model mass of the primary of 1.043 ± 0.012 M_⊙. The companion mass found in Venner et al. (2021) has a precision on their mass estimate of 1.2%, but relies on a modeled mass of the primary star, unlike the primary mass we derive in this work. The mass of HD 92987 B determined by Venner et al. (2021) of 268.4 ± 4.7 M_Jup is in good agreement with our determined mass of M_Jup, with an expected increase on the error bar as we used an uninformed log-flat prior on the primary mass when doing the orbital fit, leading to a completely model independent derived mass of HD 92987 A and HD 92987 B.

A precise dynamical mass can be well constrained for HD 92987 B thanks to the well-sampled orbital phase coverage from the CORALIE survey, as shown in Fig. 4. The relative astrometric point from the VLT/SPHERE observations is able to underpin the tight fit of the relative astrometric orbit shown in Fig. 4, which could be further constrained with additional direct imaging observations. The fitted orbital parameters are listed in Table 4. The resulting corner plots from the probed parameter space from the MCMC fit with orvara and the posterior distributions are shown in Fig. B.1.

The HD 92987 system has an almost face-on orbit with an inclination of 175.76 ± 0.08 degrees. This emphasizes the importance of not taking the mass derived from radial-velocities alone (which only gives the minimum mass, m sin i) as the true mass. The companion nature of such detections are undetermined and could be a stellar or substellar companion masquerading as an exoplanetary candidate. This is the case for HD 92987 B, which joins an increasing list of substellar companions discovered using radial velocities alone that are in fact stellar in nature (e.g., Bean et al. 2007; Sozzetti & Desidera 2010; Sahlmann et al. 2011a; Wilson et al. 2016; Moutou et al. 2017; Benedict & Harrison 2017; Kiefer et al. 2019).

5.2 HD 142234 (HIP 77931)

HD 142234 is a G5V star that has been observed with CORALIE at La Silla Observatory since January 2003. Twenty-six measurements were taken with CORALIE-98, an additional 14 measurements were taken with CORALIE-07, and 7 measurements obtained with CORALIE-14.

HD 142234 was then observed with VLT/SPHERE program 105.20SZ.001 (PI: Rickman) on 2021 May 16, with a total observation time of 1 h 50 min, with a 64-s integration time for the coronagraphic frames using IRDIS, as well as several sky, center, and flux frames taken which was necessary for the postprocessing of the data. The imaging data from the VlT/SPHERE observations are shown in Fig. 2.

From the VLT/SPHERE images we found HD 142234 to be a triple system. Given the linear trend of the radial velocities (top left plot in Fig. 5) and the closeness of the B and C components, we conclude that this is a hierarchical triple system with BC in a tight orbit around each other, and a longer orbit around HD 142234 A. As orvara does not yet have the capability to perform orbital fits of this configuration, we performed an orbital fit with the photocenter of the relative astrometry between the two components as the barycenter to estimate the combined mass of HD 142234 B+C as if it were one point source, and then used the photometric ratio of the two resolved point sources to estimate the ratio of the masses of the individual components.

Only a small fraction of the orbital phase is sampled by the radial-velocity data, and so we implemented a Gaussian prior on the primary mass derived using the Geneva stellar-evolution models (Ekström et al. 2012; Georgy et al. 2013), as described in Sect. 2. We used a value of 0.83 ± 0.02 M_⊙ for the primary (see Table 1). We ran orvara using 500000 steps in each chain with the CORALIE radial-velocity data, the relative astrometry from VLT/SPHERE taken to be a single point as the photocenter between HD 142234 B and HD 142234 C, and the HIPPARCOS-Gaia astrometry. The resulting orbital fit is shown in Fig. 5, with a burn-in phase of 700 multiplied by every 250th step that is saved on the chain with 500 randomly drawn orbits.

We note that the relative astrometric orbital fit is not as tightly constrained as the other systems in this paper. This is primarily due to the unknown true barycenter of HD 142234 B and HD 142234 C, leaving a large uncertainty on the relative astrometric position. We take the photocenter to be the point that is equidistant between HD 142234 B and HD 142234 C, giving the astrometry to be separation ρ = 1.117 ± 0.003″ and position angle θ = 258.13 ± 0.23°. Furthermore, we only sample a relatively uninformative portion of its RV orbit leading to an uncertainty on the orbital period. The orbital constraints provide a good first estimate of the companion masses, and in future versions of orvara with the updated capability to incorporate this orbital configuration, an updated orbital fit and comparison can be carried out.

From this orbital fit, we calculate a combined mass of HD 142334 B+C to be M_Jup with an orbital period of yr. Using the photometric ratio of HD 142234 B to HD 144234 C from the calculated contrasts shown in Table 3, we estimate the masses of each component from the combined orbital fit. In this way we obtain an estimate of HD 142234 B and HD 142234 C to be M_Jup and M_Jup, respectively.

The fitted orbital parameters are listed in Table 4. The resulting corner plots from the probed parameter space from the MCMC fit with orvara and the posterior distributions are shown in Fig. B.2.

As shown in Fig. B.2, the eccentricity of HD 142234 B+C is not very well constrained, with the posterior distribution spanning from the lower to the upper bounds, with a skew toward a higher eccentricity, and a resulting fit of . The large uncertainty is expected in this case, as the orbital phase coverage is not completely covered by the radial velocities. As seen in the top left panel of Fig. 5, the radial velocities mostly reveal a long-term linear trend, making it difficult to tightly constrain the orbital parameters. This is also due to the orbital fit of HD 142234 B and HD 142234 C fitted as a single point source; there are large uncertainties in the true relative astrometry of the system. However, the orbital fit does provide good first constraints on this system.

Fig. 4 Orbit fits of HD 92987 using the orbit-fitting code orvara. Top left: radial-velocity orbit induced by HD 92987 B over a full orbital period. Shown are the radial-velocity data of COR-98 (blue points), COR-07 (yellow points), and COR-14 (green points). The thick line shows the highest likelihood fit; the thin colored lines show 500 orbits drawn randomly from the posterior distribution. Top right: relative astrometric orbit of HD 92987 B. The thick black line represents the highest likelihood orbit; the thin colored lines represent 500 orbits drawn randomly from the posterior distribution. Dark purple corresponds to a low companion mass and light yellow corresponds to a high companion mass. The dotted black line shows the periastron passage, and the arrow at the periastron passage shows the direction of the orbit. The dashed line indicates the line of nodes. Predicted past and future relative astrometric points are shown by black circles with their respective years, while the observed relative astrometric point from VLT/SPHERE data is shown by the blue-filled data point, where the measurement error is smaller than the plotted symbol. Bottom: acceleration induced by the companion on the host star as measured from absolute astrometry from Hipparcos and Gaia. The thick black line represents the highest likelihood orbit; the thin colored lines are 500 orbits drawn randomly from the posterior distribution. Darker purple represents a lower companion mass and light yellow represents a higher companion mass. The residuals of the proper motions are shown in the bottom panels.

Fig. 5 Same as Fig. 4, but for HD 142234 B+C. Here the orbital fit is performed with HD 142234 B+C together.

Fig. 6 Same as Fig. 4, but for HD 143616 B.

5.3 HD 143616 (HIP 78551)

HD 143616 is a G6/8V star that has been observed with CORALIE at La Silla Observatory since April 2003. Twelve measurements were taken with CORALIE-98, an additional 12 measurements were taken with CORALIE-07, and 48 measurements were taken with CORALIE-14.

HD 143616 was then observed with VLT/SPHERE program 105.20SZ.001 (PI: Rickman) on 2021 July 17, with a total observation time of 40 min, with a 64-s integration time on the coronagraphic frames using IRDIS, as well as several sky, center, and flux frames taken which was necessary for the postprocessing of the data. The direct detection of HD 143616 B is shown in Fig. 2.

From the radial velocities alone, we already see that the orbit is extremely eccentric (see Fig. 6). Due to its extreme eccentricity, we implemented a starting position for the orbital fit, which we determined from a fit on the radial-velocity data alone.

To perform an initial fit on just the radial-velocity data, we used the Data and Analysis Center for Exoplanets (DACE; Buchschacher & Alesina 2019)⁴. The DACE online tool fits a Keplerian model using the formalism described in Delisle et al. (2016). The modeling of instrumental noise and stellar jitter is described in Díaz et al. (2016) and Delisle et al. (2018). We used the parameters derived for the semimajor axis (a), sin ω, and cos ω as the start points for the orbital fit with orvara, which can be implemented using a start file. The start points we implemented from the DACE radial-velocity fit are 9.0 ± 0.143 AU for the semimajor axis a; −0.036460 ± 0.000330 for sin ω; and −0.9790787 ± 0.0000866 for cos ω.

The choice of these start points was motivated by the wellsampled periastron passage with the radial-velocity data and the extremity of the eccentricity of the system. Using these start points we fit the orbit combining all of the CORALIE radial-velocity data, with the relative astrometry from the VLT/SPHERE imaging, as well as the proper motion anomalies from the HGCA. We ran orvara with 500000 steps in each chain, and for the orbital plots shown in Fig. 6 we defined a burn-in phase of 700 multiplied by every 250th step that is saved on the chain, with the best fit orbit shown, with 500 orbits randomly selected from the posterior distribution plotted alongside the RVs, and proper motions from Hipparcos and Gaia. From the combined orbital fit, we derive a dynamical mass of HD 143616 B of M_Jup, an orbital period of yr, and an extreme eccentricity of .

We placed an uninformed log-flat prior on the primary mass when doing the orbital fit, leading to a completely model-independent derived mass of HD 143616 A and HD 143616 B. From the orbital fit, we find the dynamical mass of the primary star to be 0.935 ± 0.017 M_⊙, which is in good agreement with the mass derived from the Geneva stellar-evolution models (Ekström et al. 2012; Georgy et al. 2013), as described in Sect. 2, where we find a value of 0.92 ± 0.02 M_⊙, as listed in Table 1.

The fitted orbital parameters are listed in Table 4. The resulting corner plots from the probed parameter space from the MCMC fit with orvara and the posterior distributions are shown in Fig. B.3 in Appendix B.

Fig. 7 Same as Fig. 4, but for HD 157338 B. Top left: radial-velocity orbit induced by HD 157338 B. The white data points show the HIRES data, the blue data points show the COR-98 data, the cyan data points show the HARPS data, the yellow points show the COR-07 data, and the green points show the COR-14 data. Top right: predicted past and future relative astrometric points are shown as black circles with their respective years, while the observed relative astrometric points from the 2006 VLT/NACO data and the 2021 VLT/SPHERE data are shown as blue filled data points, where the measurement error is smaller than the plotted symbol.

5.4 HD 157338 (HIP 85158)

HD 157338 is a G0/1V star that has been observed with CORALIE at La Silla Observatory since April 1999. Thirty-one measurements were taken with CORALIE-98, an additional 26 were taken with CORALIE-07, and 53 points with HARPS (Mayor et al. 2003). We also include 33 radial-velocity measurements from HIRES (Butler et al. 2017) that date back to April 1997, increasing the time span of the CORALIE survey sampling by an additional two years.

HD 157338 was first directly imaged by Montagnier (2008) with VLT/NACO (Lenzen et al. 2003; Rousset et al. 2003) in 2006 across the F1 (λ = 1.575 µm), F2 (λ = 1.600 µm), and F3 (λ = 1.625 bands. From those observations, Montagnier (2008) reported the detection of HD 157338 B. We used these archival NACO data in our orbital fit, as shown in Figs. 7 and 8, which helped further refine our derived orbital parameters and dynamical mass.

HD 157338 was then observed with VLT/SPHERE program 105.20SZ.001 (PI: Rickman) on 2021 June 03, with a total observation time of 30 min, with a 32-s integration time on the coronagraphic frames using IRDIS, as well as several sky, center, and flux frames taken which was necessary for the postprocessing of the data. A follow-up observation of HD 157338 was taken on 2021 June 30, providing an additional relative astrometric point for our orbital fit. The direct detection of HD 157338 B is shown in Fig. 2.

We ran orvara for the orbital fit with 500 000 steps in each chain, using the radial-velocity data from CORALIE, HARPS, and HIRES, along with relative astrometry from VLT/NACO and VLT/SPHERE, and astrometric accelerations from the HGCA. The resulting orbital fits are shown in Fig. 7, which are plotted with a burn-in phase of 700 multiplied by every 250th step that is saved on the chain with 500 orbits randomly selected from the posterior distribution plotted alongside the RVs, and proper motions from Hipparcos and Gaia. From the orbital fit we derive a companion mass of 359.3 ± 10.5 M_Jup, and an orbital period of yr

The long baseline between direct imaging observations means, despite only having a fraction of the orbital phase covered with radial velocities, that the orbit is well-constrained. The significance of the additional imaging data with VLT/SPHERE is shown in Fig. 8, where the long baseline between direct imaging observations and the gain in the precision of the relative astrometry contribute greatly to constraining the orbital parameters.

Montagnier (2008) derived a mass of 0.405 ± 0.03 M_⊙ (424.3 ± 31.4 M_Jup) for HD 157338 B. The mass we derive using orvara is less massive at 359.3 ± 10.5 M_Jup, but as seen in Fig. 8 and from an additional 16 yr of radial-velocity data, and proper motion data from the HGCA, we are able to determine a much more precise updated dynamical mass value.

We placed an uninformed log-flat prior on the primary mass when doing the orbital fit, leading to a completely model independent derived mass of HD 157338 A and HD 157338 B. From the orbital fit, we find the dynamical mass of the primary star to be M_⊙, which is in agreement with the primary mass derived using the Geneva stellar-evolution models (Ekström et al. 2012; Georgy et al. 2013), as described in Sect. 2 and listed in Table 1.

The greatest uncertainty in the orbital fit comes from the lack of full orbital phase coverage from the radial-velocity sampling due to the long orbital period (>100 yr). However, because the system will approach its periastron passage in late 2024, the change in the orbital motion will be at its greatest, which is an ideal time to gain additional radial velocity or direct imaging astrometry. This would place even tighter constraints on an already relatively well-constrained orbit.

The fitted orbital parameters are listed in Table 4. The resulting corner plots from the probed parameter space from the MCMC fit with orvara and the posterior distributions are shown in Fig. B.4.

Fig. 8 Relative separation and position angle of HD 157338 B across three epochs of data. The 2006 data corresponds to VLT/NACO data as reported by Montagnier (2008), and the two 2021 data points correspond to the VLT/SPHERE data presented in this paper. The second epoch of direct imaging observations over a long baseline is able to constrain the orbital parameters significantly, as shown in the residual plots in the corresponding bottom panels.

5.5 HD 195010 (HIP 101145)

HD 195010 is a G8/K0V star that has been observed with CORALIE at La Silla Observatory since September 2000. Forty-two measurements were taken with CORALIE-98, 15 were taken with CORALIE-07, and 51 with CORALIE-14.

HD 195010 was directly imaged by Montagnier (2008) with VLT/NACO (Lenzen et al. 2003; Rousset et al. 2003) in 2006 across the F1 (λ = 1.575 µm), F2 (λ = 1.600 µm), and F3 (λ = 1.625 µm) bands. From those observations, Montagnier (2008) reported the detection of HD 195010 B. We used these archival NACO data in our orbital fit, as shown in Figs. 9 and 10, to help further refine our derived orbital parameters and dynamical mass.

We observed HD 195010 with VLT/SPHERE program 105.20SZ.001 (PI: Rickman) on 2021 June 03 and 2021 June 30, with an observation time of 30 minfor each epoch, with a 32-sintegration time on the coronagraphic frames using IRDIS, as well as several sky, center, and flux frames taken which was necessary for the post-processing of the data. The direct detection of HD 195010 B from the VLT/SPHERE observations is shown in Fig. 2.

We ran the orbit-fitting code orvara with 500 000 steps in each chain using the radial-velocity data from the CORALIE survey, the relative astrometry from high-contrast imaging with NACO in 2006, as well as the relative astrometry from high-contrast imaging with VLT/SPHERE in 2021, and the astrometric accelerations from the Hipparcos-Gaia catalog of accelerations (HGCA; Brandt 2021).

The best orbital fits, with 500 orbits randomly selected from the posterior distribution and a burn-in phase of 700 multiplied by every 250th step that is saved on the chain, are plotted alongside the RVs, and proper motions from Hipparcos and Gaia as shown in Fig. 9. From this we derived a companion mass of M_Jup and an orbital period of yr. We compare this to the value obtained in Montagnier (2008), who find a mass from combining the VLT/NACO imaging data with CORALIE radial velocities of 0.26 ± 0.03 M_⊙ or 272.4 ± 31.4 M_Jup.

The dynamical mass that we determine for HD 195010 B is lower, but we now have far more significant constraints than the previously reported mass from Montagnier (2008). This is due to an increase in the baseline of CORALIE radial velocities of 16 yr, which is a significant increase in orbital phase coverage with the radial-velocity data, especially throughout its periastron passage, which is key in fitting such a well-constrained orbit. Furthermore, our new data increases the baseline of the relative astrometry, with constraints from both the VLT/NACO and VLT/SPHERE data, as shown in Fig. 10. In addition, the proper motion data from the HGCA helps place additional constraints on the orbital parameters, as shown in Fig. 9.

We placed a log-flat prior on the primary mass when doing the orbital fit, leading to a completely model independent derived mass of HD 195010 A and HD 195010 B. Much like in the case for HD 92987 B, the system HD 195010 has a highly inclined orbit that is close to face-on, with an inclination of i = 162.79 ± 0.43 degrees. This once again demonstrates the importance of obtaining dynamical masses by combining these mass measuring techniques to confirm the companion type of detected objects, which could indeed be stellar objects masquerading as brown dwarfs or exoplanetary candidates. From fitting the radial velocities alone, with no additional astrometric information, we obtain the minimum mass (m sin i) of HD 195010 B to be M_B ≈ 51 M_Jup. Without the additional relative astrometry or proper motion information, this would be interpreted as a brown dwarf, leading to a different understanding of this system and emphasizing the importance of deriving these dynamical masses.

The fitted orbital parameters of HD 195010 B are listed in Table 4. The resulting corner plots from the probed parameter space from the MCMC fit with orvara and the posterior distributions are shown in Fig. B.5 in Appendix B.

Fig. 9 Orbit fits of HD 195010 using the orbit-fiting code orvara. Top left: same as Fig. 4, but for HD 195010 B. Top right: predicted past and future relative astrometric points are shown as black circles with their respective years, while the observed relative astrometric points from the 2006 VLT/NACO and the 2021 VLT/SPHERE data is shown as blue filled data points, where the measurement error is smaller than the plotted symbol.

Fig. 10 Relative separation and position angle of HD 195010 B across three epochs of data. The 2006 data corresponds to the VLT/NACO data reported by Montagnier (2008), and the 2021 data point correspond to the VLT/SPHERE data presented in this paper. The residuals are shown in the bottom two panels.

5.6 HIP 22059

HIP 22059 is a K5V star that has been observed with CORALIE at La Silla Observatory since January 2001. Twenty-nine measurements were taken with CORALIE-98, 9 additional RV measurements were obtained with CORALIE-07, and 30 additional RV measurements were obtained with CORALIE-14. HIP 22059 was also observed with HARPS (Mayor et al. 2003) with 11 RV measurements from December 2003 to August 2007.

HIP 22059 was then observed with VLT/SPHERE program 0104.C-0724(A) (PI: Rickman) on 2019 November 27, with a total observation time of 1 h 30 min, with a 32-s integration time on the coronagraphic frames using IRDIS, as well as several sky, center, and flux frames taken which was necessary for the postprocessing of the data. The direct detection of HIP 22059 B is shown in Fig. 2.

The orbital fit was run using orvara with 500 000 steps in each chain by combining CORALIE and HARPS radial velocities, with the relative astrometry outlined in Table 3, and epoch astrometry from Hipparcos and Gaia eDR3. We imposed an uninformative log-flat prior on the primary and secondary mass and uniform priors on the remaining orbital parameters, meaning that the tight constraints on the primary mass come from the orbital fit alone, and leads to precise dynamical masses of the two components in the system. From this we determine the mass of the primary star to be M_⊙.

For the plots shown in Fig. 11, we implemented a burn-in phase of 700 multiplied by every 250th step that is saved on the chain, and randomly drew 500 orbits showing the posterior distribution. The highest likelihood fit is shown by the thick black line for the RVs, proper motions from Hipparcos and Gaia, as well as the relative astrometric orbit. From the orbital fit we find a very low-mass stellar companion with a companion mass of M_Jup, with an orbital period of yr.

The largest uncertainty in the mass is from only having partial orbital phase coverage from the radial-velocity data, as shown in the top left panel of Fig. 11. Full phase coverage would take more than 90 yr to collect, which is more than a scientific lifetime. However, given the dispersion of the uncertainty in the orbit in just a few years from now, an additional relative astrometric point in five to ten years would tightly constrain the orbital parameters of this system.

The fitted orbital parameters are listed in Table 4. The resulting corner plots from the probed parameter space from the MCMC fit with orvara and the posterior distributions are shown in Fig. B.6.

6 Dynamical vs. isochronal masses

We compare the companions’ dynamical masses to their luminosity masses inferred from stellar isochrones. We extracted the narrow-band H2 and H3 absolute magnitudes for each companion using the BT-NextGen models (Allard et al. 2012), and applying our measured contrast values listed in Table 3 and Gaia eDR3 parallaxes listed in Table 1. We then convert the extracted H2 and H3 magnitudes into a H_2MASS magnitude for each companion, making use of the 2MASS filter system (Skrutskie et al. 2006), to be able to directly compare with the mass-luminosity relation in the H band shown in Fig. 12. The results are summarized in Table 5.

In Fig. 12, we compare the masses and magnitudes of the companions to the mass-luminosity relation near the low-mass end of PARSEC isochrones (Bressan et al. 2012; Chen et al. 2014). We adopted three different isochrones at low (Z = 0.05, dotted orange line), intermediate (Z = 0.15, solid red line), and high (Z = 0.35, dashed orange line) metallicities, accounting for variations in the mass-luminosity relation due to compositional differences. The effect of age on the mass-luminosity relation is negligible for low-mass stars. We take an age of 10 Gyr for all three isochrones. Most of our companions agree reasonably well with the theoretical mass-luminosity relations, except HIP 22059 B.

In the case of HIP 22059, which is under-luminous for its derived dynamical mass, this could be explained by being an unresolved binary. One or more unseen components would contribute up to ~0.07 M_⊙ to the dynamical mass, but only negligibly to the flux that we measure. This is especially true if the unresolved components are brown dwarf companions, where the flux drops off dramatically with age, and would be difficult to detect such a low-mass component in the presence of a brighter component, as well as the flux from the primary star. Alternatively, we could be observing two equal-mass components, but with a blended flux. Given that the PSF of HIP 22059 B is unresolved, these two objects must have a projected separation of less than ~40 mas, which is 1.2 AU at the distance of the host star. This is compatible with observed distributions of field M-dwarf binaries (Shan et al. 2015; Winters et al. 2019). In either case, taking into account the additional mass component would place HIP 22059 nicely on the mass-luminosity relation in Fig. 12.

We also note that the large uncertainty in the derived dynamical masses of HD 142234 B and HD 142234 C are mostly due to unconstrained orbits, as described in Sect. 5.2, the result of treating HD 142234 B and HD 142234 C as a combined mass, with their photocenter taken to be their center of mass. Furthermore, we placed an informed prior on the primary star for the orbital fit, and therefore these two comparisons with the mass-luminosity relation shown in Fig. 12 should be taken with caution.

Fig. 11 Same as Fig. 4, but for HIP 22059 B. Top left: radial-velocity orbit induced by HIP 22059 B. The blue data points show the COR-98 data, the yellow points show the HARPS data, the green points show the COR-07 data, and the cyan points show the COR-14 data.

Fig. 12 Comparison of the dynamical masses and the H-band magnitudes of the companions to theoretical mass-luminosity relations. We adopted the mass-luminosity relations from PARSEC models for low-mass stars (Bressan et al. 2012; Chen et al. 2014) at low (Z = 0.05, dotted orange line), intermediate (Z = 0.15, solid red line), and high (Z = 0.35, dashed orange line) metallicities. Most of our companions (blue dots with error bars and labels) agree reasonably well with the theoretical mass-luminosity relations, except for HIP 22059 B, as explained in Sect. 6.

7 Discussion and conclusions

In this paper we present new direct detections of a number of low-mass stellar companions, obtained by combining radial velocities, relative astrometry, and astrometric accelerations to derive precise dynamical masses of these companions. The dynamical masses of stellar companions is fundamental for a number of different reasons:

• From an exoplanetary and brown dwarf perspective, it is important to understand the true nature of these massive companions that are detected from their radial velocities with companions that are masquerading as much less massive components, ultimately leading to a different scientific interpretation of these detections.

• It is vital for future direct imaging surveys that such systems are thoroughly investigated for such low-mass stellar companions in order to make sure that targeted high-contrast imaging surveys are probing the most promising candidates for exoplanet and brown dwarf studies. This is especially important with expensive telescope facilities, like the James Webb Space Telescope (JWST), and the upcoming Extremely Large Telescopes (ELTs), where the nature of the detection needs to be understood better before spending extensive observing time for a comprehensive in-depth characterization, leading to far more efficient direct imaging surveys.

• The detection and characterization of low-mass stellar companions is vital in probing and populating stellar mass-luminosity relations, as described in Sect. 6.

As we continue to measure radial velocities and astrometry of our nearby stellar neighbors, it will become increasingly important that we understand the nature of the companions that we are investigating, which all have important astrophysical aspects to learn from, but can be scientifically interpreted very differently if their masses are unknown. By combining these techniques as described in this paper, we are able to detect companions by probing the giant planet to brown dwarf to low-mass stellar regime, and ultimately it will help us define the blurred boundaries between these different types of companions.

In conclusion, we report the discovery of new directly imaged companions orbiting HD 142234, HD 143616 and HIP 22059. Furthermore, we report the first direct detection of HD 92987 B, and update the orbital parameters and dynamical masses of two previously reported companions, HD 157338 B and HD 195010 B.

For each system we perform an orbital fit using orvara (Brandt et al. 2021c) by combining long baseline radial-velocity data from the CORALIE survey for extra-solar planets (Udry et al. 2000), along with long baseline proper motion anomalies from Hipparcos data (ESA 1997; van Leeuwen 2007; van Leeuwen & Michalik 2021) and Gaia eDR3 (Gaia Collaboration 2021) reported in the Hipparcos-Gaia catalog of accelerations (HGCA; Brandt 2021), with high-contrast images taken with VLT/SPHERE, to reveal the orbital properties and precise dynamical masses of each companion. The summary of each system is presented below:

• HD92987: We present the first direct detection of HD 92987 B from high-contrast imaging with VLT/SPHERE. By doing a combined orbital fit of CORALIE radial velocities, relative astrometry from direct imaging, and astrometric accelerations, we confirm that HD 92987 B is a low-mass star with a dynamical mass of M_Jup on a yr orbital period.

• HD 142234: We present the CORALIE radial-velocity and direct imaging observations of HD 142234 and discovered that it is a triple stellar system, where HD 142234 B and HD 142234 C are likely to be in a tightly bound orbit, orbiting around the primary star HD 142234 A. For this reason, we performed an orbital fit with the photocenter of the relative astrometry between the two components treated as the barycenter to estimate the combined mass of HD 142234 B+C, and then used the ratio of the contrast values to estimate the ratio of the mass of the individual components. From this we find the mass of HD 142234 B and HD 142234 C to be M_Jup and M_Jup, respectively.

• HD 143616: We present the CORALIE radial-velocity data for HD 143616, as well as direct imaging observations from VLT/SPHERE and discover a very low-mass star, with a dynamical mass of M_Jup. From an orbital fit combining the radial velocities, relative astrometry, and astrometric accelerations, we find this system to be highly eccentric, with an extreme eccentricity of and an orbital period of yr.

• HD 157338: We present new direct images of HD 157338 B with VLT/SPHERE. In addition, we fold in relative astrometry from 2006 VLT/NACO observations reported by Montagnier (2008), combined with radial-velocity measurements from HARPS, HIRES, and CORALIE, along with absolute astrometry to determine the dynamical mass of HD 157338 B to be 359.3 ± 10.5 M_Jup with an orbital period of yr.

• HD 195010: We present new direct images of HD 195010 B with VLT/SPHERE. In addition, we fold in relative astrometry from 2006 VLT/NACO observations reported by Montagnier (2008), combined with radial-velocity measurements from CORALIE, along with absolute astrometry to determine the dynamical mass of HD 195010 B to be M_Jup with an orbital period of yr.

• HIP 22059: We present the new detection of HIP 22059 B with CORALIE and HARPS radial-velocity data, as well as high-contrast imaging from VLT/SPHERE and absolute astrometry from HGCA to determine the dynamical mass of HIP 22059 B to be M_Jup with an orbital period of yr.

As we enter the next era of instrumentation in direct imaging exoplanets (i.e., JWST and the ELTs), it is crucial that we use the orbital information available to us through radial velocities and astrometry to thoroughly vet targets for low-mass stellar companions, and to make informed choices on target selection to prioritize the most appropriate targets for direct imaging substellar companions and exoplanets. The importance of this is demonstrated through the increasing number of studies that are using this valuable precursor information from RVs and/or astrometry to select informed targets for direct imaging searches, and combining such multiple techniques to improve orbital constraints on previously detected companions (e.g., Bonavita et al. 2020, 2022; Brandt et al. 2021a; Swimmer et al. 2022; Kuzuhara et al. 2022; Mesa et al. 2022).

This work also serves as a cautionary tale to warn about the danger of using a single detection technique to discover giant planets or low-mass brown dwarfs, which can turn out to be stellar companions, like the cases presented here for HD 92987 B and HD 195010 B. While this contributes to our understanding of low-mass stars, the science case is fundamentally different. Ultimately, this work acts as a precursor observation program as we look toward the future of the direct imaging of exoplanets.

Acknowledgements

This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. The authors acknowledge the financial support of the SNSF. This publications makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of the SIMBAD database and the VizieR Catalogue access tool, both operated at the CDS, Strasbourg, France. The original descriptions of the SIMBAD and VizieR services were published in Wenger et al. (2000) and Ochsenbein et al. (2000). This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

Appendix A DACE links

The radial-velocity measurements and the additional data products discussed in this paper are available in electronic form on the Data Analysis Center for Exoplanets (DACE) web platform for each individual target:

• HD 92987:

  https://dace.unige.ch/radialVelocities/?pattern=HD92987

• HD 142234:

  https://dace.unige.ch/radialVelocities/?pattern=HD142234

• HD 143616:

  https://dace.unige.ch/radialVelocities/?pattern=HD143616

• HD 157338:

  https://dace.unige.ch/radialVelocities/?pattern=HD157338

• HD 195010:

  https://dace.unige.ch/radialVelocities/?pattern=HD195⍉1⍉

• HIP 22059:

  https://dace.unige.ch/radialVelocities/?pattern=HIP22⍉59

Appendix B Posterior distributions of the orbital fits

Here we show the corner plots for the posterior distributions of the orbital fits for each system fitted using orvara (Brandt et al. 2021c). For each system the posterior distributions for the primary stellar mass (M_⊙), the companion mass (M_sec), the semimajor axis (a), the eccentricity (e), and the orbital inclination (i) are shown. In the case of HD 142234 the posterior distributions are shown for HD 1422234 B and HD 142234 C fitted as one single companion, as described in Section 5.

Fig. B.1 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 92987 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ± 1σ values are given at the top of each 1D distribution. The 1, 2, and 3σ contour levels are overplotted on the 2D posterior distribution.

Fig. B.2 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 142234 B and HD 142234 C combined and treated as a single companion. The posterior distributions correspond to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from HIPPARCOS and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ±1σ values are given at the top of each 1D distribution. Here 1, 2, and 3σ contour levels are overplotted on the 2D posterior distribution.

Fig. B.3 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 143616 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ±1σ values are given at the top of each 1D distribution. Here 1, 2, and 3σ contour levels are overplotted on the 2D posterior distribution.

Fig. B.4 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 157338 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions, with the median ±1σ values are given at the top of each 1D distribution. The 1, 2, and 3 σ contour levels are overplotted on the 2D posterior distribution.

Fig. B.5 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 195010 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions, with the median ±1σ values are given at the top of each 1D distribution. The 1, 2, and 3 σ contour levels are overplotted on the 2D posterior distribution.

Fig. B.6 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HIP 22059 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from HIPPARCOS and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ±1σ values are given at the top of each 1D distribution. The 1, 2, and 3 σ contour levels are overplotted on the 2D posterior distribution.

References

All Tables

All Figures

Fig. 1 Contours showing the astrometric predictions of each companion relative to the host star in ascension (Δα) and declination (Δδ) using the orvara orbit fitting package (Brandt et al. 2021c), with the measured companion positions shown by the overplotted gray stars. Relative astrometric positions indicate whether a companion can be directly detected with imaging given its projected separation between the host star and the companion itself. The 1, 2, and 3σ contour levels of the predicted positions of each companion are plotted from a full orbit fit using just the radial velocities and the HGCA astrometry, with the astrometric predictions being shown for the year in which each companion was then subsequently directly imaged with VLT/SPHERE in order to visually compare the prior predictions to the direct detections. The measured position of each companion from the VLT/SPHERE imaging are shown by the overplotted gray stars. In the case of HD 142234, the left star corresponds to HD 142234 B and the right star corresponds to HD 142234 C. The measured positions of each companion from direct imaging observations agree well with the predicted positions from the orbit fitting using just the radial velocities and astrometry from the HGCA.

In the text

	Fig. 2 High-contrast images taken with VLT/SPHERE using IRDIS coronagraphic imaging in the H2 band. In each image the star is masked behind the central white circle. Some of the coronagraphic images shown here show some saturation, which is why the flux frames were used to obtain the photometry of each companion. The image for HD 142234 is shown on a wider angular scale than the other images.
In the text

Fig. 3 Measured contrast of each companion at the measured separations with stacked 5σ contrast curves corresponding to measured on-sky performance from VLT/SPHERE data for IRDIS H2 bands. The detection probability scale is shown on the right. The darker color corresponds to a lower likelihood of detectability with imaging with VLT/SPHERE IRDIS H2 band, and the lighter colors correspond to a higher likelihood of detectability at a given separation and contrast ratio. The astrometry and photometry for HD 157338 and HD 195010 are each shown as an average of the two 2021 SPHERE H2 band data points.

In the text

Fig. 4 Orbit fits of HD 92987 using the orbit-fitting code orvara. Top left: radial-velocity orbit induced by HD 92987 B over a full orbital period. Shown are the radial-velocity data of COR-98 (blue points), COR-07 (yellow points), and COR-14 (green points). The thick line shows the highest likelihood fit; the thin colored lines show 500 orbits drawn randomly from the posterior distribution. Top right: relative astrometric orbit of HD 92987 B. The thick black line represents the highest likelihood orbit; the thin colored lines represent 500 orbits drawn randomly from the posterior distribution. Dark purple corresponds to a low companion mass and light yellow corresponds to a high companion mass. The dotted black line shows the periastron passage, and the arrow at the periastron passage shows the direction of the orbit. The dashed line indicates the line of nodes. Predicted past and future relative astrometric points are shown by black circles with their respective years, while the observed relative astrometric point from VLT/SPHERE data is shown by the blue-filled data point, where the measurement error is smaller than the plotted symbol. Bottom: acceleration induced by the companion on the host star as measured from absolute astrometry from Hipparcos and Gaia. The thick black line represents the highest likelihood orbit; the thin colored lines are 500 orbits drawn randomly from the posterior distribution. Darker purple represents a lower companion mass and light yellow represents a higher companion mass. The residuals of the proper motions are shown in the bottom panels.

In the text

	Fig. 5 Same as Fig. 4, but for HD 142234 B+C. Here the orbital fit is performed with HD 142234 B+C together.
In the text

	Fig. 6 Same as Fig. 4, but for HD 143616 B.
In the text

Fig. 7 Same as Fig. 4, but for HD 157338 B. Top left: radial-velocity orbit induced by HD 157338 B. The white data points show the HIRES data, the blue data points show the COR-98 data, the cyan data points show the HARPS data, the yellow points show the COR-07 data, and the green points show the COR-14 data. Top right: predicted past and future relative astrometric points are shown as black circles with their respective years, while the observed relative astrometric points from the 2006 VLT/NACO data and the 2021 VLT/SPHERE data are shown as blue filled data points, where the measurement error is smaller than the plotted symbol.

In the text

Fig. 8 Relative separation and position angle of HD 157338 B across three epochs of data. The 2006 data corresponds to VLT/NACO data as reported by Montagnier (2008), and the two 2021 data points correspond to the VLT/SPHERE data presented in this paper. The second epoch of direct imaging observations over a long baseline is able to constrain the orbital parameters significantly, as shown in the residual plots in the corresponding bottom panels.

In the text

Fig. 9 Orbit fits of HD 195010 using the orbit-fiting code orvara. Top left: same as Fig. 4, but for HD 195010 B. Top right: predicted past and future relative astrometric points are shown as black circles with their respective years, while the observed relative astrometric points from the 2006 VLT/NACO and the 2021 VLT/SPHERE data is shown as blue filled data points, where the measurement error is smaller than the plotted symbol.

In the text

	Fig. 10 Relative separation and position angle of HD 195010 B across three epochs of data. The 2006 data corresponds to the VLT/NACO data reported by Montagnier (2008), and the 2021 data point correspond to the VLT/SPHERE data presented in this paper. The residuals are shown in the bottom two panels.
In the text

	Fig. 11 Same as Fig. 4, but for HIP 22059 B. Top left: radial-velocity orbit induced by HIP 22059 B. The blue data points show the COR-98 data, the yellow points show the HARPS data, the green points show the COR-07 data, and the cyan points show the COR-14 data.
In the text

Fig. 12 Comparison of the dynamical masses and the H-band magnitudes of the companions to theoretical mass-luminosity relations. We adopted the mass-luminosity relations from PARSEC models for low-mass stars (Bressan et al. 2012; Chen et al. 2014) at low (Z = 0.05, dotted orange line), intermediate (Z = 0.15, solid red line), and high (Z = 0.35, dashed orange line) metallicities. Most of our companions (blue dots with error bars and labels) agree reasonably well with the theoretical mass-luminosity relations, except for HIP 22059 B, as explained in Sect. 6.

In the text

Fig. B.1 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 92987 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ± 1σ values are given at the top of each 1D distribution. The 1, 2, and 3σ contour levels are overplotted on the 2D posterior distribution.

In the text

Fig. B.2 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 142234 B and HD 142234 C combined and treated as a single companion. The posterior distributions correspond to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from HIPPARCOS and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ±1σ values are given at the top of each 1D distribution. Here 1, 2, and 3σ contour levels are overplotted on the 2D posterior distribution.

In the text

Fig. B.3 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 143616 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ±1σ values are given at the top of each 1D distribution. Here 1, 2, and 3σ contour levels are overplotted on the 2D posterior distribution.

In the text

Fig. B.4 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 157338 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions, with the median ±1σ values are given at the top of each 1D distribution. The 1, 2, and 3 σ contour levels are overplotted on the 2D posterior distribution.

In the text

Fig. B.5 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HD 195010 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from Hipparcos and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions, with the median ±1σ values are given at the top of each 1D distribution. The 1, 2, and 3 σ contour levels are overplotted on the 2D posterior distribution.

In the text

Fig. B.6 Marginalized 1D and 2D posterior distributions for selected orbital parameters of HIP 22059 B corresponding to the fit of the RV, relative astrometry from direct imaging observations, and absolute astrometry from HIPPARCOS and Gaia with the use of orvara (Brandt et al. 2021c). Confidence intervals at 15.85%, 50.0%, 84.15% are overplotted on the 1D posterior distributions; the median ±1σ values are given at the top of each 1D distribution. The 1, 2, and 3 σ contour levels are overplotted on the 2D posterior distribution.

In the text