© ESO 2022

1 Introduction

Dynamical masses for components of stellar binaries are critical for benchmarking the computational models used to estimate the masses and radii of isolated objects based on their luminosities. While abundant dynamical masses for high-mass and solar-mass stellar binaries are reported in the literature, fewer dynamical masses are available for M dwarfs with masses <0.3 M_⊙ and ultracool dwarfs, especially atyoung (<1 Gyr) ages. In the last decade, empirical mass-luminosity relationships have been constructed for late K dwarf and M dwarf stars. Benedict et al. (2016) constructed a mass-luminosity relationship for M dwarf stars via Hubble Space Telescope orbital monitoring of 47 stars with masses from 0.08–0.62 M_⊙. Mann et al. (2019) derived an empirical mass-luminosity-metallicity relation for stellar masses between 0.075 and 0.7 M_⊙ by measuring total system masses from 62 nearby M dwarf binaries.

Many of the dynamical masses in the literature determined via direct imaging of visual binaries are in fact total masses for the system. Direct imaging of the motion of a companion relative to its primary must be combined with either the radial velocity (RV) or absolute astrometry of the primary in order to determine individual component masses. Combining Canada France Hawaii Telescope (CFHT) absolute astrometry with adaptive-optics-assisted resolved imaging relative to the primary, Dupuy & Liu (2017) present 38 individual masses for ultracool dwarfs (spectral type ≥M6), ranging from 30 to 115 M_Jup, the largest such sample for substellar objects and ultracool dwarfs in the literature. A number of low-mass and substellar companions have been directly imaged after their presence was indirectly noted via an RV trend (e.g., Crepp et al. 2016). For companions to relatively bright stars, the combination of Gaia Data Release 2 (DR2) or Early Third Data Release (EDR3) proper motions with HIPPARCOS proper motions can provide additional absolute astrometry orbital constraints (Kervella et al. 2019). Only in the last few years have these techniques been fully combined to determine dynamical masses for a range of stars, substellar objects, and exoplanets (Bowler et al. 2018; Calissendorff & Janson 2018; Dieterich et al. 2018; Snellen & Brown 2018; Grandjean et al. 2019; Dupuy et al. 2019; Brandt et al. 2019, 2020). Here we present dynamical mass determinations, combining RV, direct imaging, and Gaia-HIPPARCOS data for all components of the HIP 113201 and HIP 36985 binary systems.

HIP 113201AB and HIP 36985AB are two M1 + mid-M binary systems. The existence of these binaries was first suggested via their proper motion anomalies (Kervella et al. 2019) – slight differences between the proper motions of the primaries as measured by HIPPARCOS and Gaia suggested the presence of an additional component in each system. Both binaries have since been confirmed via direct imaging and RV observations. HIP 113201B was discovered by Bonavita et al. (2021) in the course of the SpHere INfrared survey for Exoplanets (SPHERE-SHINE) survey (Chauvin et al. 2017), a ~500 star survey targeting young stellar systems using the SPHERE planet-finding camera (Beuzit et al. 2019) at the Very Large Telescope (VLT) in order to directly image young giant exoplanets (Desidera et al. 2021; Langlois et al. 2021; Vigan et al. 2021). The HIP 113201AB system also has undergone extensive RV monitoring with the High Accuracy Radial velocity Planet Searcher (HARPS). For HIP 36985B, Kammerer et al. (2019) first imaged the companion via VLT/NaCo kernel phase imaging and direct VLT/NaCo imaging. Grandjean et al. (2020) confirmed that HIP 36985B was bound and covered a significant fraction of the companion’s orbit via HARPS RV observations. Bonavita et al. (2021) provided additional epochs of VLT-SPHERE imaging of this companion, and Baroch et al. (2021) has recently published additional CARMENES and FEROS RV observations. With preliminary dynamical masses <0.2 M_⊙ (Bonavita et al. 2021; Baroch et al. 2021), HIP 113201B and HIP 36985B are two of the lowest-mass stellar companions imaged as part of the SPHERE-SHINE survey.

In Sect. 2 we summarize the properties of both stellar systems and present updated ages for the systems. In Sect. 3 we describe the parallel-tempering Markov chain Monte Carlo (PT-MCMC) code used to fit orbits for both systems. Section 4 covers results from our PT-MCMC fits to both systems. We compare our derived dynamical masses for all components of both systems to existing mass estimates in Sect. 5 and summarize our conclusions in Sect. 6.

2 Properties of the primary stars

2.1 HIP 113201/GJ 4303/BPM 28050

Gaidos et al. (2014) fit the optical spectrum of HIP 113201 and find a best fit spectral type of M1, with best fit values of T_eff = 3693 ± 82 K, luminosity = 0.038 ± 0.009 L_⊙, radius = 0.49 ± 0.05 R_⊙, and mass = 0.52 ± 0.07 M_⊙. Hawley et al. (1996) find a similar but slightly earlier M0.5 spectral type from low-resolution optical spectroscopy using the Cerro Tololo Inter-American Observatory (CTIO) 1.5 m telescope. We used the M_K-mass relationship for stars with masses between 0.075 and 0.7 M_⊙ derived from orbital fits to 62 nearby binaries by Mann et al. (2019) to estimate a mass of 0.52 ± 0.01 M_⊙ for HIP 113201, matching the Gaidos et al. (2014) spectroscopic mass estimate. Stellar properties are presented in Table 1.

HIP 113201 was included in the SPHERE-SHINE survey sample because its kinematics are a reasonable match to both the Tucana-Horologium and β Pic moving groups. However, using a variety of age diagnostics, the star appears to be older than comparable spectral type members in both groups. Its weak H-alpha emission coupled with a lack of lithium absorption suggest an age equal toor greater than the Hyades – in other words, somewhat younger than typical field stars, but not <300 Myr.

To robustly determine the age via gyrochronology, we undertook a photometric monitoring campaign in October-November 2016 at the Remote Observatory Atacama Desert (ROAD) located in the Atacama Desert close to the town of San Pedro de Atacama, Chile. The telescope is a 40 cm f /6.8 Optimized Dall-Kirkham equipped with a 4K × 4K pixel FLI ML16803 charge-coupled device camera (9 μm pixel size) with a 40′   × 40′ field of view and BV I filters. We observed HIP 113201 for 36 almost consecutive nights spanning a time interval of 40 days. We collected a total of 179 frames in V filter and 179 frames in I filter (generally one telescope pointing per night with five consecutive frames in each filter that were averaged to get one average magnitude per night and corresponding σ). We used aperture photometry to extract the magnitudes and computed differential magnitudes of HIP 113201with respect to an ensemble comparison consisting of three non-variable stars (TYC 8453 806 1; CD−539254; TYC 8453 856 1). The photometric precision is 0.007 mag in both the I and V filters. We performed Lomb-Scargle (LS) and CLEAN periodogram analysis on the average data for the V, the I, and the (V + I)/2 time series. The analysis of the (V + I)/2 time series provides themost precise results, which are summarized in Fig. 1. In the periodogram of V-filter data we find a period of 19.9 ± 0.7 days and a light curve amplitude 0.03 mag. In the periodogram of I-filter data we find a period of 19.2 ± 0.7 days and a light curve amplitude 0.02 mag. In the periodogram of the combined light curve, we find a period of 19.6 ± 0.5 days, which we adopt as the photometric rotation period for this star. The rotation period is highly significant (confidence level >99%) with coverage over multiple rotation periods. The V and I magnitude variations are correlated with linear Pearson correlation coefficient 0.6 and significance 99.9%, suggesting that the photometric variability is dominated by either cool spots or hot spots (faculae).

