Updated orbital monitoring and dynamical masses for nearby M-dwarf binaries

Young M-type binaries are particularly useful for precise isochronal dating by taking advantage of their extended pre-main sequence evolution. Orbital monitoring of these low-mass objects becomes essential in constraining their fundamental properties, as dynamical masses can be extracted from their Keplerian motion. Here, we present the combined efforts of the AstraLux Large Multiplicity Survey, together with a filler sub-programme from the SpHere INfrared Exoplanet (SHINE) project and previously unpublished data from the FastCam lucky imaging camera at the Nordical Optical Telescope (NOT) and the NaCo instrument at the Very Large Telescope (VLT). Building on previous work, we use archival and new astrometric data to constrain orbital parameters for 20 M-type binaries. We identify that eight of the binaries have strong Bayesian probabilities and belong to known young moving groups (YMGs). We provide a first attempt at constraining orbital parameters for 14 of the binaries in our sample, with the remaining six having previously fitted orbits for which we provide additional astrometric data and updated Gaia parallaxes. The substantial orbital information built up here for four of the binaries allows for direct comparison between individual dynamical masses and theoretical masses from stellar evolutionary model isochrones, with an additional three binary systems with tentative individual dynamical mass estimates likely to be improved in the near future. We attained an overall agreement between the dynamical masses and the theoretical masses from the isochrones based on the assumed YMG age of the respective binary pair. The two systems with the best orbital constrains for which we obtained individual dynamical masses, J0728 and J2317, display higher dynamical masses than predicted by evolutionary models.


Introduction
The study of the multiplicity of stars is a useful diagnostic for obtaining insight into their formation and dynamical evolution, as it allows for important properties such as binary fraction, semimajor axis distribution, and mass ratios to be constrained (e.g. Burgasser et al. 2007). Since low-mass M-dwarf stars form a nat-ural link between the substellar brown dwarfs and the solar-type stars, and the multiplicity frequency tend to decline with lower masses and later spectral types (Duchêne & Kraus 2013; Moe & Di Stefano 2017;Winters et al. 2019), it becomes even more crucial to discover and characterise low-mass M-dwarf multiples. Hence, a rigorous understanding of the multiplicity characteristics and their evolution within this transitional mass-region that is made up of M dwarfs is vital for constraining formation scenarios of low-mass stars and brown dwarfs. Astrometric monitoring of such binary systems allow for dynamical masses to be derived, which become essential in the efforts of empirical calibrations of fundamental properties such as the mass-luminosity relation and evolutionary models (Dupuy & Liu 2017;Mann et al. 2019;Rizzuto et al. 2020). This becomes even more important at the lowest stellar masses for which the current theoretical models have been shown to systematically underestimate Mdwarf masses below M ≤ 0.5 M by 5−50 % (e.g. Hillenbrand & White 2004;Montet et al. 2015;Calissendorff et al. 2017;Biller et al. 2022). Since M dwarfs evolve slowly and remain in their pre-main sequence phase for ∼ 100 Myrs (Baraffe et al. 1998), M-dwarf binaries become valuable benchmark targets for astrophysical calibrations comparing dynamical masses from observational data to isochronal models (e.g. Janson et al. 2017). As such, groups and associations of stars that can be expected to have originated from the same mutual cluster or region, commonly referred to as young moving groups (YMGs), have seen an increase in interest of late (e.g. Torres et al. 2008;Malo et al. 2013). Thus, M dwarfs residing in YMGs can be isochronally dated, and binaries with estimated dynamical masses can also be used to robustly test the coevality of the YMGs.
Motivated by these arguments for binary characterisation and low-mass multiplicity studies, the AstraLux Large M-dwarf Multiplicity Survey systematically studied over 1,000 X-ray active M dwarfs with the lucky imaging technique, identifying ≈ 30 % as multiple systems, many of which are known YMG members (Bergfors et al. 2010;Janson et al. 2012Janson et al. , 2017. Although most of these binaries have separations which correspond to orbital periods of several decades to hundreds of years, some have periods short enough so that they can already be mapped out after a few years of monitoring. The SpHere INfrared Exoplanet (SHINE) project utilising the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2019) instrument at the Very Large Telescope (VLT) is surveying 500 stars with the purpose of directly detecting substellar companions to the stars in order to better understand their formation and early evolution (Chauvin et al. 2017b;Vigan et al. 2021). As an auxiliary result from the survey, several low-mass binaries that coincide with the AstraLux survey sample (Desidera et al. 2021) have been observed with highcontrast imaging which provides high quality precision measurements (Langlois et al. 2021), which are excellent for astrometry and useful for constraining orbital motion.
Here we present the latest results from the combined effort of the AstraLux M-dwarf multiplicity monitoring programme and the SHINE M-dwarf filler programme. We have identified the 20 most prominent systems for fundamental properties to be characterised from the AstraLux Large Multiplicity Survey which have sufficient orbital coverage with which first-hand constraints can be made from orbital fitting routines. Out of these 20 systems, eight have strong indicators to place them in YMGs and thereby have their ages constrained.
The paper is divided up into the following sections that cover the following areas: Section 2 where we go into detail on how the target sample was collected from the different surveys we combined, and the observations taken along with how data were reduced. In section 3 we explain the orbital fitting procedures and present the main results and discussions in Section 4 where dynamical masses are compared to theoretical isochronal models. Finally we provide a summary and conclusions in Section 5. The collected astrometric data are given in Appendix A, together with the resulting orbital fits for the binaries in Appendix B.

