Dynamical Mass of the Ophiuchus Intermediate-mass Stellar System S1 with DYNAMO-VLBA

We report dynamical mass measurements of the individual stars in the most luminous and massive stellar member of the nearby Ophiuchus star-forming region, the young tight binary system S1. We combine 28 archival data sets with seven recent proprietary Very Long Baseline Array (VLBA) observations obtained as part of the Dynamical Masses of Young Stellar Multiple Systems with the VLBA project (DYNAMO–VLBA), to constrain the astrometric and orbital parameters of the system, and recover high-accuracy dynamical masses. The primary component, S1A, is found to have a mass of 4.11 ± 0.10 M ⊙, significantly lower than the typical value ∼6 M ⊙ previously reported in the literature. We show that the spectral energy distribution (SED) of S1A can be reproduced by a reddened blackbody with a temperature between roughly 14,000 and 17,000 K. According to evolutionary models, this temperature range corresponds to stellar masses between 4 M ⊙ and 6 M ⊙, so the SED is not a priori inconsistent with the dynamical mass of S1A. The luminosity of S1 derived from SED fitting, however, is only consistent with models for stellar masses above 5 M ⊙. Thus, we cannot reconcile the evolutionary models with the dynamical mass measurement of S1A: The models consistent with the location of S1A in the Hertzsprung-Russel diagram correspond to masses higher by 25% at least than the dynamical mass. For the secondary component, S1B, a mass of 0.831 ± 0.014 M ⊙ is determined, consistent with a low-mass young star. While the radio flux of S1A remains roughly constant throughout the orbit, the flux of S1B is found to be higher near apastron.


INTRODUCTION
The Ophiuchus molecular cloud complex is one of the nearest and best-studied sites of active star formation (see review by Wilking et al. 2008a).It has, indeed, played a crucial role in our overall understanding of star formation.Most of the ongoing star formation activity occurs in the Ophiuchus core (also known as Lynds 1688; L 1688), located at 137.3 ± 1.2 pc (Ortiz-León et al. 2017, OLK17 hereafter).The most luminous and massive member of L 1688, and of the entire Ophiuchus region, is the young Herbig-Be star Oph-S1 (also known as S1).Early radio observations of S1 showed that its radio emission originates mainly from two components: a thermal extended (∼ 20 ′′ ) component associated to a compact H II region and a non-thermal compact source embedded within the H II region (Andre et al. 1988).Very Long Baseline Interferometry (VLBI) observations later showed that the non-thermal component is associated with a magnetically active star; S1 was the first young stellar object directly detected with the VLBI technique (Andre et al. 1991).These radio observations, as well as early infrared measurements by Lada & Wilking (1984), suggested that S1 was a young intermediate-mass star of spectral type B3 to B5, with a mass of about 6 M ⊙ .
While the properties of the compact radio emission in S1 are consistent with a magnetically active star, the detection of non-thermal radio emission in an intermediate-mass star is a priori unexpected.This is because the interiors of young intermediate-mass stars (3-8 M ⊙ ) are thought to be fully radiative while they move toward the main-sequence.Under these circumstances, the dynamo effect, which requires a rotating core surrounded by a convective envelope, cannot operate (Stelzer et al. 2005(Stelzer et al. , 2006;;Hubrig et al. 2009).Observational evidence, however, has shown that a significant fraction of Herbig Ae/Be stars and main sequence A/B stars show signs of magnetic activity (Hubrig et al. 2009;Wade et al. 2011;Hubrig & Schöller 2021).A plausible explanation is that the magnetic activity is associated with a lowmass companion star rather than the intermediate-mass star itself (e.g., Stelzer et al. 2005Stelzer et al. , 2006)).In the specific case of S1, lunar occultation observations have indeed revealed the existence of a lower mass companion with a separation around 20 mas (Richichi et al. 1994;Simon et al. 1995).More recently, Very Long Baseline Array (VLBA) high angular resolution observations obtained as part of the Gould's Belt Distances Survey (GOBELINS; OLK17) detected two stellar components in S1.From these data, OLK17 obtained preliminary dynamical masses of 5.78±0.15and 1.18±0.10M ⊙ for the two stars in S1.While these VLBA observations confirmed the existence of a low-mass companion, they also show that both the low-mass companion and intermediatemass star emit non-thermal radio emission.Thus, a radio emission mechanism immanent to the intermediate star re-mains necessary.Possible mechanisms include the existence of a fossil magnetic field, and the temporary apparition of a thin convective layer on the surface of the star (Palla & Stahler 1993;Tout & Pringle 1995;Skinner et al. 2004).
Regardless of the exact emission mechanism, the binarity of S1 provides us with a rare opportunity to dynamically measure the mass of an Herbig Be star.In this paper, we aim at improving on the mass determination of OLK17.We present new VLBA observations of S1, obtained as part of the Dynamical Masses of Young Stellar Multiple Systems with the VLBA (DYNAMO -VLBA) 1 project.To improve the constraints on the mass of the stars in S1, we combine these new observations with previous VLBA data, particularly those obtained as part of GOBELINS.In the rest of the paper, we will call the primary of the system S1A (this is the young Herbig Be star; OLK17 and Loinard et al. 2008, LTM08 hereafter) and the secondary component S1B.

