Substellar Companions of the Young Weak-line TTauri Star DoAr21

The most popular method to search for periodicities in data is the so-called Lomb–Scargle periodogram. However, this method performs optimally only under an important implicit assumption: all the other signals (e.g., linear trend, an average offset, etc.) can be subtracted from the data without affecting the significance of the signal. This assumption does not hold for astrometry because the proper motion and the parallax are also a significant part of the signal and they typically correlate with the periodic motion of a companion (see Black & Scargle 1982).

Following the procedure presented by Anglada-Escudé et al. (2010), we use instead a circular least-squares (CLS) periodogram. In this approach, the weighted least-squares solution is obtained by fitting all the free parameters in the model for a given period. The sum of the weighted residuals divided by N is the so-called χ² statistic, where N is the number of data points. Then, each χ_P² of a given model with k_P parameters can be compared to the ${\chi }_{0}^{2}$ of the null hypothesis with k₀-free parameters by computing the power, z, as

$$\begin{eqnarray}&&z(P)=\displaystyle \frac{({\chi }_{0}^{2}-{\chi }_{P}^{2})/({k}_{P}-{k}_{0})}{{\chi }_{P}^{2}/({N}_{\mathrm{obs}}-{k}_{P})},\end{eqnarray} \tag{ 1 }$$

where a large z is interpreted as a very significant solution. The values of z follow a Fisher F-distribution with k_P − k₀ and N_obs − k_P degrees of freedom (Scargle 1982; Cumming 2004). Even if only noise is present, a periodogram will contain several peaks (Scargle 1982, as an example) whose existence have to be considered in obtaining the probability that a peak in the periodogram has a power higher than z(P) by chance, which is the so-called false alarm probability (FAP):

$$\begin{eqnarray}&&\mathrm{FAP}=1-{\left(1-{\rm{p}}{\rm{rob}}[z\gt z(P)]\right)}^{O},\end{eqnarray} \tag{ 2 }$$

where O is the number of independent frequencies. In the case of uneven sampling, O can be quite large and is roughly the number of periodogram peaks one could expect from a data set with only Gaussian noise and the same cadence as the real observations. We adopt the recipe O ∼ 2ΔT/P_min given in Cumming (2004, Section 2.2), where ΔT is the time span of the observations and P_min is the minimum period searched. For instance, assuming that ΔT = 2300 days and P_min = 20 days, the astrometric data is expected to have O ∼ 230 peaks.

In our case, the null hypothesis is the basic kinematic model with k₀ = 5: two reference coordinates, two proper motions, and the parallax. As a first approach, our simplest non-null hypothesis considers circular orbits. For a given period, the number of free parameters is then k_P = 9: five kinematic parameters plus the four Thiele–Innes elements X1, Y1, X2, and Y2 (Green 1993).

As a second approach, we do a least-squares fit of a full Keplerian orbit to the astrometric data, and if we obtain a reasonable fit for the eccentricity of the orbit, we recalculate the periodogram with the eccentricity fixed to this value. To find the new least-squares periodogram, we find the time passage at the periastron τ using a Monte Carlo method for each given period (and the eccentricity fixed), and we perform a least-squares periodogram sampling of a grid of fixed eccentricity–period (eP; where e is fixed either to 0 or the fitted value) pairs and fitting all other parameters. For each eP pair and estimated τ, k_P is 9: the null-hypothesis parameters plus all the other Keplerian elements (using the Thiele–Innes elements): ω, Ω, a₁, and i, where a₁ is the semimajor axis of the orbit of the star due to the companion. As mentioned in Section 2, in the fitting of the data, we only use the positions of DoAr21 published in Ortiz-León et al. (2017) and those that we report here, which cover a time span of about 2375 days. Since the model is linear in all nine parameters, the power can be efficiently computed for many periods between 20 and 3500 days (i.e., a period larger than the time span of the observations) to obtain the familiar representation of the periodogram (see Figure 2).

Zoom In Zoom Out Reset image size

Here, we present a model of μas astrometric data of the sort that can be provided by VLBI, such as the VLBA, as well as GAIA (and in the future by the next generation Very Large Array (ngVLA)).

The source barycentric two-dimensional position is described as function of time, accounting for the (secular) effects of proper motions (μ_α and μ_δ), the (periodic) effect of the parallax Π, and the (Keplerian) gravitational perturbation induced on the host star by one or more companions (low-mass stars, substellar companions, or planets; mutual interactions between companions are not taken into account). The proper motions and the annual parallax terms can be expressed as follows:

$$\begin{eqnarray}\begin{array}{rcl}\alpha (t) & = & {\alpha }_{0}+({\mu }_{\alpha }\cos (\delta ))(t-{t}_{0})+0.5({a}_{\alpha }\cos (\delta )){\left(t-{t}_{0}\right)}^{2}\\ & & +\ {\rm{\Pi }}{F}_{\alpha }(t)+\displaystyle \sum {G}_{\alpha }(t),\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}\delta (t) & = & {\delta }_{0}+{\mu }_{\delta }(t-{t}_{0})+0.5{a}_{\delta }{\left(t-{t}_{0}\right)}^{2}+{\rm{\Pi }}{F}_{\delta }(t)\\ & & +\ \displaystyle \sum {G}_{\delta }(t),\end{array}\end{eqnarray} \tag{ 4 }$$

where (α₀, δ₀) is a reference position, t₀ is a reference time (usually the mean epoch of the multi-epoch observations), and a_α and a_δ are acceleration terms needed to take into account the dynamical effect induced by long-period stellar companions to the primary (e.g., when the primary is part of an unseen long-period stellar binary; these acceleration terms can also be necessary for close-by high-velocity stars). However, the acceleration terms are not always necessary. F_α and F_δ refer to the astrometric displacement due to the parallax in α and δ directions during observations, which are given by (Seidelman 1992):

$$\begin{eqnarray}&&{F}_{\alpha }(t)=(X\sin (\alpha )-Y\cos (\alpha ))/(15\cos (\delta )),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{F}_{\delta }(t)=X\cos (\alpha )\sin (\delta )+Y\sin (\alpha )\sin (\delta )-Z\cos (\delta ),\end{eqnarray} \tag{ 6 }$$

here (X, Y, Z) represent the Cartesian components in equatorial coordinates of the position of the observatory at the time of the observations, t, with respect to the solar system barycenter (in units of au when Π is in arcsec) and α and δ are the coordinates of the barycentric place of the source at each epoch. These values for the Earth are available from the NASA Jet Propulsion Laboratory Solar System ephimerides.⁵

The term (G_α(t),G_δ(t)) is the induced Keplerian orbit due to an unseen companion, and the parallax factors are defined using the classic formulation by Green (1993). The Keplerian orbit of each companion is scaled and projected onto the plane of the sky through (Green 1993)

