© The Authors 2026

1. Introduction

Binary and multiple stellar systems play a fundamental role in the star formation process and in the early dynamical evolution of stellar clusters. Their frequency, separation distribution, and mass ratios provide essential tests for theories of core fragmentation and disk instability (e.g., Duchêne & Kraus 2013). In young star-forming regions, direct measurements of orbital motions are particularly valuable because they yield dynamical masses that can be compared with predictions from pre-main-sequence (PMS) evolutionary models. Such measurements, however, are often challenging: optical spectroscopy provides only projected velocity amplitudes, adaptive-optics imaging offers limited resolution for close systems, and extinction frequently prevents optical access to the youngest members.

Radio interferometry provides a complementary avenue for studying young binaries. Compact, nonthermal radio emission from magnetically active young stellar objects (YSOs) can be monitored at milliarcsecond (mas) scales with very long baseline interferometry (VLBI). Facilities such as the Very Long Baseline Array (VLBA) allow orbital motion to be traced directly over timescales of a few years. VLBI astrometry has already yielded accurate dynamical masses for nearby young multiple systems by directly tracing their stellar motions (e.g., Loinard et al. 2007; Dzib et al. 2010; Azulay et al. 2017; Ortiz-León et al. 2017; Zhang et al. 2020), demonstrating its unique capability to resolve orbits with stellar separations on the astronomical unit scale.

The Orion molecular complex, at a mean distance of ∼400 pc (Kounkel et al. 2017; Großschedl et al. 2019), provides a rich environment to extend such studies to a larger and more diverse population. It is composed of several nearby and well-studied regions of ongoing star formation, including the Orion Nebula Cluster (ONC) and NGC 2024 (see Fig. 1). It hosts hundreds of nonthermal radio-emitting YSOs spanning a broad range of masses and evolutionary stages (Kounkel et al. 2017; Forbrich et al. 2021), making it ideal for investigating how multiplicity evolves with stellar mass and age. Within this framework, the dynamical masses of multiple YSOs with the VLBA (DYNAMO-VLBA) project was designed to obtain multi-epoch VLBI observations aimed at resolving orbital motions and determining dynamical masses in nearby star-forming regions (Ordóñez-Toro et al. 2024), including the Orion complex, as well as ρ-Ophiuchus, Taurus, and Serpens. The program also records all compact radio sources within each observed field, allowing the measurement of parallaxes and proper motions for isolated objects. In Dzib et al. (2026, hereafter Paper I), we introduce the DYNAMO-VLBA observations of Orion, and present all detected sources, the astrometric fitting procedures, the derived astrometric parameters (including distances for individual sources and per region), and a detailed comparison of VLBA and Gaia astrometric results.

Fig. 1. Orion binaries studied in DYNAMO-VLBA for which orbital motion or significant astrometric acceleration could be derived. Central panel: Wide-field infrared view of the Orion star-forming region using WISE imaging. RGB image is composed from W3 (12 μm, red), W2 (4.6 μm, green), and W1 (3.4 μm, blue) bands. Cyan regions indicate areas affected by saturation in the WISE images. Blue symbols mark the positions of 13 binaries observed in our sample, where orbital or acceleration parameters could be derived. Marker shapes are unique to each system for clarity. Zoomed views of the ONC and NGC 2024 regions are also shown. The yellow star marks the position of θ1 Ori C for reference. The surrounding subplots display the projected stellar orbital motions, plotted as apparent motion on the sky. In panels showing two ellipses, the black and red curves correspond to the best-fit astrometric models for the primary and secondary components, respectively. In panels with a single ellipse, the curve represents the orbital motion of the detected radio source. In the cases of V* V1230 Ori B and COUP 450, the arrows represent the measured acceleration vectors and errors. For V* V1230 Ori B, the black cross indicates the relative position of the associated Gaia source, likely corresponding to the primary component of the system.

In this second paper of the DYNAMO-VLBA Orion series, we focus on the subsample of 15 binaries and multiple systems (see Table 1 and Sect. 2). The observations, as well as astrometric and orbital fitting methodology adopted in this work are described in Sect. 3, following the same framework introduced in Paper I. By combining new VLBA observations with earlier VLBA epochs from the Gould’s Belt Distances Survey (GOBELINS; Ortiz-León et al. 2017 and Kounkel et al. 2017, the latter hereafter KHL17), we determined full or partial orbits for these systems and derived total and, when possible, individual stellar masses. Our results (Sect. 4) provide a homogeneous set of VLBI-based dynamical masses for young binary and multiple systems in the Orion molecular cloud complex, derived using a uniform astrometric methodology. In Sect. 5, we discuss the implications of these results, including the occurrence of compact, nonthermal radio emission in intermediate-mass systems, and in Sect. 6 we summarize our main findings and conclusions.

2. Binaries in Orion studied within DYNAMO-VLBA

The DYNAMO-VLBA survey has targeted a sample of young stars in the Orion complex to investigate multiplicity at milliarcsecond scales. This section summarizes the main binary and multiple systems identified or characterized in the course of the program. The systems span a wide range of stellar masses and evolutionary stages, from low-mass embedded YSOs to massive OB-type stars, and include both previously known binaries and new candidate multiples. In Table 1 we compile previous identifications, spectral classifications, and mass estimates. We briefly discuss evidence for multiplicity based on earlier observations.

2.1. VLBI visual binaries.

These are systems in which both stellar components are directly detected in VLBA images. They are of particular importance within the DYNAMO–VLBA survey, since for those with sufficient orbital coverage the VLBA data allow a direct determination of the individual stellar masses. Within this category we identify four systems. CXOU J054146.1–015622 is a deeply embedded young star that has only been detected in X-rays and radio wavelengths (Skinner et al. 2003; Rodríguez et al. 2003; Broos et al. 2013). KHL17 reported two compact radio components (VLBA 61 and 62) separated by ∼15 mas, and provided preliminary mass estimates for both components. The young star Brun 656 (Walker 1969) is located in the Orion Nebula Cluster (ONC) and in which KHL17 detected two radio sources (VLBA 4 and 107) and provided lower-limit mass estimates of ∼1.7 M_⊙ for both components, suggesting an intermediate-mass system. Later, Valegård et al. (2021) performed a detailed spectral energy distribution (SED) analysis and estimated a stellar mass of $2.{72}_{-0.15}^{+0.21}$ M_⊙, proposing that this system may be a precursor to a Herbig Ae/Be star. Finally, while KHL17 reported single radio sources associated with HD 294300 and 2MASS J05414134−0153326, our new VLBA observations reveal the companions to both, confirming their binary nature.