OBSERVATIONS AND DATA REDUCTION
The DYNAMO-VLBA observations (VLBA project code: BD215) used the same correlator setup as GOBELINS to record the radio continuum flux at λ = 6.0 cm (ν = 5 GHz) using the VLBA of the National Radio Astronomy Observatory (NRAO).S1 was observed six times, approximately every four months, from early 2018 to late 2019.During the sixth observation, the antennas on all sites were not pointing correctly for the first half of the observation time, and an additional observation was obtained four days later.Thus, S1 was observed a total of seven times as part of the DYNAMO-VLBA project.The observational details and the measured position are provided in Table 1.
The total duration of each observation was three hours.The first and last 25 minutes were dedicated to observing about a dozen quasars distributed over the entire sky.The rest of the observation was spent on the target source and its gain calibrator (J1627-2427) in cycles of two minutes on-target and one minute on the calibrator.This specific gain calibrator was chosen because it is located in the core of L 1688 at an angular distance of 10 ′ from S1 (Dzib et al. 2013).For all observations, the julian date reported in Table 1 corresponds to the mid-point of the three-hour observing run.
As mentioned earlier, S1 was observed with the VLBA prior to the DYNAMO-VLBA project.Multi-epoch, phasedreferenced observations were obtained as part of VLBA projects with codes BL128 (six epochs), BT093 (eight epochs) as well as the GOBELINS survey with code BL175 (14 epochs).These observations were first reported in LTM08, OLK17, and Ortiz-León et al. (2018) and we include them in our analysis.Thus, this paper considers a total of 35 VLBA observations of S1, distributed between 2005 to 2019.
The initial data calibration followed standard procedures for phase-referenced VLBI observations and was carried out with the Astronomical Image Processing System (AIPS) software (Greisen 2003).The steps have been described in detail in Loinard et al. (2007), Dzib et al. (2010) and OLK17.
The phase calibrator used in the DYNAMO-VLBA and GOBELINS observations is located only 10 ′ from the target source, but its flux density is only about 50 mJy.On the other hand, the observations of BL128 and BT093 used the stronger (∼ 1.0 Jy) quasar J1625-2527 located ∼ 1. • 1 south S1.As a consequence, larger residual calibration errors are expected to affect the final images obtained from BL128 and BT093.S1A is a compact source with a flux density of about 6.0 mJy (LTM08; OLK17) so, given the expected noise levels in our observations, it ought to be detected on individual baselines with a signal-to-noise ratio of order 5 on integrations of two minutes.Self-calibration with a solution interval as short as two minutes can therefore be applied (see Cornwell & Fomalont 1989, for further details).
The self-calibration procedure used the following steps.First, we imaged the externally calibrated visibilities with a few clean components so the image contained only the strongest peak.Second, phase corrections on time intervals of 15 minutes are obtained using this image as a model.Then, the first and second steps are repeated two more times, decreasing the solution time interval to eight and two minutes, respectively.Final corrections in amplitude and phase are calculated, with time intervals of fifteen and eight minutes, using the final image of the previous step as the model.The final image is obtained using these corrections.
After self-calibration, the data were imaged using a pixel size of 100 µas and a natural weighting scheme (robust = 5 in AIPS).The rms noise level in the images (25-80 µJy beam −1 , see Table 1) is consistent with the values expected theoretically.The fluxes and positions of the detected sources were measured using a two-dimensional Gaussian fitting procedure (task JMFIT in AIPS); the full list is presented in Table 1. Figure 1 shows the images of S1 as seen in the seven new epochs from the DYNAMO-VLBA project.
For the observations of projects BL128 and BT093 corresponding to the 14 older epochs, the measured positions were corrected to account for the most recent measured position of the gain calibrator used for these observations, J1625-2527.The positional offset was estimated by OLK17 to be ∆α = −23 µas and ∆δ = −1.2mas.These values were added to the measured positions of S1A and S1B in all epochs corresponding to projects BL128 and BT093, and are already taken into account in the values listed in Table 1.

