Dynamical masses of M-dwarf binaries in young moving groups: I - The case of TWA 22 and GJ 2060

Evolutionary models are widely used to infer the mass of stars, brown dwarfs, and giant planets. Their predictions are thought to be less reliable at young ages ($<$ 200 Myr) and in the low-mass regime ($\mathrm{<1~M_{\odot}}$). GJ 2060 AB and TWA 22 AB are two rare astrometric M-dwarf binaries respectively members of the AB Doradus and Beta Pictoris moving groups. As their dynamical mass can be measured within a few years, they can be used to calibrate the evolutionary tracks and set new constraints on the age of young moving groups. We find a total mass of $\mathrm{0.18\pm 0.02~M_\odot}$ for TWA 22. That mass is in good agreement with model predictions at the age of the Beta Pic moving group. We obtain a total mass of $\mathrm{1.09 \pm 0.10~M_{\odot}}$ for GJ 2060. We estimate a spectral type of M$1\pm0.5$, $\mathrm{L/L_{\odot}=-1.20\pm0.05}$ dex, and $\mathrm{T_{eff}=3700\pm100}$ K for GJ 2060 A. The B component is a M$3\pm0.5$ dwarf with $\mathrm{L/L_{\odot}=-1.63\pm0.05}$ dex and $\mathrm{T_{eff}=3400\pm100}$ K. The dynamical mass of GJ 2060 AB is inconsistent with the most recent models predictions (BCAH15, PARSEC) for an ABDor age in the range 50-150 Myr. It is 10 to 20\% (1-2 sigma, depending on the assumed age) above the models predictions, corresponding to an underestimation of $0.10$ to $0.20~\mathrm{M_\odot}$. Coevality suggests a young age for the system ($\sim$ 50 Myr) according to most evolutionary models. TWA 22 validates the predictions of recent evolutionary tracks at $\sim$20 Myr. On the other hand, we evidence a 1-2 sigma mismatch between the predicted and observed mass of GJ 2060 AB. This slight departure may indicate that one of the star hosts a tight companion. Alternatively, this would confirm the models tendency to underestimate the mass of young low-mass stars.