2.2. Spectroscopic binaries.

A second group within our sample consists of spectroscopic binaries previously identified at optical wavelengths and now detected as compact radio sources with the VLBA. Because the spectroscopic mass ratios and orbital parameters are known for these systems, individual stellar masses can, in principle, be derived by combining the spectroscopic solutions with the VLBA astrometry, a possibility already explored by KHL17. HD 37017 (V* V1046 Ori, VLBA 34) and Parenago 1540 (VLBA 19) are well-established spectroscopic binaries with orbital periods of 18.6561 ± 0.0002 and 33.73 ± 0.03 days, respectively (Bolton et al. 1998; Marschall & Mathieu 1988). Both were detected as radio sources by KHL17, who provided preliminary component mass estimates. θ¹ Ori E (VLBA 9) is another short-period spectroscopic binary (P ≈ 9.9 days), also eclipsing, with accurately measured masses of 2.807 and 2.797 M_⊙ (Costero et al. 2006; Morales-Calderón et al. 2012).

2.3. Hierarchical triple and higher-order systems.

These are systems where more than two stars are gravitationally bound. V* NU Ori (HD 37061) is a B0.5V star and the primary ionizing source of the M43 H II region (Simón-Díaz et al. 2011). It is a known hierarchical quadruple system. At optical wavelengths, it appears as a visual binary (components A and B) with a separation of $0\stackrel{″}{.}47$ (Preibisch et al. 1999). Component A is itself a spectroscopic triple system: the primary Aa and its close companion Ab orbit each other every 14 days, while a third component (C) orbits the Aa+Ab binary with a period of 467 days. From spectroscopic and evolutionary model analysis, Shultz et al. (2019) derived component masses of M_Aa = 14.9 ± 0.5 M_⊙, M_Ab = 3.9 ± 0.7 M_⊙, and M_C = 7.8 ± 0.7 M_⊙. No mass has yet been derived for V* NU Ori B. All known components are of intermediate or high mass. V* NU Ori was also suggested to be a magnetic active star; the magnetic fields were initially attributed to the primary (Petit et al. 2008), but Shultz et al. (2019) showed it is actually associated with component C. KHL17 detected two radio sources toward this system, VLBA 27 and 28, but with only 3 and 2 detections respectively, no astrometric solution was derived. Our reanalysis indicates that one detection previously labeled as VLBA 28 likely corresponds to VLBA 27, and the other is likely a spurious detection, leaving a single confirmed radio component associated with V* NU Ori. Finally, θ¹ Ori A is a hierarchical triple system, composed of a spectroscopic binary (θ¹ Ori A₁ and θ¹ Ori A₃) and a wide companion (θ¹ Ori A₂). While KHL17 associated their radio source VLBA 11 with the spectroscopic pair A₁+A₃, Dzib et al. (2021) showed it instead corresponds to the wider companion A₂, a ∼4 M_⊙ star at $\sim 0\stackrel{″}{.}18$ separation (Petr et al. 1998; Schertl et al. 2003; GRAVITY Collaboration 2018).

2.4. Suspected astrometric binaries.

Several embedded Orion stars show nonlinear or anomalous motions in VLBA data, indicating the presence of unresolved companions (KHL17; Dzib et al. 2021). These include COUP 450 (GMR A, VLBA 5), ALMA J053514.5010–052238.674 (COUP 639, GMR H, VLBA 7), V* MT Ori (GMR G, VLBA 8), and V* V1399 Ori (GMR V). In addition, 2MASS J05414134−0153326 (VLBA 153), V* V1230 Ori (VLBA 149), and CXOU J054145.8−015411 (VLBA 124), which had only two detections in KHL17 and were not previously analyzed astrometrically, show positional deviations in our new data consistent with binarity. Although only one component is detected at radio wavelengths, and therefore individual stellar masses cannot be determined, a description of its motion can put constraints on the system properties, as discussed in the next section.

3. Data and analysis methods

3.1. VLBA observations

The observations analyzed in this paper are part of the DYNAMO-VLBA project and were obtained with the VLBA under project code BD215. The full observational setup, calibration procedures, and astrometric analysis have been described in detail in Paper I. Here we summarize only the main aspects relevant to the binary and multiple systems.

Each of the seven DYNAMO-VLBA Orion pointings (blocks F-L, see Paper I) was observed at a frequency of 5.0 GHz with a total bandwidth of 256 MHz. Between February 2018 and December 2019, six epochs were obtained per block, separated by 2–6 months depending on the expected orbital periods of the known binaries. Additional earlier epochs from the GOBELINS project (Kounkel et al. 2017) were included when available, extending the temporal baseline back to 2016.

Data calibration and imaging were carried out in AIPS following the procedures described in Paper I. These include the standard calibration scheme, multiphase calibrator referencing, and geodetic block corrections. The mean synthesized beam was 2.5 × 1.0 mas, and the typical image rms noise was about 30 μJy beam⁻¹. Astrometric positions for all detected components were obtained through two-dimensional Gaussian fitting with JMFIT.

From the 216 compact radio sources detected in the full DYNAMO-VLBA+GOBELINS Orion dataset, 19 were identified to be related to 15 binaries or higher–order multiples showing measurable positional changes inconsistent with linear motion. These systems constitute the sample analyzed in the present paper. Their selection was based on (i) multiple detections of at least two compact components, or (ii) significant curvature or acceleration in the single–source astrometric fits from Paper I. For each system, we combined the positions of the individual components across all available epochs to model their relative and absolute motions.

3.2. Motion fitting and mass estimation

The astrometric positions (α,  δ) of all detected components were modeled following the same formalism described in Paper I. In this approach, each trajectory is represented as the combination of trigonometric parallax (ϖ) and linear proper motions (μ_α ⋅ cos(δ),  μ_δ). For binary and multiple systems, the model was extended as required by the data to include either additional acceleration terms (for long-period orbits, only partially covered by the VLBA time baseline) or full Keplerian orbital terms describing the absolute motion of the detected component(s). The equations used to describe these motions are