HIP 113201 was also monitored by the Transiting Exoplanet Survey Satellite (TESS) continuously over 25 days in both Sector 1 (2019) and Sector 28 (2020). Visual inspection of the TESS light curves shows that the Pre-search Data Conditioned Simple Aperture Photometry (PDCSAP) light curves appear to have additional trends added by the PDCSAP process, and hence we adopt the SAP light curve instead. An LS periodogram analysis of the SAP light curve from Sector 1 yields a broad double peaked periodogram with the peak covering periods from 17 to 23 days and highest individual peak powers at 17.8 and 22.7 days. The LS periodogram for Sector 28 has a power peak at the period of 18.1 days. The TESS monitoring only covers ~1 rotational period for this star, so does not have the time baseline necessary for a robust period determination. Hence we adopt the 19.6 ± 0.5 days period from our dedicated ground-based monitoring campaign.

With a comparatively long rotation period given its spectral type, HIP113201 is on the slow rotator branch. Using the Delorme et al. (2011) calibration with the 19.6 ± 0.5 days period yields an age of 1.2 ± 0.1 Gyr for this star. The stellar properties of HIP 113201 are summarized in Table 1.

Fig. 1 Summary of periodogram analysis for HIP 113201. Top-left panel: combined V + I-band magnitudes versus HJD (-2457670). Measurements on each night were taken in both the V band and the I band; the average of the two bands is plotted here. Top-middle panel: LS periodogram. The horizontal dotted lines show the 99.9 and 98% confidence levels, whereas the red dotted line shows the spectral window function. Top-right panel: clean periodogram. Bottom panel: light curve phased with the 19.6 d rotation period. Red circles are data points from the first half of the observations; blue circles are data points from the second half of the observations. The red dotted line is a sinusoidal fit to the phased light curve.

2.2 HIP 36985/GJ 282C

With a separation of >50 000 AU, HIP 36985 is one of the widest known companions, and is a triple component to the K dwarf binary GJ 282AB (Poveda et al. 2009). Gaidos et al. (2014) fit the optical spectrum of HIP 36985 and find a best fit spectral type of M1, with best fit values of T_eff = 3744 ± 82 K, luminosity = 0.044 ± 0.01 L_⊙, radius = 0.51 ± 0.05 R_⊙, and mass = 0.54 ± 0.07 M_⊙. This M1 spectral type is independently confirmed by Alonso-Floriano et al. (2015), from low-resolution optical spectroscopy with the Calar Alto Faint Object Spectrograph (CAFOS) at the 2.2 m Calar Alto telescope. Using the M_K -mass relationship for stars with masses between 0.075 and 0.7 M_⊙ derived from orbital fits to 62 nearby binaries by Mann et al. (2019), we estimate a mass of 0.57 ± 0.01 M_⊙, matching well the mass estimate from spectroscopy.