ASTROMETRIC FITTING PROCEDURE
The displacements of S1A and S1B on the celestial sphere are described by the combination of their common trigonometric parallax π, the uniform proper motion of their center of mass in right ascension and declination, µ α and µ δ , and their orbital motion.For the primary component, the complete equations describing the motions are: (1) where α 0 and δ 0 are the coordinates of the center of mass of the system at the chosen reference epoch, f α and f δ are the projections of the parallactic ellipse (Seidelmann 1992), and a 1 is the semi-major axis of the orbit of the primary around the center of mass.Q α (t) and Q δ (t) represent the projections of the orbital motions and are functions of the orbital elements: where θ is the true anomaly and r is the radius from the center of mass along a Keplerian orbit.The inclination angle is i, the angle of the line of nodes Ω, and the angle from the ascending node to periastron ω.For the secondary component, the semi major axis a 2 is used instead of a 1 in equations ( 1) and ( 2) and θ is rotated 180 • in equations ( 3) and ( 4) (see e.g.Kounkel et al. 2017 for details).
To solve the previous equations we translated the procedure described by Kounkel et al. (2017) from IDL to python2 .This routine simultaneously fits equations (1) to (4) to the data by the method of least squares using MPFIT.It is able to simultaneously process absolute positions of both stars from VLBA astrometry, as well as relative positions of the secondary with respect to the primary from optical data if they exist.The free orbital parameters in this model are the period P, the semimajor axis of primary a 1 , the eccentricity e, the inclination i, the angle of the line of nodes Ω, the time of periastron passage T , the angle from node to periastron ω, and the mass ratio m 2 /m 1 .The total mass of the system is derived from Kepler's third law and the distance -obtained from the very same fit through the trigonometric parallax.
This routine is sensitive to the initial guesses of several non-linear orbital parameters, namely P, e, ω, and Ω; with poorly chosen initial guesses the fitter may not converge, getting stuck in a local minimum instead.The remaining parameters are not affected by this, and can always be scaled.Thus, to fully probe the parameter space, we performed the fit 1000 different times with different initial guesses.The routine then iterates on these guesses, enabling optimal convergence for several dozens of these sets of iterations, which are then evaluated using χ 2 .These sets of iterations with χ 2 < 1.1χ 2 best are selected (where χ 2 best is the minimum chi squared of the entire set.The uncertainties are evaluated in two ways: from examining the intrinsic scatter from all of the valid solutions, and from calculating a weighted average error of all of the errors returned by MPFIT, weighted by χ 2 .
Typically, both methods produce comparable distributions.

RESULTS
The S1 system was detected in all seven new VLBA observations (Figure 1).S1A was detected in all these epochs, with a mean flux density of 7.59 mJy; S1B was detected five times with a mean flux density of 0.85 mJy.In the ancillary observations, thanks to the improvement obtained by applying self-calibration, we have detected S1B in more epochs than  1.Images are centered on the measured position of S1A, the primary component of the S1 binary system as listed in Table 1 and shown in the top-left corner of each plot.The white ellipse at the bottom-left of each plot represents the synthesized beam of the corresponding observed epoch.
previously reported.In addition, we found that a previous detection of S1B (BL128 GE from 2006 June 3) was spurious: the emission peak at that position in fact corresponded to a prominent side-lobe.We do detect S1B in the self-calibrated image of the same epoch but at a different position.Therefore, the measured S1B position on 2006 June 03 by OLK17 is discarded in our analysis and we use, instead, the newly measured position from the self-calibrated image.S1A was detected at all 35 epochs -a 100% detection rate was also achieved in the publications reporting on older observations (LTM08; OLK17, Ortiz-León et al. 2018).On the other hand, the companion, S1B, is detected in 14 epochs, of which only four detections were reported by OLK17.All the measured positions and flux densities for S1A and S1B are listed in Table 1.

Astrometry
The final measured positions of S1A and S1B were used as inputs to our astrometric fitting procedure.The resulting parameters from our best fit are listed in Table 2 and are plotted in Figures 2 (total sky motion) and 3 (orbit motion).The total and individual masses are derived from Kepler's third law, taking into account the obtained values for the period, semi-major axis, and the mass ratio from our model fit.As a result, we obtain masses for the primary and secondary components of 4.11±0.10M ⊙ and 0.83±0.01,respectively.It is important to note that these mass estimates are somewhat at odds with those reported by OLK17: 5.78±0.15M ⊙ for the primary component and 1.18±0.10M ⊙ for the secondary component.This discrepancy arises from two main factors.Firstly, we have significantly increased the number of S1B detections, yielding more robust and accurate measurements.Secondly, we have corrected the position of S1B at epoch 2006 June 3rd (JD 2453889.789).As mentioned earlier, the detection of S1B reported by OLK17 turned out to be spurious, and the new position is significantly different from the previous one.It is worth mentioning that OLK17 based their fits on only four detections of S1B, so the inclusion of a spurious position had an especially large detrimental effect.The addition of seven new, high-quality data, obtained as part of DYNAMO-VLBA was particularly important to constrain the orbital fit.From our study, we conclude that S1A is indeed an intermediate-mass star, but less massive than previously suggested.
We also fitted the data with the MCMC Orbitize!package (Blunt et al. 2017) to assess the accuracy of the orbital parameters obtained from the MPFIT fitting.This algorithm explores the multidimensional parameter space, effectively sampling the posterior distribution of the orbital parameters.The results obtained demonstrate remarkable consistency with the values derived from the MPFIT fit, confirming the reliability of our orbital parameter estimates.A detailed description of the MCMC analysis and its results is presented in Appendix 6.

Flux variations
Figure 4 shows the flux of S1A and S1B as a function of the orbital phase, defined as the fractional part of the equation where T 0 is the time of periastron (JD 2457163.69),and P is the orbital period -both taken from the fitting results given in Table 2.
The projects BL128 and BT093 were recorded at a frequency of 8 GHz, while the GOBELINS project used 5 GHz, except for the first epoch which was observed at 8 GHz.As mentioned in Section 2, the DYNAMO-VLBA observations were also recorded at 5 GHz.To examine flux variations, we  1).The red and blue curves show the best fits, described in the text, to the positions of S1A and S1B, respectively.have to account for the fact that the data have been observed at two different frequencies.For S1A from the individual measured fluxes, we calculated the mean flux densities at both frequencies: 5.49 ± 0.45 mJy at 5 GHz and 6.95 ± 0.57 mJy at 8 GHz.This is roughly consistent with a flat spectrum over these frequencies, as also seen in previous VLA observations (Dzib et al. 2013).Then, the flux values measured at 8 GHz were scaled to the 5 GHz values by adding the difference between the mean value at both frequencies.
For S1A (top panel of Figure 4), there are no obvious variations with phase in the individual flux measurements.To search further for evidence of systematic flux variations, we averaged the data in bins of width 0.2 (in orbital phase).The horizontal black lines and vertical gray bars in the top panel of Figure 4 indicate the mean and standard deviation in each bin, respectively.From this analysis, we do not find any evidence for variation with orbital phase.Variability at the level of about 50% is present in the data, but it appears to be stochastic, with no relation to the orbital phase.
For S1B, on the other hand, detections (particularly associated with higher fluxes) appear to be clustered around an orbital phase value of 0.5 (bottom panel of Figure 4).More quantitatively, S1B was detected in 9 out of 14 observations (64%) carried out in the range of orbital phases between 0.4 and 0.6, while it was detected in only 5 out of 21 obser-   The measured positions of S1B obtained from infrared observations reported in Richichi et al. (1994) and Cheetham et al. (2015) are also plotted, as red triangles.These infrared positions are listed in Table 4.We can examine the significance of the detection rates in each orbital phase range by modeling the detectability of S1B with a binomial distribution.In this statistical model, if the events are independent and the probability of a success in a given experiment is p, the probability of k successes in n experiments is given by: -30 -20 -10 0 10 (mas) S1B is detected in 14 out 35 experiments, so the overall probability of detection (a success) is p = 0.4.Given this probability and assuming that the chance of detection is independent of orbital phase, there is only a 4% chance that 9 detections would be obtained out of 14 observations (the results obtained near apoastron), while there is only a 6% chance that only 5 detections would be obtained out of 21 observations (the results obtained near periastron).We conclude that the higher detection fraction near apoastron than periastron is very unlikely to be the result of chance alone and that the source is intrinsically more likely to be detectable near apoastron.