$$\begin{array}{cc}\hfill \alpha (t)& ={\alpha }_0+{\mu }_{\alpha }·cos(\delta )·t+\varpi ·f_{\alpha }(t)+\mathrm{\Delta }{\alpha }_{\mathit{orb}}(t),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \delta (t)& ={\delta }_0+{\mu }_{\delta }·t+\varpi ·f_{\delta }(t)+\mathrm{\Delta }{\delta }_{\mathit{orb}}(t),\hfill \end{array}$$(2)

where (f_α, f_δ) represent the parallactic factors projected along right ascension and declination. Δα_orb(t) and Δδ_orb(t) represent the orbital contribution induced by the presence of a companion. These terms may correspond to the contribution of a constant-acceleration approximation (${\dot{\mu }}_{\alpha }·cos(\delta )$, ${\dot{\mu }}_{\delta }$) for systems with long orbital periods, or to the full Keplerian description. The latter case is defined by the orbital elements: Period (P), semi-major axis (a), eccentricity (e), orbit inclination (i), the argument of ascending node (Ω), the argument of pericenter (ω), and the time of passage of pericenter (T₀).

For systems in which both components are detected at radio wavelengths, the relative orbit and the absolute motions of each star can be jointly fitted. The combination of the orbital semi-major axes (a₁, a₂) and the period (P) then yields the total and individual stellar masses through Kepler’s third law, M₁ + M₂ = (a₁ + a₂)³/P², and M₁/M₂ = a₂/a₁. This approach provides model-independent dynamical masses when the orbital coverage is sufficient.

For systems with only one radio-detected component, useful dynamical information can still be derived depending on the type of binary. In spectroscopic binaries, where the mass ratio (q = M₁/M₂) and some orbital parameters are known from optical spectroscopy, these quantities were used as priors in the VLBA modeling, allowing individual masses to be inferred. In sources showing periodic astrometric motion but lacking spectroscopic information, the fitted orbital curvature yields an astrometric mass function, defined as

$$\begin{array}{c}\hfill f(M)=\frac{M_2^3{sin}^{3}i}{{(M_1+M_2)}^2}=\frac{a_1^3}{P^2},\end{array}$$(3)

where a₁ is the semi-major axis of the detected component (in AU) and P is the orbital period (in yr). The mass function provides a lower limit on the companion mass and, when combined with an estimate of the primary mass or inclination, allows further constraints on the total system mass.

Finally, in long-period binaries characterized only by measurable accelerations, the observed curvature in the trajectory can be used to constrain the total system mass. All results from the fits are presented in Sect. 4 and discussed in Sect. 5.

4. Results

4.1. Overview of the astrometric solutions

The parallaxes and proper motions of all compact radio sources in the DYNAMO–VLBA survey, including the targets analyzed here, are presented in Paper I. In this paper, we focus on the subset of sources that are known members of multiple systems or that exhibit clear signatures of binarity, such as orbital motion, measurable accelerations, or periodic residuals in position. For these objects, we derive and report the corresponding orbital (Table 2) or acceleration parameters. From the orbital parameters we estimate either the individual stellar masses or the astrometric mass function (Table 3).

From the 15 analyzed systems, we derived orbital parameters for 11 and acceleration parameters for two. Figure 1 shows the orbital or acceleration solutions obtained for these 13 systems. The motions of the two remaining sources are consistent with linear trajectories affected only by the trigonometric parallax. The best-fit models to the positions of all 15 systems are shown in Fig. 2, and a detailed description of each case is provided below.

Fig. 2. Best-fit astrometric motions of the radio sources associated with stellar multiple systems. Filled circles mark the measured positions relative to the first detection, solid lines show the best-fit trajectories, and open diamonds indicate the expected positions from the fitted model at the corresponding epochs. Different colors represent the individual stellar components detected within each system.

4.2. Radio sources related to multiple YSOs

Binary and multiple systems constitute a significant fraction of the radio-detected YSOs in Orion. In total, we identify 19 sources associated with multiple systems, which we classify into four categories based on the nature of the detections and orbital constraints: (i) radio visual binaries, where two radio components are resolved and full orbital solutions can be derived; (ii) spectroscopic binaries, where external spectroscopic constraints on mass ratios and/or orbital elements are combined with our astrometry; (iii) astrometric binaries, identified from sinusoidal residuals in the positional fits of single detected components; and (iv) wide binaries, for which orbital periods are much longer than the VLBA time baseline and binarity manifests primarily as constant accelerations rather than closed orbits.

For categories (i)–(iii), we apply the orbital fitting procedures described in Paper I, and the reported reference positions in that paper correspond to the system barycenter rather than the locations of individual radio components at the reference epoch. The resulting orbital parameters and dynamical masses are listed in Tables 2 and 3. For category (iv), no orbital solution is attempted; instead, as outlined in Paper I, we model the motion as linear plus constant acceleration to capture the observed curvature over the available time span. In the following subsections, we present the results for each category and discuss their astrometric signatures and implications for dynamical mass estimates.

4.2.1. Radio visual binary stars

Visual binaries resolved in our VLBA observations provide the most direct means of determining dynamical masses, as the orbital motions of both components can be traced without additional assumptions. When two radio sources are detected in the same system with separations of a few tens of mas, their relative motion around the common center of mass yields strong constraints on orbital parameters and stellar masses (e.g., Ordóñez-Toro et al. 2024).

We identified four such systems: CXOU J054146.1–015622, Brun 656, HD 294300, and 2MASS J05414134–0153326 (see also their corresponding panels in Fig. 2). In the latter two cases, we report for the first time the detection of the companion radio sources to those previously identified by KHL17, which we designate VLBA 145B and VLBA 153B, respectively. For these visual binaries, both stellar orbits are well sampled, and the resulting dynamical masses are highly robust, as they follow directly from the measured motions without requiring external constraints.

4.2.2. Spectroscopic binaries

Spectroscopic binaries provide a unique advantage for astrometric studies, as independent spectroscopic measurements constrain key orbital parameters such as mass ratios and, in some cases, orbital elements. When combined with VLBA astrometry, these constraints enable more robust orbital fits and improved stellar mass determinations, even if only one radio component is detected.