$$\begin{eqnarray}&&{G}_{\alpha }(t)=r[\cos (\nu +\omega )\sin ({\rm{\Omega }})+\sin (\nu +\omega )\cos ({\rm{\Omega }})\cos (i)],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{G}_{\delta }(t)=r[\cos (\nu +\omega )\cos ({\rm{\Omega }})-\sin (\nu +\omega )\sin ({\rm{\Omega }})\cos (i)],\end{eqnarray} \tag{ 8 }$$

where i is the inclination of the orbital plane (such that i = 0 corresponds to a face-on, anticlockwise orbit), ω is the longitude of the periastron, Ω is the position angle of the line of nodes, ν is the true anomaly, and r is the radius vector, which can be expressed in terms of the true anomaly ν, or the eccentric anomaly E, using the dynamical equations

$$\begin{eqnarray}&&r=\displaystyle \frac{a(1-{e}^{2})}{(1+e\cos (\nu ))}=a(1-e\cos (E)),\end{eqnarray} \tag{ 9 }$$

where e is the eccentricity, and a is the apparent semimajor axis of the star's orbit around the systemic barycenter, i.e., the astrometric signature. The eccentric anomaly is the solution of Kepler's equation:

$$\begin{eqnarray}&&E-e\sin (E)=M,\end{eqnarray} \tag{ 10 }$$

with the mean anomaly M, expressed in terms of the orbital period P and the epoch of the periastron passage τ:

$$\begin{eqnarray}&&M=\displaystyle \frac{2\pi (t-\tau )}{P}.\end{eqnarray} \tag{ 11 }$$

Finally, the true anomaly ν is function of the eccentricity and the eccentric anomaly:

$$\begin{eqnarray}&&\tan \left(\displaystyle \frac{\nu }{2}\right)={\left(\displaystyle \frac{1+e}{1-e}\right)}^{1/2}\tan \left(\displaystyle \frac{E}{2}\right).\end{eqnarray} \tag{ 12 }$$

Note that since a is the apparent semimajor axis of the star's orbit in units of arcsec, its relationship to the mass of the secondary depends on the astrometric method used. For astrometric perturbations due to an unseen companion, we have from Kepler's third law:

$$\begin{eqnarray}&&{a}^{3}=\displaystyle \frac{{{\rm{\Pi }}}^{3}{m}_{2}^{3}}{{\left({m}_{1}+{m}_{2}\right)}^{2}}{P}^{2},\end{eqnarray} \tag{ 13 }$$

where a is measured in arcseconds when P is measured in years, Π is the parallax in arcseconds, and m₂ is the mass of the unseen companion and m₁ is the mass of the star, both in solar masses.

The unknown parameters are a reference position, proper motions, acceleration terms, parallax, and 7 × n_p orbital elements, where n_p is the number of companions that are fitted (P, τ, e, ω, Ω, a, and i for each companion). There are a total of 14 free parameters in the case of a single companion and 21 free parameters for two companions, etc.

In the case of an astrometric binary (when the two stars in the binary are detected), it is necessary to fit the Keplerian parameters of the binary, the proper motions of the system's barycenter, and the parallax simultaneously (see above). For the secondary, ω is rotated 180°, and a₂ is used instead, which is scaled from a₁ by the mass ratio q = (m₂/m₁). In this case, two additional unknown parameters, a₂ and q, need to be taken into account to fit the secondary component, making, in this case, a total of 16 free parameters. However, in most cases only 14 free parameters are needed. The acceleration terms are only needed when there is something that perturbs the motion of this compact system, which could be the case when the binary is part of a wider multiple system.

In the case of having relative astrometry observations of the binary plus absolute observations of only one of the components in the binary (we will call it the main component), the fitting procedure is similar to the case of a binary system, but in this case the two additional parameters (ω₂ and a₂ or q = (m₂/m₁)) are used to fit the relative astrometry, and the other orbital parameters are the same as those used for the absolute astrometry: P, τ, e, ω, Ω, a, and i. In this case there are also 16 free parameters. In the case that the primary is the most massive component in the system (q ≪ 1), 7 additional parameters can be added to fit a second companion, making a total of 23 free parameters. Since the acceleration terms are only needed when the compact system is part of a wider multiple system (e.g., a triple or quadruple system, where the other stellar components are farther away from the observed compact system and have orbital periods much larger than the orbital period in the compact system), in general only 14 free parameters (in the case of a single companion) or 21 free parameters (in the case of two companions) will be fitted.

In this case, we follow the procedure described in the previous section, but here the orbital elements (ω, Ω, a₁, and i) are obtained using the Thiele–Innes elements X1, Y1, X2, and Y2 (Green 1993) instead of Equations (7) and (8). In this case,

$$\begin{eqnarray}&&{G}_{\alpha }(t)=X1x(t)+X2y(t),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{G}_{\delta }(t)=Y1x(t)+Y2y(t),\end{eqnarray} \tag{ 15 }$$

where the Thiele–Innes constants are expressed as

$$\begin{eqnarray}&&X1=a[\cos (\omega )\sin ({\rm{\Omega }})+\sin (\omega )\cos ({\rm{\Omega }})\cos (i)],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&Y1=a[\cos (\omega )\cos ({\rm{\Omega }})-\sin (\omega )\sin ({\rm{\Omega }})\cos (i)],\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&X2=a[-\sin (\omega )\sin ({\rm{\Omega }})+\cos (\omega )\cos ({\rm{\Omega }})\cos (i)],\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&Y2=a[-\sin (\omega )\cos ({\rm{\Omega }})-\cos (\omega )\sin ({\rm{\Omega }})\cos (i)].\end{eqnarray} \tag{ 19 }$$

The elliptical rectangular coordinates x(t) and y(t) are given in terms of the dynamical equations (Equation (9)) and the true anomaly ν (Equation (12)) by

$$\begin{eqnarray}&&x(t)=\left(\displaystyle \frac{r}{a}\right)\cos (\nu ),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&y(t)=\left(\displaystyle \frac{r}{a}\right)\sin (\nu ).\end{eqnarray} \tag{ 21 }$$

The linear parameters (reference position, proper motions, acceleration terms, parallax, and 4 × n_p orbital elements) are fitted by matrix inversion, while the nonlinear parameters (P, τ, and e) are found using a Monte Carlo method. We use the correlation matrix to obtain the uncertainty of the fitted parameters.