Sample selection and observations
The target list was created from a sample of known binaries and higher hierarchical-systems from the AstraLux M-dwarf multiplicity survey (Janson et al. 2012), consisting of over 200 Mdwarf multiple systems with separations within 0.08 − 6.0 , and the extended AstraLux sample (Janson et al. 2014a) of ≈ 60 multiples with spectral types M5 and later. We selected 20 systems for which we had identified to either have sufficient orbital coverage that dynamical masses could be robustly constrained, or undergone enough monitoring that ≥ 25 % of the orbit can be mapped out to provide some useful information. We present here new observations from our survey and some previously unpublished observations of the targets.
We compared the space velocities and positions of the targets with those of known YMGs and associations using the BANYAN Σ-online tool (Gagné et al. 2018), the LACEwING code (Riedel et al. 2017) and the GALEX convergence tool (Rodriguez et al. 2013). Unless otherwise specified, parallaxes and proper motions were obtained from the Gaia archive (Gaia Collaboration et al. 2016), both Gaia Data Release 2 (DR2; Gaia Collaboration et al. 2018) and Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021). Spectral types presented in Table 1 were derived by Janson et al. (2012), using the (i − z ) photometry obtained from the AstraLux observations and following the methods by Daemgen et al. (2007). Some of the systems have additional information from the resolved near-IR medium resolution spectra from SINFONI which Calissendorff et al. (2020) used to derive near-IR spectral types from the JHK bands and surface-gravity sensitive emission lines.
As the systems we present in this survey are binaries, we refer to the distance to the systems from us as d in pc, and the separation between the binary components as s, either in projected separations of milliarcseconds (mas) or physical as AU. The target binaries all have designated Two Micron All-Sky Survey (2MASS) identifiers, and we abbreviate them by their first four to six digits as Jhhmm(ss). The target systems are presented in Table 1. The YMGs of interest and their adopted ages are shown in Table 2, and the target parallaxes and space velocities shown in Table 3 which were used to derive membership probabilities for each source. Not all YMGs and associations are included in each of the YMG membership probability tools, and we did not introduce additional YMGs to the existing code.

AstraLux
The AstraLux Large Multiplicity Survey has been ongoing for over a decade, collecting data of numerous visual binaries by applying the lucky imaging technique. The survey primarily employs two principle instruments; AstraLux Norte on the 2.2m telescope in Calar Alto, Spain (Hormuth et al. 2008), as well as AstraLux Sur at the 3.5m New Technology Telecope (NTT) at La Silla, Chile (Hippler et al. 2009). The full frame field of view for the respective AstraLux instrument is ≈ 24 × 24 for  Norte and ≈ 15.7 ×15.7 for Sur, although typical observations utilise subarray redouts in order to minimise readout times. As-traLux observations are mainly carried out in the SDSS z -and i -bands, with a preference towards the z -band due to its smaller susceptibility to atmospheric refraction compared to the i -band (Bergfors et al. 2010).
Our observations typically consisted of 10 000 -20 000 short exposures of just 15 − 30 ms each, adding up to a total of 300 s integration time. Both AstraLux Norte and AstraLux Sur data were reduced with the real-time pipeline at the time of the observations, producing a final image from each observation where a subset of 1 − 20% of the best frames taken were kept. Generally images where 10% of the frames were kept provided a decent trade-off between sensitivity and resolution. Occasionally for closely separated binaries of similar magnitudes the pipeline would centre the frames on the secondary instead of the primary star, leading to a false stellar ghost to appear at the same separation but shifted at a 180 • from the real secondary (Bergfors et al. 2010).
We performed calibrations for the astrometric measurements with AstraLux by comparing observations of the Orion Trapezium Cluster and M15 to reference observations of the same fields from McCaughrean & Stauffer (1994) and van der Marel et al. (2002). The calibrations were performed by measuring the positions of bright stars within the same field that were recognisable and easily identified. We employed between 5-14 reference stars for the calibrations depending on the quality of the point spread function (PSF) and brightness. We assigned the brightest star in the field of view as the main reference, for which we calculated the relative separation and positional angle to for all other reference stars. We then compared the separations and positional angles for our AstraLux measurements to those of the reference observations, taking the average ratio of the separation as the plate scale and the standard deviation as its uncertainty. Correction for True North was performed in similar manner where the average difference in positional angle was used and standard deviation from the average assigned as the uncertainty. The final AstraLux astrometric calibrations calculated here and those obtained from earlier literature are listed in Table. 4. For two epochs, February and April of 2015, we did not have proper reference fields to calibrate the astrometry to with AstraLux Norte, and we assumed a mean pixel scale and correction for True North from the other AstraLux Norte epochs with proper calibration. This only affected two observations of the J1036BC binary and we assumed that the instrument had not changed significantly at this time compared to our other observed epochs. An alteration in plate scale from our smallest to largest calibration values would only change the resulting projected separation for the binary by ∼ 1 mas.

NaCo
We downloaded the NaCo data from the ESO archive together with their associated calibration files and performed basic reductions using custom scripts with Python. These basic reductions included corrections for bias, dark, flatfield division and combination of multiple frames.
We applied the astrometric corrections from Chauvin et al.
(2010) using plate scales 27.01 ± 0.05 mas/pxl for observations in the L -band and 13.25 ± 0.05 mas/pxl for shorter wavelengths. We do not correct the observing frames here for True North, but instead add a factor of ±0.20 deg to the uncertainty for the positional angle of each astrometric data point given by NaCo data, which is in line with the True North corrections obtained by Chauvin et al. (2010).