The X-ray emission of GJ282A and B (individual components from ROSAT pointed observations, Schmitt & Liefke 2004) yields log (L_x∕L_bol) = −4.62 and −4.57, respectively. These values are within the distribution of the Hyades members but slightly above the median values and well below the values observed for Pleiades and AB Dor MG members (Desidera et al. 2015). The log $R\text{^{'}\_HK}$ of GJ282A (−4.39 Wright et al. 2004) is at the upper envelope of the Hyades members. The activity indicators then suggest an age slightly younger than the Hyades.

Lithium can be used as additional diagnostic: Ramírez et al. (2012) find A(Li) = 0.13 ± 0.04 for HIP 36985, while lithium non-detection has been reported by Mishenina et al. (2012) and Luck (2017). At the color of the star, a lithium non-detection implies an age older than about 550 Myr. The spectroscopic analysis in the literature supports a metallicity close to solar ([Fe/H] = −0.09: Tabernero et al. 2017; −0.12: Mishenina et al. 2012; +0.01: Ramírez et al. 2012; 0.00: Luck 2017; [M/H] = −0.05, Valenti & Fischer 2005).

HIP 36985 has been monitored photometrically in I and Z by the All Sky Automated Survey (ASAS; Pojmanski 2002), with inhomogeneous sampling over a time baseline of >2000 days. Using these data, Kiraga (2012) and Díez Alonso et al. (2019) derive consistent rotation periods of 12.2 ± 0.1 days and 12.16 days, respectively. HIP 36985 was also monitored by the TESS continuously over 25 days. Similarly as for the case of HIP 113201, visual inspection of the TESS light curves shows that the PDCSAP light curve appear to have additional trends added by the PDCSAP process. Thus, we again adopt the simple aperture photometry (SAP) light curve instead. An LS periodogram analysis of the Simple SAP light curve from TESS yields a periodogram power peak at the period of ~12.1 days.

With a best rotational period of 12.2 ± 0.1 days from literature values, HIP 36985 is a relatively fast rotator for its spectral type. Since the best adopted period is shorter than the period at the convergence time on the slow rotator branch for this spectral type, gyrochronological relations such as those from Delorme et al. (2011) cannot be applied. In this case, gyrochronology can only indicate that the age of the star is compatible with any age between 0 and the age of the Hyades. Current best estimates for the age of the Hyades range from 625 ± 50 Myr (Delorme et al. 2011), 650 ± 70 Myr (Martín et al. 2018), and 750 ± 100 Myr (Brandt & Huang 2015). Given that the rotation period of HIP 36985 is only slightly shorter than the period at the convergence time for an M1 spectral type, its age is likely similar to or slighter younger than the age of the Hyades, consistent with the age estimates from coronal and chromospheric emission.

The two higher-mass components of this system have K spectral types and are thus actually better targets for gyrochronological age dating as they are older than the convergence time for these earlier spectral types. Both were observed as part of ASAS, but no significant periodicities were found in these data. GJ 282A and B were covered in TESS Sector 7. An LS periodogram analysis yielded the following periods for GJ 282 A and B: P_A = 13 ± 3.5d and P_B = 14 ± 4.0d, of insufficient precision for an accurate gyrochronological age dating. Given the inconclusive results from the ASAS and TESS light curves,we conducted a photometric monitoring campaign for GJ 282A and B, from January 15 until April 20, 2021, for a total of 75 nights at the ROAD observatory (see Sect. 2.1 for details). Observations were collected in the V and B filters and consisted of five consecutive frames per filter on each night, totaling 750 frames. The magnitudes of GJ282A and B were extracted using aperture photometry together with two nearby stars that served as comparison (C; BD-03 2003) and check (CK: HD 61723) stars. The comparison and check stars turned out to be constant in flux during the monitoring with a standard deviation σ_V = 0.017 mag and σ_B = 0.012 mag. Both A and B components show a low level of variability with a standard deviation of their differential magnitude time series similar to that of the C−CK differential magnitudes with ${\sigma }_{V_{\text{A}}}$ = 0.017 mag; ${\sigma }_{B_{\text{A}}}$ = 0.013 mag and ${\sigma }_{V_{\text{B}}}$ = 0.015 mag; ${\sigma }_{B_{\text{B}}}$ = 0.015 mag. The periodogram analysis, carried out with LS and CLEAN, did not reveal any significant periodicity for the A component. In contrast, the periodogram analysis of the B component showed a significant power peak in both V and B band at P = 12.10 ± 0.77 days. In Fig. 2, we summarize the results of our periodogram analysis. To make the phase rotational modulation more easily visible, we averaged all magnitudes within bins of 0.05 in phase.

Using the gyrochronology relationship from Delorme et al. (2011) and adopting an age for the Hyades of 625 Myr, the period of 12.10 ± 0.77 days for GJ 282B implies an age of 740 ± 100 Myr. This iscompatible with the gyrochronological age limit of the Hyades age or younger that we found for HIP 36985. Thus, combining this with the agelimits placed by lithium non-detection and activity indicators, we adopted a system age similar to that of the Hyades for GJ 282AB/HIP 36985.

The system has been considered for membership in the ~300 Myr Ursa Major (UMa) moving group in several works (e.g., Montes et al. 2001). Most recently, Baroch et al. (2021) recalculated the galactocentric velocity for HIP 36985, finding 25.1, −2.5, and −7.8 km s⁻¹ for the U, V, and W components, respectively, and suggest that it may be a member of the UMa moving group. These values are very similar to those of GJ 282A (U,V,W = 25.20 ± 0.26, −3.07 ± 0.25, −7.65 ± 0.09 km s⁻¹) listed by Tabernero et al. (2017), further supporting the physical association of HIP 36985 with GJ 282AB. Montes et al. (2001) find a mean velocity vector of the UMa nucleus of U, V, W = 14.9, 1.0, −10.7 km s⁻¹; King et al. (2003) and Mamajek et al. (2010) find updated mean velocity vectors of U, V, W = 14.2 ± 0.7, 2.8 ± 1.3, −8.7 ± 1.8 km s⁻¹ and U, V, W = 15.0 ± 0.4, 2.8 ± 0.7, −8.1 ± 1.0 km s⁻¹, respectively. Both HIP 36985 and GJ 272A have velocity vectors close to but significantly outside the UMa nucleus. Combined with the higher system age we find from gyrochronology, and consistent with estimates from coronal and chromospheric emission, we conclude that this stellar system is not part of the nucleus of the UMa group. The stellar properties of HIP 36985 are summarized in Table 1.

3 Description of the parallel-tempering Markov chain Monte Carlo orbital fitting code

We used the PT-MCMC ensemble sampler as implemented in emcee v2.1.0 (Foreman-Mackey et al. 2013) to fit the combination of RV, direct imaging, and Gaia-HIPPARCOS astrometry available for each system. The PT-MCMC method runs chains at different “temperatures” (Earl & Deem 2005), where higher temperature chains sample large volumes of phase space, while lower temperature chains sample more precisely a local region of phase space. At the end of the MCMC run, only the results from the lowest temperature chain are retained; however, in the PT-MCMC implementation, the higher and lower temperature chains swap information, preventing the lower temperature chains from becoming trapped in local minima. Thus, the PT-MCMC method is particularly useful for evaluating multimodal or complex posteriors.

Fig. 2 Summary of periodogram analysis for GJ 282B. Top-left panel: V-band magnitudes versus HJD (-2459000). Top-middle panel: LS periodogram. The horizontal dotted lines show the 99.9 and 98% confidence levels, whereas the red dotted line is the spectral window function. Top-right panel: clean periodogram. Bottom panel: light curve phased with the 12.15 d rotation period. The red line is a sinusoidal fit to the phased light curve.

3.1 Simulating trial orbits

At each step of each MCMC chain, we simulate trial orbits using the following parameters (including the standard orbital elements): the orbital period (P) in years; the orbital inclination (i) in degrees; the orbital eccentricity (e); the RV measurement offset (γ) in km s⁻¹; the longitude of the ascending node (Ω) in degrees; the argument of periastron (ω) in degrees; the time of periastron passage (T₀) in Julian days (JD); the distance to the system (d) in pc; the total system mass (M_tot) in M_⊙; and the companion mass (M₂) in M_⊙.

We constructed simulated orbits following the approach given in Argyle (2012). Time (t) in JD at each epoch is converted tomean anomaly (M) via $$M=\frac{2\pi }{P}(t-T_0).$$(1)

Mean anomaly (M) is converted to eccentric anomaly (E) by using theNewton-Raphson method to iteratively solve Kepler’s equation: $$M=E-\sin E.$$(2)

Radial velocity variations were modeled at each epoch t as $$V(t)=K[\cos (\omega )+\nu (t))+e\cos (\omega )],$$(3)

where ω is the argument of periastron, ν(t) is the true anomaly at epoch t, e is the eccentricity, and the true anomaly ν(t) is related to the eccentric anomaly E by $$\nu =2 \text{arctan}\left(\left(\sqrt{\frac{1+e}{1-e}}\right)\text{tan}\left(\frac{E}{2}\right)\right).$$(4)

As the arctan function in python only returns angles from $-\frac{\pi }{2}$ to $\frac{\pi }{2}$, we adopted the following form for ν to cover the full range of possible true anomaly values: $$\nu =2 \arctan 2\left(\left(\sqrt{1+e}\right)\sin \left(\frac{E}{2}\right),\left(\sqrt{1-e}\right)\cos \left(\frac{E}{2}\right)\right).$$(5)

The K is the RV semi-amplitude, given by $$K={\left(\frac{2\pi G}{P(\text{sec})}\right)}^{\frac{1}{3}}\frac{M_2 \sin (i)}{{(1-e)}^2{M}_{\text{tot}}^{\frac{2}{3}}},$$(6)

where P(sec) is the orbital period in seconds, M₂ is the companion mass, M_tot is the total system mass, e is eccentricity and i is the inclination of the system. We adopt an additional offset parameter γ for each RV data set, to take into account the systematic errors between instruments.

Direct imaging relative companion positions (with respect to the stellar position) were modeled using the Thiele-Innes elements. The semimajor axis a of the system in AU is calculated from Kepler’s third law: $$a={\left(P^2 M_{\text{tot}}\right)}^{1/3}.$$(7)

This is projected onto the sky using the Thiele-Innes elements: $$A=a (\cos (\Omega ) \cos (\omega )-\sin (\Omega ) \sin (\omega ) \cos (i))$$(8) $$B=a (\sin (\Omega ) \cos (\omega )+\cos (\Omega ) \sin (\omega ) \cos (i))$$(9) $$F=a (-\cos (\Omega ) \sin (\omega )-\sin (\Omega ) \cos (\omega ) \cos (i))$$(10) $$G=a (-\sin (\Omega ) \sin (\omega )+\cos (\Omega ) \cos (\omega ) \cos (i)).$$(11)

At a given time (t), mean anomaly (M), and eccentric anomaly (E), the projected separations (in AU) of the companion from the primary star in the x direction, y direction, or radially are given by $$X=\cos (E)-e$$(12) $$Y=\sqrt{1-e^2} \sin (E)$$(13) $$R=a \sqrt{X^2+Y^2}.$$(14)

Following a similar approach as Grandjean et al. (2019), the tangential velocity of the companion at time t is given by $$V_X=-\frac{2 \pi a^2 \sin (E)}{P R}$$(15) $$V_Y=\frac{2 \pi a^2 \cos (E)}{P R} \sqrt{1-e^2}.$$(16)

Projecting onto the sky to obtain separations and velocities in the right ascension (RA) and declination (Dec) directions then provides the following companion positions and tangential velocities relative to the primary star: $$\Delta \text{Dec}=A X+F Y \text{in} \text{AU}$$(17) $$\Delta \text{RA}=B X+G Y \text{in} \text{AU}$$(18) $$\Delta v_{\text{Dec}}=\frac{A V_X+F V_Y}{a} \text{in} \text{AU/yr}$$(19) $$\Delta v_{\text{RA}}=\frac{B V_X+G V_Y}{a} \text{in} \text{AU/yr}.$$(20)

To convert onto the on-sky projection requires a distance measurement, we included a fit to the Gaia DR2 distance in our likelihood function to help pin the distance of the system to an appropriate value.

We fit the RV and direct imaging points directly; as epoch astrometry is not yet available for Gaia, we fit the average model tangential velocity in RA and Dec to the HIPPARCOS and Gaia proper motion anomaly measurements. To correctly compare Gaia and HIPPARCOS proper motion measurements to our model orbits then requires also taking the (model-dependent) barycentric motion of the system into account.

3.2 Simulating HIPPARCOS-Gaia proper motion anomalies

For each primary star, the combination of HIPPARCOS and Gaia data yields three measurements: μ_H, the proper motion as measured over the HIPPARCOS mission, μ_G, the proper motion measured over the portion of the Gaia mission corresponding to a Gaia data release (e.g., Gaia DR2 or EDR3), and μ_HG, the scaled positional difference between the HIPPARCOS and Gaia positional measurements over the full ~25 yr baseline between both missions (Brandt 2018, 2021; Brandt et al. 2019). In the case of a single star, all three of these measurements will be the same. A binary companion tugging on the primary star will cause deviations in proper motion on the timescale of the orbit of the companion (see, e.g., Fig. 1 of Kervella et al. 2019). Thus, the “proper motion anomaly”, the difference between the measured proper motion at the HIPPARCOS or Gaia epoch and the long-term proper motion of the photocenter of the system, as defined in Kervella et al. (2019), $$\Delta {\mu }_{\text{H/G2}}={\mu }_{\text{H/G2}}-{\mu }_{\text{HG}},$$(21)

can be an important indicator of binarity and can provide additional orbital constraints. For both HIP 113201 and HIP 36985, the significant proper motion anomaly in both cases (Kervella et al. 2019; Bonavita et al. 2021) strongly pointed to the existence of the binary companion. For the purpose of orbit fitting, the proper motion anomaly (an average quantity over a given mission time span) is compared with either the instantaneous proper motion from the model at the midpoint of that mission time span or the mean value of proper motion from the model over the mission time span. A shorter proper motion monitoring time will thus provide a proper motion anomaly measurement that is a better approximation to the instantaneous or mean model proper motion. If a given mission time span catches a part of the orbit with particularly high tangential velocity, a believable fit to the proper motion anomaly for that mission may not be possible. Thus, in the future, for orbit fitting, it would be preferable to directly fit epoch-by-epoch measurements. However, epoch astrometry will not be available until Gaia DR4 and later. For now, average measurements from astrometry over 1.5–2 yr time spans, while not entirely accurate, can place valuable additional constraints on the orbit. Because of this issue, we chose to fit Gaia DR2 proper motions instead of the more recent Gaia EDR3 proper motions, as Gaia DR2 covers a shorter on-sky period and hence provides a more “instantaneous” measurement.

The proper motion anomaly is measured relative to the photocenter of the system (assumed to be the primary star), whereas model quantities are calculated relative to the barycenter of the system. Thus, to compare observed and model quantities, we must account for the barycenter position of the system relative to the photocenter. At a given time, t, the instantaneous on-sky tangential velocity (in mas) of the primary (V_primary) due to the pull of the secondary (v_secondary) relative to the barycenter of the system in RA and Dec, respectively, is $$V_{\text{primary}, \text{RA}}=-\frac{M_2}{M_{\text{tot}}} v_{\text{secondary}, \text{RA}} \frac{1000}{d}$$(22) $$V_{\text{primary}, \text{Dec}}=-\frac{M_2}{M_{\text{tot}}} v_{\text{secondary}, \text{Dec}} \frac{1000}{d}.$$(23)

However, as epoch astrometry is not available for intermediate Gaia data releases, what we measure instead is the average over all the Gaia or HIPPARCOS individual measurement epochs for each data release. For each model orbit, we calculate this as the average of the tangential velocities for each of the observed epochs: $$\langle V_{\text{primary}, \text{RA}}\rangle =-\frac{M_2}{M_{\text{tot}}} \langle v_{\text{secondary}, \text{RA}}(t)\rangle \frac{1000}{d}$$(24) $$\langle V_{\text{primary}, \text{Dec}}\rangle =-\frac{M_2}{M_{\text{tot}}} \langle v_{\text{secondary}, \text{Dec}}(t)\rangle \frac{1000}{d}.$$(25)

For HIPPARCOS, we retrieved the observation dates using the HIPPARCOS intermediate observation app¹ and for Gaia, we used the Gaia Observation forecast tool (GOST)² to select observations within the Gaia DR2 observation period, calculated the instantaneous tangential velocity for each individual observation, and then adopted the average of these values to compare against the (averaged) HIPPARCOS or Gaia proper motion.

The tangential velocity calculated from the model gives the instantaneous motion on the sky due to the orbit of the secondary, whereas the proper motion anomaly calculated via Eq. (21) subtracts the long-term proper motion trend from the “instantaneous” Gaia or HIPPARCOS proper motion to approximate this same value. However, all Gaia and HIPPARCOS proper motion values are calculated relative to the photocenter of the system (assumed to be the primary in this case, as the photometric contribution of the secondary is negligible), while the average model tangential velocity calculated via Eq. (25) is relative to the barycenter of the system. Thus, we must correct the proper motion anomaly measurement to the barycenter of the system. This correction depends on the masses of both components of the binary, and hence, is model dependent. Again, following Grandjean et al. (2019), the position of the star relative to the barycenter of the system due to reflex motion from the influence of the companion is given by $$\Delta {\text{Dec}}_{*}=-\Delta \text{Dec} \frac{M_2}{M_{\text{tot}} d}$$(26) $$\Delta {\text{RA}}_{*}=-\Delta \text{RA} \frac{M_2}{M_{\text{tot}} d}.$$(27)

The barycenter proper motion, as calculated from the scaled HIPPARCOS-Gaia positional difference, is then given by $${\mu }_{\text{bary} \text{RA}}=-\frac{M_2}{M_{\text{tot}} d}\frac{{\langle \Delta {\text{RA}}_{*}\rangle }_{\text{G}}-{\langle \Delta {\text{RA}}_{*}\rangle }_{\text{H}}}{\langle t_{\text{H}}\rangle -\langle t_{\text{G}}\rangle }$$(28) $${\mu }_{\text{bary} \text{Dec}}=-\frac{M_2}{M_{\text{tot}} d}\frac{{\langle \Delta {\text{Dec}}_{*}\rangle }_{\text{G}}-{\langle \Delta {\text{Dec}}_{*}\rangle }_{\text{H}}}{\langle t_{\text{H}}\rangle -\langle t_{\text{G}}\rangle }.$$(29)

To compare between model and data quantities, we then subtracted out the model barycenter proper motion from the observed proper motion anomaly (in RA and Dec, respectively): $$\Delta {\mu }_{\text{bary}}={\mu }_{\text{H/G2}}-{\mu }_{\text{HG}}-{\mu }_{\text{bary}}.$$(30)

We can then compare Δμ_ref in RA and Dec with V_primary, RA and V_{primary, Dec}. Errors for Δμ_ref are added in quadrature.

3.3 Likelihood function and priors

Combining direct imaging, RV, and HIPPARCOS-Gaia astrometry, our final combined likelihood function for the PT-MCMC runs is $$\begin{array}{l}\mathcal{L}=-\frac{1}{2}[{\sum }_{i=1}^{N_{DI}}{\left(\frac{x_i^{DI}-f_i^{DI}}{{\sigma }_i^{DI}}\right)}^2+{\sum }_{i=1}^{N_{DI}}{\left(\frac{y_i^{DI}-f_i^{DI}}{{\sigma }_i^{DI}}\right)}^2\\ +{\sum }_{i=1}^{N_{\text{inst}}}{\sum }_{j=1}^{N_{\text{RV}}}{\left(\frac{{\text{RV}}_j-(f_j^{\text{RV}}+{\gamma }_i)}{{\sigma }_j^{\text{RV}}}\right)}^2+{\left(\frac{\Delta {\mu }_{\text{bary},\text{H},\text{RA}}-V_{\text{primary},\text{H},\text{RA}}}{{\sigma }_{\Delta {\mu }_{\text{bary},\text{H},\text{RA}}}}\right)}^2\\ +{\left(\frac{\Delta {\mu }_{\text{bary},\text{H},\text{Dec}}-V_{\text{primary},\text{H},\text{Dec}}}{{\sigma }_{\Delta {\mu }_{\text{bary},\text{HIP},\text{Dec}}}}\right)}^2+{\left(\frac{\Delta {\mu }_{\text{bary},\text{G},\text{RA}}-V_{\text{primary},\text{G},\text{RA}}}{{\sigma }_{\Delta {\mu }_{\text{bary},\text{G},\text{RA}}}}\right)}^2\\ +{\left(\frac{\Delta {\mu }_{\text{bary},\text{G},\text{Dec}}-V_{\text{primary},\text{G},\text{Dec}}}{{\sigma }_{\Delta {\mu }_{\text{bary},\text{G},\text{Dec}}}}\right)}^2+{\left(\frac{{\varpi }_{DR2}-(1000/d)}{{\sigma }_{{\varpi }_{\text{DR2}}}}\right)}^2],\end{array}$$

where $x_i^{DI}$ and $y_i^{DI}$ are the direct imaging offsets of the secondary from the primary position in RA and Dec (ΔRA and ΔDec from Sect. 3.1), respectively, for each data epoch, $f_i^{DI}$ are the model predictions for the secondary position relative to the primary at each data epoch, RV_j is the measured RV at each data epoch, $f_j^{\text{RV}}$ is the model prediction for RV at that data epoch, γ_i are the instrumental offsets between the model RV and the measured RV for each instrument i, Δμ_bary values are the barycenter corrected proper motion anomaly in RA and Dec for HIPPARCOS and Gaia, respectively, ϖ_DR2 is the measured Gaia DR2 parallax in mas, and (1000∕d) is the model prediction for the parallax of the system (with distance d given in pc). The values labeled as σ are the respective errors on each measured property. For properties that combine multiple measurements (e.g., various Δμ_bary values), ${\sigma }_{\Delta {\mu }_{\text{bary}}}$ has been calculated as the sum in quadrature of all the constituent measured values. The PT-MCMC fits do not actually serve to put useful constraints on the distance to the system; the fit to the Gaia DR2 parallax is included instead as an additional prior to weight the PT-MCMC runs to distance values that are consistent with the measured Gaia parallax.

For all PT-MCMC runs, we adopted the following uniform priors: for P, 0 to 70 yr; for i,Ω,ω, 0 to 2π radians; for e, 0 to 1; for γ, 20 km s⁻¹ to 20 km s⁻¹, fit separately for each instrument used; for T₀, uniform between 2 450 000 days < T₀ < 2 464 500 days for HIP 113201 and between 2 454 000 days < T₀ < 2 466 000 days for HIP 36 985; for M_tot, 0–2 M_⊙; and for M₂, 0–0.5 M_⊙.

4 Orbital fits

4.1 HIP 113201

For direct imaging points, we fit the VLT-SPHERE astrometric data points presented in Bonavita et al. (2021). HIP 113201 also has multiple years (spanning 2008–2017) of RV monitoring using the HARPS instrument on the 3.6 m telescope at the La Silla observatory. We derived RV measurements from all spectra available on the ESO archive using the SAFIR (Spectroscopic data via Analysis of the Fourier Interspectrum Radial velocities) pipeline described in Galland et al. (2005). These measurements are presented in Appendix A.

We performed PT-MCMC fits to (1) the combination of direct imaging and RV data by itself, (2) the combination of direct imaging, RV, and both HIP-Gaia proper motion anomalies, (3) the combination of direct imaging, RV, and only the HIP proper motion anomaly (still using Gaia astrometry to correct for barycentric motion and fit for the system distance), and (4) the combination of direct imaging, RV, and only the Gaia proper motion anomaly. For most runs, we ran the PT-MCMC sampler described in Sect. 3 with 15 temperatures, 50 walkers, and 40 000 steps, and used the final 20 000 steps of the coldest temperature walkers for figures shown here. For the fit to just the direct imaging and RV points, the chains did not reach convergence within 40 000 steps, so we reran this fit with 100 walkers. For all MCMC runs, we inspected the walkers from the coldest temperature chain by eye to determine when the chain had converged.

Best parameters derived from each MCMC run are presented in Table 2. We adopt the DI+RV+HIPPARCOS fit as the best overall fit to the system. Corner plots for all parameters of this fit are presented in Fig. 3. The best fit orbit, as well as 100 orbits randomly selected from the posterior probability distribution function (pdf) are plotted alongside the RV and direct imaging data in Fig. 4. The best fit orbit and 100 randomly selected orbits from the posterior of the model tangential motion on the sky compared to the barycenter-corrected HIPPARCOS and DR2 proper motion anomalies are presented in Fig. 5. Similar plots for the other fits are presented in Appendix B for the fit to just direct imaging and RV data, Appendix C for the fit to direct imaging data, RV data, and both Gaia and HIPPARCOS proper motion anomalies, and Appendix D for the fit to direct imaging data, RV data, and only the Gaia proper motion anomaly. Unfortunately, the Gaia DR2 epoch caught the orbit in a period with extreme tangential motion – the averaging technique we used here to compare the tangential motion of the companion to the Gaia-HIP astrometry is not appropriate for this epoch, given the extreme motion on the sky, especially in RA. Forcing a fit of the high-precision Gaia measurement to the uncertain model astrometry over the full 1.5 yr Gaia DR2 observation period drives the period to larger values. Consequently, to compensate for the longer period from the Gaia constraint, the fit to direct imaging, RV, and Gaia-HIP astrometry produces an anomalously high mass for the primary. This dynamical mass is incompatible with both spectroscopic and photometric mass estimates for this star (Gaidos et al. 2014; Bonavita et al. 2021; Mann et al. 2019) – with an M1 spectral type, we do not expect a mass beyond 0.6 M_⊙ for HIP 113201A. Hence, we adopt the fit using only the HIPPARCOS constraint for our final values, as the HIPPARCOS epochs of observationscover a much more gradual, well-behaved portion of the orbital motion.

4.2 HIP 36985

HIP 36985 has extensive multi-instrument monitoring available for both direct imaging and RV observations. For our fits, we included the3 VLT-NaCo direct imaging observations presented in Table 3 of Baroch et al. (2021) and the VLT-SPHERE astrometric data points presented in Table 7 of Bonavita et al. (2021). HIP 36985 has multiple years (spanning 2014–2019) of RV monitoring using the HARPS instrument on the 3.6 m telescope at the La Silla observatory (Grandjean et al. 2020) and a significant proper motion anomaly between HIPPARCOS and Gaia DR2 measurements (Kervella et al. 2019). We again draw astrometry values from the HIPPARCOS-Gaia Catalog of Accelerations (Brandt 2018, 2021) to account for systematics on the sky and ensure all astrometry is in the Gaia DR2 reference frame. Additional RV monitoring with FEROS andCARMENES is presented in Baroch et al. (2021). We fit all HARPS, FEROS and CARMENES data points, taken from Table A.1 of Baroch et al. (2021). To account for the different RV offsets, we fit 3 offset terms (γ₁, γ₂, γ₃), one per instrument.

We performed two PT-MCMC fits to combinations of these data: (1) only direct imaging and RV points and (2) direct imaging, RV, and the HIP-Gaia proper motion anomalies. In both cases, we ran the PT-MCMC sampler for 15 temperatures, 50 walkers and 50 000 steps, and used the final 20 000 steps of the coldest temperature walker for figures shown here. We inspected 100 walkers from the coldest temperature chain by eye to determine when the chain had converged. The Gaia DR2 observationsfor HIP 36985 cover a less extreme part of the orbital motion than for HIP 113201; thus, we found that including both HIPPARCOS and Gaia proper motion anomalies in our fits improved the overall fit and produced physically consistent masses for both theprimary and the secondary.

Best parameters derived from each PT-MCMC run, as well as from the fit from Baroch et al. (2021), are presented inTable 3. We find comparable fits from fitting only the direct imaging and the RV data and from the fit to all direct imaging data, RV data, and HIPPARCOS-Gaia proper motion anomalies. Thus, we adopt the DI+RV+HIPPARCOS-Gaia fit as the best overall fit to the system. This fit matches well with the model parameters found using a similar method and the orvara fitting code from Baroch et al. (2021). This is unsurprising, as we fit nearly the same data points as Baroch et al. (2021), with a very similar method. Our fit additionally incorporates three SPHERE-SHINE observations, providing slightly stronger constraints on the orbital motion.

Corner plots for all parameters of the DI+RV+HIPPARCOS-Gaia fit are presented in Fig. 6. The best fit orbit, as well as 100 orbits randomly selected from the posterior pdf are plotted alongside the RV and direct imaging data in Fig. 7. The best fit orbit and 100 randomly selected orbits from the posterior of the model tangential motion on the sky compared to the barycenter-corrected HIPPARCOS and DR2 proper motion anomalies are presented in Fig. 8.

Fig. 3 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging, RV, and the HIPPARCOS proper motion anomaly in the orbital fit.

5 Discussion – Comparison of companion dynamical masses with existing mass estimates

Before full orbital fitting was possible for these objects, a number of authors estimated the mass of all components of these systems via model isochrones. Our dynamical mass measurements thus provide valuable tests both to photometric estimates of very low-mass star masses as well as to empirical mass-luminosity relationships from orbital fitting. Here we compare our dynamical masses to both model estimates based on photometry and other dynamical mass measurements in the literature.

Bonavita et al. (2021) estimate component masses using BT-Settl pre-main-sequence isochrones (Allard 2014) with SPHERE photometry and Gaia EDR3 distances. Their age estimates for HIP 113201 and HIP 36985 were derived from a mix of kinematic and gyrochronological age analysis. Adopting ages of 750 ± 250 Myr and 500 ± 200 Myr, respectively, Bonavita et al. (2021) find primary and secondary masses of 0.53 M_⊙ and 0.10 M_⊙ for HIP 113201AB, and 0.58 M_⊙ and 0.19 M_⊙ for HIP 36985AB.

Mann et al. (2019) derive an empirical M_K-mass relation for stellar masses between 0.075 and 0.7 M_⊙ by finding total system masses from 62 nearby M dwarf binaries. No K-band photometric measurements are available for HIP 36985B and HIP 113201B, but using the reddest value Δmag measurements from Bonavita et al. (2021) and assuming that the companions are uniformly redder than their primaries (i.e., that Δmag(K) will be less than the Δmag value we adopt from bluer SPHERE measurements), we can estimate minimum masses using the M_K-mass expression from Mann et al. (2019). For HIP 113201B, Bonavita et al. (2021) report 5 ΔJ measurements at different epochs ranging from 3.54 to 4.09 mag. We adopt the median value of these, and the standard deviation as the error: ΔJ = 3.93±0.31, yielding a minimum mass using the Mann et al. (2019) relationship of 0.099 ± 0.010 M_⊙. Adopting ΔH₃ = 2.848 ± 0.245 for HIP 36985B from Bonavita et al. (2021), we estimate a minimum mass using the Mann et al. (2019) relationship of 0.174 ± 0.021 M_⊙.

Combining existing HARPS RV measurements with SPHERE astrometry and the HIPPARCOS-Gaia proper motion anomaly measurement, Bonavita et al. (2021) find a dynamical mass for HIP 36985B of 0.180 ± 0.007 M_⊙. Using the orvara orbit-fitting code (Brandt et al. 2020), Baroch et al. (2021) fit HARPS, CARMENES, and FEROS RV measurements, VLT-NaCo direct imaging measurements and the HIPPARCOS-Gaia proper motion anomaly, finding dynamical masses of 0.554${}_{-0.049}^{+0.058}{M}_{\odot }$ and 0.1881${}_{-0.0047}^{+0.0048}{M}_{\odot }$ for the primary and secondary, respectively.

Our best fit dynamical masses for both primaries agree well with photometric estimates, the empirical M_K-mass relation from Mann et al. (2019), and also existing dynamical mass measurements in the literature. Both primaries are young main-sequence stars, and thus their model masses derived from their luminosity should be insensitive to the age of the star.

A comparison of dynamical masses and model mass estimates for HIP 113201B and HIP 36985B is presented in Fig. 9. All mass estimates and measurements for HIP 36985B agree on a mass of around 0.19 M_⊙. There is considerably more divergence between mass estimates for the lower-mass HIP 113201B companion. Estimates combining photometry and the system age with model isochrones find a mass of ~0.1 M_⊙; we find a considerably higher dynamical mass of 0.13–0.15 M_⊙. We only find a lower limit to mass with the relationship from Mann et al. (2019), but this is consistent with our dynamical mass estimate. However, the 750 ± 250 Myr age adopted for HIP 113201B in Bonavita et al. (2021) is considerably lower than the age of 1.2 ± 0.1 Gyr adopted here. To evaluate whether this is the source of the divergence between the dynamical mass and the model mass estimate from Bonavita et al. (2021), in Fig. 10, we overplot luminosity as a function of age using the Baraffe et al. (2015) models for objects with masses between 0.1 and 0.25 M_⊙ with the bolometric luminosities and age ranges of the two companions considered in this paper. We estimated bolometric luminosities for HIP 113201B and HIP 36985B using the J-band bolometric correction from Filippazzo et al. (2015), with the 2MASS J measurements for each primary stars and the Δ(mag) values in theSPHERE J2 band reported in Bonavita et al. (2021). We adopted a mid-M spectral type for each companion and used the M6 J-band bolometric correction from Filippazzo et al. (2015), as this was the earliest spectral type covered in this paper. The bolometric corrections from Filippazzo et al. (2015) do not change much between M6–M8 spectral types and we expect a roughly similar bolometric correction value for M4–M5 spectral types. At the younger ages reported for both stars in Bonavita et al. (2021) and the older age ranges adopted here, both companions have clearly reached the main sequence; therefore, the mass estimate based on luminosity will be largely age independent. We then firmly conclude that while all mass estimates for HIP 36985B agree well, HIP 113201B is anomalously faint given its observed dynamical impacts on its system. An undetected brown dwarf companion to HIP 113201B could be a natural explanation for this apparent discrepancy. Adding the luminosity from a 30 M_Jup brown dwarf to the luminosity of HIP 113201B as modeled by the 0.1 M_⊙ model track from Baraffe et al. (2015) would make almost no difference to the total luminosity. From the Baraffe et al. (2015) models, the Δmagnitude in H band at 1 Gyr between a 0.1 M_⊙ star and a 30 M_Jup companion would be 5.47 mag, so the total system brightness would change by less than 1% in H-band. Meanwhile, the presence of such an unseen companion would bring the total system mass up to 0.13 M_⊙, in line with the dynamical data. Such a companion at a close separation could easily elude detection in all of the existing data.

Fig. 4 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

Fig. 5 Model tangential velocity compared to HIPPARCOS and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. HIPPARCOS and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The HIPPARCOS mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which HIPPARCOS and Gaia measurements were acquired.

6 Conclusions

We present updated ages, orbital fits, and dynamical masses for HIP 113201AB and HIP 36985AB, two M1 + mid-M dwarf binary systems monitored as part of the SPHERE-SHINE survey (Desidera et al. 2021; Langlois et al. 2021; Vigan et al. 2021). The main results of this work are as follows:

• To robustly determine the age of both systems via gyrochronology, we undertook a photometric monitoring campaign for HIP 113201 and for GJ 282AB, the two wide K star companions to HIP 36985, using the 40 cm telescope ROAD located in the Atacama Desert. HIP 113201 is a slow rotator with a period of 19.6 ± 0.5 days, corresponding to a gyrochronological age of 1.2 ± 0.1 Gyr using the relationship from Delorme et al. (2011). HIP 36985 is a wide companion to the K star binary system GJ 282AB (Poveda et al. 2009). The M1 star HIP 36985 has a rotation period shorter than the period at the convergence time on the slow rotator branch for this spectral type, so only an upper limit on age for this star can be placed from gyrochronological relations such as those from Delorme et al. (2011). However, its K star companions have reached convergence; hence, we use the rotational periods for the K star components of the system to set the gyrochronological age of the system as a whole. Using the gyrochronology relationship from Delorme et al. (2011) and adopting an age for the Hyades of 625 Myr, the period of 12.10 ± 0.77 days for GJ 282B implies an age of 750 ± 100 Myr. This is compatible with the gyrochronological age limit of 650 Myr that we found for HIP 36985. Thus, combining this withthe age limits placed by lithium non-detection and activity indicators, we adopt a system age similar to that of the Hyades for GJ 282AB-HIP 36985. These system ages are sufficiently old that we expect that all components of these binaries have reached the main sequence.

• To derive dynamical masses for all components of the HIP 113201AB and HIP 36985AB systems, we used the PT-MCMC sampler implemented in EMCEE (Foreman-Mackey et al. 2013) to fit a combination of RV, direct imaging, and Gaia and HIPPARCOS astrometry.

• As epoch astrometry is not available yet for Gaia, we did not fit individual HIPPARCOS and Gaia data points but rather the proper motion anomaly (Kervella et al. 2019) at both the Gaia and HIPPARCOS epochs.

• For HIP 113201AB, combining direct imaging, RV, and both HIPPARCOS and Gaia proper motion anomalies yields an unphysically large primary mass of ~0.8 M_⊙ given the M1 spectral type of HIP 113201A, as well as a considerably longer period compared to fits of just the direct imaging and RV data. The proper motion anomaly method averages over 1.5–2 yr of measurement – we found in the case of HIP 113201AB that the Gaia proper motion anomaly did not serve as an appropriate constraint due to the extreme orbital motion of the companion in RA during the Gaia DR2 observation epochs. Fitting the  HIPPARCOS proper motion anomaly (which still relies on Gaia DR2 measurements to correct the barycenter motion) alongside the direct imaging and RV data yields a primary mass of ~0.6 M_⊙, marginally consistent with the M1 spectral type of the primary, and a secondary mass of 0.13 M_⊙. Fitting only the direct imaging and RV data yields a primary mass of 0.54 ± 0.03 M_⊙, fully consistent with its M1 spectral type, and a secondary mass of ~0.145 M_⊙. The secondary masses derived with and without including HIPPARCOS and Gaia data are all considerably more massive than the 0.1 M_⊙ estimated mass from the photometry of the companion (Bonavita et al. 2021). Thus, the dynamical impacts of this companion suggest that there is more mass in the system than expected from its photometry. An undetected brown dwarf companion to HIP 113201B could be a natural explanation for this apparent discrepancy. At an age >1 Gyr, a 30 M_Jup companion to HIP 113201B would make a negligible (<1%) contribution to the system luminosity but could have strong dynamical impacts.

• For HIP 36985AB, the dynamical masses found in this work for both the primary and secondary agree well with the photometric estimates of component masses, the masses estimated from the Mann et al. (2019) M_K – mass relationship, and previous dynamical masses in the literature (Baroch et al. 2021).

In the end, the two systems studied here are rare and not exemplary of M dwarf multiple systems. HIP 36985AB + GJ 282AB is an uncommon K+M dwarf quaternary system. In this paper we have shown preliminary evidence that HIP 113201AB may comprise an equally rare M dwarf star + brown dwarf triple system. Both systems present interesting challenges to formation mechanisms for low-mass stellar multiple systems.

Fig. 6 Corner plot for the PT-MCMC fit to the orbit of HIP 36985AB, incorporating direct imaging, RV, and both HIPPARCOS and Gaia proper motion anomalies in the orbital fit.

Fig. 7 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 36985AB. FEROS data points are plotted as red squares, CARMENES data points are plotted as green circles, and HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data, RV data, and both the HIPPARCOS and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: direct imaging data versus model comparison for HIP 36985AB. NaCo data points are plotted as orange circles, and SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data, RV data, and both the HIPPARCOS and Gaia proper motion anomalies is plotted as asolid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

Fig. 8 Model tangential velocity compared to HIPPARCOS and Gaia barycenter-corrected proper motion anomalies for HIP 36985AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. HIPPARCOS and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The HIPPARCOS mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which HIPPARCOS and Gaia measurements were acquired.

Fig. 9 Comparison of companion dynamical masses from this paper, dynamical mass from Baroch et al. (2021), and model estimate and dynamical masses from Bonavita et al. (2021). Mass values for HIP 36985B have been slightly offset from one other to improve readability. While all dynamical mass measurements and mass estimates agree well for HIP 36985B, thedynamical mass measurements and model mass estimate diverge considerably for HIP 113201B, which is extremely faint given the observed dynamical effects it produces.

Fig. 10 Best age and Lbol ranges for HIP 113201B (yellow rectangle) and HIP 36985B (cyan rectangle) overplotted on model evolutionary tracks for very low-mass stars from Baraffe et al. (2015). Both companions have reached the main sequence, and hence model mass estimates for these objects should be age independent for system ages greater than ~0.6 Gyr.

Acknowledgements

B.B. acknowledges funding by the UK Science and Technology Facilities Council (STFC) grant no. ST/M001229/1. B.B. would like to acknowledge valuable conversations with Trent Dupuy and Mariangela Bonavita which improved this manuscript. This work has been supported by the PRIN-INAF 2019 “Planetary systems at young ages (PLATEA)” and the ASI-INAF agreement n.2018-16-HH.0. A.Z. acknowledges support from the FONDECYT Iniciación en investigación project number 11190837. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF–Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). We also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). It has also been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). M.R.M., H.M.S., and S.D. are pleased to acknowledge this financial support of the SNSF. Finally, this work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon, also supported by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process.