The AGA fitting procedure that we use here is similar to that described by Cantó et al. (2009) and Curiel et al. (2011), where the fitting procedure was used to find Keplerian orbits, using the RVs of the host star. This method can be extended to the problem of fitting the Keplerian elements of an orbit to astrometric data, such as those that can be obtained with the VLBA, GAIA, and in the future, with the ngVLA. As described by Cantó et al. (2009), the curve or model fitting is essentially an optimization problem. Given a discrete set of N data points (α_i, δ_i) with associated measurement errors σ_i, one seeks for the best possible model (in other words, the closest fit) for these data using a specific form of the fitting function, (α(t), δ(t)). This function has, in general, several adjustable parameters, whose values are obtained by minimizing a "merit function," which measures the agreement between the data (α_i, δ_i) and the model function (α(t), δ(t)). The maximum likelihood estimate of the model parameters (c₁, ..., c_k) is obtained by minimizing the χ² function:

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\min }^{2} & = & \displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\alpha }_{i}-\alpha ({t}_{i};{c}_{1},\; ...,{c}_{k})}{{\sigma }_{i}}\right)}^{2}\\ & & +\ \displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\delta }_{i}-\delta ({t}_{i};{c}_{1},\; ...,{c}_{k})}{{\sigma }_{i}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 22 }$$

where each data point (α_i, δ_i) has a measurement error that is independently random and distributed as a nominal distribution about the "true" model with standard deviation σ_i. We then apply the AGA method to search for the best solution using this initial guess. Here we have made two important improvements to the original algorithm. First, we do an initial search in the hypercube space of possible solutions to find an initial guess for the parameters that are fitted. Second, the error estimates, σ(c_j), for the fitted parameters c_j can be estimated as the projection of the confidence region of the m-dimensional space parameter for which χ² does not exceed the minimum value by an amount Δ(m, α), where α is the significance level (0 < α < 1). Following Avni (1976) and Wall & Jenkins (2003), the probability,

$$\begin{eqnarray}&&{\rm{p}}\mathrm{rob}[{\chi }^{2}-{\chi }_{\min }^{2}\leqslant {\rm{\Delta }}(m,\alpha )]=\alpha ,\end{eqnarray} \tag{ 23 }$$

is that of a chi-square distribution with m degrees of freedom. Thus, Δ(m, α) is the increment of χ² such that if the observation is repeated a large number of times, a fraction α of times the values of the parameters fitted will be inside the confidence region, i.e., in the interval c_j ± σ(c_j) (see Estalella 2017, and references therein).

In this case, we also follow the procedure previously described in Section 3.2. This fitting procedure allows us to fit simultaneously all free parameters. However, in the case of a single star, since we do not know a priori if it has a single or multiple companions, we first fit the proper motions and the parallax of the star using the least-squares method and the AGA method, and if necessary, we include the acceleration terms. We then analyze the residuals in order to find if there may be an astrometric companion orbiting around the star. To do this, we compute the least-squares periodogram of the astrometric data comparing the null solution and the Keplerian solution (see Section 3.1). If the periodogram shows possible companions (e.g., significant peaks with FAPs smaller than 1%), we then include a possible companion in the fitting procedure of the data, using again both fitting methods: least-squares and AGA algorithms. We apply this procedure for each signal present in the periodogram. In the case that there is more than one signal in the periodogram, and there are enough observations (enough data points) to fit simultaneously all the required free parameters, we fit simultaneously more than one astrometric companion. In the case that more than one component is detected with the VLBA, we fit simultaneously the two components, either by doing full astrometry to the main component, or by combining the relative astrometry of the system and the absolute astrometry of one of the main components.

Here we follow the standard nomenclature to name the host star and the companions that were detected: we use a capital letter for the host star and low-mass stellar companions, and a lowercase letter for substellar companions. Thus, in what follows, DoAr21 or DoAr21A is the host star, DoAr21B is the low-mass stellar companion, DoAr21b is the inner substellar companion, and DoAr21c is the outer substellar companion.

Before fitting the astrometric data that we present here, we used both algorithms, as well as the least-squares periodogram, to fit the data sets of several sources previously investigated by our team in the Ophiuchus region (Ortiz-León et al. 2017). We fitted single sources (when only one source was detected with the VLBA) and binary systems (when two sources were detected with the VLBA). We found that our results are in good agreement with those reported by Ortiz-León et al. (2017). As an example, we present in Table 2 the fit of the binary system LFAM15. The main difference between our fit and that obtained by Ortiz-León et al. (2017) is the time of the periastro of the orbit, which is offset by two orbital periods of the binary system. We obtained a periastro time near the center of the observing interval, while Ortiz-León et al. (2017) obtained a periastro time near the beginning of the observing interval.

Table 2.  Binary Astrometry Fits^a

	LFAM15
Parameter	This work	Ortiz-León et al. (2017)
	Parameters Fitted
μx (mas yr−1)	−6.303 ± 0.020	−6.31 ± 0.02
μy (mas yr−1)	−26.964 ± 0.050	−26.95 ± 0.05
Π (mas)	7.259 ± 0.079	7.253 ± 0.054
P (days)	1308.75 ± 10.39	1311.61 ± 6.68
T0 (yr)	2014.1931 ± 0.076	2007.008 ± 0.039
e	0.5292 ± 0.0075	0.528 ± 0.005
ω (deg)	56.09 ± 0.52	55.54 ± 1.02
Ω (deg)	338.01 ± 0.40	337.93 ± 0.81
a1 (mas)	7.760 ± 0.070	...
i (deg)	110.24 ± 0.27	110.30 ± 0.49
ω2 (deg)	236.09 ± 0.52	235.54 ± 1.02
a2 (mas)	8.637 ± 0.091	...
a (mas)	16.40 ± 0.16	16.40 ± 0.13
	Other Parameters
D (pc)	137.77 ± 1.48	137.9 ± 1.0
m (M⊙)	0.898 ± 0.042	0.89 ± 0.01
m1 (M⊙)	0.473 ± 0.022	0.469 ± 0.015
m2 (M⊙)	0.425 ± 0.019	0.421 ± 0.010
a1 (au)	1.069 ± 0.021	...
a2 (au)	1.190 ± 0.025	...
χ2, ${\chi }_{\mathrm{red}}^{2}$	75.53, 3.43	...

Note.

^aThe parameters presented here were obtained with AGA. Very similar results were obtained with the least-squares fitting method. The subindex 1 corresponds to the main component and the subindex 2 corresponds to the secondary component.

Download table as:  ASCIITypeset image