NOT FastCam
The Lucky Imaging FastCam is an instrument at the Roque de los Muchachos Observatory on La Palma in the Canaries, Spain, designed and capable of obtaining high-resolution images in optical wavelengths from medium-sized ground-based telescopes at the observatory (Oscoz et al. 2008). The instrument features a 512 × 512 pixels L3CCD from Andor Technology, and a special software package that reduces images in parallel with the data acquisition at the telescope, so that a small fraction of images with minimal atmospheric turbulence can be evaluated in real-time. For our observations, the FastCam instrument was mounted at the 2.56-m Nordic Optical Telescope (NOT), and the observations were carried out in August and November of 2016 using the I-band at 820 nm.
The astrometric calibrations were made in a similar way to our AstraLux astrometric calibrations, comparing reference fields of the M15 stellar cluster with images taken by the Hubble Space Telescope. We obtained a platescale of 30.6 ± 0.1 mas/pixel for the 2016.63 epoch in August, and a platescale of 30.5 ± 0.1 mas/pixel for the November epoch of 2016.87. The corresponding corrections for True North were +3.64 ± 0.01 • and −1.54 ± 0.1 • respectively.

SPHERE
The SPHERE data were collected as part of a sub-programme for the SHINE survey (Chauvin et al. 2017a). The filler programme from which our observations was taken were devoted to astrometric monitoring of tight visual binaries, many of which had been discovered in the AstraLux survey.
All SPHERE data, both IRDIS and IFS modes, were downloaded and reduced using the SPHERE data centre (Pavlov et al. 2008;Delorme et al. 2017). The reductions carried out by the automated pipeline included basic corrections for bad pixels, dark current, flat field, as well as corrections for the instrument distortion (Maire et al. 2016a) and rotation. Calibrations of the platescale and for the True North angle were performed in accordance to Maire et al. (2016b), with typical corrections of ≈ 12.267 mas/pixel for IRDIS and ≈ 7.46 mas/pixel for IFS observations, with a True North of ≈ −1.75 • . The specific plate scale and True North corrections handled by the pipeline are stated for each individual observing data point.
From the astrometric measurements described in Section 2.2 we also obtained accurate flux-ratios for the components in each system. We summarised the resulting contrast magnitudes in the SPHERE dual-band images in Table 5. Not all observed epochs had separate dual-band images available and thus missing from the table.
We did not use the IFS data to perform any spectral analysis of the targets here, some of which have superseding spectral information from the SINFONI observations instead (Calissendorff et al. 2020). Instead, we collapsed the data cubes and performed astrometry on a single frame from the IFS data.