Appendix A HARPS RV measurements for HIP 113201

Appendix B HIP 113201 orbital fits – Direct imaging and RV-only fit

Corner plots for all parameters of this fit are presented in Fig. B.1. The best fit orbit, as well as 100 orbits randomly selected from the posterior pdf are plotted alongside the RV and direct imaging data in Fig. B.2. The best fit orbit and 100 randomly selected orbits from the posterior of the model tangential motion on the sky compared to the barycenter-corrected Hipparcos and DR2 proper motion anomalies are presented in Fig. B.3.

Fig. B.1 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging data and RV data in the orbital fit.

Fig. B.2 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data and RV data is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: Direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data and RV data is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a blackstar symbol.

Fig. B.3 Model tangential velocity compared to Hipparcos and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit to the direct imaging data and RV data is plotted as a solid black line; blue lines depict 100 randomorbits taken from the final converged PT-MCMC posterior pdf. Hipparcos and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The Hipparcos mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which Hipparcos and Gaia measurements were acquired.

Appendix C HIP 113201 orbital fits – Direct imaging, RV, and Hipparcos -Gaia proper motion anomalies

Corner plots for all parameters of this fit are presented in Fig. C.1. The best fit orbit, as well as 100 orbits randomly selected from the posterior pdf are plotted alongside the RV and direct imaging data in Fig. C.2. The best fit orbit and 100 randomly selected orbits from the posterior of the model tangential motion on the sky compared to the barycenter-corrected Hipparcos and DR2 proper motion anomalies are presented in Fig. C.3.