The least-squares periodogram with a circular orbit (CLSP) of the astrometric data (see Figure 2) shows a prominent peak (with 999 days period) with high power and very low FAP. There is another significant peak (with 87 days period) with lower power and higher FAPs. There is also a third weaker peak that appears in the CLSP periodogram (see Figure 2, upper panel) with a period of 151 days. However, when we fix the eccentricity to 0.368 (see discussion below), the peak with a period of about 999 days becomes more prominent, while this weak peak remains constant and a somewhat stronger peak appears with a period of 129 days (see Figure 2). This new peak is quite weak, as shown in the upper panel of Figure 2. Thus these two weaker peaks are dubious, while the peak with a period of about 87 days remains consistent in both periodograms. Thus, we investigate here only the two candidates with the strongest peaks in the periodograms (with periods of about 999 and 87 days) to obtain their "significances."

We use the two methods described above (least-squares and AGA) to fit the astrometric data. First we use both methods to fit the 11 astrometric observations to obtain the proper motions and the parallax of this source without taking into account any companion (single star solution). The results of a single star solution are shown in Table 3 and Figure 3. We did not find it necessary to include an acceleration term for the fitting of the external companion DoAr21c. However, in the case of the internal companion DoAr21B, we included acceleration terms to take into account the contribution of the external companion. The acceleration terms are small and produce a small effect in the fitting of DoAr21B.

Zoom In Zoom Out Reset image size

Table 3.  Absolute Astrometry Fits: One Companion^a

Parameter	Single Star Solution	DoAr21c	DoAr21B
		Parameters Fitted
μx (mas yr−1)	−19.6953 ± 0.0020	−19.7089 ± 0.0028	−19.5565 ± 0.0028
μy (mas yr−1)	−26.9477 ± 0.0050	−26.9637 ± 0.0069	−26.8091 ± 0.0069
accx (mas yr−2)	...	...	0.0335 ± 0.0020
accy (mas yr−2)	...	...	−0.0015 ± 0.0048
Π (mas)	7.5330 ± 0.0068	7.5001 ± 0.0095	7.4482 ± 0.0095
P (days)	...	1018.43 ± 3.98	92.892 ± 0.030
T0 (days)	...	2457163.20 ± 3.89	2457293.95 ± 0.31
e	...	0.368 ± 0.036	0.3272 ± 0.0095
ω (deg)	...	236.90 ± 1.91	32.24 ± 1.05
Ω (deg)	...	40.41 ± 2.82	41.43 ± 1.09
a1 (mas)	...	0.377 ± 0.024	0.706 ± 0.014
i (deg)	...	102.25 ± 1.83	89.41 ± 0.96
		Other Parameters
D (pc)	132.74 ± 0.12	133.33 ± 0.17	134.26 ± 0.17
m (M⊙)b	...	2.9441	2.8745
m1 (M⊙)	...	2.8920 ± 0.0032	2.397 ± 0.010
m2 (M⊙)	...	0.0521 ± 0.0032	0.477 ± 0.010
a1 (au)	...	0.0503 ± 0.0032	0.0947 ± 0.0021
a2 (au)	...	2.7890 ± 0.0041	0.4760 ± 0.0019
χ2, ${\chi }_{\mathrm{red}}^{2}$	2309.29, 128.29	91.02, 8.27	70.80, 7.87

Notes.

^aThe parameters presented here were obtained with AGA. Very similar results were obtained with the least-squares fitting method. The subindex 1 corresponds to the main component (i.e., the star) and the subindex 2 corresponds to the secondary component (i.e., companion DoAr21c or DoAr21B). ^bThe mass of the system is fixed and corresponds to the masses obtained from the relative plus absolute, simultaneous astrometry fit of components b and c (see Section 4.3). In the case of component c we use the total mass obtained from the simultaneous fit of components b and c. In the case of component B we used the inner mass of the system, after removing the masses of components b and c.

Download table as:  ASCIITypeset image

We find that after the fitting the residuals are large compared with the noise present in the data and the astrometric precision obtained with the VLBA. Figure 3 shows residuals up to about 0.3 mas and the astrometric precision obtained with the VLBA is about 30 μas. Furthermore, the residuals have a temporal trend that suggests the presence of at least one companion. We then use the least-squares method and the AGA algorithms to fit the astrometric observations of this source, including a single companion (i.e., we performed an independent fit for each companion candidate). In both cases we obtain consistent solutions. The parallax and the orbital fits for the two candidates are presented in Figure 4. Table 3 summarizes the best fits and their ${\chi }_{\mathrm{red}}^{2}$ per degree of freedom (${\chi }_{\mathrm{red}}^{2}$ = χ²/(N_data − N_par − 1), where N_data = 2 × N_points and N_par is the number of fitted parameters). The fit of the astrometric data clearly improves when including a companion, as seen by the ${\chi }_{\mathrm{red}}^{2}$.

Zoom In Zoom Out Reset image size

The two components have a similar position angle of the line of nodes (Ω), within the errors. The eccentricity of the orbits seem to be reasonably well constrained: component c has an orbital eccentricity ∼0.37, while component DoAr21B has an eccentricity of ∼0.33. The astrometric signature of the source (a_*) due to the gravitational pull of the companion is larger in the case of component DoAr21B, by almost a factor of 2. The inclination of the orbit (i) is somewhat different for the two companions. The inclination angle is larger than 90° for component DoAr21c (∼102°, which indicates a retrograde orbit), while component B has an orbital inclination of about 90° (indicating an edge-on orbit), which suggests that the orbits of these two companions are not coplanar. This result is somewhat surprising since one would expect a similar inclination for the orbits of all the companions. We do not have an explanation for this result. However, it is important to point out that the combined fit (relative astrometry plus absolute astrometry) of two components (components DoAr21b and DoAr21c; see Section 4.3) also give retrograde orbits (i > 90°). With these fits we cannot estimate the dynamical mass of the system, thus to estimate the mass of the companions we use the mass of the system obtained with the combined fit of the astrometric data (a fixed mass of M ∼ 2.94 M_⊙ for the fit of component DoAr21c and ∼2.87 M_⊙ for the fit of component DoAr21B). These masses were obtained with the combined fits of the orbits of components DoAr21b and DoAr21c (see Section 4.3). With this assumption, we obtain that the component with a compact orbit (component DoAr21B) has a mass consistent with a low-mass star (m_B ∼ 0.48 m_⊙), while component DoAr21c has a mass consistent with a brown dwarf (m_c ∼ 54.57 m_jup). The orbits of the two companions have semimajor axis a_c ∼ 2.79 and a_B ∼ 0.48 au, respectively. This suggests that the mass of the companion decreases with the distance to the source. The brown dwarf is the component with the wider orbit. The estimated distance to the source is similar in both cases, however, the estimated distance is slightly larger when including a companion than when fitting only the proper motions and the parallax of the host star.