Four of our targets (V* NU Ori, θ₁ Ori E, Parenago 1540, and HD 37017) are known spectroscopic binaries. In all cases, only one radio source was detected in our images. We therefore adopted the mass ratios derived from spectroscopy and, where available, fixed additional orbital parameters to their known values, following the procedure outlined by KHL17. This approach allows us to estimate the masses of the detected radio star and its unseen companion.

In the case of V* NU Ori, the best astrometric fit was obtained by associating the detected radio source with component C, which orbits the close pair Aa and Ab. The identification is also consistent with previous suggestions that component C is the magnetically active star. The mass ratio of V* NU Ori was fixed to M_Aa + Ab/M_C = (14.9 + 3.9)/7.8 = 3.2. For Parenago 1540 or HD 37017, no additional VLBA observations were obtained, but we carried out new fits to determine their barycentric positions at the reference epoch 2016.0, which had not been reported previously. For Parenago 1540 our orbital and mass estimates are in good agreement with those of KHL17, whereas for HD 37017 we obtained lower masses. In both cases, our error bars are also smaller than those of KHL17; however, this likely reflects the fact that we fixed orbital parameters from the literature to their most probable values, without propagating their uncertainties. Finally, the tightness of θ₁ Ori E system and its short orbital period (∼9.9 days), prevents us from deriving reliable orbital and mass parameters, and we therefore treat it as a single star in our astrometric analysis and following analysis.

4.2.3. Astrometric binaries

Astrometric binaries are systems in which only one component is detected, but the orbital motion can be inferred from deviations in the measured position relative to a simple linear motion plus parallax model. In such cases, periodic residuals provide strong evidence for binarity, even when the companion is not directly observed.

In our sample, six radio sources (ALMA J053514.5010–052238.674, V* MT Ori, V* V1399 Ori, and CXOU J054145.8–015411) display significantly larger astrometric residuals than typical single stars, consistent with earlier classifications of these objects as binary candidates. The residuals exhibit sinusoidal patterns, suggesting orbital periods comparable to the observational time baseline. For these systems, we therefore modeled the motion by including full orbital parameters, despite the absence of a directly detected companion. The orbital period itself remains uncertain, but we constrained our fits to periods shorter than 10 years, guided by the variability observed in the astrometric residuals.

It is important to emphasize that the orbital solutions derived for these astrometric binaries should be treated with caution. Since only one component is detected, the fits rely on indirect signatures of binarity, which can lead to underestimated uncertainties in the derived parameters. As a result, while the fitted orbits are consistent with the observed residuals, the true parameter uncertainties are likely larger than those formally reported by the model.

4.2.4. Wide binaries

Wide binary systems, with orbital periods much longer than the VLBA time baseline, require a different treatment than close binaries. In these cases, only partial orbital coverage is available, and binarity signatures appear primarily as constant accelerations in the astrometric fits rather than complete orbital solutions.

In the case of source θ₁ Ori A₂ the companion at $\sim 0\stackrel{″}{.}18$ (∼70 AU) from θ₁ Ori A₁ we proceed differently. As the total mass of the system is ∼20 M_⊙, the period is estimated to be ∼130 years, significantly longer than the ∼5 years covered by the VLBA observations used in this work. For this system, we performed an astrometric fit that includes linear motion with a constant acceleration (see Paper I). The resulting acceleration terms are (${\dot{\mu }}_{\alpha }^{\ast }$, ${\dot{\mu }}_{\delta }$) = (−0.01 ± 0.03, −0.01 ± 0.06) mas yr⁻², both consistent with zero within the uncertainties. Thus, from an astrometric perspective, the system can be treated as a single star with the currently available data, and our final fit to this source does not consider the constant acceleration components.

Two additional sources, V* V1230 Ori and COUP 450, are classified as suspected astrometric binaries, likely with long orbital periods. Both display significantly larger astrometric residuals after fitting linear motion and parallax compared to typical single stars. Including acceleration terms improves the fits, and significant detections are obtained. For V* V1230 Ori, the accelerations (${\dot{\mu }}_{\alpha }^{\ast }$, ${\dot{\mu }}_{\delta }$) = (−0.60 ± 0.03,−1.83 ± 0.25) mas yr⁻² are detected at the 20σ and 7σ levels in RA and Dec, respectively. For COUP 450, we measure (${\dot{\mu }}_{\alpha }^{\ast }$, ${\dot{\mu }}_{\delta }$) = (−0.14 ± 0.04, 0.08 ± 0.05) mas yr⁻², with the RA component detected at a significance level of 3.5σ. The implications of these acceleration detections are discussed further in Sect. 5.

The radio source associated with V* V1230 Ori is likely a companion to the massive B1 star seen at optical wavelengths. The radio-optical position difference of 35.80 mas (see Table 6 in Paper I) supports this view. This will, a posteriori, explain the necessity of adding constant-acceleration terms to the motion fits.

5. Discussion

The VLBA astrometry presented in this work enables the first homogeneous dynamical characterization of multiple young stellar systems across the Orion complex. Beyond providing parallaxes and proper motions, the DYNAMO–VLBA survey allows the direct determination of stellar masses and orbital parameters, thereby offering crucial tests of PMS stellar evolution models and a benchmark for multi-wavelength studies of young binaries. The discussion below places our results in the context of these broader astrophysical questions.

We begin by comparing the derived dynamical masses with previous estimates obtained from spectroscopy, SED fitting, or earlier interferometric campaigns. These comparisons assess the internal consistency of our VLBA modeling and highlight cases in which discrepancies may reflect additional complexity, such as orbital inclination effects, limited temporal sampling, or unresolved multiplicity. For two well-characterized systems, Brun 656 and HD 294300, we further perform SED modeling to test the predictions of PMS evolutionary tracks against directly measured stellar masses.

For systems where only one radio component is detected, we discuss how the astrometric mass function and measured accelerations can still provide meaningful dynamical constraints on unseen companions. Finally, we examine the implications of our findings for the occurrence of magnetic activity in intermediate-mass young stars, particularly those exhibiting nonthermal radio emission, which may trace the onset of magnetism in a regime where standard convective dynamos are not expected to operate.

5.1. Mass comparisons