Fig. C.1 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging, RV, and both Hipparcos and Gaia proper motion anomalies in the orbital fit.

Fig. C.2 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data, RV data, and both Hipparcos and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: Direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data, RV data, and both Hipparcos and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

Fig. C.3 Model tangential velocity compared to Hipparcos and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit to the direct imaging data, RV data, and both Hipparcos and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Hipparcos and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The Hipparcos mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which Hipparcos and Gaia measurements were acquired.

Appendix D HIP 113201 orbital fits – Direct imaging, RV, and Gaia proper motion anomaly

Corner plots for all parameters of this fit are presented in Fig. D.1. The best fit orbit, as well as 100 orbits randomly selected from the posterior pdf are plotted alongside the RV and direct imaging data in Fig. D.2. The best fit orbit and 100 randomly selected orbits from the posterior of the model tangential motion on the sky compared to the barycenter-corrected Hipparcos and DR2 proper motion anomalies are presented in Fig. D.3.

Fig. D.1 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging data, RV data, and the Gaia proper motion anomaly in the orbital fit.

Fig. D.2 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data, RV data, and Gaia proper motion anomaly is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: Direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data, RV data, and Gaia proper motion anomaly is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

Fig. D.3 Model tangential velocity compared to Hipparcos and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit to the direct imaging data, RV data, and Gaia proper motion anomaly is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Hipparcos and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The Hipparcos mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which Hipparcos and Gaia measurements were acquired.