The orbital periods of the orbits of components DoAr21c and DoAr21B are consistent with those found in the CLSP (see Figure 2). The difference in the orbital periods obtained with the astrometric fits and the CLSP is due to the fact that in the case of CLSP we assume circular orbits, while in the astrometric fits we include the eccentricity of the orbits as a free parameter to be fitted. In Figure 2, we also show the least-squares periodogram of the observed astrometric data fixing the eccentricity obtained for companion DoAr21c (e_c ∼ 0.368). In this case we obtain an orbital period for component DoAr21c that is consistent with that obtained with the astrometry fit. This figure shows that the power of component c increases substantially when fixing the eccentricity of this companion. On the other hand, the power of component DoAr21B does not change substantially when using this eccentricity in the periodogram.

Figure 4 shows the orbits of both companions. As was mentioned before, a companion was detected with the VLBA in three different epochs. The relative position of the companion with respect to the host star is included in Figure 4. This figure shows that the relative position of none of the detected companions coincide with the fitted orbits. Bellow we attempt to fit simultaneously the orbits of the two inferred components (DoAr21c and DoArt21B). We also attempt to fit simultaneously the relative astrometry of the companion detected with the VLBA (DoAr21b) and the absolute astrometry of the host star (combined fit; Section 4.3).

The 11 observed epochs are in principle enough to simultaneously fit the orbits of the two companions (DoAr21c and DoAr21B), plus the parallax and the proper motions of the host star. We find that the astrometry fit improves when simultaneously fitting both companions. Table 4 summarizes the parameters obtained with the best two-companion astrometry fit of the data. Figure 5 shows the parallax of the host star and the orbital motion of the host star due to the gravitational pull of both companions. This figure also shows the orbits of both companions. For this fit we have also fixed the total mass of the system (m = 2.9441 m_⊙; see Section 4.3 and Table 6). The fitted parameters are similar to those found from the single-companion fit. However, some of the fitted parameters change substantially. For instance, the eccentricity, semimajor axis, and inclination of the orbits are somewhat larger, and the estimated masses of the companions are somewhat smaller in this case than in the case of a single-companion astrometry fit. Furthermore, the position angle of the line of nodes of both companions differs by more than 20° (Ω ∼ 32.7 and 59.7 for DoAr21B and DoArt21c, respectively), while in the case of a single-companion fit we obtained a similar Ω for both companions (Ω ≃ 41°). These new results suggest that the orbits of these two companions are not coplanar.

Zoom In Zoom Out Reset image size

Table 4.  Absolute Astrometry Fit: Two Companions^a

Parameter	DoAr21c and DoAr21B
	Fitted Parameters
μx (mas yr−1)	−19.6442 ± 0.0035
μx (mas yr−1)	−26.8756 ± 0.0084
Π (mas)	7.474 ± 0.012
P (days)	1034.50 ± 7.74
T0 (days)	2457157.64 ± 6.94
e	0.404 ± 0.065
ω (deg)	240.43 ± 3.52
Ω (deg)	59.72 ± 6.98
a1 (mas)	0.203 ± 0.024
i (deg)	109.99 ± 5.62
P (days)	92.921 ± 0.091
T0 (days)	2457295.80 ± 0.94
e	0.391 ± 0.029
ω (deg)	39.61 ± 2.78
Ω (deg)	32.73 ± 2.50
a1 (mas)	0.346 ± 0.020
i (deg)	93.03 ± 2.03
	Other Parameters
D (pc)	133.80 ± 0.21
m (M⊙)b	2.9441
m1 (M⊙)	2.681 ± 0.014
m2 (M⊙)	0.0279 ± 0.0031
m3 (M⊙)	0.236 ± 0.014
a1–2 (au)	2.869 ± 0.014
a1 (au)	0.0271 ± 0.0032
a2 (au)	2.842 ± 0.011
a1–3 (au)	0.5736 ± 0.0004
a1 (au)	0.0463 ± 0.0027
a3 (au)	0.5273 ± 0.0024
χ2, ${\chi }_{\mathrm{red}}^{2}$	31.72, 7.93

Notes.

^aThe parameters presented here were obtained with AGA. In this case, two components were fitted simultaneously (components B and c) using the absolute astrometry model. The subindex 1 corresponds to the main component (i.e., the star) and the subindices 2 and 3 correspond to the companions (DoAr21c and DoAr21B). ^bThe mass of the system is fixed and corresponds to the mass obtained from the relative plus absolute, simultaneous astrometry fit of components DoAr21b and DoAr21c (see Section 4.3).

Download table as:  ASCIITypeset image

In addition, Figure 5 shows that the companion detected with the VLBA does not coincide with the estimated orbits of these two companions. We obtained a similar result from the single-companion astrometry fit. This suggests that there is probably another companion that does not appear in the periodogram, but that was detected with the VLBA.

Bellow we attempt to fit simultaneously the relative astrometry of the companion observed with the VLBA and the absolute astrometry of the host star (combined fit).

We investigate here the possibility that the relative positions of the detected companion are associated with another companion, different from the two companions that we have already found. We find that the three secondary detections with the VLBA can be fitted by combining relative astrometry and absolute astrometry (combined model) using both the least-squares and the AGA algorithms. These three detections correspond to a new component that we call DoAr21b. Table 5 summarizes the parameters that were obtained with the best fit of the data and the ${\chi }_{\mathrm{red}}^{2}$ per degree of freedom. The parallax and the astrometric signature of the host star due to the gravitational pull of this companion are presented in Figure 6. The orbital fit of this companion is also presented in Figure 6. Since DoAr21b is an internal companion, we also included acceleration terms to take into account the contribution of the external companion DoAr21c. The acceleration terms are small and produce a small effect in the fitting of DoAr21b.

Zoom In Zoom Out Reset image size

Table 5.  Relative Plus Absolute Astrometry Fit^a

Parameter	DoAr21b
	Fitted Parameters
μx (mas yr−1)	−19.6682 ± 0.0030
μx (mas yr−1)	−26.9746 ± 0.0074
accx (mas yr−2)	0.0057 ± 0.0021
accy (mas yr−2)	−0.0091 ± 0.0052
Π (mas)	7.350 ± 0.010
P (days)	452.97 ± 0.15
T0 (days)	2457507.45 ± 1.07
e	0.076 ± 0.012
ω (deg)	220.90 ± 0.91
Ω (deg)	53.20 ± 1.78
a1 (mas)	0.286 ± 0.015
i (deg)	103.57 ± 1.17
q (m1/m2)	39.86 ± 2.67
ω2 (deg)	40.90 ± 0.91
a (mas)	11.68 ± 0.44
	Other Parameters