We have determined dynamical masses for seven systems (Table 3) and compare them here with previously reported values listed in Table 1. This comparison provides a first benchmark for assessing the reliability of our results and identifying systems where discrepancies may indicate additional astrophysical complexity. Masses reported in the literature mainly were derived from spectroscopic orbital fitting or from placement in the Hertzsprung–Russell diagram and comparison with evolutionary models. Only in the case of CXOU J054146.1−015622 were the previously estimated masses also obtained from VLBA astrometry, albeit with fewer detections.

For several systems, such as Brun 656, V* NU Ori, and HD 294300, the VLBA-derived dynamical masses are broadly consistent with values obtained from spectroscopy or SED fitting, lending confidence to the robustness of our orbital modeling. In particular, Brun 656 shows excellent agreement between our measurement (2.64 ± 0.13 M_⊙ for the primary) and the SED-based estimate of $2.{72}_{-0.15}^{+0.21}$ M_⊙ by Valegård et al. (2021), confirming its intermediate-mass nature. Likewise, for V* NU Ori, our dynamical mass estimate for component C (6.96 ± 0.32 M_⊙) is consistent with the spectroscopic mass estimate of 7.8 ± 0.7 M_⊙ from Shultz et al. (2019) for the same component, supporting its association of nonthermal radio emission and magnetic activity.

In contrast, for the close spectroscopic binaries Parenago 1540 and HD 37017, our dynamical masses are significantly lower than those previously reported. In the case of Parenago 1540 our estimated masses are M₁ = 0.34 ± 0.04 M_⊙ and M₂ = 0.44 ± 0.05 M_⊙, whereas earlier estimates were 1.7 M_⊙ and 1.25 M_⊙, respectively. Similarly, for HD 37017 we derive M₁ = 2.28 ± 0.07 M_⊙ and M₂ = 1.20 ± 0.03 M_⊙, while previous estimates range from 3.5 to 8.5 M_⊙ for the main star and from 1.8 to 4.5 M_⊙ for its companion. This discrepancy most likely reflects the sparse temporal sampling of the available VLBA data (with cadences set by the GOBELINS survey) and the short orbital periods of the systems, which together limit the accuracy of orbital fitting. These results suggest that the formal uncertainties in these cases are underestimated, and that additional VLBA monitoring with improved cadence will be essential to obtain more reliable dynamical masses.

Among the four systems in which both stellar components were detected with the VLBA, the derived dynamical masses offer an exceptional opportunity for direct comparison with independent estimates from SED modeling. Such a comparison is particularly valuable because it allows the empirical calibration of PMS evolutionary tracks using dynamically measured stellar masses. However, for two of these systems, CXOU J054146.1−015622 and 2MASS J05414134−0153326, the available photometric data are too limited to perform a reliable SED analysis. Therefore, in the following we focus on Brun 656 and HD 294300, the two systems with sufficient multiwavelength coverage to enable detailed SED modeling and comparison with PMS evolutionary predictions.

5.2. Spectral energy distribution (SED) modeling

The precise determination of dynamical masses in binary systems provides a fundamental reference for calibrating PMS evolutionary models, which still show significant discrepancies, particularly in the intermediate–mass regime. A clear example is the previous study of the Oph S1 system, where PMS models overpredicted the mass by 25–50% relative to the dynamical value (Ordóñez-Toro et al. 2024, 2025). To extend such tests to the Orion sample, we analyze the SEDs of Brun 656 and HD 294300, the two binaries in which both components are resolved with the VLBA and have broad multi-wavelength photometric coverage.

Modeling the SEDs allows us to derive the stellar effective temperatures (T_eff), bolometric luminosities (L_★), and extinction values (A_V), thereby enabling a direct comparison between dynamical and model-predicted stellar masses. The photometric characterization and SED fitting were performed with the ARIADNE code (Vines & Jenkins 2022), which implements a fully Bayesian approach with Bayesian Model Averaging (Fragoso et al. 2018), combining several stellar atmosphere libraries according to their Bayesian evidence. Multiband photometry was automatically compiled through astroquery (Ginsburg et al. 2019) from Gaia DR3 (Gaia Collaboration 2023), 2MASS (Skrutskie et al. 2006), Strömgren-Paunzen (Paunzen 2015), Strömgren–Hauck (Hauck & Mermilliod 1998; Mermilliod 1986), TYCHO-2 (Høg et al. 2000), Johnson UBV (Mermilliod 1986), SDSS DR12 (Alam et al. 2015), and TESS (Stassun et al. 2019); low-quality SkyMapper data were discarded. The SEDs were modeled using the PHOENIX v2 (Husser et al. 2013), BT-Settl, BT-NextGen, and BT-Cond (Allard et al. 2011, 2012; Hauschildt et al. 1999), and Kurucz/CK04 (Castelli & Kurucz 2003) atmosphere grids. Extinction was treated as a free parameter using the law of Fitzpatrick (1999) with R_V = 3.1.

Distances were fixed to 403 ± 9 pc for Brun 656 and 405 ± 10 pc for HD 294300 based on our VLBA parallaxes (Paper I), and T_eff priors were set according to the literature spectral types (Walker 1969; Valegård et al. 2021; Pinzón et al. 2021). Posterior sampling was performed with dynesty (Speagle 2020). The resulting parameters are: for Brun 656, $L_{★}=25.0_{-2.2}^{+2.7}{L}_{\odot }$ and A_V = 1.65 ± 0.10 mag; for HD 294300, $L_{★}=25.7_{-2.0}^{+2.3}{L}_{\odot }$ and A_V = 0.67 ± 0.06 mag. In both systems, the SEDs are well reproduced from the optical to mid-infrared wavelengths without significant excess and residuals < 2.5σ (see Fig. 3), indicating that the emission is dominated by stellar photospheres.

Fig. 3. Spectral energy distributions of Brun 656 and HD 294300 obtained with ARIADNE. The solid blue line represents the best-fit model, and black markers denote the observed photometry used in the fit.

5.3. Comparison with PMS evolutionary models

Because the observed SED represents the combined emission of the system, it was necessary to estimate the individual contribution of each component to the total luminosity before comparing with PMS models. For each system, we used the total luminosity derived from the SED fit obtained with ARIADNE.