References

All Tables

All Figures

Fig. 1 Summary of periodogram analysis for HIP 113201. Top-left panel: combined V + I-band magnitudes versus HJD (-2457670). Measurements on each night were taken in both the V band and the I band; the average of the two bands is plotted here. Top-middle panel: LS periodogram. The horizontal dotted lines show the 99.9 and 98% confidence levels, whereas the red dotted line shows the spectral window function. Top-right panel: clean periodogram. Bottom panel: light curve phased with the 19.6 d rotation period. Red circles are data points from the first half of the observations; blue circles are data points from the second half of the observations. The red dotted line is a sinusoidal fit to the phased light curve.

In the text

Fig. 2 Summary of periodogram analysis for GJ 282B. Top-left panel: V-band magnitudes versus HJD (-2459000). Top-middle panel: LS periodogram. The horizontal dotted lines show the 99.9 and 98% confidence levels, whereas the red dotted line is the spectral window function. Top-right panel: clean periodogram. Bottom panel: light curve phased with the 12.15 d rotation period. The red line is a sinusoidal fit to the phased light curve.

In the text

	Fig. 3 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging, RV, and the HIPPARCOS proper motion anomaly in the orbital fit.
In the text

Fig. 4 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

In the text

Fig. 5 Model tangential velocity compared to HIPPARCOS and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. HIPPARCOS and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The HIPPARCOS mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which HIPPARCOS and Gaia measurements were acquired.