DISCUSSION
In this section, we mostly discuss the implications of our mass measurements for the properties of the S1 system.It is worthwhile mentioning at the outset that, to the best of our knowledge, no evidence for cicumstellar (or circumbinary) disks has ever been reported for this system.

Spectral Energy Distribution of S1
Before dynamical measurements became possible, the mass of S1 had been estimated from its photometry, and some spectroscopy, to be about 6 M ⊙ .Let us briefly review these works.Based on early infrared observations, Grasdalen et al. (1973) classified S1 as a late B star.It was later classified as B3V by Elias (1978) due to the surrounding HII region causing strong CO absorption in the 2.3 µm band.Combining spectroscopy and photometry, Cohen & Kuhi (1979) obtained a B2 classification.Andre et al. (1988) estimated a B3.5 spectral type on the basis of the radio flux of the surrounding HII region.From spectroscopic observations, Bouvier & Appenzeller (1992) classified it as a B4 star since it exhibits strong H α absorption but no lithium absorption; this same classification is adopted in Andre et al. (1991).Analysis by Lada & Wilking (1984) of the spectral energy distribution implied that S1A has a luminosity of 1,500 L ⊙ and an effective temperature of 16,000 K consistent with a main-sequence star of spectral type B3-B5.Finally, spectroscopy by Wilking et al. (2005a) implied a spectral type earlier than B8 -they adopted a B3 classification based on the effective temperature from Lada & Wilking (1984).The classification we just described implies a mass between 5 and 6 M ⊙ -this is significantly higher than the value we derive dynamically.
Recently, Mookerjea et al. (2018) compiled photometric observations of S1A available through the astronomical catalog service VizieR to plot the SED of S1A and fit it with a reddened blackbody.They adopted a blackbody temperature of 17,000 K, which corresponds to the effective temperature of a B3/B4 ZAMS star, as suggested by the studies mentioned above.They assume a stellar radius R = 4.35 R ⊙ which, combined with the temperature mentioned above, results in a luminosity L = 1, 420 L ⊙ .Their best fit implies a reddening E(B − V ) = 3.8 mag, a selective extinction ratio R V ≡ A V /E(B − V ) = 3.35, and a total V band extinction A V = 12.7 mag.
In order to examine the compatibility between our dynamical mass of S1A and its observational properties, We performed an analysis similar to that of Mookerjea et al. (2018).We included the photometric points currently available in VizieR; which cover the optical and infrared bands from 0.4 µm to 8 µm.The complete list of included measurements is provided in table 6 of Appendix 6.We adopted a standard selective extinction ratio R V = 3.1, and the extinction curves of Fitzpatrick (1999) for the range between 0.4 µm and 1.25 µm and Rieke et al. (1989) for the range between 1.25 µm and 8 µm.Figure 6 shows that the data points can be reproduced by a range of effective temperatures, from about 14,000 K to 17,000 K.The specific values T eff = 14,000 K, 15,700, and 17,000 K shown in Figure 6 were chosen because they correspond to spectral types B5, B4 and B3 according to Pecaut & Mamajek (2013).Values in the lower part of this range are consistent with the dynamical mass of S1 determined here.The extinction derived for the models shown in Figure 6 vary from 10.5 (for T eff = 14, 000 K) to 11.8 (T eff = 17, 000 K). Plots of the residuals between the models and the photometric data are shown in the bottom panels.6) while the fits (color lines) are based on the extinction curves of Rieke et al. (1989) and Fitzpatrick (1999) [See text].The figure also includes residual plots for each model (bottom panels).