Introduction
Our understanding of stellar evolution has made a lot of progress since the introduction of the Hertzsprung-Russell diagram (HRD) a hundred years ago.The beginning of a star life, before it reaches the zero age main sequence, has been in particular deeply investigated through the development of evolutionary models.The latter rely on equations of state describing the stellar interior structure, and can make use of atmospheric models to define boundary conditions and predict emergent spectra.Different families of models exist (D'Antona & Mazzitelli 1997;Siess et al. 2000;Tognelli et al. 2012;Bressan et al. 2012;Feiden 2015; Baraffe et al. 2015), and their physical and chemical ingredients (e.g., nuclear rates, opacity, atmospheric parameters) have been updated in the recent years (e.g., Baraffe et al. 2015).The models can predict the age and mass of stellar and substellar objects from the measured broad band photometry, surface gravity, radius, luminosity, and effective temperature.The mass is the fundamental parameter which allows to comprehend the object nature and formation pathways.
The models predictions remain to be calibrated in various mass and age regimes (e.g., Hillenbrand & White 2004;Mathieu et al. 2007).Uncertainties related to the object formation process (formation mechanism, early accretion history, etc.) exist in the pre-main sequence (PMS) regime (e.g., Baraffe et al. 2002;Baraffe & Chabrier 2010).Further uncertainties may be added for low-mass stars, which have strong convection, rotation and magnetic activity (Mathieu et al. 2007).About 50 low-mass (below 1 M ) pre-main sequence stars had their mass determined thus far (e.g., Simon et al. 2000;Gennaro et al. 2012;Stassun et al. 2014;Mizuki et al. 2018).Most of these systems have been studied through their disk kinematics, and are thus younger than 10 million years (e.g., Guilloteau et al. 2014;Simon et al. 2017).Moreover, this method only allows to determine the total mass of the system, disk included.The disk mass can be a non negligible fraction of the total mass (e.g., Andrews et al. 2013, Fig. 9), so that uncertainties remain on the stellar mass.A dozen of the stars with dynamical mass are SB2 eclipsing binaries, for which the orbital inclination can be strongly constrained and the mass determined from the orbit.However, eclipsing binaries are very tight stellar pairs (orbital periods 1-10 d) so that each star strongly influences the other one (tides, high rotation speed, convection inhibition).Thus, their evolution may not be representative of typical stars (Chabrier et al. 2007;Kraus et al. 2011;Stassun et al. 2014).Consequently, evolution models remain poorly constrained for low-mass stars for most of the premain sequence stellar evolution.This can induce systematic offsets and disparate mass predictions (Hillenbrand & White 2004;Mathieu et al. 2007).
Some rare young (age<200 Myr) and nearby (d<100 pc) binaries resolved with high resolution imaging techniques (adaptive optics, speckle interferometry, lucky imaging, sparse aperture masking) have orbital periods shorter than a decade.Combined with a precise parallax, astrometric follow-up of the two components relative orbit gives the total dynamical mass of the system.Knowledge about the individual masses can then be gained from additional radial velocity measurements.These systems offer a good prospect for calibrating the PMS tracks and the underlying models physics.To date and to our knowledge, only 9 such systems in the intermediate PMS regime (10-100 Myr) have dynamical mass estimates below 1 M , with various model agreements: HD 98800 B (Boden et al. 2005), TWA 22 (Bonnefoy et al. 2009), HD 160934 (Azulay et al. 2014), AB-Dor (Azulay et al. 2015;Close et al. 2007), GJ 3305 (Montet et al. 2015), V343 Nor A (Nielsen et al. 2016), NLTT 33370 (Dupuy et al. 2016), GJ 2060 (this work) and GJ 1108 (Mizuki et al. 2018).We will here provide a refined dynamical mass for TWA 22 and a first determination for GJ 2060.
The calibration is nonetheless often limited by uncertainties on the age and distance of these benchmarking systems.These uncertainties are mitigated for systems belonging to known young nearby associations and moving groups (YMG).The age of the YMG can be inferred from several approaches (lithium depletion boundary, kinematics,...) and parallaxes can be measured for individual members (Gaia, Arenou et al. 2017, Hipparcos, Van Leeuwen 2007).Moreover, those systems share the same age  as the substellar companions resolved during direct imaging surveys (planets and brown dwarfs; e.g., Chauvin et al. 2004;Lagrange et al. 2010;Marois et al. 2008Marois et al. , 2010;;Rameau et al. 2013a,b) and whose mass determination also depends on PMS evolutionary models.
TWA 22 and GJ 2060 are two precious astrometric M-dwarf binaries with orbital period of a few years.They are proposed members of the young βPic and AB Dor moving groups, respectively.Both systems have a well-measured parallax.We initiated their follow-up in 2004 with various ground-based facilities in order to measure their dynamical masses and characterize their components.This paper presents an in-depth study of these two systems using published and additional observations and discusses the agreement between their orbit, their atmospheric properties, the age of their moving group and the PMS evolutionary models.We first review the observations and membership studies previously performed (section 2), and then present new imaging and spectroscopic data (section 3).We analyze the spectroscopic properties of GJ 2060 (section 4).We derive in section 5 the dynamical masses from orbital fits, and use them to probe the evolutionary models (section 6).The agreement between models and data is finally discussed in section 7.

Age and Membership of TWA 22 and GJ 2060
2. 1. TWA 22 TWA 22 (2MASS J10172689-5354265), located at d = 17.5 ± 0.2 pc (Teixeira et al. 2009), was originally proposed as a member of the ∼10 Myr old (Bell et al. 2015) TW Hydrae association (TWA) by Song et al. (2003).This classification was based on its strong Li 6708 Å absorption and Hα emission lines and sky position near other TWA members.A subsequent kinematic analysis of all TWA members proposed at the time by Mamajek (2005) revealed that the available kinematics of TWA 22 were largely inconsistent with the bulk of other TWA members and provided a low probability of membership.Possible membership in either TWA or the older β Pictoris moving group (∼ 25 Myr, Bell et al. 2015) was then proposed by Song et al. (2006).
TWA 22 was included in 2003 as a target in an adaptive optics (AO) imaging survey to search for low-mass companions (Chauvin et al. 2010).It was resolved into a ∼100 mas, equal luminosity binary and was as a potential benchmark target for dynamical mass measurements and model calibration.For this purpose, Teixeira et al. (2009) measured the parallax, provided revised proper motion and radial velocity measurements, and performed a detailed kinematic analysis of TWA 22 and found further evidence for membership in the β Pic group, but were unable to fully rule out TWA membership.Then, Bonnefoy et al. (2009) presented resolved spectra of the components, measured spectral types (later refined by Bonnefoy et al. 2014a to M5 ± 1 for TWA 22 A and M5.5 ± 1 for TWA 22 A), and performed an astrometric orbit fit to the available observations.This revealed that the total mass of the system was incompatible with model predictions, considering an age range consistent with the age of TWA.The authors however noted that the models may simply be underpredicting the system mass at such a young age.
TWA 22 has now been adopted as a bona fide member of the β Pic group on the basis of Bayesian methods to determine membership to kinematic moving groups (BANYAN I, Malo et al. 2013;BANYAN II, Gagné et al. 2014).TWA 22's kinematics were used to develop the β Pic group kinematic model implemented in the BANYAN Bayesian estimator (it has >99% probability of membership).The amount of lithium observed in TWA 22 is consistent with the age of the TWA association, but we know now that it is also compatible with its membership to the β Pic group, as Li may still subsist in the components at the age of β Pic.The age of the β Pic group has been revised multiple times in recent years using isochronal methods that rely on all group members (Malo et al. 2014;Bell et al. 2015), the lithium depletion boundary of the group (Binks & Jeffries 2014;Malo et al. 2014;Messina et al. 2016;Shkolnik et al. 2017), the rotation distribution of known members (Messina et al. 2016), and model comparisons to dynamical masses of binaries in the group (Montet et al. 2015;Nielsen et al. 2016).This large variety of age determination methods converge toward a group age of ∼25 Myr (see Table 1).
In this work we adopt the β Pic group age for TWA 22, provide new astrometric measurements of the binary components, combine these data with previous data to perform an updated orbital fit and measure the system mass, and compare the derived mass to estimates from the latest stellar evolution models.The binary period is relatively short (∼ 5 yr) and TWA 22 has been regularly observed from 2004 to 2007 and later in 2013 and

GJ 2060
GJ 2060 (2MASS J07285137-3014490) is an early M-dwarf at d = 15.69 ± 0.45 pc (Van Leeuwen 2007) that was first identified as a small separation binary by the Hipparcos satellite (Dommanget & Nys 2000).The star was subsequently identified as a nearby young star in the paper presenting the discovery of the AB Doradus moving group (Zuckerman et al. 2004).This work presented GJ 2060 and ∼30 other stars as having both Galactic kinematics consistent with the well studied young system AB Dor and independent indicators of youth (X-ray and H-alpha emission, large vsin i, etc.).Along with AB Dor itself and six other nearby stars within a ∼5 pc radius, GJ 2060 is a member of the AB Dor moving group nucleus.The system has since been verified as a bona fide member of the AB Dor moving group using revised group kinematic distributions and Bayesian methods with an estimated membership probability of >99% (Malo et al. 2013;Gagné et al. 2014).GJ 2060 was first resolved into an 0.175" multiple system by Daemgen et al. (2007) using adaptative optics imaging.The system has been observed multiple times since with high resolution imaging and exhibited rapid orbital motion (see Janson et al. 2014).The age of the AB Dor moving group, and thereby GJ 2060, was first proposed to be ∼50 Myr by Zuckerman et al. (2004).Yet, the age of the group remains relatively poorly constrained and ages ranging from the original ∼50 Myr to ∼150 Myr have been proposed over the last decade (e.g., Close et al. 2005;Nielsen et al. 2005;Luhman et al. 2005;Lopez-Santiago et al. 2006;Ortega et al. 2007;Torres et al. 2008;Barenfeld et al. 2013;Bell et al. 2015).The system components of the groups' namesake quadruple system AB Dor have been studied in detail (Close et al. 2005;Nielsen et al. 2005;Guirado et al. 2011;Azulay et al. 2015) and comparisons to stellar evolution models indicate discrepancies between the measured masses of the components and point toward ages <100 Myr.This is in conflict with group ages estimated both from the individual components of AB Dor and from the ensemble of stars using HR diagram placement (Luhman et al. 2005;Bell et al. 2015), rotation periods (Messina et al. 2010), and Li depletion (Barenfeld et al. 2013).These works find that individual members and the ensemble of proposed AB Dor members have properties consistent with the Pleiades open cluster and likely have a comparable age (∼120 Myr; Stauffer et al. 1998;Barrado y Navascués et al. 2004;Dahm 2015).Here we use new astrometric and radial velocity measurements of the GJ 2060 system to derive component masses and perform similar comparisons to stellar evolution models.No orbital fit had been performed on this system yet, so that its orbital elements and dynamical mass are first determined in the present article.

Observation and data processing
A summary of the new observations of TWA 22 and GJ 2060 is given in Table 3.We describe the datasets and related reduction processes in more details below.).The S13 camera was associated to the H-band filter (λ c = 1.66 µm, ∆λ = 0.33 µm), yielding a square field of view of 13.5 arcsec.The wavefront sensing was achieved in the near-infrared on the pair (seen as a whole).We acquired 32 frames (NEXPO) of the binaries consisting of 0.345 s × 30 (DIT × NDIT ) averaged exposures each.Small (±3") dithers were applied every four frames to allow for an efficient sky and bias subtraction at the data processing step.We observed immediately after TWA 22 AB the M6 star GSC08612-01565 to calibrate the point-spread function (PSF) of the instrument using the same adaptive-optics setup and the same DIT, NDIT, and NEXPO as for TWA 22 AB.We observed the following night the crowded field of stars around Θ Ori C to calibrate the platescale and field orientation with the same filter and camera and the visible wavefront sensor.That astrometric field was already used in Bonnefoy et al. (2009) for the previous observations of TWA 22.All the data were reduced with the eclipse sofware (Devillard 1997).The eclipse routines carried out the basics cosmetic steps: bad pixel flagging and interpolation, flat field calibration, sky subtraction, and cross-correlation and shift of the dithered frames.We extracted the position of the Θ Ori stars and compared them to those reported in McCaughrean & Stauffer   (1994) to infer a platescale of 13.19 ± 0.08 mas/pixel and a True North of −0.90 ± 0.15 • for those observations.We used a deconvolution algorithm dedicated to the stellar field blurred by the adaptive-optics corrected point spread functions to deblend the overlaping point-spread-functions of TWA22 A and B in the final NaCo image (Veran & Rigaut 1998) and measure the position and the photometry of each components.The same tool was used in Bonnefoy et al. (2009).The algorithm is based on the minimization in the Fourier domain of a regularized least square objective function using the Levenberg-Marquardt method.It is well suited to our data which are Nyquist-Sampled.We crosschecked our results using the IDL Starfinder PSF fitting pack-age (Diolaiti et al. 2000) which implements a custom version of the CLEAN algorithm to build a flux distribution model of the binary but does not perform any spatial deconvolution.We find a contrast ∆H = 0.52 ± 0.05 mag consistent with the values derived at previous epochs (Bonnefoy et al. 2009).The binary is found at a PA = 1.15 ± 0.15 • and separation ρ = 100 ± 3 mas.

SPHERE observations
The binary was observed on February 3, 2015 as part of the SHINE (SpHere INfrared survey for Exoplanets) survey (Chauvin et al. 2017) with the high-contrast instrument SPHERE at UT3/VLT (Beuzit et al. 2008).The observations were scheduled as part of a sub-program (filler) of SHINE devoted to the astrometric monitoring of tight binaries.SPHERE was operated in field-tracking mode.No coronagraph was inserted into the light path.The IRDIFS_EXT mode enabled for simultaneous observations with the dual-band imaging sub-instrument IRDIS (Dohlen et al. 2008;Vigan et al. 2010) in the K1 (λ c = 2.110 μm; ∆λ = 0.102 μm) and K2 (λ c = 2.251 μm; ∆λ = 0.109 μm) filters in parallel with the lenslet-based integral field spectrograph (IFS, Claudi et al. 2008;Mesa et al. 2015) in the Y to H band (0.96 − 1.64 μm).Only the IRDIS data were exploited because the low-resolution (R∼30) IFS observations are superseeded by the SINFONI spectra (R∼ 1500−2000) of the binary exploited in Bonnefoy et al. (2009) and Bonnefoy et al. (2014a).
We acquired 240×4s IRDIS frames of the binary.The IRDIS dataset was reduced at the SPHERE Data Center1 (DC) using the SPHERE Data Reduction and Handling (DRH) automated pipeline (Pavlov et al. 2008;Delorme et al. 2017).The DC carried out the basic corrections for bad pixels, dark current, and flat field.It also includes correction for the instrument distortion (Maire et al. 2016a).
The wavefront-sensing of the adaptive optics system SAXO (Fusco et al. 2006;Petit et al. 2014) could operate on the target in spite of its faintness at optical wavelengths (V=13.8mag; Zacharias et al. 2005) and of the adverse observing conditions (Table 3).But the tip-tilt mirror occasionally produced strong undesired offset of TWA 22 in the field of view and part of the sequence was affected by low Strehl ratio.We then selected by eye 71 frames with the best angular resolution.We measured the relative position of the binary in the remaining frames using a custom cross-correlation IDL script.The frames were then realigned using sub-pixel shifts with a tanh interpolation kernel.The registered frames were averaged to produce a final frame using the Specal pipeline (Galicher et al., in prep).
TWA 22 A and B are well resolved into the final K1 and K2 images (see Fig. 1).We did not observe any reference star to calibrate the point-spread-function so that we could not use deconvolution algorithm for that epoch.But the high Strehls of the SPHERE observations mitigate the cross-contamination of the binary components.We measured their position in the K1 image (offering the best angular resolution) fitting a Moffat function within an aperture mask (4 pixel radii) centered on the guessed position of the stars.We varied the aperture size (±1 pixels in radius) and considered alternative fitting function (Gaussian, Lorentzian) to estimate an error on the astrometry.We used a True North value of 1.72 ± 0.06 • and a platescale of 12.267 ± 0.009 mas/pixel derived from the observations of Θ Orionis C as part of the long term analysis of the SHINE astrometric calibration (same field as the one observed with NaCo; Maire et al. 2016a,b).This leads to a position angle PA = 114.90± 0.10 • and a separation ρ = 103 ± 1 mas between the two components of TWA 22.

NaCo data
We observed GJ 2060 with NaCo (Program 090.C-0698; PI Delorme) in the H-band in the course of a direct imaging survey of M-dwarfs (Delorme et al. 2012;Lannier et al. 2016).The observations were performed in field tracking mode with the de- tector cube mode enabling for short integration time (0.15s).We also observe the astrometric calibrator Θ Ori C with the same setup.The data were all reduced with the eclipse tool.We find a platescale of 13.19 ± 0.06 mas/pixel and a True North value of −0.60 ± 0.33 • for those observations.GJ 2060 AB was tight (69 mas) in the images.This required the use of a deconvolution algorithm to deblend the binary components.We reduced for that purpose the data of GJ 3305 observed the same night (Table 3).GJ 3305 is itself a tight pair of M-dwarfs and is a member of the Beta Pictoris moving group (Zuckerman et al. 2001).The separation of GJ 3305 in November 2012 (290 mas) and the Strehl ratio were sufficient to mitigate the self-contamination of the binary component.We could then extract a subfield centered on GJ 3305 A that could serve as a reference PSF.We nonetheless used in addition three isolated bright stars from the θ Orionis field observed four nights before GJ 2060 at a close airmass as to evaluate the dependency of the results related to the PSF choice.The deconvolution algorithm of Veran & Rigaut (1998) yields a PA = 232.2± 2.3 • and ρ = 69 ± 2 mas for GJ 2060 AB.This measurement is confirmed by the Starfinder tool.
Both sequences contain two sets of frames which correspond each to a position of the source on the detector.We averaged each set of three frames to produce two resulting frames.Those resulting frames were then used to subtract the sky and bias contributions into the 6 original frames.We registered the frames on a common origin, applied a rotation to re-align them to the North, and averaged them to produce the final frames.The last step enabled to filter out part of the bad pixels.
We fitted a Moffat function on each star flux distribution to retrieve their relative position.For both epochs, we considered the platescale (9.971 ± 0.005 mas/pixel) and the absolute orientation on the sky (0.262 ± 0.022 • ) reported in Service et al. (2016).The estimated contrasts and astrometry are reported in Tables 5 and 6, respectively.

SPHERE observations
The binary was observed as part of the SHINE survey on February 2015 and February 2017 in field and pupil tracking mode, respectively.For both nights, the IRDIFS_EXT mode was used.
The reduction of the IFS data was performed following the procedure described in Mesa et al. (2015) and Zurlo et al. (2014).The calibrated spectral datacubes are made of 39 narrow band images.We rotated the datacubes corresponding to each exposures to align them to the North and averaged them.We extracted from the resulting cube the flux ratio between each component of the binary for both epochs (74 mas wide circular aperture).
We made use of the IRDIS data for the astrometric monitoring.We followed the procedure described in Section 3.1.2to reduce those data.We used the True North and platescale values reported in Section 3.1.2for the 2015 observations.We adopted a true North on sky of 1.702 ± 0.058 • and a platescale of 12.250 ± 0.009 mas/pix from the observations of NGC3603 obtained on February 7, 2017 as part of the long-term astrometric calibration of the instrument (Maire et al. 2016b).The binary position in the final images was measured with a Moffat function and is reported in Table 6.Fig. 1 displays the 2017 epoch.

AstraLux observations
Three of the AstraLux data points presented here are previously unpublished.These were obtained as a continuation of the As-traLux orbital monitoring campaign for young M-dwarf binaries, with a particular focus on young moving group members (Janson et al. 2014(Janson et al. , 2017)).The new data were acquired in March and December of 2015 with the lucky imaging camera AstraLux Sur (Hippler et al. 2009) at the ESO NTT telescope (programs 094.D-0609(A) and 096.C-0243(B)).They were reduced in an identical way as previously in the survey (e.g., Janson et al. 2014).For the March run, the cluster NGC 3603 was used as astrometric calibrator, giving a pixel scale of 15.23 mas/pixel and a true North angle of 2.9 • .In the December run, the Trapezium cluster was used for astrometric calibration, yielding a pixel scale of 15.20 mas/pixel, and a true North angle of 2.4  (Bonnet et al. 2003) to the integral field spectrograph SPIFFI (Eisenhauer et al. 2003) operating in the near-infrared (1.10-2.45μm).SPIFFI slices the field of view into 32 horizontal slitlets that sample the horizontal spatial direction and rearranges them to form a pseudo long slit.That pseudo-slit is dispersed by the grating on the 2048 × 2048 SPIFFI detector.GJ 2060 A was bright enough at R band to allow for an efficient adaptive optics correction.We used the pre-optics providing 12.5 mas × 25 mas rectangular spaxels on sky and a square field of view of 0.8" side.The target was observed during two consecutive sequences with the J and H+K gratings, covering the 1.10 − 1.40 and 1.45 − 2.45 μm wavelength range at R∼2000 and 1500 resolving powers, respectively.We obtained 11 frames with the binary in the field of view.In-between each frame, the binary was dithered to increase the final field of view and filter out residual non-linear and hot pixels.We also obtained an exposure on the sky at the end of each sequence to efficiently subtract the sky emission lines, detector bias, and residual detector defects.The observatory obtained observations of HIP 036092 immediately after GJ 2060.HIP 036092 is a B8V star that was used to evaluate and remove the telluric absorption lines.
We used the version 3.0.0 of the ESO data handling pipeline (Abuter et al. 2006) through the workflow engine Reflex (Freudling et al. 2013) which allowed for an end-to-end automatized reduction.Reflex performed the usual cosmetic steps on the bi-dimensional raw frames (flat field removal, bad-pixel flagging and interpolation).Those steps rely on calibration frames taken the days following our observations.The distortion and wavelength scale were calibrated on the entire detector.The positions of the slitlets on the detector were measured and used to build the datacubes containing the spatial (X,Y) and spectral dimensions (Z).In the final step, the cubes corresponding to individual exposures were merged into a master cube.
GJ 2060 is well resolved into the J and H+K mastercubes but the sources contaminate each other.We applied the CLEAN3D tool described in Bonnefoy et al. (2017) to deblend the sources at each wavelength.The PSF at each wavelength is built from the duplication of the profile of GJ 2060 A following a PA=0 • .The tool produced two datacubes where one of the two components of the system is removed.We extracted the J and H+K band spectra of each component integrating the flux within circular apertures of radius 147 and 110 mas at each wavelength in the datacubes, respectively.We extracted the telluric standard star spectrum using the same aperture sizes and corrected its continuum with a 12120 K black body (Theodossiou & Danezis 1991).The hydrogen and helium lines were interpolated using a third order Legendre polynomial.GJ 2060 A and B spectra could then be divided by the telluric standard star spectrum to correct for atmospheric absorptions.
We computed 2MASS J, H and K-band contrasts as well as the K1 and K2 SPHERE contrasts from GJ 2060 A and B spectra prior to the telluric line correction (Table 5).The H, K1 and K2 synthetic contrasts match those derived from the SPHERE and NaCo data within error bars.We therefore used the synthetic 2MASS H and K-band contrasts and the 2MASS magnitude of the system (Cutri et al. 2003) to retrieve the individual magni- tudes of GJ 2060 A and B. J-band contrasts could be extracted from the SPHERE IFS data.They agree with the one derived with SINFONI.We used the contrast value of the 2015 SPHERE data to derive the J-band magnitude of the system components.
The 2MASS J magnitudes could then be used to fluxcalibrate the J-band spectra using the 2MASS filter response curves and tabulated zero points2 .We used the K1 magnitude measured with VLT/SPHERE and a spectrum of Vega (Mountain et al. 1985;Hayes 1985) to flux-calibrate the H+K spectra.

HARPS data
High S/N spectra have been acquired with HARPS (Mayor et al. 2003): 1 night in April 2014 (JDB = 2456774.493808)and 5 nights in October 2016 (between JDB = 2457666.881702and 2457671.850830).Each spectrum contains 72 spectral orders, covering the spectral window [3800 Å, 6900Å].The spectral resolution is approximately 100 000.The S/N of the spectra is ≈ 100 at 550 nm.The number of spectra per night is 2 (consecutive), except for the first night, for which only one has been taken.The data as provided by HARPS's Data Reduction Software 3.5 (DRS) were first processed with SAFIR, a home-built tool that uses the Fourier interspectrum method described in Chelli (2000) and in Galland et al. (2005) to measure radial velocities of stars with high v sin i. SAFIR also estimates other observables, such as the cross-correlation functions, as defined in Queloz et al. (2001), and the bisector velocity spans (BVS), R0HK indexes, etc.For a detailed description of SAFIR, see Galland et al. (2005).
The values obtained in October 2016 show a very strong dispersion, probably due to the high magnetic activity of the star.Indeed, the orbit of the binary is ∼ 8 years long, so that we do not expect the radial velocities to vary more than ∼ 0.01 km/s within a few consecutive days, very different from the 0.40 km/s variation we observed.Moreover, we note a strong correlation between the star bisector and the radial velocity measurements.We will therefore add this noise to the instrument uncertainty.The jitter evidenced in the HARPS data (section 3.2.6)must be taken into account in the FEROS set.Thus, we combined quadratically this estimated activity-related noise (0.40 km/s) to each FEROS uncertainty.

Spectrophotometric analysis
We compared the SINFONI spectra of GJ2060 A and B to the medium-resolution (R∼2000) spectra of K and M-dwarfs from the IRTF library (Cushing et al. 2005;Rayner et al. 2009).The 1.1-2.5 μm spectral slopes of GJ 2060 A is best reproduced by the Gl 229 A spectrum (M1V; Fig. 2).The detailed absorptions and slopes of the J-band and K-band spectra are also reproduced Fig. 2: 1.1-2.45µm SINFONI spectra of GJ 2060 A and B renormalized at 1.55 µm.
by that template (Fig. 13).The lack of water band absorption from 1.3 to 1.4 µm in the spectrum of GJ 2060 A confirms that the object has a spectral type earlier than M2.The M0.5 and M1.5 dwarfs Gl 846 and Gl 205 fit equally well the K and H band spectrum of GJ 2060 A, respectively.Therefore, we estimate that GJ 2060 A is a M1±0.5 dwarf.
The spectral slope of GJ 2060 B is reproduced by the spectrum of the M3 dwarf Gl 388.The comparison at J-band evidences departures from 1.1 to 1.2 μm and 1.24 to 1.33 μm between our object spectrum and the templates (Fig. 13).These departures are also evidenced in the SINFONI spectra of GJ 3305A and B obtained as part of our observation program (Durkan et al. 2018).It likely arises from the SINFONI instrument.The multiple atomic lines (K I, Na I, Fe I, Al I) and the water band absorption from 1.3 to 1.35 μm indicate that the object has a spectral type later than M2.The H band spectrum is best represented by the one of the M3.5 dwarf Gl 273 while the K-band is perfectly reproduced by the spectrum of the M3 template (Fig. 13).We conclude that GJ 2060 B is a M3±0.5 dwarf.
Those spectral types confirm the estimates made in Bergfors et al. (2010) from the optical colors.We used them together with the bolometric corrections of Pecaut & Mamajek (2013) and the J-band magnitude (Tab.5) of each component to infer a log(L/L ) = −1.20 ± 0.05 dex and log(L/L ) = −1.63 ± 0.05 dex for GJ 2060 A and B, respectively.
We performed a χ 2 comparison of GJ 2060 A and B spectra to a grid of BT-SETTL atmosphere models (Baraffe et al. 2015) and show the best fitting solutions in Fig. 3.The grid covers 1500 ≤ T eff (K) ≤ 5500 (in steps of 100 K), 2.5 ≤ log g(dex) ≤ 5.5 (in steps of 0.5 dex) and considers solar abundances.We find T eff = 3700 ± 100 K and log g>4.0 dex for GJ 2060 A. Similarly, we find T eff = 3400 ± 100 K and log g≥3.5 dex for GJ 2060 B. The T eff are in good agreement with the estimates ( T eff = 3615 − 3775 K for A and T eff = 3300 − 3475 K for B) derived from Table 5 of Pecaut & Mamajek (2013) for the estimated spectral type of the binary components.As these results rely on atmosphere models, they do not depend on the system age.Both the T eff and bolometric luminosities are used as input of evolutionary models for the calibration of their mass predictions in Section 6.

Orbital fit and dynamical mass
Both systems orbits have been observed on several occasions covering a timespan longer than their periods, so that we are now able to derive precise estimates of their orbital elements.In both cases, we fit the relative orbit of the B component with respect to the A component, assuming a Keplerian orbit projected on the plane of the sky.In this formalism, the astrometric position of the companion can be written as: where Ω is the longitude of the ascending node (measured counterclockwise from North), ω is the argument of the periastron, i is the inclination, θ is the true anomaly, and r = a(1 − e 2 )/(1 + e cos θ) is the radius, where a stands for the semi-major axis and e for the eccentricity.
The orbital fit we performed uses the observed astrometries depicted in Table 4 and 6 to derive probability distributions for elements a, P (period), e, i, Ω, ω, and time for periastron passage t p .Elements a and P are probed separately, so that we can deduce the probability distribution of the total mass as a by-product, as a function of the distance of the star d.
We used two complementary fitting methods, as described in Chauvin et al. (2012), (i) a least squares Levenberg-Marquardt (LSLM) algorithm to search for the model with the minimal reduced χ 2 , and (ii) a more robust statistical approach using the Markov-Chain Monte Carlo (MCMC) Bayesian analysis technique (Ford 2005(Ford , 2006) ) to probe the distribution of the orbital elements.Ten chains of orbital solutions were conducted in parallel, and we used the Gelman-Rubin statistics as convergence criterion (see the details in Ford 2006).We picked randomly a sample of 500,000 orbits into those chains following the convergence.This sample is assumed to be representative of the probability (posterior) distribution of the orbital elements, for the given priors.We chose the priors to be uniform in x = (ln a, ln P, e, cos i, Ω + ω, ω − Ω, t p ) following Ford (2006).For any orbital solution, the couples (ω,Ω) and (ω + π,Ω + π) yield the same astrometric data, this is why the algorithm fits Ω + ω and ω − Ω, that are not affected by this degeneracy.The system distance has to be given to the algorithm.No input on the mass    is needed, as it can be derived directly from a and P by Kepler's third law.The resulting MCMC distributions are well peaked when the data sample adequately the orbits, as is the case in this study.The complete set of posterior distributions and correlations are given in the appendix.

TWA 22
Bonnefoy et al. ( 2009) already performed an orbital fit of TWA 22 based on astrometric data from 2004 to 2007.The data covered at that time about three-quarters of a period.The authors used a pure Levenberg-Marquardt algorithm, that finds local minima and estimates the uncertainties from the resulting covariance matrix.We intend here to improve the orbital fit by using the new astrometric data (two periods are now covered) and a refined algorithm, described above, that allow a fine sampling of the phase parameters and a robust determination of the probability distributions.
The astrometric measurements gathered with NaCo on the system are particularly homogeneous and sample well the orbit.Therefore, we excluded the SPHERE point from the fit at first in order to avoid the possible bias associated with the change of instrument.We then checked the agreement between the results and the SPHERE point afterwards.
The MCMC algorithm gives an estimate of the orbital elements (see Table 8), with a precision of 0.02 on the eccentricity, 0.05 au on the semi-major axis, 0.04 yr for the period or 6 • on the inclination (see appendix A).The portrayed orbit has a low eccentricity (∼ 0.1) and inclination (∼ 22 • ), as can be hinted from its on-sky representation in figure 4.This figure shows the best fit obtained with the LSLM algorithm together with a hundred orbits picked up randomly within the 500,000 total sample used to derive the posteriors.The orbital elements derived by Bonnefoy et al. (2009) are all retrieved within 1σ.
The total system mass was computed from the semi-major axis and period corresponding to each orbit explored by the MCMC chains.For any distance d, we find a resulting total mass of m tot = a 3 /T 2 = 0.179 ± 0.018 M d 17.5 pc 3 .Using the parallax distance and propagating its uncertainty, we finally obtain a dynamical mass of m tot = 0.18 ± 0.02 M for the pair.We checked the consistency between the fitted orbit and the SPHERE point that we did not consider: the astrometry falls within the 68% of confidence interval of the orbital fit, between 0.4 and 0.9 σ from the probability peak (2-3 • in position angle, 2-3 mas in radius).Running the algorithm with this extra data point give very similar orbital elements (all well inside the 68% confidence interval).It yields to the same dynamical mass, with smaller error-bars (0.18 ± 0.01 M ).
A dynamical mass of m tot = 0.21 ± 0.02 M was obtained in Bonnefoy et al. (2009) with less than 4 years coverage via a LSLM algorithm.This value is close to the one we obtain, but our peak value is outside the 1 σ confidence interval.However, the error bars on the orbital elements in Bonnefoy et al. (2009) may be slightly underestimated, as they are roughly estimated from the covariance matrix.The present determination should therefore be more robust.

GJ 2060
Radial velocity measurements (RVs) from HARPS and FEROS (Durkan et al. 2018) help here refining the orbital fit.The binary is not resolved by the spectrometers (SB1).We only considered the FEROS data to get homogeneous measurements.This is le- gitimate, as taking into account HARPS data would not bring significant constraints.Indeed, HARPS data come down to two epochs (April 2014 and October 2016) that are close to FEROS epochs (see Fig. 6), and we have to fit an additional RV offset if we want to include data from another instrument.The code we use is a slightly modified version of the code used for TWA 22, similar to the code used in Bonnefoy et al. (2014b) for β Pic b.In addition to the orbital elements, it evaluates the probability distributions of the offset velocity v 0 and amplitude K of the radial velocity, with a prior uniform in (v 0 , ln K) assumed for these extra variables (Ford 2006).In the formalism described previously, assuming a Keplerian orbit, the radial velocity is If the binary is a pure SB1, the amplitude derives from the fractional secondary mass m B /m tot as The introduction of the radial velocity breaks the degeneracy of the couple (Ω, ω) and unique values can thus be derived for these two variables.
The astrometric data are more numerous than in the case of TWA 22, but less homogeneous.Therefore, small systematic errors may bias the orbital fit (Table 6).Such errors are discussed in section 6.1.These 15 years of data cover approximately twice the relative orbit, but the passages near the periastron are not very well constrained and suggest a very quick displacement in that zone, hinting for a high eccentricity.The results of the MCMC algorithm are displayed in Table 9.The distribution of orbital elements are very peaked, especially the one on the eccentricity (see appendix B).Indeed, we obtain a precision of 0.01 on the eccentricity, 0.04 au on the semi-major axis, 0.04 yr on the period and 3 • on the inclination.Noticeably, the eccentric distribution peaks at e = 0.89, but does not extend up to e = 1: the components are bound.This orbital elements, and in particular the eccentricity, are very robust, and we obtain the same constraint when we fit only the astrometry.A hundred orbits, selected randomly within the 500,000 orbits used in to derive the posteriors, are plotted on Fig. 5. Fig. 6 displays the radial velocity data.The portrayed orbit confirms the very high < 1 eccentricity.for the pair.The fractional mass could be computed for each orbit thanks to the fit of the radial velocity amplitude (Eq.3).Considering our system as a pure SB1 (Eq.4), we obtain a fractional mass of m B /m tot = 0.26 ± 0.10 15.69 pc d .However, this naive approach is questionable giving the flux ratio of the two components at FEROS wavelengths (∼ 0.25).Thus, we used the method proposed by Montet et al. (2015) and considered our RVs to be the flux-weighted sum of the two individual RVs.The amplitude K fitted by the orbital fit could then be written as

FEROS HARPS
Fig. 6: Plots of a hundred radial velocity evolution obtained with the MCMC algorithm for system GJ 2060.Radial velocity measurements are displayed in their instrument color, along with their position on the fit.Only FEROS data are used in the orbital fit.The orbit in black, obtained with the LSLM algorithm, corresponds to the lower χ 2 .
where The uncertainty on the total mass mainly comes from the uncertainty on the distance d, as ∆m tot /m tot = 3∆d/d.d and ∆d derive from the parallax released within the new reduction of Hipparcos data (Van Leeuwen 2007).The binarity of the system was taken into account in Hipparcos reduction through two additional variables in the proper motion fit.Moreover, the high eccentricity of the orbit prevents a good sampling of the radial velocity, in particular close to the periastron passage.That leads to high error bars in the inclination and velocity amplitude K.These errors propagate on the fractional secondary mass (see equation 4).A higher accuracy on the orbital elements determination would certainly be achieved if the periastron passage was sampled in the available astrometric and spectroscopic data.This is unfortunately not the case yet.
On the other hand, the uncertainty on the fractional mass mainly comes from the very low constraints provided by the RVs.During most of the orbital revolution, the RV variation has a similar magnitude than the noise evidenced by the HARPS measurements.Only the sampling of the periastron passage could provide meaningful points that can refine the fractional mass.The next passage corresponds to October 2020.We can also notice that the difference is very significant between the naive (SB1) and corrected (flux-weighted) approach: the fractional mass nearly doubles.Averaging the flux-weighted RVs is a first order method, and is probably not precise enough to dis-entangle the two lines in our case, where the luminosity of the secondary is not negligible compared to the primary.A more refined method (e.g., Czekala et al. 2017) would be necessary to trace back the individual RVs from our measurements and compute a robust fractional mass.Thus, we use only the total mass in the next sections.

Comparison to the models
Both our systems now have a dynamical mass and an estimated age given by their membership to moving groups, as well as a robust estimate of their bolometric luminosities L and effective temperatures T eff .Thus, we are able to probe the accuracy of the PMS evolutionary models at these mass ranges.
There exist several evolutionary models for PMS stars, relying on different physics (e.g., atmospheric models, convection efficiency).Two of them are suitable for 0.1 M objects, the DM97 model (D'Antona & Mazzitelli 1997) and the BHAC15 model (Baraffe et al. 2015).Four more models becomes suitable for higher-mass PMS stars: the SDF00 model (Siess et al. 2000), the PISA model (Tognelli et al. 2011(Tognelli et al. , 2012)), the PARSEC model (Bressan et al. 2012) with Chen et al. ( 2014)'s corrections for low-mass stars, and the Darmouth model (Dotter et al. 2008) with Feiden (2016)'s integrations of magnetic field.Testing the predictions of different models enables to compare the relevance of their approach, and thus to achieve a better understanding of the underlying physics.

TWA 22
According to the previous sections, TWA 22 has a total dynamical mass of 0.18 ± 0.02 M and an age of ∼25 Myr.We first considered the isochrones and iso-masses predicted by evolutionary models in a (T eff , L) plane.We used the bolometric luminosities and effective temperatures derived by Bonnefoy et al. (2014a).Fig. 9 compares our observed T eff and L to the BHAC15 tracks.The two components, A and B, are not located on the same isochrone, the primary on 10 Myr and the secondary 20 Myr, but their positions are consistent with coevality between 10 and 25 Myr within 1 sigma.On the other hand, the predicted masses are respectively around 0.06 and 0.07 M for A and B, at the lower end of the stellar regime, which gives a total mass of 0.13 M .Nevertheless, when we impose coevality at the moving group age and allow for a shift of T eff within the 1σ interval, we retrieve the total dynamical mass, with masses of about 0.08 and 0.1 M for A and B. The corresponding diagram is shown in the appendix for the DM97 model.Underprediction of the total mass and non coevality are again retrieved, but once again the discrepancy disappears when we impose coevality at the moving group's age and allow for a shift of T eff .
As opposed to the bolometric luminosity and dynamical mass, the effective temperature predicted by the models is not robust, as it depends strongly on the atmosphere model.For each component, we thus used the measured luminosity to compute the predicted mass for a range of ages with the BHAC15 and DM97 models.The corresponding plot is displayed in Fig. 8, data have been linearly interpolated where necessary to provide predictions at the required ages.The prediction at the moving group age is consistent with the dynamical mass.Fig. 8: Comparison of TWA 22 direct mass measurement for two different βPic-MG age estimations with the predicted masses of the BHAC15 and DM97 tracks derived from the bolometric luminosity.Errors on the photometry are propagated on predictions (colored zone).The error on the distance is taken here into account.

GJ 2060
According to the previous sections, GJ 2060 has a total dynamical mass of m tot = 1.09 ± 0.10 M , and its age estimate can go from 30 to 200 Myr.
We first considered the isochrones and iso-masses predicted by evolutionary models in a (T eff , L) plane.Fig. 9 compares our observed T eff and L to the BHAC15 tracks.The two components, A and B, are located on the same isochrone, around 40 Myr, which is consistent with the younger estimations of the ABDor-MG age.On the other hand, the predicted masses are respectively around 0.55 and 0.3 M for A and B, which gives a total mass of 0.85 M , far (2σ) from the 1.09 M obtained by the orbital fit.We tried to impose a total mass of 0.85 M in the orbital fit in order to evaluate how this would change the distribution of χ 2 .In that case, this leads to orbits with χ 2 red > 8.5 for a distance of 15.69 pc, and χ 2 red > 2.5 for the 1 σ distance 15.24 pc, compared to 0.5 when the mass is set free.Therefore, the predicted mass does not account for the astrometry of the system.
The corresponding diagrams are shown in the appendix for all the other models.Underprediction of the total mass are retrieved in each case.Coevality is sometimes only marginally achieved (PARSEC), and a very young age can be predicted (20 Myr, DM97).
As in the TWA 22 case, we then used each component measured luminosity to compute the predicted mass for a range of ages with the six models (BHAC15, DM97, PARSEC, PISA, Darmouth and SDF00) and we deduce a plot linking the mass and age for the observed luminosity.These plots are displayed in Fig. 10.We retrieve the ∼ 20% underestimation of the total mass (2σ deviation) if a young age is assumed.Conversely, an old age (> 150 M ) gives a mass marginally compatible (1σ) with the orbital fit.
From the plots, we computed the predicted mass for each model and different ages of the AB Dor moving group.The results are displayed in Table 10.In order to avoid summing correlated errors, we drew the mass-age relation for several distances, and determined the system mass in each case.We computed the spread and deduced the uncertainty due to the distance σ d .For the most probable distance, we then derived the age uncertainty due to the luminosities σ L .The final ages uncertainties are then the quadratic sum of the independent errors σ L and σ d .Only the >100 Myr case fits marginally within the 68 % interval of confidence of the MCMC probability distribution of the dynamical mass.That age is inconsistent with the positions of both stars on the temperature-luminosity diagram for all models, except for PARSEC.Table 10: Predicted mass (in solar mass) for system GJ 2060 depending on the evolutionary model, from its luminosity and for several assumed ages.The error on the distance is taken here into account.The dynamical mass is 1.09 ± 0.10 M .

Discussion
The discrepancy evidenced on the dynamical mass of GJ 2060 AB ranges from 1 to 2 σ, depending on the system age.Such a disparity is not statistically impossible, as it represents respectively the edge of the 68% and 95% confidence interval.As an example, a 1 to 2 σ overestimation of the parallax could solve the issue while being a legit statistical realization.We present below some alternative hypotheses.

Data derivation and interpretation
The SINFONI spectra of GJ 2060 A and B are well fitted by the BT-SETTL model (Fig. 3).Consistent estimates of the bolometric luminosities can be inferred from the multi-epoch observations of the pair (Table 5).Thus, we can focus on vetting the dynamical mass estimate.
GJ 2060 astrometry has been measured with many different instruments, so that systematic errors can lead to important biases if they are not accounted for in the error bars.However, we performed the MCMC fit to each instrument astrometry, with and without the radial velocity measurements, and found very close values to the one given in section 5, for both orbital elements and total mass peak values, that can definitely not account for the 0.1 or 0.2 M difference.We find a total mass of m tot = 1.10 ± 0.14 M when we only consider the largest homogeneous sample of astrometric epochs (AstraLux), and m tot = 1.05 ± 0.12 M when we consider the three main sets of astrometric epochs (AstraLux, SPHERE, Keck).We estimate 1.08 ± 0.12 M when we consider all the astrometry and exclude the radial velocity measurements.In all these cases, as the orbit is less constrained, the mass can agree within 1σ with the model predictions for the old age ranges of the AB Dor moving group.
The absolute orientation of the field are usually inferred from the observations of different reference astrometric fields (clusters).This orientation could not be derived in a homogeneous way for all our astrometric epochs and instruments.Therefore, it may introduce a bias on the orbital parameter determination of GJ 2060 AB.A systematic on the pixel scale of the instruments is less likely to change our results given the short separation of GJ 2060 AB.Therefore, we added parameters to the MCMC algorithm in order to estimate and account for systematic angular offsets in the astrometry.Our astrometry consists of 5 different samples, but only three of them contain more than 2 data points, AstraLux, SPHERE and Keck.We performed an orbital fit with only these samples along with the RVs, with AstraLux data (more numerous) taken as reference.Two offset parameters α 1 and α 2 had thus to be added to the original MCMC algorithm : the AstraLux data are fitted as they are, with a model corresponding to Eqs. ( 1) and ( 2), while the SPHERE data are first rotated through the angle α 1 and the Keck data through the angle α 2 .We modified the algorithm to allow any number N of samples (that is N different instruments) as astrometry input.One sample has to be designated as the reference and the algorithm fits N-1 angular offsets, assuming a flat prior (e.g., Montet et al. 2015).The results are then displayed with posterior distributions and correlations to the other parameters.We found a distribution centered around -0.18 • for the offset between AstraLux and SPHERE data, with standard deviation 0.2 • .For the offset between AstraLux and Keck, the distribution is centered around 0.10 • , with standard deviation 0.3 • .Near apoastron, an angular offset of 0.3 • corresponds to 2-3 mas of offset on the right ascension.However, the total and fractional mass and remain unaffected.We are then confident that our dynamical mass estimate is not strongly affected by these angular offsets.

Models imprecision at the moving group ages
Pre-main sequence models have a well known tendency to underestimate significantly the mass of low-mass stars (Hillenbrand & White 2004;Mathieu et al. 2007).Mathieu et al. (2007) studied the 23 PMS stars for which a dynamical mass had been derived and compared these masses to the predictions of the evolutionary models, given the bolometric luminosity and effective temperature of the stars.They highlighted a mean underestimation of 20-30% for low-mass stars (< 1 M ), similar to the underestimation of GJ 2060 mass.Since then, new evolutionary models have been designed.Moreover, dozens of new dynamical masses have been obtained for PMS stars in the mean time (most of them for stars younger than 10 Myr).These studies often con-firm the previously noticed mass discrepancy (e.g., Simon et al. 2017;Mizuki et al. 2018).
Among these systems, some are comparable to our objects.In the ABDor moving group, the systems ABDor Bab (Azulay et al. 2015;Janson 2018) and C (Close et al. 2005;Luhman & Potter 2006;Close et al. 2007;Boccaletti et al. 2008;Azulay et al. 2017) have been deeply analyzed in relation to the discussion about the age of the moving group.ABDor Ba and Bb's dynamical masses and luminosities are not consistent with the PMS isochrones (Janson 2018, Paper II).At the given luminosities and for a moving group age of 150 Myr, the predicted masses are ∼ 25% below the dynamical masses, which are similar to GJ 2060 B. The study of AB Dor C is consistent with any age between 20 and 120 Myr, and the mass derived from the models is slightly underestimated (10%) but still consistent with the dynamical mass, which is similar to that of TWA 22 components (Azulay et al. 2017).
On the other hand, NTT 33370 AB is a 80 Myr low-mass binary very similar to TWA 22 in terms of masses (Schlieder et al. 2014;Dupuy et al. 2016).Both individual masses are strongly underpredicted by the BHAC15 model (2σ, 46 +16 +19 %), which contrasts with the perfect agreement we found for TWA 22.
This issue does not disappear for older ages of the PMS regime.System 2M1036 is a triple M-dwarf stellar system in the Ursa Major moving group whose age is estimated at 400-500 Myr (Brandt & Huang 2015;Jones et al. 2015).Calissendorff et al. (2017) evidenced a 1σ underestimation on each component mass.Conversely, in the same moving group, the K binary system NO UMa dynamical masses are in good agreement with the model predictions (Schlieder et al. 2016).
The evolutionary models have not yet entirely mastered the physics of PMS stars, as can be seen from the frequent mass underestimation.It is in particular surprising that two very similar systems can encounter very different prediction agreements.Facing up to the mass over-estimation of system GJ 1108 A, Mizuki et al. (2018) compared a dozen PMS stars dynamical masses with the predictions of the BHAC15 model, and report a ≤ 10% offset towards underestimation, and ∼ 20% scatter.Their results also confirm that the tendency to underpredict the mass is neither associated with a mass range nor with an age.
Magnetic activity is also often brought up as a cause of discrepancy in low-mass stars, as it affects greatly the convection and induces large spot coverage fractions, that may lead to displacements on the HR diagram (eg., Feiden 2015;Somers & Pinsonneault 2015).The high jitter in HARPS RV measurements (∼ 400 m/s) indicates that GJ 2060 has a strong magnetic activity.Somers & Pinsonneault (2015) studied the influence of spots on PMS stars and showed that it could lead to non negligible radius inflation, that would then lower the stars effective temperatures and luminosities.The gap bewteen the normal and spotted case depends a little on the star's age and strongly on the star's mass.Following Fig. 1 B, we assumed a ∆L of -10% and a ∆T eff of -5% for the primary, and ∆L of -20% and a ∆T eff of -8% for the secondary.We then plot the new position of GJ 2060 A and B on the BHAC15 isochrones on Fig. 18 (f) in the appendix.The positions are shifted of 0.04 dex and 0.1 dex towards the brighter luminosities, and 200 K and 300 K towards the hotter temperatures.The diagram is now consistent with coevality at 150 Myr, and the total mass that is derived matches the dynamical estimate.These corrections are computed with a spot surface coverage of 50%.The intense activity of the stars could thus account for the disagreement with the models.A shift in temperature could also resolve the slight mismatch of TWA 22 that appears in the HR diagram.However, a higher luminosity would lead to an excessive mass.This hints for a lesser activity-induced effect for TWA 22 components.In order to test this hypothesis, it would be worth comparing the activity indicators of different PMS binaries with their predicted mass discrepancy.Such study was done recently by Stassun et al. (2014) for eclipsing binaries, where they evidenced that activity was not the only cause of the disparity.Finally, in our case where the two components are regularly (at each periastron passage) very close to each other (< 1 au), tidal interactions may also affect the evolution of the stars, although the effects are expected to be weaker than in the eclipsing binary cases, that are constantly undergoing strong interactions.An in-depth study would be needed however to determine the effect of tidal forces on the interiors of eccentric binaries.
On the other hand, Simon et al. (2017) suggested that all underestimations come from hidden components within the systems.If it is unlikely that such explanation accounts for all the observed discrepancies, in particular within tight binaries, it is nevertheless a suggestion worth studying for GJ 2060, especially given its unusually high eccentricity (e ∼ 0.9).

Missing mass: existence of GJ 2060 Ab or Bb
Hidden mass close to the primary could explain the strong disagreement between models and data for GJ 2060.Indeed, a 0.1 to 0.2 M (depending on the system age) additional companion could account for the mass underestimation (see table 12) and could have been missed in the SPHERE datasets if close enough from one of the two components.Quick dynamical simulations with a symplectic integrator (SWIFT_ HJS Beust2003, see appendix) evidenced that a additional companion should orbit closer than 0.1 au (6 mas) from either of the companion to remain bound for the system's lifetime.
The PSFs of GJ 2060 A and B are not elongated, even in the SPHERE images (FWHM∼30 mas).In the latter, we injected models of putative companions with different fluxes and separations and checked whether they would induce a PSF lengthening that could be noticed by eye.The models of the putative companion were built using a flux-normalized PSF.The PSF is the other component of the system (e.g., B if the binarity of A is investigated, and vice versa).Using the BHAC15 models, we estimate that we would have just missed a 0.25 M companion at a projected separation 0.45 au.Thus, 0.2 solar mass object at 0.1 au would have gone unnoticed from the imagers.
As for the spectrograph, the available RV data are too sparse to resolve in frequency an additional orbit, and the flux ratio in the optical prevents the detection of any spectral signature.However, the closer the object, the stronger the radial velocity perturbation amplitude.A simple comparison between the perturbation amplitude on GJ 2060 A and B radial velocities and our measurements standard deviation σ is summarized in Fig. 11 for the circular case, in a (semi-major axis, mass) and a (semimajor axis, inclination) diagrams.We used the predicted mass of GJ 2060 A and B from Table 11 at an age of 100 Myr for that purpose.We chose 3σ dispersion of the RVs as a detection threshold, and represented the corresponding frontier on the plots.The limit of the dynamical stability has been set to 0.1 au; the accurate stability limit depends on the third companion mass and inclination, but is in all cases 0.1 au.In the coplanar case, a mass higher than 0.1 M could have been unnoticed around the secondary.That is not the case for a putative companion around the primary, because the light we observe comes mostly from the primary, so that most signals would be easily spotted.However, a 0.1 or even 0.2 M at 0.1 au could be compatible with our deviation in both situations, primary or secondary, for small enough inclinations, respectively 10 and 5 • around the primary, and 45 and 25 • around the secondary.If there is indeed a hidden companion, its luminosity would add to the luminosity of the nearest component, so that the latter measured flux would be biased.According to the BHAC15, a 50-75 Myr 0.2 M companion have log(L/L ) = −1.9,−2 dex and a 0.1 M companion have log(L/L ) ∼ −2.3 dex.The component hosting a hidden companion would appear overluminous for its temperature (slightly for the primary, significantly for the secondary), shifting its position on Fig. 9 towards the younger isochrones and straining coevality.In the PARSEC isochrones in the appendix, a significant luminosity shift (corresponding to a 0.2 M companion) of the primary could achieve coevality.Conversely, the same companion around the secondary would induce a luminosity shift that would break coevality in all models.
Finally, the high eccentricity (0.90 ± 0.01) of the visual orbit is noticeable, and one can wonder if it could indicate strong dynamical interactions.From the ORB6 catalog3 , we computed the periods and eccentricities of visual binary stars with reliable orbital elements (according to the grades given in the catalog).From the SB9 catalog (Pourbaix et al. 2004), we computed the periods and eccentricities of spectroscopic binary stars.Our two binaries fall near the limit of each catalog's period coverage, so that while it gives an interesting overview, more binaries would be needed to draw robust statistical conclusions.Gaia next data releases will significantly contribute to overcome this lack.The resulting diagrams are shown in Fig. 12.While not so common, GJ 2060's eccentricity seems not so rare at this range of periods and age (before circularization).Moreover, no mechanism are known to enhance the eccentricity of an outer companion at such period ratio (greater than 200, that excludes any meaningful mean-motion resonance).Only close encounters could dynamically raise the eccentricity, but the configuration would then not be stable.All in all, the high eccentricity is likely uncorrelated to the potential existence of a third companion.

Conclusion
We considered two systems of young astrometric M-dwarf binaries, TWA 22 and GJ 2060, and used existing astrometric and spectroscopic data along with new Keck, SPHERE, NaCo, HARPS and FEROS data to derive the total mass of these systems.We consolidated the total dynamical mass estimate of TWA 22: 0.18 ± 0.02 M .We derived the first estimate of the total mass of GJ 2060: 1.08 ± 0.10 M .The orbits of the two systems are well constrained thanks to the use of our additional data and the errors are carefully estimated through the MCMC approach.The orbit of GJ 2060 has an unusually high eccentricity, around 0.9.The cross-contamination of GJ 2060 primary and secondary spectra into the FEROS and HARPS data prevents us from deriving accurate dynamical masses of the individual components.
The study of the photometry and spectroscopy of the two systems, along with their membership to moving groups and accurate distances, allow to test the PMS evolutionary models predictions.The dynamical mass of TWA22 AB is correctly predicted by the models at the age of the β Pictoris moving group.The placement of GJ 2060 A and B on evolutionary tracks confirms the system coevality at an age compatible with the AB Doradus moving group (∼ 50 Myr).However, all models underpredict the total mass of GJ 2060 AB, by 10 to 20% (0.1 to 0.2 M , 1-2 σ).A new precise parallax (likely to come in Gaia DR3 release) would strongly decrease the uncertainty on the dynamical mass and could improve the statistical relevance of the discrepancy.
GJ 2060 AB's underpredicted mass is consistent with a trend found on other systems in the same mass range.It could be explained by luminosity and temperature drop caused by high starspots coverage.In that case, we would retrieve coevality at 150 Myr.We also discussed the potential existence of a third companion close to one component of GJ 2060 that could account for this disagreement.Dynamical modeling shows that such a companion would have to be very close to one of the stars, less than 0.1 au (6 mas), in order to remain stable for millions of years.Such a close companion could have gone unnoticed, although the RVs are putting some constraints on its mass and inclination.Astrometric and spectroscopic data at periastron as well as the use of RV disentanglement techniques could help clarifying the origin of the discrepancy, and in particular if only one of GJ 2060 AB's component has an underpredicted mass.
A dozen of new PMS stellar mass measurements have become available in the last decade.A complete reassessment of the dynamical mass determinations of sub-solar mass stars and a homogeneous comparison of theses measurements to the latest PMS models would help to conclude on the models reliability, On the other hand, the upcoming data releases of the Gaia mission should yield a statistical sample of dynamical mass determination of low-mass stars (Pourbaix 2011).Additional studies should in any case be needed to infer the luminosity and temperatures of theses many systems and allow for a detailed comparison of the masses to evolutionary models predictions.
was observed in field-tracking mode on February 11, 2013 with the NAOS-CONICA (NaCo) adaptive-optics instrument mounted on the VLT/UT4 (Lenzen et al. 2003; Rousset et al. 2003) as part of a program dedicated to the orbit monitoring of young binaries (PI Bonnefoy; program ID 090.C-0819

FFig. 3 :
Fig.3: Spectra (apparent flux) of GJ 2060 A and B compared to the best fitting BT-SETTL synthetic spectra.

Fig. 4 :
Fig.4: Plots of a hundred orbits obtained with the MCMC algorithm for system TWA 22. Astrometric measurements are displayed in their instrument color, along with their position on the fit.Only NaCo data are used in the orbital fit.The orbit in black, obtained with the LSLM algorithm, corresponds to the lower χ 2 .

Fig. 5 :
Fig.5: Plots of a hundred orbits obtained with the MCMC algorithm for system GJ 2060.Astrometric measurements are displayed in their instrument color, along with their position on the fit.The orbit in black, obtained with the LSLM algorithm, corresponds to the lower χ 2 .
the components luminosities in the visible spectrum and K A and K B are respectively the amplitudes from A and B. From this relation, for any distance d, we obtain a fractional mass of m B /m tot = 0.46 ± 0.10 15.69 pc d .Using the parallax distance and propagating its uncertainty, we finally obtain m tot = 1.09 ± 0.10 M and m B m tot = 0.46 ± 0.10.

Fig. 7 :
Fig. 7: Isochrones and iso-mass curves predicted by BHAC15.The 1, 10, 20, 25 and 50 Myr isochrones have been drawn in dotted lines (1 Myr is up), while the iso-mass curves are in solid line .The 25 Myr isochrone correspond to the age of the βPic-MG, and is drawn in red.The blue patterns correspond to the observed values and their error bars for each component of system TWA 22, A and B.

Fig. 9 :
Fig. 9: Isochrones and iso-mass curves predicted by BHAC15.The 10, 20, 50, 75, 150 and 600 Myr isochrones have been drawn with dashed lines (1 Myr is up), while one iso-mass is drawn in solid line every 0.1 M from 0.1 (left) to 1 M (right).The 50, 75 and 150 Myr isochrones correspond to possible ages for ABDor-MG, and are drawn in red.The blue patterns correspond to the observed values and their error bars for each component of system GJ 2060, A and B.

LFig. 10 :
Fig. 10: Mass age relations according according to the six different evolutionary models, for GJ 2060 observed luminosities.A distance of 15.69 pc is assumed.The dynamical mass and its uncertainty are depicted in red.

Fig. 11 :
Fig. 11: Radial velocity detection limits around GJ 2060 A and B as a function of the inner orbit radius (coplanar case), in terms of the inner companion (a) mass (b) inclination.For (a), the orbits have been chosen coplanar (i = 36 • ) and for (b), the inner mass has been set to 0.1 M .

Fig. 13 :
Fig. 13: Comparison of the J and K band spectra of GJ 2060 A (red) and GJ 2060 B (blue) to M-dwarf spectra.

Fig. 14 :
Fig. 14: Comparison of the H band spectra of GJ 2060 A (red) and GJ 2060 B (blue) to M-dwarf spectra.

Fig. 15 :Fig. 16 :
Fig.15: Distribution and correlations of each of the orbital element fitted by the MCMC algorithm for system GJ 2060.The black lines and points depict the best fitting orbit (better χ 2 ), obtained with the LSLM algorithm.The color scale is logarithmic, blue corresponds to 1 orbit and red to 1000.

Table 1 :
Age estimates of the Beta Pictoris moving group.The two components are the least massive stars in the β Pic group for which a dynamical mass has been computed.They complete the mass sampling between the giant planet β Pictoris b (Lagrange et al. 2010) and the higher mass binaries GJ 3305 (total mass 1.1 M ; Montet et al. 2015) and V343 Nor (total mass 1.4 M ; Nielsen et al. 2016).TWA 22 is thus an essential benchmark to test the predictions of the evolutionary models in the young group in the 0.1 M mass range.

Table 2 :
Age estimates of the AB Doradus moving group.

Table 3 :
Observing log.The field rotation θ is given when the observations are performed in pupil-tracking mode.

Table 5 :
Contrasts and apparent magnitude of GJ 2060 A and B .2.5.SINFONI integral field spectroscopy GJ 2060 AB was finally observed on November 22, 2013 with the SINFONI instrument mounted on the VLT/UT4 as part of our dedicated program for the orbital characterization of dynamical calibrators (PI Bonnefoy; program ID 090.C-0819).SINFONI (Spectrograph for INtegral Field Observations in the Near Infrared) couples a modified version of the adaptive optics module MACAO • .3

Table 7 :
Durkan et al. (2018)ial velocity measurements of the SB1 GJ 2060.Ten radial velocity measurements have been obtained using the Fiberfed Extended Range Optical Spectrograph (FEROS;Kaufer et al. 1999) mounted at the ESO-MPG 2.2 m telescope at La Silla Observatory.The data reduction process is described therein.FEROS is an echelle spectrograph covering the wavelength range 3500 -9200 Å across 39 orders with R ≈ 48000.The measurements are reported inDurkan et al. (2018), as part of a radial velocity monitoring survey of young, low-mass binaries.They cover a 12 year span, from 2005 to 2017.

Table 8 :
Bonnefoy et al. (2009)the MCMC fit of TWA 22 relative orbit, compared to the last orbit determination byBonnefoy et al. (2009).The uncertainties on the fitted parameters correspond to the 68% interval of confidence of the distribution probabilities (see appendix A).The astrometric data only allow determination of the couple (Ω, ω) modulo π.

Table 9 :
Orbital elements from the MCMC fit of GJ 2060 AB relative orbit.The uncertainties on the fitted parameters correspond to the 68% interval of confidence of the distribution probabilities (see appendix B).

Table 11 :
Mean predicted mass (in solar mass) for GJ 2060 A and B, from their luminosities.

Table 12 :
GJ 2060 missing mass (in solar mass) depending on models and age.Errors propagation have been obtained from the MCMC posteriors dispersion, the luminosities uncertainty at given distances and the errors on the distance, assuming independancy.