In the text

	Fig. 6 Corner plot for the PT-MCMC fit to the orbit of HIP 36985AB, incorporating direct imaging, RV, and both HIPPARCOS and Gaia proper motion anomalies in the orbital fit.
In the text

Fig. 7 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 36985AB. FEROS data points are plotted as red squares, CARMENES data points are plotted as green circles, and HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data, RV data, and both the HIPPARCOS and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: direct imaging data versus model comparison for HIP 36985AB. NaCo data points are plotted as orange circles, and SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data, RV data, and both the HIPPARCOS and Gaia proper motion anomalies is plotted as asolid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

In the text

Fig. 8 Model tangential velocity compared to HIPPARCOS and Gaia barycenter-corrected proper motion anomalies for HIP 36985AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit from the direct imaging, RV, and the HIPPARCOS proper motion anomaly PT-MCMC run is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. HIPPARCOS and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The HIPPARCOS mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which HIPPARCOS and Gaia measurements were acquired.

In the text

Fig. 9 Comparison of companion dynamical masses from this paper, dynamical mass from Baroch et al. (2021), and model estimate and dynamical masses from Bonavita et al. (2021). Mass values for HIP 36985B have been slightly offset from one other to improve readability. While all dynamical mass measurements and mass estimates agree well for HIP 36985B, thedynamical mass measurements and model mass estimate diverge considerably for HIP 113201B, which is extremely faint given the observed dynamical effects it produces.