Stellar evolution models with no rotation for S1A
Dynamical masses obtained in binary systems are essential to constrain stellar evolution models.Nowadays, this is particularly important for young stars because pre-main sequence evolution models are less reliable and less mature than their main sequence counterparts (e.g., Hillenbrand & White 2004;Stassun et al. 2014).Models for intermediate mass young stars are particularly uncertain given the paucity of existing observational constraints.It is, therefore, worthwhile to compare the observational properties of S1A with the predictions of theoretical models for the mass determined dynamically in the present paper.(Tognelli et al. 2011).These models have a metalicity Z = 0.02, helium abundance Y = 0.2880, mixing-length α = 1.90 and deuterium abundance XD = 4×10 −5 .The blue, red, and green data squares indicate, respectively, stars with Teff = 17,000 K and L = 2,120 L⊙, Teff = 15,700 K and L = 1,430 L⊙, and Teff = 14,000 K and L = 730 L⊙.These correspond to the three SED models, of the same colors, in Figure 6.The grey square shows the locus of a star with Teff = 14,000 K, L = 300 L⊙, which would be consistent with a 4 M⊙ star shortly before reaching the ZAMS.
In Figure 7, we show evolutionary tracks for masses between 4 and 6 M ⊙ , based on the PISA pre-main sequence stellar evolution models presented in Tognelli et al. (2011).The loci of S1A in the HR diagram for the three SED models shown in Figure 6 are indicated in Figure 7 as color points.The luminosity and effective temperature of S1A appear to be consistent with masses between about 5 and 6 M ⊙ , significantly higher than the measured dynamical mass of 4.1 M ⊙ .An effective temperature of 14,000 K could be consistent with a mass of 4 M ⊙ , but only provided that the luminosity was about 300 L ⊙ -instead of the 730 L ⊙ obtained from modelling SED.This could only be achieved if the radius of the star were of order 3 R ⊙ -instead of 4.6 R ⊙ obtained from modelling SED.
For completeness, we have repeated the same analysis with the Y 2 (Yonsei-Yale) stellar models (Yi et al. 2001), the Palla-Stahler models (Palla & Stahler 1999), the PARSEC models (Nguyen et al. 2022, see below), and the YaPSI YALE-POTSDAM stellar models (Spada et al. 2017).We obtained similar conclusions to those reached with the PISA models: the theoretical predictions for a mass of 4 M ⊙ are inconsistent with the locus of S1A in the HR diagram.Specifically, the luminosity of S1A is at least twice larger than predicted by the models.Conversely, given the locus of S1A in the HR diagram, the theoretical models would predict a mass for S1A larger than measured dynamically, by 20 to 50%.Inter-estingly, Stassun et al. (2014) also reported systematic discrepancies between masses predicted by evolutionary models and measured in binary systems, but for stars with masses below 1 M ⊙ .For masses above 1 M ⊙ , the model predictions and observed masses were found to be in agreement at the level of about 10%.We note, also, that previous dynamical mass measurements of intermediate-mass stars (e.g., Johnston et al. 2019) did not find a significant discrepancy with evolutionary models.
To finish with this section, it is worth emphasizing that the location of S1A above and to the right of the main sequence on the HR diagram implies that it has not yet reached the main sequence.Since intermediate mass stars have a very short pre-main sequence lifetime, this places stringent constraints on the age of the S1 system.For instance, according to the evolutionary models of Tognelli et al. (2011) shown in Figure 7, a 5 M ⊙ star reaches the main sequence in just about 1 Myr.The possible location of S1A along the track of a 5 M ⊙ star shown by the green square of Figure 7 would indicate an age of 0.7 Myr according to the same models.Of course, we should proceed cautiously with obtaining age estimates from the models since we have just shown that the models cannot reproduce the location of S1A on the HR diagram given its measured dynamical mass.Nevertheless, it seems reasonable to conclude that the age of S1A is likely between 0.5 and 1 Myr.

Stellar evolution models with rotation for S1A and magnetic activity
The evolutionary models considered in the preceding section do not include rotation.To investigate the possibility that rotation might help resolve the discrepancy between the measured dynamical mass of S1A and the evolutionary tracks, we considered the PARSEC V2.0 evolutionary models presented by Nguyen et al. (2022).These models include recent re-calibrations of mixing by overshooting and rotation.They parametrize rotation via the ω parameter defined as the ratio between the stellar angular velocity Ω and the breakup angular velocity Ω c (Ω c is the angular velocity for which the centrifugal force becomes equal to the effective gravity at the equator).
Figure 8 shows the PARSEC pre-main sequence tracks for a 4 M ⊙ star of Solar metallicity, appropriate for S1A, for ω values between zero (no rotation) and 0.99 (near breakup).As can be seen, the effect of rotation is imperceptible until the star approaches the ZAMS.There, the net effect of rotation is to make the star redder and slightly less luminous.Thus, rotation does not help reconcile the models with the dynamical mass of S1A since the models without rotation produced a luminosity that was already lower than expected for the dynamical mass of S1A.The PARSEC models account for rotation but do not include magnetic fields and this may be an important missing element for the interpretation of S1A since the very existence of non-thermal radio emission requires the presence of a substantial superficial magnetic field (10 2 -10 3 G; Andre et al. 1988).The origin of this field is unclear.In lowmass stars, superficial magnetic fields of this magnitude are expected from dynamo theory, but this requires the presence of deep convective envelopes (Parker 1955).In intermediate mass stars such as S1A, however, the energy transport from the core to the surface is radiative, including during the premain sequence phase (Kippenhahn & Weigert 1990).Several scenarios have been proposed to account for the presence of strong superficial magnetic fields in intermediatemass stars.Perhaps the most accepted possibility is a fossil origin (Schleicher et al. 2023).An alternative is the generation of magnetic field through the so-called Tayler-Spruit dynamo mechanism in a rotating star (Nandal et al. 2023).Regardless of the origin of the magnetic field, its effect on the early evolution of intermediate-mass stars remains poorly investigated.