D (pc)	136.05 ± 0.19
m (M⊙)	2.6064 ± 0.0018
m1 (M⊙)	2.5426 ± 0.0029
m2 (M⊙)	0.0638 ± 0.0035
a1–2 (au)	1.5886 ± 0.0007
a1 (au)	0.0389 ± 0.0021
a2 (au)	1.5497 ± 0.0014
χ2, ${\chi }_{\mathrm{red}}^{2}$	903.11, 64.51

Note.

^aThe parameters presented here were obtained with AGA. Very similar results were obtained with the least-squares fitting method. The subindex 1 corresponds to the main component (i.e., the star) and the subindex 2 corresponds to the secondary component (i.e., companion DoAr21b).

Download table as:  ASCIITypeset image

The combined fitted orbit of DoAr21b has a period of P_b ∼ 453 days, a semimajor axis of a₁ ∼ 0.29 mas (or ∼0.04 au), a position angle of the line of nodes Ω ∼ 53°, an eccentricity e ∼ 0.08, an inclination angle i ∼ 104°, and a longitude of the periastron ω ∼ 221°. Component DoAr21b has an orbit with a semimajor axis of a_b ∼ 11.40 mas (or ∼1.55 au) and a longitude of the periastron ω_b ∼ 41°. Since we are doing a combined fit (relative plus absolute orbital fit), we obtain the orbital fit around the barycenter position of the system for both the main source and the companion, thus we obtain the dynamical mass of the system, as well as the masses of each individual component. The dynamical mass of the system is M_⋆b ∼ 2.606 M_⊙, the mass of the main source is M_⋆ ∼ 2.542 M_⊙, and the mass of the companion is M_b ∼ 0.064 M_⊙ (or about 66.8 M_jup). The orbit of this new companion, DoAr21b, lies between the orbits of the other two companions, and its estimated mass is consistent with being a brown dwarf.

We find that the position angle of the line of nodes (Ω) and the inclination angle (i) of the orbit are similar to those found in DoAr21c obtained with the two-companion simultaneous fit (see Table 4). The eccentricity of the orbit (e) is smaller than those found for the other companions, and seems to be well constrained. The astrometric signature of the source (a_⋆) due to the gravitational pull of this companion is somewhat smaller than those found for the other companions (see Tables 3–5). Since the astrometric signature of the source is relatively small, the residuals of the fit (see Figure 6) are larger than those obtained from the fits of the other two companions (see Figure 5). In addition, since the astrometric signature of DoAr21b is quite small (<0.1 mas; see Tables 5 and 6) the least-squares fit of the orbit of this companion basically fails. In other words, the χ² of the fit is similar to the null solution, and thus the estimated power is small, within the expected "noise" in the periodogram (see discussion in Section 3.1). These results are consistent with the fact that companion DoAr21b does not appear in the periodogram (see Figure 2). In other words, we have found DoAr21b only because this companion was detected at several epochs with the VLBA, otherwise the astrometric signature due to this companion would be embedded in the residuals of the fits of the other two companions.

Table 6.  Relative Plus Absolute, Simultaneous Astrometry Fit^a

Parameter	DoAr21b and c	DoAr21b and B
	Fitted Parameters
μx (mas yr−1)	−19.6986 ± 0.0036	−19.5478 ± 0.0036
μx (mas yr−1)	−26.9768 ± 0.0088	−26.8373 ± 0.0088
accx (mas yr−2)	...	0.0232 ± 0.0025
accy (mas yr−2)	...	−0.0012 ± 0.0061
Π (mas)	7.428 ± 0.012	7.385 ± 0.012
P (days)	446.52 ± 0.19	453.22 ± 0.22
T0 (days)	2457045.94 ± 1.22	2457063.16 ± 1.38
e	0.090 ± 0.036	0.095 ± 0.018
ω (deg)	250.02 ± 1.11	225.86 ± 1.20
Ω (deg)	55.80 ± 2.39	55.41 ± 2.43
a1 (mas)	0.099 ± 0.015	0.073 ± 0.018
i (deg)	104.19 ± 1.52	105.43 ± 1.67
q (m1/m2)	121.91 ± 19.16	151.89 ± 37.59
ω2 (deg)	70.02 ± 1.11	45.86 ± 1.20
a1–2 (mas)	12.11 ± 0.54	11.13 ± 0.50
P (days)	984.93 ± 6.39	92.935 ± 0.041
T0 (days)	2457158.61 ± 6.34	2457295.83 ± 0.43
e	0.234 ± 0.060	0.355 ± 0.013
ω (deg)	243.17 ± 2.78	41.35 ± 1.31
Ω (deg)	44.12 ± 4.36	46.06 ± 1.60
a1 (mas)	0.321 ± 0.026	0.650 ± 0.018
i (deg)	99.53 ± 2.80	90.07 ± 1.34
	Other Parameters
D (pc)	134.62 ± 0.22	135.40 ± 0.22
m (M⊙)	2.9477 ± 0.0025	2.2246 ± 0.0021
mA, mA (M⊙)	2.8782 ± 0.0027	1.838 ± 0.018
mb, mb (M⊙)	0.0236 ± 0.0036	0.0146 ± 0.0036
mc, mB (M⊙)	0.0459 ± 0.0032	0.372 ± 0.012
aA−b, aA−b (au)	1.6308 ± 0.0009	1.5074 ± 0.0010
aA, aA (au)	0.0133 ± 0.0020	0.0099 ± 0.0024
ab, ab (au)	1.6175 ± 0.0011	1.4976 ± 0.0014
aA−c, aA−B (au)	2.7778 ± 0.0092	0.3231 ± 0.0004
aA, aA (au)	0.0433 ± 0.0036	0.0880 ± 0.0025
ac, aB (au)	2.7346 ± 0.0056	0.4350 ± 0.0021
χ2, ${\chi }_{\mathrm{red}}^{2}$	56.13, 6.24	60.32, 8.62

Note.

^aThe parameters presented here were obtained with AGA. In this case, two components were fitted simultaneously (components b and c, and components b and B) using the combined model. The subindex 1 corresponds to the main component (i.e., the star) and subindices 2 and 3 correspond to the two companions (i.e., DoAr21b and DoAr21c, or DoAr21b and DoAr21B).

Download table as:  ASCIITypeset image