To determine the contribution of each component, we employed PMS models from which theoretical luminosities L₁ and L₂ were extracted over an age range consistent with previous estimates: 1–3 Myr (Brun 656) and 1.5–3.5 Myr (HD 294300) (e.g., Valegård et al. 2021; Pinzón et al. 2021). For each model, we computed the luminosity ratio between components and analyzed its evolution with age in order to quantify the relative contribution of each star to the system’s total luminosity.

Individual luminosities were estimated directly from the luminosity ratio predicted by the evolutionary models. For each evolutionary grid and within the adopted age range of each system, we computed the luminosities expected for the primary and secondary components at the dynamically determined masses and derived the corresponding luminosity ratio, R_L = L₁/L₂, on a fine age grid. The full set of R_L values obtained by sampling all models and ages was then combined into a single distribution, which was summarized by its median, with percentile ranges used to quantify the associated dispersion. Assuming that the observed total luminosity satisfies L_tot = L₁ + L₂, the median luminosity ratio was applied to the total luminosity derived from the SED fit using ARIADNE to obtain the individual luminosities of each component. In this way, the luminosity partition reflects the theoretical proportions predicted by PMS models for the adopted dynamical masses and age ranges, while avoiding the introduction of multiple, model-dependent luminosity values.

The PMS models considered throughout this work (shown in Fig. 4) are: Siess (Siess et al. 2000), Pisa (Tognelli et al. 2011), MIST-MESA (Choi et al. 2016; Paxton et al. 2011), YaPSI (Spada et al. 2017), BaSTI (Hidalgo et al. 2018) and PARSEC (Nguyen et al. 2022).

Fig. 4. Location of the primary components of Brun 656 and HD 294300 on the HR diagram and evolutionary tracks based on pre–main-sequence stellar models.

The results of the luminosity partition are: Brun 656 with $L_1=20.{43}_{-3.26}^{+3.31}{L}_{\odot }$ and $L_2=4.{58}_{-2.65}^{+2.63}{L}_{\odot }$; HD 294300 with $L_1=13.{84}_{-1.40}^{+1.44}{L}_{\odot }$ and $L_2=11.{89}_{-1.25}^{+1.21}{L}_{\odot }$. The effective temperatures used correspond to spectroscopic estimates representative of the optically dominant primary component. Figure 4 shows the location of the primary components together with isochrones of 1, 3, and 5 Myr.

Anchoring the analysis to the VLBA dynamical masses (M_Brun 656 = 2.64 ± 0.13 M_⊙; M_HD 294300 = 2.01 ± 0.12 M_⊙), we derived the ages implied by each PMS grid. For Brun 656, all tracks place the star between the 2.5 and 3.0 M_⊙ sequences, yielding ages of 2–3 Myr. For HD 294300, most models reproduce the observed position with a 2.0–2.5 M_⊙ track at ∼3 Myr, except Siess and YaPSI, which require a slightly older age (∼5 Myr). Overall, the scatter among models is modest (≲0.3–0.5 M_⊙ in mass or ≲2 Myr in age), and all are broadly consistent with the independent age estimates from optical studies.

Conversely, by adopting the literature age priors (1–3 Myr for Brun 656; 1.5–3.5 Myr for HD 294300) and reading off the corresponding model masses, a systematic difference emerges between the two systems. For Brun 656, the MIST–MESA, BaSTI, Pisa, and YaPSI grids reproduce the VLBA dynamical mass very well, predicting M ≃ 2.5–2.7 M_⊙, in close agreement with the measured value of 2.64 ± 0.13 M_⊙. The Siess models yield slightly lower masses, while PARSEC predicts somewhat higher values (∼3.0 M_⊙). In contrast, for HD 294300 (M_dyn = 2.01 ± 0.12 M_⊙), nearly all models (MIST–MESA, BaSTI, Pisa, and PARSEC) predict systematically larger masses (2.2–2.6 M_⊙) at the nominal 2–3 Myr age, whereas only the Siess grid yields a value consistent with the dynamical mass, but at an older age of ∼5 Myr. This behavior suggests that while current PMS models reproduce the evolutionary locus of Brun 656 remarkably well, they may overestimate luminosities (or, equivalently, underpredict contraction rates) for slightly less massive stars such as HD 294300.

Among the six grids explored, the MIST–MESA and BaSTI models achieve the best overall consistency between dynamical and theoretical predictions for both systems, followed closely by Pisa. On the contrary, the PARSEC models show the larger disagreement by overpredicting the masses of both stars by ∼0.5 M_⊙. The remaining models display small but systematic offsets of ∼0.2–0.4 M_⊙, likely reflecting differences in input physics (e.g., convection efficiency, boundary conditions, or magnetic inhibition of convection). Overall, the DYNAMO–VLBA results highlight Brun 656 and HD 294300 as powerful empirical anchors for PMS model calibration in the intermediate–mass regime, providing direct, model-independent benchmarks to refine future evolutionary tracks.

5.4. Masses of remaining stars

While the dynamical masses of the remaining systems cannot be directly determined, the derived astrometric mass functions and acceleration vectors provide valuable constraints on the stellar masses and orbital configurations of these sources. Although the inferred parameters are subject to projection effects and limited time coverage, they nonetheless provide useful lower bounds on companion masses and orbital periods for embedded Orion systems where only one component is detected.

5.4.1. The astrometric mass function

The astrometric mass function, through Eq. (3), provides a relation between the companion mass, the primary mass, and the orbital inclination. When the primary mass and inclination are constrained, the only remaining unknown in Eq. (3) is the mass of the companion, M₂. By numerically solving Eq. (3), the mass function can therefore be used to estimate the mass of the companion.

From our orbital fitting, we derived the astrometric mass function f(M) (Table 3) and the orbital inclinations (Table 2) for four systems. Two of these systems, V* MT Ori and V* V1399 Ori, have known spectral types (SpTs) from the literature (see Table 1), which allow us to estimate the masses of their primary components independently. Their SpTs were converted to effective temperatures (T_eff) using the intrinsic scale for young stars from Pecaut & Mamajek (2013) and assuming an age of 3 Myr for both sources, consistent with typical ages of ONC members (e.g., Da Rio et al. 2010; Kounkel et al. 2018)¹.