The nature of S1B
Our mass measurement and analysis of the SED of S1A implies that S1B is a young (age of order 0.5 to 1 Myr, assuming co-evality with S1A) Solar type star (M ∼ 0.8 M ⊙ ).It is interesting to see if this is consistent with the few existing resolved infrared observations of S1A.The fluxes used in Figure 6 and listed in Table 4 correspond to a K band apparent magnitude of about 6.3 -this is for the entire system, but we can also take this value for S1A alone since S1B is much weaker.Richichi et al. (1994) find a magnitude dif-ference ∆K = 1.6 between S1A and S1B, while Cheetham et al. (2015) report ∆K = 1.2 and 2.3 in their two separate observations.This results in an apparent K-band magnitude for S1B between 7.5 and 8.6.The distance modulus for S1 is µ = 5.7 and the A V extinction of 12.7 reported from our fit to the SED of S1A corresponds to A K = 1.4.Thus, we can estimate an absolute K band magnitude for S1B between 0.4 and 1.5.
According to the pre-main sequence models of Baraffe et al. ( 2015), a 0.8 M ⊙ star with an age between 0.5 and 1 Myr should have an absolute K band magnitude between 1.5 and 2. We conclude that, within the errors, the observed K band magnitude of S1A is consistent with its dynamical mass.It is important to mention, in addition, that our age and magnitude estimates are quite uncertain.The measurements of Richichi et al. (1994) and Cheetham et al. (2015) suggest significant variability (either intrinsic to S1B or caused by changing foreground absorption).The age of S1B is based on the assumption of co-evality of S1A and S1B and models for S1A that we have shown to be unable to reproduce the actual location of S1A in the HR diagram (Section 5.2).The extinction is assumed to be the same for S1A and S1B, but the presence of circumstellar material probably results in different absorptions for the two stars.Future dedicated infrared resolved observations ought to enable a more detailed comparison between the models and the properties of S1B.

Comparison with Gaia Astrometry
The Gaia satellite is the foremost astrometric mission at optical wavelengths (Gaia Collaboration et al. 2016a).Its most recent data releases (Gaia DR3; Babusiaux et al. 2022 and Gaia EDR3;Gaia Collaboration et al. 2021;Lindegren et al. 2021) were published in 2022 and 2021.The angular resolution of Gaia, 0. ′′ 18 (Lindegren et al. 2021), is insufficient to resolve the S1 system3 since the two components have a maximum separation of about 25 mas (see Figure 3).The Gaia astrometric solution for S1 is listed in Table 3 (Gaia source ID: Gaia DR3 6049161961233192576) where, for easy comparison, we also show the astrometric results derived from our VLBA observations.The Gaia position does not coincide with the VLBA positions for S1A and S1B, as shown graphically in Figure 9.It falls near the VLBA position for S1A, but is slightly displaced from it in the direction of S1B (or, equivalently, in the direction of the center of mass of the system).This is not unexpected: S1A is about 5 times more massive than S1B, and therefore presumably about at least 500 times more luminous.Thus, the photocenter of the S1 system traced by Gaia should indeed be located near S1A.
It is interesting to consider the consequences of the binarity of S1 on the Gaia results.The Gaia DR3 astrometric fitting assumes that S1 is a single star that moves in a linear and uniform fashion.Such a fit does not adequately describe the motion of S1 given the binary nature of the system.The Renormalised Unit Weight Error (RUWE) is a parameter provided by Gaia DR3 and used to assess if the single-star model provides a good fit to the astrometric observations (Lindegren et al. 2018).A RUWE value above 1.4 is considered an indication that a single-star model does not properly fit the astrometric data.In the case of S1, the RUWE = 2.75, confirming that a single-star model does not provide a satisfactory description of the data.The VLBA, on the other hand, resolves the components of the S1 system, and the astrometric fitting takes into account the orbital motion of the system.Thus, for S1, the VLBA astrometric results are certainly superior to those of Gaia.In addition, the VLBA covers a significantly longer time span (14.39 years against 2.76 for Gaia), further improving the quality of the VLBA results (e.g., Dzib et al. 2017).The superiority of the VLBA astrometry is clearly seen in the accuracy of the parallax measurements (0.02 mas for the VLBA vs. 0.14 mas for Gaia).We also note that the two parallax values are barely consistent at the 3σ level.

Comparison with infrared astrometry
The companion of S1A has also been inferred from infrared observations.The S1B component was first identified using the Lunar occultation (LO) method by Richichi et al. (1994).Later, Cheetham et al. (2015) also detected the companion using Sparse Aperture Masking (SAM).Both IR methods provide measurements of the total angular separation and the position angle (P.A.) between the components of the binary.The separations found in the IR observations are between 20 and 30 mas, consistent with the values we have measured at the VLBA (Figure 3).From our best fit, we determined the orbital phase at the epochs of the IR observations, and the corresponding expected separations and P.A.The results are listed in Table 4 and shown graphically in Figure 3.We find that the IR observations were made close to apastron; within the errors, the expected separations and P.A. are in reasonable agreement with the IR results.