The 11 observed epochs and the three direct detections of DoAr21b are in principle enough to simultaneously fit the orbits of two companions (DoAr21b and DoAr21c, or DoAr21b and DoAr21B), plus the parallax and the proper motions of the host star. Table 6 summarizes the parameters that were obtained with the best fit of the data and the ${\chi }_{\mathrm{red}}^{2}$ per degree of freedom. The parallax and the astrometric signature of the host star due to the gravitational pull of each companion are presented in Figures 7 and 8. The orbital fit for each pair of companions is also presented in these figures. For the simultaneous fitting of DoAr21b and DoAr21B, we included acceleration terms to take into account the contribution of the external companion DoAr21c. The acceleration terms are small and produce a small effect in the fitting of these two companions. We find that the fitted parameters are in general consistent with the previous fits. In particular, the astrometric fit of DoAr21b is quite consistent in all the fits of this component. This is because the relative astrometry of this companion, using the three detected epochs, constrain the orbit of this companion. The astrometric fits of the other two companions are similar to those obtained previously, however, there are some differences. For instance, the position angle of the line of nodes (Ω), the inclination angle (i), and the semimajor axis (a₁) of the orbit are somewhat different from those found previously (see Table 4). It is important to mention that the estimated mass of the system is somewhat different in both fits, even if we take into account the mass of the companion DoAr21c in the simultaneous fit of DoAr21b and DoAr21B. The difference in the estimated total mass is ∼0.72 M_⊙ (or ∼0.68 M_⊙, taking into account the mass of DoAr21c in both fits), which is much higher than the estimated mass of DoAr21c (∼0.046 M_⊙). A better estimate of the masses in this multiple system can be obtained by simultaneously fitting the orbits of all the components in the system. Thus, further observations of this multiple system will be needed in order to obtain a better estimate of the total mass of the system and the individual masses of all the companions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows the astrometric signature of the host star due to the companions DoAr21b and DoAr21c, and the parallax of the system as a function of time, and Figure 8 shows the astrometry of DoAr21b and DoAr21B. Figure 7 suggests that the orbits of DoAr21b and DoAr21c cross each other, however, the results presented in Table 6 indicate that this is a projection effect due to the large inclination of the orbits (i ∼ 104° and ∼100°, respectively) and the difference in the position angle of the line of nodes of both orbits (Ω ∼ 56° and ∼44°, respectively). The semimajor axis of the orbit of DoAr21c (a_c ∼ 2.73 au) is nearly a factor of 2 larger than that of DoAr21b (a_b ∼ 1.62 au), and their orbits have low eccentricities (e_b ∼ 0.09 and e_c ∼ 0.23). Thus the orbits of these two companions do not cross each other. In addition, the orbit of DoAr21B (a_B ∼ 0.44 au) is about a factor of 3 smaller than that of DoAr21b (a_b ∼ 1.50 au). Therefore the orbits of this multiple system seem to scale with a relationship close to 3:1 between DoAr21B and DoAr21b, and close to 2:1 between DoAr21b and DoAr21c.

The inclination angle of the three companions is not the same. The orbit of DoAr21B appears to be edge-on (i ∼ 901), while the orbits of DoAr21b and DoAr21c are retrograde, with inclination angles ∼104° and 100°, respectively. Thus the inclination angle of the system seem to be about 97° ± 7°. However, this large difference in the inclination angle will have to be confirmed with further observations.

Since a companion (DoAr21b) was detected in three epochs, we were able to simultaneously fit the relative position of the companion and the absolute position of the host star. Furthermore, we were able to simultaneously fit the orbits of the pairs of companions DoAr21b−DoAr21c and DoAr21b−DoAr21B. With this relative plus absolute simultaneous combined fit, we have obtained the dynamical mass of the system, as well as the mass of the individual companions. The combined fit of the pairs of components give the dynamical masses for the system of ∼2.948 M_⊙ for the pair DoAr21b−DoAr21c and ∼2.225 M_⊙ for the pair DoAr21b−DoAr21B. There is a significant difference in the estimated mass of the system of about 0.68 M_⊙ (taking into account the mass of DoAr21c in both fits). This discrepancy is probably due to the astrometric contribution of the other companion (which is not taken into account in the fit), and the fact that DoAr21b was detected at only three epochs. We consider that at a first approximation of the total mass of the system is probably somewhere between these two values.

There is also a large difference in the estimates mass of the host star between 1.84 and 2.51 M_⊙ (we have taken into account the mass of DoAr21B in both estimates; see Table 6). The difference of about 0.67 M_⊙ is similar to the difference of about 0.68 M_⊙ found for the total mass of the system. Thus, we also consider that the mass of the star is somewhere between these two estimated masses. The estimated mass of DoAr21 is similar to the previously estimated masses of ∼2.2 and ∼1.8 M_⊙ (Jensen et al. 2009), which were obtained by assuming that all the luminosity of the source is associated with a single star and by splitting the luminosity of this source in two equal stars, respectively. The dynamical mass that we obtained here is, at present, the best estimated mass of the WTTS DoAr21.

Tables 3–6 show that the best fitted parameters of the three companions of DoAr21 do not change considerably when including one or more companions in the astrometric fit. However, taking into account that the different fits give slightly different values for all the orbital parameters, and that the orbits of the companions are not yet fully constrained, we have calculated the weighted average of the orbital parameters for all the companions, as well as the weighted average values for the mass of the host star and the companions. The weighted average parameters are presented in Table 7.

The estimated masses of the host star and the three companions are: 2.04 ± 0.70 M_⊙ for the host star, 0.35 ± 0.12 M_⊙ for DoAr21B, 0.034 ± 0.026 M_⊙ (35.6 ± 27.2 M_jup) for DoAr21b, and 0.042 ± 0.013 M_⊙ (44.0 ± 13.6 M_jup) for DoAr21c (see Table 7). The inner companion, DoAr21B, has an estimated mass consistent with a low-mass star, while the other two companions have masses consistent with being brown dwarfs.

The masses of DoAr21b and DoAr21c are probably consistent with a spectral types M7-M8 (e.g., Luhman et al. 2009, and references therein), however, the classification of brown dwarfs is based on their spectral type, in particular based on the elements seen in their optical spectrum. These two brown dwarfs cannot be classified this way since they, in principle, cannot be separated from the host star in the optical nor in the infrared. In addition, it is expected that brown dwarfs change in class type as they burn their deuterium (and their lithium in the case of the most massive brown dwarfs). Since DoAr21 is estimated to have an age of about 0.4–0.8 Myr (Jensen et al. 2009), the two brown dwarfs orbiting this WTTS must be extremely young. In addition, these two brown dwarfs have broad orbits, with semimajor axis of ∼1.6 and ∼2.8 au (see Table 7). These results suggest that these brown dwarfs were formed far away from the host star, probably close to their current orbits. This is consistent with the idea that these brown dwarfs were formed by fragmentation of the disk where this young star was formed. It is not clear if DoAr21B was formed by fragmentation, and it is also not clear if it was formed close to its current orbit.