The variable star MT Ori (SpT ∼ K2) corresponds to a 1.0 ± 0.1 M_⊙ young star, while V* V1399 Ori (SpT ∼ G8) corresponds to a 1.3 ± 0.1 M_⊙ star. Using these primary masses, together with orbital inclinations listed in Table 2, and the corresponding values of f(M), we solve Eq. (3) to infer companion masses. The resulting values are $0.{27}_{-0.14}^{+2.88}$ M_⊙ and $0.{63}_{-0.06}^{+0.07}$ M_⊙ for MT Ori and V1399 Ori, respectively.

5.4.2. Acceleration

Constant angular acceleration parameters were derived for two additional stars, COUP 450 (GMR A) and V* V1230 Ori. The measured angular accelerations can be used to constrain the lower limit of the mass of the radio source companion if the projected separation is known. Moreover, the acceleration vector is expected to point toward the companion.

COUP 450 is a deeply embedded object previously detected only in radio and X-rays. As discussed by Dzib et al. (2021), its measured proper motion varies between VLBA epochs, and our results confirm this behavior. However, because no secondary component has been detected or reported, further characterization of the system is not possible at this stage.

V* V1230 Ori, on the other hand, is associated with a B1-type star. As noted by Paper I, the separation between the VLBA radio source and the position reported in Gaia DR3 is 35.8 mas, indicating that the radio emission likely arises from a lower-mass companion rather than from the massive B1 star itself. As shown in Fig. 1, the measured acceleration vector points toward the optical primary, as expected if it dominates the gravitational potential of the system. In this context, the derived angular acceleration of the radio source can be used to estimate a lower limit for the mass of the massive optical star.

For simplicity, and in the absence of additional orbital information, we assume a circular and face-on orbit. From Newton’s law,

$$\begin{array}{c}\hfill m{\mathtt{a}}_{\mathrm{tot}}=\frac{GmM}{r^2}\Rightarrow M=\frac{{\mathtt{a}}_{\mathrm{tot}}{r}^2}{G},\end{array}$$

where r is the projected separation between the stars, a_tot is the physical acceleration of the radio source, and M and m denote the masses of the primary and secondary components, respectively. We adopt r ≃ 14 AU (2.1 × 10¹⁴ cm) and a total angular acceleration of ${\dot{\mu }}_{\mathrm{tot}}=1.9$ mas yr⁻² which corresponds to physical acceleration of a_tot = 0.011 cm s⁻². These values yield a minimum mass of M ≃ 3.7 M_⊙ for the optical star, consistent with its identification as a B1-type star, whose expected mass is ∼12 M_⊙.

Assuming the radio-emitting component is less massive than the optical star, we can also estimate a rough orbital period. Taking a total system mass between 12 and 20 M_⊙ and the observed projected separation (∼14 AU) as the semi-major axis, we obtain a period of P ≃ 12–15 yr. We note, however, that this represents only a lower limit, as the true period depends on the (unknown) orbital inclination and eccentricity.

5.5. Intermediate mass stars with nonthermal radio emission

In VLBI observations, the detection of radio emission at milliarcsecond scales provides strong evidence for a nonthermal emission mechanism, as it implies sources with brightness temperatures T_b > 10⁶ K, well above those attainable by thermal free–free emission (Thompson et al. 2017). While not a primary objective of the DYNAMO-VLBA project, the sensitivity and angular resolution of our observations allow us to identify nonthermal emission associated with some intermediate-mass systems.

Nonthermal radio emission is commonly interpreted as a signature of magnetic activity. However, magnetic activity is generally not expected in intermediate- and high-mass stars, since their interiors are dominated by radiative energy transport and they lack the convective envelopes required to sustain a solar-type dynamo (Kraft 1967; Mestel 1968). In particular, stars more massive than ∼2 M_⊙ are predicted to evolve with radiative outer layers throughout most of their PMS and main-sequence lifetimes, preventing the magnetic field amplification mechanisms observed in lower-mass T Tauri stars (Gregory et al. 2012). Nevertheless, several observational studies have revealed signatures of strong magnetic fields and activity in some young intermediate-mass stars. An initial explanation was that these signatures arise not from the intermediate-mass star itself, but from an unresolved lower-mass companion. Observationally, in the ONC, hundreds of X-ray–emitting A- and late-B–type stars have been detected whose temporal variability (e.g., flare-like light curves and timescales) is statistically indistinguishable from that of magnetically active T Tauri stars, strongly suggesting the contribution of hidden low-mass companions (Montmerle et al. 2005). Similarly, studies of Herbig Ae/Be systems indicate that X-ray emission is often due to chance or physical companions rather than intrinsic stellar activity (Stelzer et al. 2006). However, recent results have found evidence of magnetic activity directly associated with young intermediate-mass stars. These are often attributed either to fossil magnetic fields retained from the star formation process (Donati & Landstreet 2009; Bagnulo et al. 2020), or to magnetic interactions within close binary systems where tidal forces or accretion may trigger dynamo-like activity (Stelzer et al. 2009; Alecian et al. 2013). The detection of nonthermal radio emission in intermediate-mass YSOs thus provides critical constraints on the origin and evolution of magnetism in this transitional stellar mass regime.

Using the dynamical masses derived in this work together with the compact VLBI detections, we identified the following intermediate-mass stars associated with nonthermal radio emission: Brun 656, HD 294300, and V* NU Ori C. In the case of HD 37017, our determined dynamical mass is lower than spectroscopic estimates, likely due to underestimated errors (see Sect. 4). However, Hohle et al. (2010) pointed out that the luminosity of this system is consistent with a star with a mass of 8.5 M_⊙ and, consequently, a companion of 4.5 M_⊙. As the radio source is associated with both or any of these stars, our result is, in any case, consistent with nonthermal radio emission arising from an intermediate-mass star.

Additionally, other observed radio sources fall into the category of intermediate-mass stars with nonthermal radio emission. Although no dynamical masses were estimated for them in this work, their classification is supported by their spectral type and by the small angular separation between our radio and Gaia results. These stars are HD 290862, [SSC75] M78 11, θ₁ Ori E, HD 37150, TYC 5346-538-1, HD 294264, and V* KO Ori. Finally, θ₁ Ori A₂ is a 4 M_⊙ star, and its radio emission is well known to be nonthermal, and it also falls within the present category.