Flux variations
In at least one young stellar system, V 773 Tau, variations of the non-thermal radio flux synchronized with the orbital phase of the system have been well documented (Massi et al. 2002;Torres et al. 2012).In that specific case, the radio flux of both stars appears to increase substantially near the periastron of the system.In contrast, the radio emission of S1A shows no systematic variation with the orbital phase, while the flux of the low-mass companion S1B increases near the apastron (Section 4.2; Figure 5).
The origin of this behavior of S1B is unclear.In the case of V 773 Tau, it is plausible that the flux increase near periastron might be caused by increased flaring due to interactions between the individual magnetospheres when the stars approach each other (Massi et al. 2002).Since flaring is thought to be initiated by magnetic reconnection events, it is indeed reasonable to expect an increase in the number of events when the two magnetospheres interact.To interpret the reverse effect observed here with the same class of model would require an unlikely scenario where reconnection events are more frequent when the magnetospheres are further apart.An alternative would be to invoke the presence of an optically thick region surrounding S1A with a radius of order 20 mas.Each time S1B plunged into this hypothetical region, its emission would be absorbed by the optically thick material explaining the lack of detection when the sources are closer than about 20 mas.The main difficulty with this scenario is to explain why this region does not absorb the emission from S1A itself.This requires a somewhat ad hoc morphology, with an inner hole, for the hypothetical absorbing structure -for instance a toroid seen close to face-on.

CONCLUSIONS
We have used seven new and 28 archival VLBA observations to study the S1 system in Ophiuchus, and constrain the dynamical masses of its components.We detected the primary, S1A, in all 35 epochs, and the secondary, S1B, in 14 epochs.The best astrometric fit gives dynamical masses of M S1A = 4.11 ± 0.10 M ⊙ and M S1B = 0.83 ± 0.01 M ⊙ .This is the most accurate mass estimation of S1 to date.The mass of S1A appears to be significantly smaller than previously suggested (of order 6 M ⊙ ).Reanalyzing the spectral energy distribution of S1, we conclude that it is compatible with effective temperatures in the range of 14,000 to 17,000 K.The location of S1 in the HR diagram suggests that it has not yet reached its zero-age main sequence, but we find that pre-main sequence evolutionary models are unable to explain its locus in the HR diagram for a mass of 4.1 M ⊙ .For the models to pass through the allowed range of luminosity and effective temperature of S1A, a mass at least 20 to 50% higher would be needed.Thus, evolutionary models appear to over-predict the mass of S1A by that amount.
A comparison between the VLBA radio astrometry and Gaia DR3 shows that the Gaia position falls close to, but is not coincident with the position of S1A measured by the VLBA.This is not unexpected since Gaia surely traces the motion of the photocenter of the system which should be located close to the primary.However, the fact that Gaia cannot resolve the system nor account for its orbital motion results in a parallax measurement that is significantly less accurate than, and barely compatible with, the VLBA result.Finally, we analyzed the variations of the flux density of S1A and S1B as a function of the orbital phase.We find that the flux of S1A remains stable over the entire orbit while the flux of S1B appears to be higher when the system is close to apastron.The origin of this behavior is unclear.Richichi et al. (1994), we have used the correction of 180 • to the scan angle reported by them.
cessed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.This document was prepared using the collaborative tool Overleaf available at https://www.overleaf.com/.
Facilities: VLBA (Napier et al. 1994) As an additional step to rigorously quantify the key astrometric and orbital parameters, we implemented an additional approach that combines the power of two techniques: Orbits For The Impatient (OFTI) and Markov Chain Monte Carlo (MCMC).Both algorithms are integrated into the Orbitize!library8 , which is a specialized, open-source Python-based tool designed for characterizing orbits of binary astronomical systems and operates within a Bayesian analysis framework (Blunt et al. 2020).For our analysis, we exclusively utilized 14 observations in which both components of the binary system were simultaneously detected, adhering to the necessary criteria for orbit analysis using Orbitize!OFTI works by generating a multitude of possible orbits based on the observed astrometric data.It effectively explores the parameter space, laying the groundwork for subsequent analysis with MCMC (Blunt et al. 2017).Based on the results of the OFTI fit, we defined priors for the orbital elements, including the semimajor axis (a), eccentricity (e), inclination (i), argument of periastron (ω), position angle of nodes (Ω), and epoch of periastron passage (τ ), which represents the fraction of the orbital period that has elapsed from a given reference epoch to periastron, within a range of [0, 1], and then proceeded to the MCMC phase.To enhance the MCMC exploration, we employed a substantial ensemble of 10,000 walkers.These walkers independently traversed the parameter space, ensuring comprehensive sampling.To ensure the convergence of the MCMC chains, we executed a burn-in phase with 500 steps before extracting meaningful results.The results are presented in Table 5, and the corner plot in Figure 10.This corner plot visually demonstrates the well-behaved nature of parameter sampling, further reinforcing the reliability of our results compared to MPFIT.We note, in particular that the total mass of the system derived from Orbitize!(5.14±0.13M ⊙ ) is fully compatible with our MPIFIT result (4.94±0.11M ⊙ ).