In the text

	Fig. 10 Best age and Lbol ranges for HIP 113201B (yellow rectangle) and HIP 36985B (cyan rectangle) overplotted on model evolutionary tracks for very low-mass stars from Baraffe et al. (2015). Both companions have reached the main sequence, and hence model mass estimates for these objects should be age independent for system ages greater than ~0.6 Gyr.
In the text

	Fig. B.1 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging data and RV data in the orbital fit.
In the text

Fig. B.2 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data and RV data is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: Direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data and RV data is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a blackstar symbol.

In the text

Fig. B.3 Model tangential velocity compared to Hipparcos and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit to the direct imaging data and RV data is plotted as a solid black line; blue lines depict 100 randomorbits taken from the final converged PT-MCMC posterior pdf. Hipparcos and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The Hipparcos mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which Hipparcos and Gaia measurements were acquired.

In the text

	Fig. C.1 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging, RV, and both Hipparcos and Gaia proper motion anomalies in the orbital fit.
In the text

Fig. C.2 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data, RV data, and both Hipparcos and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: Direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data, RV data, and both Hipparcos and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

In the text

Fig. C.3 Model tangential velocity compared to Hipparcos and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit to the direct imaging data, RV data, and both Hipparcos and Gaia proper motion anomalies is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Hipparcos and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The Hipparcos mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which Hipparcos and Gaia measurements were acquired.

In the text

	Fig. D.1 Corner plot for the PT-MCMC fit to the orbit of HIP 113201AB, incorporating direct imaging data, RV data, and the Gaia proper motion anomaly in the orbital fit.
In the text

Fig. D.2 Comparison of data to model orbits from the PT-MCMC posterior pdf. Left: RV data versus model comparison for HIP 113201AB. HARPS data points are plotted as cyan triangles. The best fit orbit to the direct imaging data, RV data, and Gaia proper motion anomaly is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Right: Direct imaging data versus model comparison for HIP 113201AB. SPHERE data points are plotted as lavender triangles. The best fit orbit to the direct imaging data, RV data, and Gaia proper motion anomaly is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. The position of the primary is depicted with a black star symbol.

In the text

Fig. D.3 Model tangential velocity compared to Hipparcos and Gaia barycenter-corrected proper motion anomalies for HIP 113201AB. The tangential velocity in RA is plotted in the left panel, and the tangential velocity in Dec is plotted in the right panel. The best fit orbit to the direct imaging data, RV data, and Gaia proper motion anomaly is plotted as a solid black line; blue lines depict 100 random orbits taken from the final converged PT-MCMC posterior pdf. Hipparcos and Gaia barycenter-corrected proper motion anomalies are plotted as green points for the same 100 random orbits; because the barycenter correction depends on primary and secondary mass, these vary slightly depending on the orbit selected. The Hipparcos mission lifetime for the 100 random orbits is highlighted in green; the Gaia DR2 observation period is highlighted in cyan. The small circle points depict the dates at which Hipparcos and Gaia measurements were acquired.

In the text