Altogether, intermediate-mass stars with nonthermal radio emission are 12 of the 54 stellar systems with VLBA fitted astrometry (see also Dzib et al. 2026), corresponding to ∼22% of our sample. This fraction suggests that radio emission from magnetic activity is relatively common in this mass regime.

Of particular interest is the radio source associated with star V* NU Ori C. The derived dynamical mass M_C,  VLBA = 6.96 ± 0.32 M_⊙ is in good agreement with the result of M_C = 7.8 ± 0.7 M_⊙ by Shultz et al. (2019). This directly associates a nonthermal radio source with an intermediate-mass star lying just below the high-mass-star limit (M ≥ 8 M_⊙). Previously, it was suggested that Oph S1 could represent another such case, as it was commonly assumed to have a mass between 6 and 8 M_⊙. However, recent work by Ordóñez-Toro et al. (2024, 2025) has shown that Oph S1 instead has a mass of 4.12 ± 0.04 M_⊙. Since the other known nonthermal sources discussed above are consistent with spectral types corresponding to masses < 5 M_⊙, V* NU Ori C remains a unique example of a nonthermal radio emitter with a mass close to the high-mass-star boundary.

Because the radio sources in our sample are all compact, with sizes of only a few mas, their emission is most plausibly linked to magnetic activity, as suggested by previous X-ray studies. A detailed exploration of the physical origin of this activity is beyond the scope of the present paper, but it clearly requires a coordinated multi-wavelength approach to fully disentangle the mechanisms at work.

6. Conclusions

We have presented the results of the DYNAMO–VLBA survey focused on young binary and multiple systems within the Orion complex. In our multi–epoch VLBA observations, we resolved multiple radio components in the HD 294300 and 2MASS J05414134−0153326 systems, confirming their binary nature. More generally, we obtained precise astrometric measurements that allowed us to resolve or detect orbital motions in 15 systems and to derive dynamical masses for seven of them. These results represent one of the most comprehensive VLBI determinations of stellar masses for young binaries in a single star–forming complex.

For four systems in which both components are detected at radio wavelengths, the derived orbital solutions yield fully model-independent dynamical masses. In particular, Brun 656 and HD 294300 exhibit excellent agreement between their VLBA–based masses and values inferred from spectral energy distribution modeling, supporting the reliability of our orbital fits and providing valuable empirical benchmarks for pre–main–sequence evolutionary models. For the spectroscopic binaries V* NU Ori, Parenago 1540, and HD 37017, our combination of VLBA astrometry and published spectroscopic constraints provides updated dynamical estimates, although the limited temporal sampling of the VLBA data leads to increased uncertainties for the shortest–period systems.

In addition to these binaries with complete orbits, we identified several systems exhibiting clear astrometric signatures of multiplicity, as evidenced by acceleration or periodic residuals. For these, we derived astrometric mass functions or lower limits to companion masses, extending dynamical constraints to systems in which only one component is detected. The detection of significant accelerations in sources such as COUP 450 and V* V1230 Ori further demonstrates the capability of VLBI to uncover companions at separations of a few astronomical units, even in embedded or optically obscured regions.

Our results also reveal that a substantial fraction of the radio–emitting young stars in Orion are intermediate-mass objects exhibiting nonthermal emission. Systems such as V* NU Ori C, Brun 656, and HD 294300 provide direct evidence of magnetic activity in stars near or above the limit where standard convective dynamos are expected to cease, offering key observational constraints on the origin of magnetism in this transitional mass regime.

Altogether, the DYNAMO–VLBA observations establish VLBI astrometry as a uniquely powerful tool for measuring stellar masses and orbital dynamics in young multiple systems. The derived dynamical masses, combined with SED–based luminosities and future spectroscopic monitoring, will provide essential benchmarks for refining PMS evolutionary models and for advancing our understanding of magnetic activity and multiplicity among young intermediate-mass stars.

Acknowledgments

We thank the anonymous referee for a careful and constructive review that helped improve the clarity and presentation of this work. S.A.D. acknowledges the M2FINDERS project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant No 101018682). G.N.O.L. acknowledges the financial support provided by Secretaría de Ciencia, Humanidades, Tecnología e Innovación (Secihti) through grant CBF-2025-I-201. P.A.B.G. acknowledges financial support from the São Paulo Research Foundation (FAPESP, grant: 2020/12518-8) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, grant: 303659/2024-6). J.M.M.S. acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grants PID2022-136828NB-C41 and EX2024-001451-M funded by MICIU/AEI/10.13039/501100011033/ERDF/EU. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

All Tables

All Figures

Fig. 1. Orion binaries studied in DYNAMO-VLBA for which orbital motion or significant astrometric acceleration could be derived. Central panel: Wide-field infrared view of the Orion star-forming region using WISE imaging. RGB image is composed from W3 (12 μm, red), W2 (4.6 μm, green), and W1 (3.4 μm, blue) bands. Cyan regions indicate areas affected by saturation in the WISE images. Blue symbols mark the positions of 13 binaries observed in our sample, where orbital or acceleration parameters could be derived. Marker shapes are unique to each system for clarity. Zoomed views of the ONC and NGC 2024 regions are also shown. The yellow star marks the position of θ1 Ori C for reference. The surrounding subplots display the projected stellar orbital motions, plotted as apparent motion on the sky. In panels showing two ellipses, the black and red curves correspond to the best-fit astrometric models for the primary and secondary components, respectively. In panels with a single ellipse, the curve represents the orbital motion of the detected radio source. In the cases of V* V1230 Ori B and COUP 450, the arrows represent the measured acceleration vectors and errors. For V* V1230 Ori B, the black cross indicates the relative position of the associated Gaia source, likely corresponding to the primary component of the system.

In the text

Fig. 2. Best-fit astrometric motions of the radio sources associated with stellar multiple systems. Filled circles mark the measured positions relative to the first detection, solid lines show the best-fit trajectories, and open diamonds indicate the expected positions from the fitted model at the corresponding epochs. Different colors represent the individual stellar components detected within each system.

In the text

	Fig. 3. Spectral energy distributions of Brun 656 and HD 294300 obtained with ARIADNE. The solid blue line represents the best-fit model, and black markers denote the observed photometry used in the fit.
In the text

	Fig. 4. Location of the primary components of Brun 656 and HD 294300 on the HR diagram and evolutionary tracks based on pre–main-sequence stellar models.
In the text