Figure 1 .
Figure 1.Final radio images of S1 at 4.5 GHz corresponding to each epoch observed during the DYNAMO-VLBA project.Intensity background images are clipped to intensities between -0.1 to 0.6 mJy beam −1 .Contour levels are at -4, 4, 6, 10, 15, 30, and 60 times the noise levels of the images, listed in Table1.Images are centered on the measured position of S1A, the primary component of the S1 binary system as listed in Table1and shown in the top-left corner of each plot.The white ellipse at the bottom-left of each plot represents the synthesized beam of the corresponding observed epoch.

Figure 2 .
Figure 2. Measured positions of S1A (red dots) and S1B (blue dots) shown as offsets from the position of S1A in the first detected epoch used in the present paper (2005 June 24; see Table1).The red and blue curves show the best fits, described in the text, to the positions of S1A and S1B, respectively.

Figure 3 .
Figure3.Orbital solution of S1.Top-panel: Stellar orbits relative to the central of mass (black cross).The orbits from best fit for S1A (red) and S1B (blue), are accompanied by a subset of 100 orbit pairs (light red for S1A, and grey for S1B) to show the range of solutions allowed by the errors in our fit.The subset of orbits shown corresponds to the best 100 of 5,000 orbits generated by randomly obtaining orbital parameters from normal distributions with mean values equal to the values obtained from our best fit and with standard deviations equal to their errors.Measured positions of S1A (red squares) and S1B (blue circles) after correction of parallax and proper motions are also shown.Position errors are on the order of 0.1 mas, barely seen at the plot scale.Bottom-panel: Stellar relative positions and orbital fit model of S1.The blue dots indicate the relative positions of S1B with respect to S1A, and the errorbars consider the position errors of both components which are added in quadrature.The dashed black line traces the line of nodes from the model, and the black cross indicates the position of the primary.The measured positions of S1B obtained from infrared observations reported inRichichi et al. (1994) andCheetham et al. (2015) are also plotted, as red triangles.These infrared positions are listed in Table4.

Figure 4 .
Figure 4. Flux as a function of orbital phase for S1A (top) and S1B (bottom) at 5 GHz (blue dots) and 8 GHz (red squares).For S1A the bins (gray boxes) show the mean and standard deviation values of fluxes in widths of 0.2 in the orbital phase.For S1B the green arrows indicate the flux 3 σrms upper limits when the source is not detected.The error bars are shown and sometimes are smaller than the marker size.

Figure 5 .
Figure 5. Relative positions and orbital fit model.The green stars indicate the expected positions of S1B in undetected epochs.Orbital phases (ϕ) are indicated with vertical black lines and labeled accordingly.

TFigure 6 .
Figure6.Spectral energy distribution of S1A.The data points (black squares) are taken from VizieR (see Table6) while the fits (color lines) are based on the extinction curves ofRieke et al. (1989) andFitzpatrick (1999) [See text].The figure also includes residual plots for each model (bottom panels).

Figure 9 .
Figure9.VLBA positions at epoch 2016.0 (the reference epoch of the Gaia DR3 release) of S1A (yellow triangle), S1B (purple diamond), and the center of mass of the S1 system (cyan circle).The pink square indicates the Gaia DR3 position and is taken as the origin of the coordinate system.A black cross is shown associated with each position; the sizes of these black crosses are ten times the formal position errors in each case.
J.O., G.O.L., J.M.M. and L.F.R. acknowledge the financial support of CONACyT, México.L.L. acknowledges the support of DGAPA PAPIIT grants IN108324 and IN112820 and CONACyT-CF grant 263356.S.A.D. acknowledge the M2FINDERS project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant No 101018682).P.A.B.G. acknowledges financial support from the São Paulo Research Foundation (FAPESP) under grants 2020/12518-8 and 2021/11778-9.The Long Baseline Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),pro- , Gaia (Gaia Collaboration et al. 2016a) Software: AIPS (Greisen 2003), Astropy 4 (Astropy Collaboration et al. 2013), NumPy 5 (van der Walt et al. 2011), SciPy 6 (Jones et al. 2001), Matplotlib 7 (Hunter 2007), and APLpy (Robitaille & Bressert 2012) APPENDIX A. ORBITAL PARAMETERS ANALYSIS OF S1 USING MCMC

Figure 10 .
Figure10.Corner plot generated from the MCMC analysis using the Orbitize!library.The diagonal panels display 1D marginalized posterior distributions for each of the orbital parameters, including a, e, i, ω, Ω, τ and total mass.The off-diagonal panels illustrate 2D covariances between these parameters.

Figure 11 .
Figure 11.Allowed orbital configurations for the binary system S1, derived from the MCMC analysis conducted using the Orbitize!library.The black marker (star) designates the position of the primary component, and the red points represent radio observations of the secondary component.The upper right panel depicts the evolution of the angular separation over time, while the lower right panel illustrates the variation in the position angle throughout the orbits.

Table 1 .
S1A and S1B measured positions and flux densities from the 35 VLBA observations.
a Upper-limits are three times the noise level in the image.

Table 2 .
Parameters obtained for the best-fit astrometric model.

Table 3 .
Comparison of DYNAMO-VLBA project astrometry with Gaia Astrometry

Table 4 .
Comparison of Infrared and VLBA Astrometry a As pointed by

Table 5 .
Results obtained through the Orbitize!method.
TABLE FOR SED ANALYSIS