Astrometry
Astrometric positions were calculated with the same procedure as described in Calissendorff et al. (2017Calissendorff et al. ( , 2019Calissendorff et al. ( , 2020. Concisely, a grid in x and y positions was constructed where we scaled the brightness of a reference PSF, placing two of them on the grid which were sequentially shifted in positions to match the observed data. A residual was then calculated by subtracting the constructed model from the observed data, and the procedure was iterated while scaling the brightnesses and shifting the positions of the model until a minimum residual could be found. The basic workflow of the astrometry extraction procedure is illustrated in Figure 1 where we used the J1036BC binary and our AstraLux data from April 2015 as an example.
The astrometric extraction was performed in the same manner for all observations and instruments considered here. Generally we would try to obtain a good PSF reference from the same or close to the same epoch as the observation we were extracting astrometry from. In previous AstraLux observing campaigns, designated PSF reference in the form of single stars have been procured (Janson et al. 2012(Janson et al. , 2014a. However, that was not the case for later epochs. Instead, we identified which binaries or higher hierarchical systems had large relative separations or isolated components, which we then used as PSF references. For AstraLux, FastCam and SPHERE we had access to the primary in the triplet 2MASS J10364483+1521394 system for several epochs, which was close to an ideal PSF reference for our intended purposes given the similar M-dwarf spectral type of the component to the rest of our target sample and that it was available for most out of our observed epochs. The primary component in the system has a projected separation of ≈ 1 to the outer binary pair and can be viewed as a single star-proxy in this con- The upper plots display the observed AstraLux data for the binary pair and a PSF reference. The lower plots show the constructed model using two brightness-scaled and position-shifted PSF models, and the residual after subtracting the model from the observed data. The plots here have been normalised to the peak value of the observed data, are scaled linearly in the plots. The colour-bar markers indicate the peak and bottom values for the observed data. The residual for the AstraLux is typically of the order of ±15 %, whereas with SPHERE the residual is more commonly around ±1 %. text. Another benefit from using the primary of J1036 was that whatever aberrations afflicted the observations, altering the PSFs of the binary, would also be seen in the reference PSF so that they could be accounted for. Nevertheless, to increase the statistical certainty of the astrometric measurements we typically used between 3-10 different reference PSFs depending on target quality and instrument applied. We then calculated the mean separation and positional angles, using the standard deviation as the error which we added quadratically together with the instrumental errors (plate scale and True North errors) to the uncertainty. We did not try to mix PSF references from observations taken with different instruments or settings. Due to the larger field of view and smaller plate scales for the VLT instruments NaCo and SPHERE, we could with ease isolate single components and use them as PSF references. For the As-traLux observations we had the advantage of having a plethora of multiple systems to choose from in the AstraLux Multiplicity Survey. The FastCam observations however did not have quite the same luxury, as the coarser plate scale made it more cumbersome to isolate single components. As such, we mainly used the primaries from the J0111 observations taken in August and the J1036 observations taken in November as PSF references. We also included three additional references for our FastCam astrometry; J0103, J0916 and J1641, which are known tightly bound binaries but appeared as unresolved single sources in the FastCam observations.
Since SINFONI was not calibrated for astrometry we added an extra uncertainty term. We checked the consistency between the SINFONI astrometry presented in Calissendorff et al. (2020) and our SPHERE astrometry for the binaries which were observed at similar epochs and found no large discrepancy between the two. We therefore included the SINFONI data points into the fitting procedure when they were believed to aid in constraining orbital parameters.

Radial velocities
The orbital motion from the two components in the binaries make them subject to Doppler shifts which can be measured and useful for constraining the orbital motion further. We searched the literature and uncovered unresolved radial velocity (RV) observations for 16 binaries in our sample, which we included into our MCMC fitting to aid the orbital fitting procedure. However, the RV data in the literature is mostly compromised of unresolved measurements in which the two lines from the two components in the binary are blended together. As the strength of these lines are dependent on the spectral template used and fitting method for deriving the RV measurement, which differs from authors and instruments, most RV data were deemed unusable for our the orbital fitting. Hence, only the RV data for seven systems was used in the final orbital fits, listed here in Table 6 and shown in their respective fit in Figure 2, where targets with too few RV observations or lack of baseline were omitted.
For the seven instances where RV data were available and useful, we included two additional parameters to the MCMC code which evaluated the probability density functions (PDFs) of the offset velocity v 0 and RV amplitude K. In the adopted formalism which assumed a Keplerian orbit, the radial velocity can be described as which amplitude for pure SB1 binaries is deduced from the mass fraction of the secondary component m B /m tot as In principle, the RV data allows for fractional mass and thereby individual masses for the binary components to be derived. However, that is when considering pure SB1 binaries, which is a questionable assumption for the relatively high fluxratios (and mass-ratios) in our target sample. Therefore, we applied the same method as in Rodet et al. (2018), proposed by Montet et al. (2015), and assumed the sum of two flux-weighted individual RVs to be the considered RV measured. The orbital fit could then fit the RV amplitude as being the fractional flux, L V A and L V B being the luminosities in the visible spectrum for each component, and K A and K B the respective RV amplitude.
The majority of the RV data were obtained from the Fiberfed Extended Range Optical Spectrograph (FEROS) instrument at the ESO-2.20m telescope (Kaufer et al. 1999), with a wavelength coverage of λ 3500 − 9200 Å and resolving power of R = 48, 000. We did not perform the RV observations or any reanalysis of the data here, using only the values stated from the given literature cited in Table 6. We estimated the flux ratios from the magnitude difference in the i -band from the AstraLux observations, as the i span the wavelength λ ∼ 6700 − 8400 Å, thereby encompassing the most similar wavelength range as the FEROS RV data. For the J2016 system we did not posses any flux-ratio in the i -band and applied the flux-ratio in the I-band from the FastCam/NOT observations instead, which has a reference wavelength of λ ref ∼ 8200 Å. The flux ratios used in our calculations are shown in Table 7. The reported uncertainties of the flux-ratios are likely underestimated due to the different wavelength coverage by the photometric bands and that of FEROS.
This approach of weighing RV signals by the flux-ratio proved to have some limitations, where some estimated fluxratios resulted in fractional-masses with higher dynamical mass for the secondary component B compared to the primary A component. We kept the results we obtained from our calculations, but highlight the caveat of the method not being fully reliable for our target sample, mainly serving as a first-order method. In order to disentangle the two lines for more precise estimates require more refined methods, for example tracing back individual RVs (e.g. Czekala et al. 2017).

Orbital fitting
For the orbital fitting procedure we fit the relative orbit of the fainter component in the binary, typically denoted as B here, with respect to the brighter A component. We assumed Keplerian orbits projected on the plane of the sky, so that in the chosen formalism the astrometric position of the companion B could be written as: with ω being the argument of the periastron, θ the true anomaly, Ω the longitude of the ascending node and i the inclination. Here, r = a(1 − e 2 )/(1 + e cos θ) is the radius with a being the semimajor axis and e the eccentricity. The orbital fits were then performed using the observed astrometries to derive the most likely seven orbital parameters; a, e, ω, Ω, i, period P and time of periastron t p .
In order to derive and constrain orbital parameters for our target binaries we applied two complementary approaches and codes. Initially, we applied a grid-search of the seven orbital parameters, the same as described in Calissendorff et al. (2017) and (Köhler et al. 2008(Köhler et al. , 2012(Köhler et al. , 2013(Köhler et al. , 2016. The procedure determined Thiele-Innes elements for points in a grid by solving a linear fit to the astrometric data utilising singular value decomposition. The grid-search was then repeated until a minimum was found and refined for a smaller grid step size. The best-fitted parameters were then determined by comparing the reduced χ 2 from the resulting orbit obtained from the fitted parameters and the relative astrometry for the binary components as where 7 is the number of orbital parameters fitted (9 parameters for systems where RV measurements were applied), s and PA are   Table 6 used to constrain dynamical masses for the binaries. The binaries are listed from left to right in the figure as J0437 and J0459 on the first top row; J0532 and J0613 on the second row from the top; J0728 and J0916 on the third row from the top; and J2317 on the bottom row. The shades of grey represent the first, second and third sigma intervals of the values predicted by the fit for each epoch, where sigma is the standard deviation. the separation and positional angles respectively, and σ their uncertainty. The algorithm is based on a Levenberg-Marquardt χ 2 minimisation (Press et al. 1992), and relies heavily on the starting values which may bias certain orbital parameters if given insufficient orbital coverage or poor initial conditions. We stress that a low χ 2 ν -value is not necessarily a good indicator for a good orbital fit by itself, rather that the measured astrometric data points are well fitted to the calculated orbit. As such, a high χ 2 ν value is not inevitably a bad orbital fit, but could indicate for the uncertainty in the measured astrometric data points to be underestimated, and therefore also the uncertainty in the derived orbital parameters and resulting dynamical mass estimate. In order to address the underestimation of the uncertainty we scaled the astrometric errors for orbits in the grid by χ 2 ν and refitted the orbits while ensuring that the χ 2 ν was equal to 1. For the second approach for constraining the orbital parameters we employed an Monte-Carlo Markov Chain (MCMC) Bayesian analysis technique (Ford 2005(Ford , 2006, the same as used in Rodet et al. (2018). From the MCMC code we obtained the PDFs for the parameters. A sample of 500,000 orbits were randomly picked following the convergence criterion of the applied Gelman-Rubin statistics in the fitting. The sample was assumed as representative of the PDFs of the orbital elements given the initial priors, which were chosen to be uniform in p = (ln a, ln P, e, cos i, Ω + ω, ω − Ω, t p ). In this way, any orbital solution with the couples (ω, Ω) and (ω + π, Ω + π) yield the same astrometries, and thus the algorithm fits Ω + ω and ω − Ω to avoid this degeneracy. Nevertheless, due to how the (Ω, ω) pair is defined, some degeneracy remains and the MCMC occasionally found two families of solutions for some orbits, which is interpreted as a ±180 • ambiguity by the routine. The actual uncertainty is more centred upon the probability peak for each family of solutions. Therefore, for systems which lack RV and are subjected to this degeneracy we cut 90 • around the most probable peak and computed the error around that single interval, allowing for a well-defined uncertainty when the distribution is clearly peaked. The introduction of RV breaks the degeneracy of the (Ω, ω) couple, so that unique values for these variables could be derived for the systems which had sufficient RV data.
The observations with associated astrometric measurements gathered in this work are listed at the end of the paper in the Appendix A, where s is the separation between the binaries in mas, PA the positional angle in degrees. In the table we also listed the deviation between the orbital fit from the grid-method and the observations as |∆s|/σ s and |∆PA|/σ PA , which were calculated as (obs − fit) 2 /σ obs . The χ 2 is then related as the sum of the squares of the deviations.
The resulting orbits from the grid-search orbital fitting procedure are shown in Figures B.1 -B.20 together with their associated best-fit parameters, and the 68 % confidence interval around the probability peak for the MCMC. The astrometric measurements are included in the figures as black dots, with grey ellipses representing their associated uncertainty at the 1-σ level before scaling, and blue lines connect their expected positions from the fit. Most epochs are labelled with their date of observations, but some plots feature less explicitly spelled out dates to avoid cluttering the figure. The semi-major axis a and total system mass M s are listed in both AU and solar masses, as well as in units of mas and mas 2 /yr 3 in order to give the values without distance measurements and uncertainties incorporated. The (ω, Ω) couple which have confidence intervals defined with a ±180 • degeneracy for systems that lack RV measurements are marked with an asterisk ( * ). The results from the MCMC with associated proba-bility peaks and orbits are given in Appendix B.1 for the previously unpublished constraints of the J0613, J0916 and J232617 systems.
We noted that the minimised χ 2 ν -value was not sufficient alone to objectively disclose the robustness of a given orbital fit. Instead, we adopted a custom-made grading system loosely based on the orbital fit grading criteria from Worley & Heintz (1983) and Hartkopf et al. (2001) in order to quantitatively assess how well-constrained each orbital fit was. A direct comparison to the grading criteria from Worley & Heintz (1983) could not be made, as for example it does not account for the additional RV data we possessed for some of our systems, nor did we weight the observations or literature epochs when estimating our grades. For each orbit we calculated a value from a linear combination of the largest gap for the positional angle and phase 1 coverage for the observed epochs of the orbit, divided by the number of orbital revolutions and observed epochs. The values for all systems were sequentially scaled between 1 and 5 from lowest (best, J0008) to highest (worst, J0611) value, effectively creating five equally sized bins. The details of the calculations are outlined in Appendix C, with final grades presented in Table C.1. Because the grid-search method included more astrometric data points and to be consistent with the method used for the main results, we employed this strategy for the orbits obtained from the gridsearch approach only, and adopted the following grading scale: Reliable (1), the orbit is well-constrained; good (2), only minor changes to the orbital elements are expected; tentative (3), no major changes are expected for the orbit; preliminary (4), substantial revisions for the orbit are likely; indetermined (5), the orbital elements are not reliable or necessarily approximately correct.
One of the important factors allowing for constrained dynamical masses in our survey is due to to the updated parallaxes from the Gaia mission, as the distance to the system is typically the main contribution to the uncertainty. With that in mind, it is worth noting that the current Gaia data releases are not yet optimised for handling binaries as they are not photometrically resolved, and further improvements are expected to follow in future releases. The impact of the Gaia parallax measurements is more thoroughly explained for the individual systems in Appendix D.

Results and discussion
We were able to derive individual dynamical masses for the binary components for seven systems in our target sample. Furthermore, we were able to procure luminosities from the resolved observations and age estimates from their adopted YMG membership for the respective systems, which together with their dynamical masses we compared to pre-main sequence (PMS) evolutionary models, probing their accuracy.
Although many evolutionary models exist in the literature, here we adopted the evolutionary models from Baraffe et al. (2015, hereafter BHAC15), as they are well-suited for lower mass stars and younger ages, reflecting our target sample. The dynamical mass estimates for the binaries with individual masses are plotted against stellar isochrones from the BHAC15 models in Figure 3 in a mass-luminosity diagram. The individual mass data-points with distances to the respective system, corresponding absolute magnitudes, their approximate associate age, and estimated theoretical mass from the BHAC15 models are listed in Table 8. The age-ranges listed in the Table 8 show the ages of the isochrones used to calculate the theoretical mass, and not the given age-range of the respective YMG shown in Table 1. The absolute magnitudes were calculated from the unresolved 2MASS K-band magnitudes of the systems, together with the Kband flux-ratios from our SPHERE observations or from previous SINFONI observations (Calissendorff et al. 2020). We therefore prescribe K name convention for the absolute magnitude in the fourth column of Table 8 in in order to highlight this difference.
Overall we found a good consistency between the dynamical mass estimates and the theoretical mass from the models in Figure 3. Most systems have dynamical masses which correspond well with the ages from their respective YMG according to the isochrone tracks. We found two outliers from the prediction of the isochrones in the mass-magnitude diagram, J0459 and J0532. The space velocities for J0459 does not suggest it to belong to any known YMG or association, albeit the binary components are placed amongst the younger isochrones at around ≈ 40 Myrs in the mass-magnitude diagram in Figure 3. Nevertheless, the orbit for the system is not constrained to such a degree that stringent dynamical masses could be procured, and it is likely that our estimate is low. A higher mass would bring the system more in line with an older age. For J0532 we found a tentative orbital fit, displaying some degeneracy in the period-separation space which caused large uncertainties on the dynamical mass. Even with the large formal error bars we found a higher dynamical mass compared to what was expected by the model isochrones for the given young age. This discrepancy is likely attributed to the uncertain method of using the flux-ratio for the RV signals to derive the mass-ratio, but could also potentially be explained by the the individual components being unresolved binaries themselves, causing the observed source to appear underluminous for its mass. Nevertheless, such unseen companions would have to be in close-in orbits of less than one AU to avoid detection from our SPHERE observations. The RV data for the J0916 system aided in constraining the orbit, and was consistent with a Keplerian orbit. However, the flux-weighted approximation for inferring mass-fraction did not seem to apply for this particular system, where we obtained a mass-fraction 0.5 for the fainter companion. As such, we omitted the system from the mass-magnitude altogether. The flux-weighted approach also seemed dubious in other instances as well, including J0437, and to some extent J0728 where the mass-fraction suggested a slightly higher companion mass than primary, but well within the 68 % confidence interval and consistent with the results from Rodet et al. (2018) as well as with the isochrones, confirming the inferred age-range for the system belonging to the AB Doradus moving group of 120-200 Myrs. We did not find the same mass-discrepancy as Rodet et al. (2018), however, our method of testing dynamical against theoretical mass only consisted of the K−band and the age range 120 − 200 Myrs for a single theoretical model. For a younger age of 50 Myrs we obtain the same discrepancy of ≈ 15 % missing mass, with the models underpredicting the total mass of the system compared to the dynamical mass. Rodet et al. (2018) explored the possibility for an unseen companion to explain the missing mass, which could remain hidden and stable if it were closer in than 0.1 AU from one of the other components. From our flux-weighted RV analysis of the system we obtain a slightly higher mass-fraction for the secondary component than the primary, which is also something that would be expected from a higher order hierarchical system such as a triplet. Resolved spectroscopic data could potentially discover such a companion.
No RV data were available for the J1036 binary that could aid the orbital fit. Regardless, since the system is a known triplet we could take advantage of the orbit of the outer binary pair around their common centre of mass along their path on their common orbit around the primary A component in the system. This was previously done in Calissendorff et al. (2017) and we adopt their results of the outer binary being of equal mass. Given the new parallax measurements for the system provided by the Gaia EDR3, the uncertainty in the distance to the system was reduced. In turn, our dynamical mass estimate for the system is the most robust in our sample. The new dynamical mass from the improved orbit and distance is also well aligned with the prediction from the theoretical model isochrones, thus also abolishing the former discrepancy reported by Calissendorff et al. (2017).
The recent efforts by the Gaia mission made it possible to provide updated dynamical masses and absolute magnitudes for some of these systems which previous distance measurements remained somewhat uncertain. In principle, Gaia astrometry could also help to further constrain the orbits of binaries from exploiting the instantaneous acceleration in proper motions, for example by comparing to Hipparcos measurements (e.g. Calissendorff & Janson 2018; Brandt 2018; Brandt et al. 2019). However, only one system in our sample, J0728, exists in both the Gaia and Hipparcos catalogues, and the baseline of 24.5 years between Hipparcos and Gaia far exceeds the orbital period of the system of ≈ 7.8 years, causing some degeneracy when including proper motion acceleration into the orbital fit. Nonetheless, future data releases from Gaia may provide more advantageous information for binary systems with tentative orbital constraints that are unresolved by the space telescope itself. At the same time, one of the troubles for Gaia is to properly identify the photo-centre for unresolved binaries (Lindegren et al. 2018), which the ground-based astrometry could remedy to some extent.
Some systems exhibit a discrepancy in mass when compared to the evolutionary models. A portion of this may be attributed to undetected close-in companions, and given our sample of 20 Mdwarf binaries with the expected multiplicity frequency of ≈ 25 % and companion rate ≥ 30 % (Winters et al. 2019), it is likely that some of these systems contain yet unknown companions. Furthermore, systems which contain more mass, and thereby have faster orbits, are of notable interest for orbital monitoring programmes as the orbits can be mapped out in relatively shorter period of time, and the sample could be afflicted by a selection bias due to this. Observations with for example interferometry that can achieve smaller angular resolutions could place more stringent constraints on the parameter space for a potential visual companions, while spectroscopic monitoring could rule out spectroscopic companions at very low separations.

Summary and conclusions
We considered 20 systems of astrometric M-dwarf binaries for which we present over 75 previously unpublished astrometric data points from AstraLux Norte/CAHA, AstraLux Sur/NTT, FastCam/NOT, NaCo/VLT and SPHERE/VLT. The new astrometric data allowed us to constrain Keplerian orbits and derive dynamical masses for the binaries. We constructed our own relative scale for grading the orbital fits in an attempt to obtain a more objective assessment of the results, with improving grades based on the positional angle and phase coverage as well as revolutions made and number of epochs observed.  modest sample of 20 binaries used for our grading criteria, is likely to be skewed and not directly comparable to the orbital grading system demonstrated on the 900 binaries by Worley & Heintz (1983). The grades provided some quantitative measurement of the orbital fit as a whole, but did not always reflect the uncertainties of the individual orbital parameters or dynamical mass estimate. For example, the orbit for J0459 was labelled as grade 4 (preliminary) but showed low errors for the dynamical mass. Thus, the grades can be interpreted as how well the orbital fit can be trusted, and whether we expect small or large improvements to be made for the orbit, not how stringent the uncertainties are. We summarised the information on the orbital period, semi-major axis, total dynamical and theoretical masses from both orbital fitting methods, along with grades for each orbit in Appendix E. This was the first time orbital constraints were attempted and reported for 14 of our targeted systems. We found good orbital constraints for three systems that had not previously been published, J0613, J0916 and J232617. The PDFs for the orbital parameters obtained from the MCMC fitting procedure for three of these systems together with their associated orbit predictions and mass-distributions are presented in Appendix B.1.
Six systems, J0008, J0437, J0532, J0728, J1036 and J2317 all had previous orbital parameters, and our additional data mainly confirmed the earlier results, with a slight improvement for J0437 and J0728 with updated Gaia parallaxes. J1036 and J2317 had previously more uncertain dynamical mass estimates due to the lack of good distance measurements, which we redressed here in this work. The improved dynamical mass estimates for J1036 and J2317 now stand among the most robust in our sample. The remaining 11 systems still require further monitoring to provide reliable dynamical masses, and our new data pave the way for future orbital constraints for these systems.
Out of the 20 binaries in our sample, the four systems J0008, J0728, J1036 and J2317 received reliable orbital constraints. We determined tentative orbits for six systems; J0111, J0245, J0907, J0916, J1014 and J232611, which are expected to yield more ro-bust orbital parameter constraints in the near future if additional astrometric measurements can be procured. Six systems, J0459, J0532, J2016, J2137, J232611 and J2349, only had their orbits constrained to a preliminary level, which may help to indicate an approximate period to some degree, while most of the orbital parameters are yet too uncertain to place robust dynamical masses. For the J0225 and J0611 systems the orbits were completely undetermined.
We searched the literature for RV data for the target sample presented here, finding RV measurements that were consistent with Keplerian orbits for seven of our binaries. Since the binaries were not of pure SB1 single lined spectroscopic binaries we assumed a flux-weighted method to infer individual dynamical masses. The method proved unreliable for the J0916 binary, and dubious at most for J0437. We therefore argue that only their total dynamical masses are reliable, and that a different approach is necessary to disentangle the individual lines and masses.
When we compared the derived individual dynamical masses with theoretical masses from the BHAC15 model isochrones we found an overall good consistency between the empirical and theoretical masses. The largest discrepancy was found for J0459, for which the isochrones predict ages between 10-50 Myrs, but Bayesian membership probabilities suggest the system more likely to be belonging to the field and no known young association or group, thus expected to be much older. We could explain this discrepancy from the degeneracy in the orbital fit, and it is possible that the period is overestimated while the semi-major axis underestimated. The near-IR spectral type for the primary being an earlier M1.7 ± 0.6 does also suggest the system to be older and more massive than what our derived dynamical-mass advocates. However, spectral analysis by Calissendorff et al. (2020) shows a discrepancy between the near-IR bands, where the J−band suggest a lower surface-gravity, and thereby younger age, compared to the other bands. Our gridsearch method suggested a greater mass and confidence interval than the MCMC orbital fit for this system, and we are likely to see further improvements to the orbit with just a single more as-trometric data point in the future, which may also aid to constrain the age of the system through isochronal dating.
Out of our sample of 20 low-mass binaries, eight have strong indicators of being young and members of YMGs and associations (excluding J1036 which has uncertain Ursa-Majoris affiliation and it is questionable whether the age estimate of ∼ 400 Myrs can be considered young in this context). These systems will continue to prove to be important calibrators for evolutionary models, and as orbital monitoring continues, better estimates for important formation diagnostics such as semi-major axis and eccentricity distributions will be acquired. With the increasing sample size of young binaries with stringent orbital parameter constraints, we will soon achieve a set of empirical isochrones, which can then be utilised to evaluate more precise ages of nearby young moving groups.     The age ranges listed were used to calculate theoretical masses in the last column predicted by the models for the given absolute magnitude in 2MASS K-band. The magnitudes listed in the fourth column are derived from the unresolved 2MASS K-band magnitudes of the system and the flux ratios in the SPHERE or SINFONI K-bands, and the name convention of K is applied to mark this difference. For J0437 we used the flux ratio from the Keck/NIRC2 K−band observations from Montet et al. (2015).
Article number, page 14 of 38           and theoretical masses closer together, as well as be more in line with the expected mass from the derived spectral types.
2MASS J02451431-4344102 is a likely field binary according to all three YMG tools applied. We mainly sampled the orbit around the periastron and despite the small uncertainty in dynamical mass we obtained a grade 3 (tentative) for the orbit, which has more than half of its path yet uncharted. The MCMC method did not include the astrometric data points from 2019-2020 which explains the difference in orbital parameters determined by the two methods. The RV data in the literature (Durkan et al. 2018) had too short baseline to aid the orbital fitting and does not match the astrometry in this case and was therefore excluded from the fitting procedure. The spectral types are similar to that of the J0008 system which had a reliable orbital fit, also belonging to the field and having similar dynamical mass. This could suggest that the orbit we obtained for J0245 is actually better constrained than expected from our grading criteria.
2MASS J04373746-0229282 has have its orbit constrained previously by Montet et al. (2015), obtaining a dynamical mass of M = 1.11 ± 0.04 M . We found a good orbital fit for the system, obtaining a grade 2 on our scale, and a slightly lower mass of M = 0.92 +0.06 −0.04 M . This discrepancy is explained by the different distance measurements adopted for the system, where Montet et al. (2015) resorted to using the Hipparcos distance to the comoving system 51 Eri of 29.43 ± 0.30 pc (van Leeuwen 2007) while we had access to the Gaia EDR3 parallax corresponding to a distance of 27.77 ± 0.37 pc. The RV measurements for the system allow for individual dynamical masses to be derived for the system, where we scaled the flux according to the relative brightness of the components. However, the flux ratio between the two components show additional discrepancies. Indeed, we obtainad a ratio of F B /F A = 0.03 in the i -band (Janson et al. 2014a) corresponding to a mass-ratio of M B /M tot = 0.44. Given the magnitude difference displayed in other bands (e.g. Figure 2 in Montet et al. 2015), the flux ratio is more likely to be F B /F A = 0.2 ± 0.1, which translated to a mass-ratio of M B /M tot = 0.62. As a comparison, Montet et al. (2015) found a mass-ratio of ≈ 0.4. However, they find that the secondary B component is missing mass when compared to evolutionary models, which could be attributed to an unseen companion. Such a companion could also be the cause for the oddity in massratio being higher for the secondary than the primary component that we see. While all three of the YMG tools suggest different groups for the system, shown in Table 3, the LACEwING tool also produced a ≈ 37 % probability for BPMG membership, with most of the literature agree that the system is a member of the BPMG.
2MASS J04595855-0333123 received a grade 4 (preliminary) from our orbital fitting. The MCMC orbital fitting did not exclude the low probability of a high-mass system, and we cut the mass-distribution at ≤ 2 M . The BANYAN Σ-tool suggests the system to be in the field while LACEwING proposes the system to be a HYA member, which regardless places the system as old in comparison to the younger systems in our sample. Although the system is likely belonging to the field, the isochrones in the theoretical models predict a much higher mass for the system, which further alludes to the poor orbital constraints and that more information is necessary to obtain a better mass estimate. The models could however explain the lower dynamical mass if the age of the system was ≤ 50 Myrs. Nevertheless, the optical spectral types are similar to that of J0008 and J0245, which should indicate for similar masses given their approximate ages. The near-IR spectra on the other hand implies an earlier spec-tral type for the primary state, and its mass could potentially be heavily underestimated.
2MASS J05320450-0305291 is part of a higher hierarchical sextuple system (Tokovinin 2022) and a strong candidate of being a member of the BPMG according to the BANYNA Σ-tool, but also a potential member of ABDMG according to the convergence point tool. The orbital fit is still preliminary according to our grading scale, receiving a grade of 4 (preliminary), with most of the orbital phase not yet observed. Despite the poor orbital constraints, we obtained low minimised χ 2 ν, grid = 1.3 and χ 2 ν, MCMC = 1.6, suggesting that the astrometry is a good for the orbital fit. However, we find some degeneracy and the estimated period ranges from 20 to 100 years, where the shorter periods would suggest dynamical masses above 2 M , and thus dubious. A longer period of ∼ 80 is more consistent with the results from Tokovinin (2022), where they provide two potential orbital solutions for either P = 80 or P = 143 years. We did not include the astrometric epochs from (Tokovinin 2022) in the MCMC fit which explains the discrepancy compared to the grid model.
2MASS J06112997-7213388 received the worst grade of the orbital fits in our sample, and is assigned as undetermined. The system is however likely young, with its highest probability being a CAR member according to the BANYAN Σ-tool, with LACEwING suggesting either a COL or CAR member with the same probabilities. The convergence point tool however suggests the system to be a CARN member, which is likely just because the tool does not include CAR in its calculations. The spectral types derived from optical photometry and near-IR spectra are consistent with each other. The secondary component is estimated to move close to its apastron in its coming years and expected to show limited motion in its orbit.
2MASS J06134539-2352077 is a likely ARG member, supported by both the BANYAN Σ and LACEwING tools. The astrometric data from the epochs between 2019-2020 were only included in the grid-search approach, which explains the discrepancy in the orbital parameters obtained for the two methods and the different masses obtained. However, both procedures obtained masses greater than that of the evolutionary model, suggesting that perhaps there is some missing mass and an unseen companion in the system. The SPHERE observations should have been able to detect massive companions of ∼ 0.2 M , or at least a non-uniform PSF, down to ∼ 0.5 AU, which is not apparent from the only epoch taken with SPHERE for the system. The orbit from the grid-search method obtained a grade 2 (good) on our scale.
2MASS J07285137-3014490 is a well-studied system in the ABD moving group, for which we mainly contributed by adding additional astrometric data points and an updated distance parallax compared to the previous results in Rodet et al. (2018). The updated Gaia EDR3 parallax measurement helped to constrain the distance uncertainty to the system by a factor of ≈ 4, resulting in a more precise mass estimate. Our results were consistent with those of Rodet et al. (2018) for most part, with the exception of a slightly higher mass-fraction in our case of M B M tot = 0.51 +0.10 −0.08 compared to M B M tot = 0.46 ± 0.10 inRodet et al. (2018), where we applied the same flux ratio of 0.2 ± 0.01. The missing mass in the secondary B component that Rodet et al. (2018) found is accounted for in our case, which potentially could be caused by the different distances adopted, but the discrepancy in massfractions we found instead could also be explained by the same argument of an unseen companion. The system also possesses a notably high eccentricity of e = 0.90 which points towards significant dynamical interactions. Nevertheless, the existence of a