DoAr21 presents a large variability at radio wavelengths (see Table 8). Figure 9 shows that the radio flux of this WTTS changes in more than two orders of magnitude, from a fraction of a mJy to several tens of mJy. This figure shows that there are abrupt changes in the flux in intervals of several weeks. However, this change could occur in much shorter times, and the cadence of the observations is not adequate to find how steep are the outbursts observed in this source, nor if there is a frequency in the outbursts. It has been speculated that this variability is probably due to the interaction between the magnetosphere of the star an the magnetosphere of a close stellar companion. We have found three astrometric companions to DoAr21. DoAr21B is the most massive companion (m_B ∼ 0.35 M_⊙) and with the most compact and eccentric orbit (a_B ∼ 0.482 au, e_B ∼ 0.37) in this multiple system. Thus, DoAr21B could be responsible for the large flux variation observed in DoArt21. Given the large eccentricity of the orbit of this companion, and its proximity to the star, the interaction of their magnetospheres probably occurs when the companion is closer to the star, which occurs close to the periastron of the orbit.

Zoom In Zoom Out Reset image size

To investigate this possibility, we have first calculated the closest distance between DoAr21A and DoAr21B, which occurs close to the periastron of their orbital motions. We obtain that the closest distance between them is a_p = a(1 − e) = 0.34 au, where a = a_⋆ + a_B ∼ 0.536 au is the relative semimajor axis of the orbit of DoAr21B around the host star DoAr21, and e ∼ 0.37 is the eccentricity of the orbit (see Table 7). The closest separations between the two stars is much larger than the typical magnetosphere sizes of the stars, which are in the range of a few stellar radii (e.g., Feigelson & Montmerle 1999). Second, we investigated the possibility of a correlation between the period of the orbit of DoAr21B and the variability of the star. We have calculated the expected position of DoAr21B in its orbit for each observed epoch. In Figure 9, we present the observed fluxes of DoAr21 folded with the period of DoAr21B (P_B ∼ 93.0 days). We do not find a clear correlation between the outbursts observed in DoAr21 and the orbital period of this component. Figure 9 shows that the two strongest peak fluxes nearly coincide in the orbital phase of this companion, about halfway between the periastron and the apoastron in the orbit of DoAr21B. However, one would expect that the outbursts would occur when the low-mass stellar companion is closest to the star (near the periastron of its orbit). If DoAr21B were the responsible of the outbursts, then the outburst occurs about 20 days after this low-mass star has passed its periastron. These results indicate that DoAr21B is probably not responsible for the variability nor the X-ray outburst observed in DoAr21. This suggests that there may be another companion with a more compact orbit than those of the companions we report here that perturbs the magnetosphere of the host star. Further observations will be needed to search for this putative companion in a very compact orbit around DoAr21.

The substellar companion DoAr21b was detected at three epochs at radio wavelengths with the VLBA. This suggests that this brown dwarf is highly variable at radio wavelengths. The other brown dwarf in this system, DoAr21c, which is in a more extended orbit than DoAr21b, was not detected at any of the observed epochs. It is not clear why DoAr21b was directly detected with the VLBA, while DoAr21c was not.

This is the first time that a brown dwarf orbiting a young star has been directly detected. Several brown dwarfs have been detected with the VLA and with the VLBA, but most of them are single stars (they are not part of an stellar system) and are located close by (a few tens of parsecs away), and have ages of more that 100 Myr (e.g., Berger 2002; Sahlmann et al. 2013 and Forbrich et al. 2016). A few brown dwarf candidates have been found at radio wavelengths (e.g., Dzib et al. 2013; Rodríguez et al. 2017), and they are also isolated sources. A putative brown dwarf was recently found to be orbiting a TTauri star (Ginski et al. 2018), however, its orbit lies outside the circumbinary disk, with a projected orbit of about 210 au. The estimated mass of this companion is quite uncertain, and it could be a very low-mass star.

Tables 3–6 show that the estimated distance to DoAr21 does not change substantially when including one or more companions in the astrometric fit. Taking into account that the different fits give slightly different values for the distance, and that the orbits of the companions are not yet fully constrained, we have calculated the weighted average of the estimated distances. We obtained that the weighted distance is d = 134.6 ± 1.0 pc, where the estimated error corresponds to the standard deviation of the fitted values, which better reflects the dispersion seen in the different astrometric fits. This estimate is an improvement to the distance to this source of 135.76 ± 2.27 pc, previously obtained with the VLBA (Ortiz-León et al. 2018). The distance to DoAr21 that we obtain is in agreement, within the estimated errors, with that recently obtained by Ortiz-León et al. (2017). The error that we obtain here, however, it is a few times smaller than those obtained previously for this source, and the typical errors for the other sources in the Ophiucus Complex, obtained by Ortiz-León et al. (2017). This is probably due to the larger number of observations used for the present astrometric fit, the accuracy of the observations we present here (see Table 1) and a better coverage of the parallax, as well as the inclusion of the astrometric signal of one or more companions, which was not taken into account by Ortiz-León et al. (2017).

We have obtained the astrometric best fits for the independent Keplerian orbits of the three DoAr21 companions. These solutions can be used to estimate an expected induced RV of the star due to the gravitational pull of each companion. Assuming a simple model of totally independent companions (e.g., Cantó et al. 2009), the expected induced RV would be

$$\begin{eqnarray}&&{K}_{j}={\left(\displaystyle \frac{2\pi G}{{T}_{j}}\right)}^{1/3}\displaystyle \frac{{m}_{j}\sin ({i}_{j})}{{\left({M}_{\star }+{m}_{j}\right)}^{2/3}}\displaystyle \frac{1}{\sqrt{1-{e}_{j}^{2}}},\end{eqnarray} \tag{ 24 }$$

where G is the gravitational constant, and T_j, M_⋆, m_j, and e_j are the orbital period, the star and companion masses, and the eccentricity of the orbit of the companion. Using the parameters presented in Sections 5.1 and 5.3 and those presented in Table 7, we obtain that the maximum induced velocities on DoAr21 by DoAr21B, DoAr21b, and DoAr21c are K ∼ 9.902, 0.564, and 0.557 km s⁻¹, respectively. These RVs could in principle be observed with high-spectral resolution spectrographs. However, TTauri stars are magnetically active, may have broad (molecular) spectral features, and present ubiquitous variability, that would make very difficult to observe such RVs. Recent optical spectroscopic observations of DoAr21 have shown possible velocity variations of ∼4.9 km s⁻¹ over the course of about 2 hr, but with a quite low precision (∼1–2 km s⁻¹) due to the high rotation velocity of this star (James et al. 2016). This velocity variation, if real, may suggest a much shorter orbital period than that estimated for DoAr21B. Future short- and long-term, high-resolution spectroscopic observations of DoAr21 may show whether this short period signal is real or not, and furthermore, may be able to detect the ∼9.9 km s⁻¹ RV signature that we find that DoAr21B probably induces on DoAr21.