3D Orbital Architecture of a Dwarf Binary System and Its Planetary Companion

GJ 896AB (a.k.a. EQ Peg, BD+19 5116, J23318+199, HIP 116132) is a nearby low-mass M dwarf binary system with an estimated age of 950 Myr (Parsamyan 1995), at a distance of 6.25 pc, and a separation of 54 with a PA of 78° (e.g., Heintz 1984; Liefke et al. 2008; Bower 2011; Pearce et al. 2020). However, some observations also suggest that the age of this binary system is ≲100 Myr (Riedel et al. 2011; Zuckerman et al. 2013). If this latter age is confirmed, these stars would be the closest pre-main-sequence stars known. This binary system is part of a quadruple system, where the other two stellar companions are further away from the main binary. Recent observations suggest that both stars have companions, but theirs orbits have not been characterized (e.g., Delfosse et al. 1999; Winters et al. 2021).

Both stars have been observed to have radio outbursts (e.g., Pallavicini et al. 1985; Jackson et al. 1989; Benz et al. 1995; Gagne et al. 1998; Crosley & Osten 2018a, 2018b; Villadsen & Hallinan 2019; Davis et al. 2020) and are flaring X-ray emitters (e.g., Robrade et al. 2014; Liefke et al. 2008). GJ 896A has been detected at several epochs with the VLBA showing compact and variable flux emission (Bower et al. 2009; Bower 2011).

An early determination of the orbital motion of the binary GJ 896AB suggested an orbital period of about 359 yr and a semimajor axis of 687 (e.g., Heintz 1984). More recent, and more accurate, relative astrometric fits indicate that the orbital parameters of this binary system are P ∼ 234 yr, a ∼5.3 arcsec, i ∼ 126°, Ω ∼ 77°, e ∼ 0.11, and ω ∼ 97° (Mason et al. 2001). However, as the time baseline used for these orbital determinations is ≲34% the estimated period, the orbital parameters were not well constrained. The more massive star, GJ 896A, is an M3.5 star with an estimated mass of 0.39 M_⊙, a radius of 0.35 R_⊙, a rotation period of 1.061 days, and a rotational inclination angle of 60°, while the less massive star, GJ 896B, is an M4.5 star with an estimated mass of 0.25 M_⊙, a radius of 0.25 R_⊙, a rotation period of 0.404 days, and a rotational inclination angle of about 60° ± 20° (e.g., Delfosse et al. 1999; Morin et al. 2008; Davison et al. 2015; Pearce et al. 2020).

Here we present the discovery of a Jovian-mass planetary companion to the young close-by (6.25 pc away from the Sun) M3.5 dwarf GJ 896A, which is the more massive star in the low-mass M dwarf binary system GJ 896AB with a separation of about 31.6 au. Section 2 describes the observations and data reduction. In Section 3 we present the methodology to fit the astrometric data. Section 4 presents the results, including the fitted 3D orbital architecture of this system. The results are discussed in Section 5, and our conclusions are given in Section 6. In Appendix A we discuss the possible contributions of the variability of the main star to the expected "jitter" of the star, and in Appendix B we present the posterior sampling of the combined astrometric fit solution using an MCMC sampler.

First, both the least-squares and the AGA algorithms were used to fit the proper motions and the parallax of GJ 896A, without taking into account any possible companion. The results of the absolute astrometric fit are shown in column 1 of Table 2 and in Figure 2. We find that the residuals are quite large and present an extended temporal trend that indicates the presence of a companion with an orbital period larger than the time span of the observations.

Zoom In Zoom Out Reset image size

Table 2. Absolute Astrometric Fits ^a

Parameters	Single Source	Single Source+Acceleration	Single Companion
	(1)	(2)	(3)
		Fitted Parameters
μα (mas yr−1)	576.707 ± 0.014	574.171 ± 0.014	574.142 ± 0.019
μδ (mas yr−1)	−59.973 ± 0.014	−60.333 ± 0.014	−60.354 ± 0.020
aα (mas yr−2)	...	0.8893 ± 0.0034	0.8871 ± 0.0047
aδ (mas yr−2)	...	0.1402 ± 0.0035	0.1400 ± 0.0049
Π (mas)	161.136 ± 0.094	160.027 ± 0.094	159.83 ± 0.13
P (days)	...	...	281.56 ± 1.67
T0 (Julian day)	...	...	2,455,696.5 ± 10.4
e	...	...	0.30 ± 0.11
ω (deg)	...	...	344.1 ± 13.1
Ω (deg)	...	...	47.7 ± 12.1
aA (mas)	...	...	0.52 ± 0.11
i (deg)	...	...	66.0 ± 15.0
		Other Parameters
D (pc)	6.2059 ± 0.0036	6.2489 ± 0.0037	6.2565 ± 0.0051
mAb (M⊙) b	...	...	0.43814 (fixed)
mA (M⊙)	...	...	0.43589 ± 0.00047
mb (M⊙).	...	...	0.00225 ± 0.00047
mb (MJ )	...	...	2.35 ± 0.49
aAb (au)	...	...	0.6386 ± 0.0025
aA (au)	...	...	0.00328 ± 0.00070
ab (au)	...	...	0.6353 ± 0.0018
ab (mas)	...	...	101.53 ± 0.42
Δα(mas) c	8.51	0.24	0.14
Δδ(mas) c	1.51	0.24	0.13
χ2, ${\chi }_{\mathrm{red}}^{2}$ d	105728.53, 4066.48	190.91, 7.95	57.43, 3.38
BIC e	105,745.86	215.17	105.93
AIC e	105,738.53	204.91	85.41

Notes.

^a The parameters presented here were obtained with AGA. Very similar results were obtained with the least-squares fitting method. Column (1) contain the astrometric fit of the VLBA data without acceleration terms. Column (2) includes acceleration terms. Column (3) corresponds to the full astrometric fit, where a single planetary companion is included. ^b The mass of GJ 896A, including possible companions, was taken from the combined fit presented in Table 3. ^c The rms dispersion of the residual. ^d χ² and reduced χ² of the astrometric fit. The residuals clearly improve when we include acceleration terms and the companion. ^e The Bayesian information criteria (BIC) and Akaike information criterion (AIC) are the criteria used to choose the best-fitted model (e.g., Liddle 2007). The smaller the BIC and AIC values, the better the model.

Download table as:  ASCIITypeset image

The recursive least-squares periodogram with a circular orbit (RLSCP) of the astrometric data (see Figure 3, top panel) shows a very strong signal that extents beyond the extend of the plot. A blind search for the orbital period of the signal indicates that the orbital period of this astrometric signal could not be constrained. This suggests that the signal may be due to a companion with an orbital period much larger than the time span of the observations (≫16 yr). As we will see below (see Sections 4.4 and 4.5), this long-term signal is due to the gravitational interaction of GJ 896A with its very-low-mass stellar companion GJ 896B. As the orbital period of the binary system (>200 yrs) is much larger than the time span of the VLBA observations (∼16 yrs), the variation trend in the residuals can be taken into account by using accelerations terms in the astrometric fit.

Zoom In Zoom Out Reset image size

The astrometric data were then fitted with the least-squares and AGA algorithms, including acceleration terms, which take into account the astrometric signature due to the low-mass stellar companion (single-source solution). The results of this single-source solution are shown in column 2 of Table 2 and Figure 4(a). The astrometric solution shows that GJ 896A has proper motions of μ_α = 574.171 ± 0.014 mas yr⁻¹ and μ_δ =−60.333 ± 0.014 mas yr⁻¹. The fitted parallax is Π =160.027 ± 0.094 mas, which corresponds to a distance of 6.2489 ± 0.0037 pc. These values are similar to those obtained from the astrometric fit where no acceleration terms were taken into account. The fitted acceleration terms are relatively large (a_α = 0.8893 ± 0.0034 mas yr⁻² and a_δ = 0.1402 ± 0.0035 mas yr⁻²). We obtain now a better fit to the astrometric data, with smaller residuals (rms ∼0.24 mas in both R.A. and decl.; see column 2 of Table 2). The total residuals (rms ∼0.34) are a factor of 25 smaller than those obtained with the astrometric fit without acceleration terms (rms ∼8.64). However, the residuals are large compared to the mean noise present in the data (rms ∼0.11 and 0.10 mas for R.A. and decl., respectively) and the astrometric precision expected with the VLBA (<70 μas). Below we investigate the possibility that these residuals are due to a planetary companion orbiting around GJ 896A.

Zoom In Zoom Out Reset image size

We carried out a RLSCP of the astrometric data, including accelerations terms. The RLSCP shows now at least two significant peaks (see Figure 3, middle panel), the most prominent one with a period of about 280 days, and the weaker signal with a period of about 237 days. The stronger signal appears to be well constrained, with a false-alarm probability (FAP) of about 0.90%. This low FAP suggests that the main signal in the periodogram is real and due to the presence of a companion in a compact orbit. We then used both the least-squares and AGA algorithms to fit the astrometric observations of GJ 896A, now including acceleration terms and a possible single companion in a compact orbit (single-companion solution). The solution of the fit is shown in column 3 of Table 2 and in Figures 4(b) and 5. The fit of the astrometric data clearly improves when including a companion, as seen by the ${\chi }_{\mathrm{red}}^{2}$. The ${\chi }_{\mathrm{red}}^{2}$ is now a factor 2.4 smaller than that obtained with the single-source solution. Table 2 and Figure 4 show that the residuals of the single-companion solution (rms ∼0.14) are a factor of 1.7 smaller than in the case of the single-source solution (rms ∼0.24).

Column 3 in Table 2 summarizes the best fit of the VLBA astrometric data, including this new companion in a compact orbit. The single-companion solution shows that GJ 896A has proper motions of μ_α = 574.142 ± 0.019 mas yr⁻¹ and μ_δ = −60.354 ± 0.020 mas yr⁻¹. The fitted parallax is Π =159.83 ± 0.13 mas, which corresponds to a distance of 6.2565 ± 0.0051 pc. The fitted acceleration terms are a_α = 0.8871 ± 0.0047 mas yr⁻² and a_δ = 0.1400 ± 0.0049 mas yr⁻². These values are similar to those obtained from the single-source solution plus acceleration.

Figure 5 shows the orbital motion of the star GJ 896A due to the gravitational pull of the companion. The orbit of the main star around the barycenter has an orbital period P = 281.56 ± 1.67 days, an eccentricity e = 0.30 ± 0.11, a longitude of the periastron ω = 3440 ± 131, a position angle of the line of nodes Ω = 477 ± 121, a semimajor axis a_A = 0.52 ± 0.11 mas, and an inclination angle i = 660 ± 150, which indicates that the orbit is prograde (i < 90°). The orbit of the star around the barycenter due to the companion is well constrained, which is consistent with the narrow signal observed in the periodogram (see Figure 3, middle panel). With this astrometric fit alone we cannot estimate the mass of the companion. However, we used the estimated mass of the planetary system GJ 896A that we obtain below, using the full combined astrometric fit of this binary system (see Section 4.5). Column 3 in Table 2 summarizes the parameters of the new companion, hereafter GJ 896Ab. We find that its mass is 0.00223 ± 0.00047 M_⊙, which is consistent with a planetary companion with a mass of 2.35 ± 0.49 M_J . The orbit of the planetary companion around the barycenter has a semimajor axis a_b = 0.6352 ± 0.0018 au (or 101.53 ± 0.42 mas).

Zoom In Zoom Out Reset image size

There is an ambiguity with the position angle of the line of nodes between Ω and Ω + 180°. This ambiguity can be solved by radial velocity (RV) observations. However, as this planetary companion has not been detected with RV observations, we will leave the fitted Ω angle in Table 2. We further discuss this ambiguity below (see Section 5.5).

The orbit of the planetary companion is well constrained. However, the rms of the residuals (∼0.14 mas) is still large compared to the mean noise present in the data and the astrometric precision expected with the VLBA. To investigate the possibility of a second possible planetary companion, we obtained a new RLSCP, including acceleration terms and a signal with a period of 281.6 days. The new RLSCP does not show a significant signal (see Figure 3, bottom panel). We also carried out a blind search for a second planetary companion using the least-squares and AGA algorithms. We did not find a possible candidate. We will need further observations to investigate whether or not this source has more planetary companions.

GJ 896AB is a visual binary that contains two M dwarf stars that has been observed in the optical and near-infrared for more than 80 yr. The Washington Double Star (WDS; Mason et al. 2001) catalog provides separation and position angle measurements of this binary system spanning these decades and go back as far as 1941. As the stellar companion GJ 896B was detected in two of our recent VLBA observations of this binary system, we include the separation and the position angle of these two epochs in the relative astrometric fit. We performed an astrometric fit to the relative astrometric measurements with both, the least-squares and the AGA algorithms (binary-system solution). As the WDS catalog does not provide estimated error bars for most of the observed epochs, we applied a null weight to the astrometric fit, which translates into a uniform weight for all the observed epochs.

The results of the binary-system solution are presented in column 1 of Table 3 and Figure 6. We find that the relative astrometric fit of this binary system is well constrained when assuming a circular orbit. When the eccentricity of the binary system is taken into account, the solution does not converge to a stable solution. The orbit of the low-mass stellar companion GJ 896B around the main star GJ 896A has an orbital period P_AB = 96606.13 ± 2.24 days (264.5 yr), a position angle of the line of nodes Ω_AB = 797083 ± 00054, a semimajor axis a_AB = 5378.79 ± 0.51 mas, and an inclination angle i_AB =1306592 ± 00087, which indicates that the orbit of the stellar companion is retrograde (i > 90°). Using the distance to the source that we obtain from the full combined astrometric fit (see Section 4.5), we obtain that the mass and the semimajor axis of this binary system are m_AB = 0.544304 M_⊙ and a_AB = 33.64265 au, respectively.

Zoom In Zoom Out Reset image size

Table 3. Relative and Combined Astrometry Fits ^a

Parameter	Relative b	Combined c	Full Combined d
	(1)	(2)	(3)
μα (mas yr−1)	...	571.515 ± 0.019	571.467 ± 0.023
μδ (mas yr−1)	. ...	−37.750 ± 0.019	−37.715 ± 0.023
Π (mas)	159.88 (fixed)	159.98 ± 0.14	159.88 ± 0.17
		Binary System
PAB (days)	96606.13 ± 2.24	83665.80 ± 1.64	83664.63 ± 1.98
T0AB (Julian day) e	2,504,877.56 ± 1.48	2,401,894.57 ± 1.00	2,401,891.34 ± 1.19
eAB	0.0 (fixed)	0.108047 ± 0.000044	0.108047 ± 0.000053
ωAB (A) (deg)	...	127.1531 ± 0.0037	127.1416 ± 0.0045
ωAB (B) (deg)	0.0 (fixed)	307.1531 ± 0.0037	307.1416 ± 0.0045
ΩAB (deg)	259.7083 ± 0.0054	255.0891 ± 0.0028	255.0919 ± 0.0034
aAB (A) (mas)	...	1381.015 ± 0.069	1385.328 ± 0.083
aAB (B) (mas)	...	3676.95 ± 0.35	3672.64 ± 0.42
aAB (mas)	5378.79 ± 0.51	5057.96 ± 0.36	5057.97 ± 0.43
iAB (deg)	130.6592 ± 0.0087	130.0664 ± 0.0085	130.065 ± 0.010
QB/A	...	0.375588 ± 0.000030	0.377202 ± 0.000037
		Companion
PAb (days)	...	...	284.39 ± 1.47
T0Ab (Julian day)	...	...	2,455,702.65 ± 17.26
eAb	...	...	0.35 ± 0.19
ωA (deg)	...	...	353.11 ± 11.81
ΩAb (deg)	...	...	45.62 ± 9.60
aA (mas)	...	...	0.51 ± 0.15
ab (mas)	...	...	102.27 ± 0.15
iAb (deg)	...	...	69.20 ± 25.61
		Other Parameters
D (pc)	...	6.2506 ± 0.0055	6.2545 ± 0.0066
mAB (M⊙)	0.544304 ± 0.000023	0.602253 ± 0.000020	0.603410 ± 0.000025
mAB (A)(M⊙) f	...	0.43782 ± 0.00054	0.43814 ± 0.00065
mAB (B)(M⊙)	...	0.16444 ± 0.00020	0.16527 ± 0.00025
aAB (au)	33.64265 ± 0.00052	31.615 ± 0.028	31.635 ± 0.033
aAB (A)(au)	...	8.6322 ± 0.0075	8.6646 ± 0.0091
aAB (B)(au)	...	22.983 ± 0.020	22.971 ± 0.024
mA (M⊙) g	...	...	0.43599 ± 0.00092
mb (M⊙)	...	...	0.00216 ± 0.00064
mb (J)	...	...	2.26 ± 0.57
aAb (au)	...	...	0.64282 ± 0.00068
aA (au)	...	...	0.00317 ± 0.00093
ab (au)	...	...	0.63965 ± 0.00067
ΔαAB (mas) h	123.11	89.60	89.60
ΔδAB (mas) h	89.14	74.28	74.28
ΔαA (mas) h	...	0.27	0.21
ΔδA (mas) h	...	0.47	0.45
ΔαB (mas) h	...	0.20	0.11
ΔδB (mas) h	...	0.12	0.042
χ2, ${\chi }_{\mathrm{red}}^{2}$ i	980846.81, 6956.36	977729.63, 5785.38	977566.99, 6034.36

Notes.

^a The parameters presented here were obtained with AGA. Very similar results were obtained with the least-squares fitting method. The subindexes A, B, and AB correspond to the main star (GJ 896A), the low-mass stellar companion (GJ 896B), and the binary system (GJ 896AB), respectively. ^b Relative astrometric fit of the secondary star GJ 896B around the main star GJ 896A. In this case the parallax is fixed using the solution of the full combined astrometric fit (column 3). ^c The combined astrometric fit is obtained by fitting simultaneously the relative astrometry of the binary system and the absolute astrometry of the main star GJ 896A and the secondary star GJ 896B (see text). All free parameters are fitted simultaneously. ^d The full combined astrometric fit is obtained by fitting simultaneously the relative astrometry of the binary system and the absolute astrometry of the two stars, including a companion associated to GJ 896A (see text). All free parameters are fitted simultaneously. ^e Time of periastron passage. In the case of the relative astrometry (column 1), T_0AB corresponds to the time of passage through the ascending line of nodes of the orbit of the secondary star around the main star. In the case of the combined astrometry (columns 2 and 3), T_0AB corresponds to the time of periastron passage of the primary star. ^f Dynamical mass of the star GJ 896A including any possible companions. ^g Dynamical mass of the star GJ 896A removing the mass of the planetary companion. ^h rms dispersion of the residuals. The first two terms correspond to the residuals of the binary system part of the fit, the next two terms correspond to the residuals of the absolute astrometry part of the fit of the main star, and the last two terms correspond to the residuals of the secondary star. ⁱ χ² and reduced χ² of the astrometric fit. In all three cases the residuals of the relative astrometry dominate over the residuals of the fit.

Download table as:  ASCIITypeset image

As it was mentioned before, there is an ambiguity with the position angle of the line of nodes between Ω and Ω + 180°. In this case, the observed RV of the two M dwarfs can be used to find the correct angle. From the GAIA catalog (Lindegren et al.2018), the observed mean radial velocity of GJ 896A is −0.02 ± 0.31 km s⁻¹. The GAIA catalog does not contain the radial velocity of GJ 896B. However, recent RV observations of this source indicate that its RV is about 3.34 km s⁻¹ (Morin et al. 2008). This means that GJ 896B is redshifted with respect to GJ 896A. Thus, the position angle of the line of nodes is consistent with these RVs is Ω_AB = 2597083 ± 00054. This correction to the position angle of the line of nodes is applied to all the astrometric fits presented in Table 3.

The optical/infrared relative astrometry of the binary system GJ 896AB (Mason et al. 2001), obtained in a time interval of about 80 yr, was also combined with the absolute radio astrometry to simultaneously fit: the orbital motion of the two stars in GJ 896AB; the orbital motion of the planetary companion around GJ 896A (see Section 4.5); and the parallax and proper motion of the entire system. However, as most of the relative astrometric observations do not have estimated errors, for the relative astrometry part we carried out a uniform weighted fit (without errors) of the observed data (73 epochs, including the two relative positions of GJ 896B obtained with the VLBA), while for the absolute astrometry of the VLBA observations of the main star GJ 896A (16 epochs) and the secondary star GJ 896B (2 epochs), we applied a weighted fit using the estimated error of the observed position at each epoch. The results presented here and in Section 4.5 were obtained by fitting simultaneously the absolute astrometric observations of both stars, and the relative astrometry of this binary system. The results are limited by: (a) the lack of error bars of the relative astrometric data, (b) the fact that the relative astrometric data cover about 35% of the orbital period of the binary system (∼229.06 yrs; see Section 4.5), and (c) the fact that the absolute astrometry covers only a time span of about 16.56 yr. Thus, all errors reported below are statistical only.

Fitting the proper motions and parallax of a binary system is complex due to the orbital motions of each component around the center of mass, especially because both stars have different masses (m_B /m_A < 1). The proper motion and parallax of a binary system can be obtained with high precision when the orbital period of the system is small compared with the time span of the astrometric observations. In contrast, when the time span of the astrometric observations is small compared with the orbital period of the binary system, it is difficult to separate the orbital motion of each component and the proper motion of the system as both spatial movements are blended. The best way to separate both movements is to simultaneously fit the proper motions, the parallax, and the orbital motion of the binary system. Our combined astrometric fit includes all these components, and thus gives an accurate estimate of the proper motions and the parallax of this binary system.

We carried out a combined fit of the relative and absolute astrometric data and found that the astrometric fit improves substantially. The solution of the combined astrometric fit is similar to that obtained from the relative astrometric fit. However, using the combined astrometric fit, we are able to obtain the parallax and the proper motions of the binary system, as well as improved orbits of the two stars around the center of mass and the masses of the system and the individual stars. The results of this combined fit are shown in Figures 7(a) and 8(a) and summarized in column 2 of Table 3. Although, Figure 8(a) shows the full combined astrometric solution (including a planetary companion; see Section 4.5), the solution for the relative astrometric part of the binary system is basically the same as that obtained from the combined astrometric solution (see columns 2 and 3 of Table 3). Thus, we present only one figure showing both fits. In contrast to the case of the relative astrometric fit (see Section 4.3), we find here that the combined astrometric fit is well constrained when considering that the orbit of the binary system could be eccentric (e ≥ 0). The proper motions of the binary system are μ_α = 571.515 ± 0.019 mas yr⁻¹ and μ_δ = −37.750 ± 0.019 mas yr⁻¹. The parallax of the binary system is Π = 159.98 ± 0.14 mas, which corresponds to a distance of 6.2506 ± 0.0055 pc. The orbit of the binary system has an orbital period P_AB = 83665.80 ± 1.64 days (229.06 yr), a position angle of the line of nodes Ω_AB = 2550891 ± 00028, a semimajor axis of the binary system a_AB = 5057.96 ± 0.36 mas, a semimajor axis of the main source a_AB (A) = 1381.015 ± 0.069 mas, and an inclination angle i_AB = 1300664 ± 00085. The estimated masses of the binary system and the two components are m_AB = 0.602253 ± 0.000020 M_⊙, m_AB (A) = 0.43782 ±0.00054 M_⊙, and m_AB (B) = 0.16444 ± 0.00020 M_⊙, respectively. The semimajor axis of the binary system and the two stellar components are a_AB = 31.615 ± 0.028, a_AB (A) = 8.6322 ±0.0075, and a_AB (B) = 22.983 ± 0.020 au, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Column 2 of Table 3 indicates that there is a relatively small improvement of about 17% in the ${\chi }_{\mathrm{red}}^{2}$. This is because the residuals of the relative data of the binary system are much larger that those of the absolute astrometry of the main star, and this dominates the residuals. However, Table 3 and Figure 7 show that the residuals of the relative astrometric data of the binary system from the combined astrometric solution (rms ∼116.38 mas) are a factor of 1.3 smaller than in the case of the relative astrometry solution (rms ∼152.0 mas). Thus, the combined astrometric solution is an improvement to the solution obtained from the pure relative astrometric fit. In addition, the residuals of the absolute astrometric data of the main source GJ 896A (rms ∼0.27 and 0.47 mas for R.A. and decl.) are large compared to the mean noise present in the data.

Comparing the residuals that were obtained from the single-source fit (see Section 4.1) of the astrometric observations of the main star GJ 896A, including acceleration terms, and the residuals of the combined astrometric fit (see column 2 of Table 2 and column 2 of Table 3, and Figures 4 and 7), we find that the rms of the residuals of the absolute astrometry are similar. In addition, we find that the temporal distribution of the residuals are also very similar. However, the residuals of the last three epochs from the combined astrometric fit are considerable larger that those obtained from the single-source plus acceleration fit. We do not find an explanation for this discrepancy. However, we notice that two out of these three epochs are the ones where the secondary star was detected, and the coordinates of this star were used in the fit of the relative astrometry of the binary system and in the absolute astrometry fit of the secondary star GJ 896B. As the secondary star seems to be a close binary (see Sections 1 and 5.1) and it was only detected in two epochs, this probably affects the astrometric fit (see Figure 7). Further observations will be needed to improve the astrometric fit of this binary system.

As an independent check of the derived parameters, we performed a combined fit of the astrometric data with the lmfit package (Newville et al. 2020), which uses a nonlinear least-squares minimization algorithm to search for the best fit of the observed data (see Appendix B for further details). We find that the astrometric fit is well constrained, and that the solution and the residuals of the fit are in good agreement with those derived with the combined astrometric fit (see Section 4.4).

These results suggest that the main star may have at least one companion. Below we investigate further the possibility that these residuals are due to a planetary companion orbiting around the main star GJ 896A (as in the case of the single-companion fit; see Section 4.2).

We carried out a combined astrometric fit of the relative orbit of the binary and the absolute astrometric data of the main and the secondary stars, now including the orbital motion of a possible companion in a compact orbit (full combined astrometric fit), and found that the astrometric fit improves substantially. The results of this full combined astrometric fit are presented in Figures 7(b), 8, and 9. We find that the full combined astrometric fit is well constrained. Column 3 of Table 3 summarizes the best-fit parameters of the full combined astrometric fit.

By combining relative and absolute astrometric data of the binary system, we are able to obtain the full dynamical motion of the system, including the close companion. From the best fit of the full combined astrometric data, we obtained proper motions μ_α = 571.467 ± 0.023 mas yr⁻¹ and μ_δ = −37.715 ± 0.023 mas yr⁻¹, and a parallax Π = 159.88 ± 0.17 mas for the binary system (see column 3 of Table 3). The estimated proper motions and parallax are very similar to those obtained from the Combined astrometric fit. However, as expected, the proper motions differ from those obtained by GAIA for both stars GJ 896A and GJ 896B, which where obtained from independent linear fits of both stars (Lindegren et al. 2018). The parallax, which corresponds to a distance of d = 6.2545 ± 0.0066 pc, is an improvement to that obtained by GAIA for both stars (Π_A = 159.663 ± 0.034 mas, and Π_B = 159.908 ± 0.051 mas). The parallax estimated by us is for the binary system, while GAIA's parallaxes are for the individual stars.

The part of the solution that corresponds to the astrometric fit of the binary system is very similar to that obtained from the combined astrometric fit (see columns 2 and 3 of Table 3). The fit of the binary system does not improve in this case. However, the best full combined fit of the astrometric data indicates that the M3.5 dwarf GJ 896A has at least one planetary companion, GJ 896Ab. The orbit of the star around the barycenter, due to this planetary companion, has an orbital period P_Ab =284.39 ± 1.47 days, an eccentricity e_Ab = 0.35 ± 0.19, a longitude of the periastron ω_A = 353.11° ± 11.81°, a position angle of the line of nodes Ω_Ab = 45.62° ± 9.60°, a semimajor axis a_A = 0.51 ± 0.15 mas, and an inclination angle i_Ab = 69.20° ± 25.61°, which indicates that the orbit is well constrained and prograde (i_Ab < 90°). However, the rms of the residuals (∼0.21 and ∼0.45 mas in R.A. and decl., respectively) are still large compared to the mean noise present in the data (∼0.15 mas) and the astrometric precision expected with the VLBA (<70 μas). This may explain the large errors of the orbital parameters. The large residuals may also indicate the presence of a second planetary companion. We carried out a blind search for a possible astrometric signal due a second planetary companion using the least-squares and the AGA algorithms; however, we did not find a candidate.

Using the fitted solution, we obtain the total mass of the system (m_AB = 0.603410 ± 0.000025 M_⊙), the masses of the two stars (m_AB (A) = 0.43814 ± 0.00065 M_⊙ and m_AB (B) = 0.16527 ± 0.00025 M_⊙). In addition, we find that the mass of the star GJ 896A is m_A = 0.43599 ± 0.00092 M_⊙, and that the mass of the planetary companion is m_b = 0.00216 ± 0.00064 M_⊙, which is consistent with being planetary in origin with an estimated mass of 2.26 ± 0.57 M_J . The semimajor axis of the orbit of the binary system is a_AB = 31.635 ± 0.033 au. The semimajor axis of the two stars around the center of mass are a_AB (A) = 8.6646 ± 0.0091 au and a_AB (B) = 22.971 ± 0.024 au, and the semimajor axis of the orbit of the planetary companion is a_b = 0.63965 ± 0.00067 au (or 102.27 ± 0.15 mas).

This is the first time that a planetary companion of one of the stars in a binary system has been found using the astrometry technique. The solution of the full combined astrometric fit of the planetary companion is similar to that obtained from the single-companion solution (see Section 4.2 and Tables 2 and 3). The fact that the astrometric signal appears in the periodogram of the absolute astrometric observations of GJ 896A, as well as in both the single-companion astrometric fit and the full combined astrometric fit obtained with two different algorithms (least-squares and AGA), supports the detection of an astrometric signal due to a companion. Furthermore, the similar astrometric fit solution obtained from the single-companion astrometric fit and the full combined astrometric fit further supports the planetary origin of the astrometric signal.

As mentioned before, the estimation of the proper motions of a binary system is complex due to the orbital motions of each component around the center of mass, in particular when both stars have different masses (m_B /m_A < 1) and the time span of the observations cover only a fraction of the orbital period of the binary system. In such a case, it is difficult to separate the orbital motion of each component and the proper motion of the system, both movements are blended. The best way to separate both movements is to simultaneously fit the proper motions, the parallax, and the orbital motion of the binary system. The full combined astrometric fit that we carried out (see Section 4.5) includes all these components and thus gives good estimates of the proper motion and the orbital motion of the binary system (see Table 3).

The full combined astrometric solution (see column 3 of Table 3) shows that the orbital motions of the two stars around the center of mass and of the planetary companion GJ 896Ab around the main star GJ 896A are well constrained. In addition, we found a similar solution for the orbital motion of the planetary companion using the full combined astrometric fit (see Section 4.5 and column 3 of Table 3) and the single-companion astrometric fit (see Section 4.2 and column 3 of Table 2). In the former case, the astrometric fit includes intrinsically the orbital motion of GJ 896A around the center of mass of the binary system, while in the latter case the included acceleration terms take into account this orbital motion.

Tables 2 and 3 show that the proper motions of the system obtained with the single-source, single-companion, and full combined astrometric fits differ significantly. Here we discuss the origin of this difference. The single-source astrometric solution gives μ_α = 576.707 mas yr⁻¹ and μ_δ = −59.973 mas yr⁻¹ (see Section 4.1 and column 1 of Table 2), where we include only the proper motions and the parallax to the absolute astrometric fit of the main star GJ 896A. Including acceleration terms (single-source plus acceleration astrometric solution), to take into account the orbital motion of the binary system, the estimated proper motions of the system are μ_α = 574.171 mas yr⁻¹ and μ_δ = −60.333 mas yr⁻¹ (see Section 4.2 and column 2 of Table 2). On the other hand, from the full combined astrometric solution, the proper motions of the system are μ_α = 571.467 mas yr⁻¹ and μ_δ = −37.715 mas yr⁻¹ (see Section 4.5 and column 3 of Table 3). The absolute astrometric solutions presented in Table 2 show that by including the acceleration terms in the fit we improve the estimated proper motions of the system and that the inclusion of an astrometric signal due to a companion in the fit does not affect the estimation of the proper motions and the estimated acceleration terms. Comparing the single-source and the full combined astrometric solutions, we find a significant difference in the estimated proper motions, particularly in the north–south direction (Δμ_α = 5.24 mas yr⁻¹ and Δμ_δ = −22.26 mas yr⁻¹). This is consistent with the fact that the VLBA astrometric observations of GJ 896A were obtained when the source was crossing the ascending node of the binary orbit and mainly in the north–south direction (see discussion below and Figures 8 and 9). These results indicate that the orbital motion of GJ 896A around the center of mass of the system blends with the proper motions of the binary system and in second order with the acceleration terms included in the absolute astrometric fit. Thus, we conclude that the best estimates of the proper motions of this binary system are those obtained with the full combined astrometric fit (see Table 3).

Zoom In Zoom Out Reset image size

We can also obtain an estimation of the velocity on the plane of the sky of the center of mass μ(CM) of this binary system, as follows:

$$\begin{eqnarray}&&{\mu }_{\alpha }(A)={\mu }_{\alpha }({CM})+{v}_{\alpha }(A),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\mu }_{\delta }(A)={\mu }_{\delta }({CM})+{v}_{\delta }(A),\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{v}_{\alpha }(A)=(2\pi {a}_{A}/{P}_{{AB}})\times {\sin }({\theta }_{{rot}})\times {\cos }(\phi ),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{v}_{\delta }(A)=(2\pi {a}_{A}/{P}_{{AB}})\times {\cos }({\theta }_{{rot}})\times {\cos }(\phi ),\end{eqnarray} \tag{ 4 }$$

where μ_α (A) and μ_δ (A) are the observed proper motion of the main source GJ 896A (single-source solution; see column 1 of Table 2), μ_α (CM) and μ_δ (CM) are the proper motions of the center of mass of the binary system, v_α (A) and v_δ (A) are the projected rotational velocity of GJ 896A on its orbital motion around the center of mass, a_A is semimajor axis of the orbital motion of GJ 896A around the center of mass, P_AB is the orbital period of the binary system, θ_rot ≈Ω_AB − 90° is the position angle of the projected rotation velocity (tangential velocity) of GJ 896A around the center of mass of the system (see Figures 8(b) and 9), ϕ = 180° − i_AB is the complementary angle of the orbital plane of the binary system, and i_AB is the fitted inclination angle of the orbital plane of the binary system (see column 3 of Table 3). Using the fitted parameters presented in column 1 of Table 2 and column 3 of Table 3, μ_α (A) = 576.707 mas yr⁻¹ and μ_δ (A) = −59.973 mas yr⁻¹ (μ_A = 579.8170 mas yr⁻¹ with PA = 9594), we obtain v_α (A) = 6.2925 mas yr⁻¹ and v_δ (A) = −23.6355 mas yr⁻¹ (v_rot =24.4588 mas yr⁻¹ with PA = 16509), and thus μ_α (CM) = 570.4145 mas yr⁻¹ and μ_δ (CM) = −36.3375 mas yr⁻¹ (μ(CM) = 571.5707 mas yr⁻¹ with PA = 9365). We find that the estimated proper motions of the center of mass are very similar to those obtained with the full combined astrometric fit μ_α =571.467 mas yr⁻¹ and μ_δ = −37.715 mas yr⁻¹ (μ = 572.72 mas yr⁻¹ with PA = 93.79°; see Table 3). The difference in the estimated proper motions of the center of mass is about 1 mas yr⁻¹. This further confirms that the full combined astrometric fit gives the best estimate for the proper motions of the center of mass of a binary system.

We can estimate the acceleration, acc(A), of GJ 896A due to the gravitational pull of the low-mass stellar companion GJ 896B, as follows:

$$\begin{eqnarray}&&\mathrm{acc}{(A)=(2\pi /{P}_{{AB}})}^{2}{a}_{A},\end{eqnarray} \tag{ 5 }$$

where acc(A) is the mean acceleration of GJ 896A, P_AB is the orbital period of the binary system, and a_A is the semimajor axis of the orbital motion of GJ 896A around the center of mass. Using the full combined solution (see column 3 of Table 3), we obtain that acc(A) ≃1.04234 mas yr⁻². In the case of the single-companion solution (see column 3 of Table 2), the fitted acceleration terms are a_α = 0.8873 mas yr⁻² and a_δ = 0.1402 mas yr⁻² , and thus the acceleration of the main star in the plane of the sky is acc(A) = $\sqrt{{a}_{\alpha }^{2}+{a}_{\delta }^{2}}$= 0.89831 mas yr⁻². Thus, the estimated acceleration of GJ 896A due to GJ 896B is consistent with the acceleration found in the single-companion astrometric fit.

The acceleration acc(A) obtained from the single-source astrometric fit of the VLBA data allows us to place a mass estimate for the stellar companion, m_B , using

$$\begin{eqnarray}&&\left(\displaystyle \frac{{m}_{B}}{{M}_{\odot }}\right)=0.02533\left(\displaystyle \frac{{acc}(A)}{{{auyr}}^{-2}}\right){\left(\displaystyle \frac{{a}_{{AB}}}{{au}}\right)}^{2},\end{eqnarray} \tag{ 6 }$$

or

$$\begin{eqnarray}&&\left(\displaystyle \frac{{m}_{B}}{{M}_{\odot }}\right)=0.29368{\left[\left(\displaystyle \frac{{acc}(A)}{{{auyr}}^{-2}}\right){\left(\displaystyle \frac{{a}_{A}}{{au}}\right)}^{2}{\left(\displaystyle \frac{{M}_{{AB}}}{{M}_{\odot }}\right)}^{2}\right]}^{1/3},\end{eqnarray} \tag{ 7 }$$

where acc(A) = $\sqrt{{a}_{\alpha }^{2}+{a}_{\delta }^{2}}$ is the estimated acceleration needed to fit the absolute astrometric data (see column 3 of Table 2), a_AB = a_A + a_B is the semimajor axis of the orbit of GJ 896B around GJ 896A, and m_AB is the total mass of the binary system (see column 3 of Table 3). From Equation (7) we find that the estimated mass of the stellar companion GJ 896B is m_B ≈0.15728 M ⊙ . This estimated mass is consistent with the mass of the stellar companion (m_B = 0.16527 M ⊙ ) obtained from the full combined astrometric fit (see column 3 of Table 3).

The solution of the full combined astrometric fit can be used to estimate an expected induced maximum RV of the star due to the gravitational pull of a companion as follows (e.g., Cantó et al. 2009; Curiel et al. 2020):

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

where G is the gravitational constant, and T, M_*, m_c , and e are the estimated orbital period, star and companion masses, and the eccentricity of the orbit of the companion. Using the full combined astrometric solution (see column 3 of Table 3), the maximum RV of GJ 896A induced by the planetary companion GJ 896Ab is K_A (b) ∼ 121 m s⁻¹, and the maximum RV induced by the stellar companion GJ 896B is K_A (B) ∼ 867 m s⁻¹. These RVs can in principle be observed with modern high-spectral resolution spectrographs.

The maximum radial velocity of both stars occurred in 2013 October, when GJ 896A and GJ 896B passed through the ascending and descending nodes, respectively, of their orbits around the center of mass of the binary system (see Figure 8). The RV signal due to GJ 896Ab can be measured with modern high-spectral resolution spectrographs on a time span shorter than one year. Recent radial velocity observations of GJ 896A show a radial velocity variability of ΔV ∼ 175 ± 37 m s⁻¹ (e.g., Gagne et al. 2016), which is consistent with the expected radial velocity of this source, induced by the planetary companion GJ 896Ab. However, these observations were taken in just seven epochs within a time span of about 4 yr. Further observations will be needed to confirm whether the RV variability observed on GJ 896A is due to this planetary companion.

The RV signal of GJ 896A (and GJ 896B) due to the stellar companion will be difficult to measure due to the binary's very long orbital period of 229.06 yr. It will be very difficult to separate this signal from the systemic velocity V₀ of the binary system. However, we can use the observed mean RV of GJ 896A (Lindegren et al. 2018), and the estimated maximum RV of this star due to its stellar companion to obtain a raw estimate of the systemic radial velocity of the binary system, as follows:

$$\begin{eqnarray}&&{V}_{\mathrm{obs}}(A)={K}_{A}(B)+{V}_{0},\end{eqnarray} \tag{ 9 }$$

where V_obs is the observed radial velocity, K_A (B) is the expected maximum RV of GJ 896A due to GJ 896B, and V₀ is the barycentric RV of the binary system. From the GAIA catalog, the observed radial velocity of GJ 896A is V_obs(A) = −0.02 ± 0.31 km s⁻¹ (Lindegren et al. 2018). As these RV observations of GJ 896A were take in a time span of about 2 yr, the RV signal from GJ 896Ab is averaged out in this mean radial velocity. The reference epoch of the GAIA DR2 observations is J2015.5, which is close to orbital ___location where the radial velocity of GJ 896A is maximum and negative (see Figure 8). Thus, K_A (B) = −867 m s⁻¹, and the resulting barycentric RV of the binary system is V₀ ∼847 m s⁻¹.

A similar calculation can be carried out for GJ 896B. In this case, the observed RV is V_obs(B) = 3340 m s⁻¹ (Morin et al. 2008), which was obtained also close to the maximum RV of this star. The reference epoch of the RV observations is 2006.0, which is close to orbital ___location where the radial velocity of GJ 896B is maximum and positive (see Figure 8). We estimate that the maximum RV of GJ 896B induced by the main star GJ 896A is K_B (A) ∼ 2299 m s⁻¹. Thus, in this case we obtain a barycenter RV of V₀ ∼ 1041 m s⁻¹. This barycentric RV is different to that obtained from the RV observation of GJ 896A (847 m s⁻¹). There is a difference of about 194 m s⁻¹. The RV measurements of GJ 896B were obtained in a time period of only a few days, and therefore the observed RV would include the contribution due to possible companions associated to this M4.5 dwarf. This suggests that GJ 896B may have at least one companion. Such a radial velocity signature should be easily detected using modern high-spectral-resolution spectrographs. Recent observations also suggest that the M4.5 dwarf GJ 896B may be an unresolved binary system (Winters et al. 2021). Further observations will show whether this very-low-mass M dwarf star has a companion.

A simple estimate of the habitable zone (HZ) can be obtained as follows. The inner and outer boundaries of the HZ around a star depends mainly on the stellar luminosity. Thus, combining the distance dependence of the HZ as a function of the luminosity of the star and the mass–luminosity relation, the inner a_i and outer a_o radius of the HZ are:

$$\begin{eqnarray}&&{a}_{i}=\sqrt{\left(\displaystyle \frac{{L}_{\star }}{{L}_{\odot }}\right)\displaystyle \frac{1}{1.1}},\,\,\,\,\,{a}_{o}=\sqrt{\left(\displaystyle \frac{{L}_{\star }}{{L}_{\odot }}\right)\displaystyle \frac{1}{0.53}}\end{eqnarray} \tag{ 10 }$$

(e.g., Selsis et al. 2007; Kopparapu et al. 2013), where L_⋆ is the luminosity of the star in solar luminosities. In the case of M dwarfs with masses below 0.43 M_⊙, the mass–stellar luminosity relation can be approximated as:

$$\begin{eqnarray}&&\left(\displaystyle \frac{{L}_{\star }}{{L}_{\odot }}\right)=0.23\times {\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{2.3}\end{eqnarray} \tag{ 11 }$$

(e.g., Kuiper 1938; Duric 2012), where M_⋆ is the mass of the star in solar masses. Using the estimated mass of GJ 896A (m_A = 0.436 M ⊙; see Table 3), the limits of the habitable zone around this M3.5 dwarf are a_i ∼ 0.18 au and a_o ∼0.26 au. This is considerably smaller than the semimajor axis of the orbit of the planetary companion GJ 896Ab (a_b = 0.6397 au). As the planetary companion is in an eccentric orbit (e_b = 0.35), the minimum and maximum distances between the planet and the star are 0.42 and 0.86 au, respectively. Thus, the orbit of the planet is located outside the habitable zone of this M3.5 dwarf star.

We can also estimate if the planetary companion GJ 896Ab is located inside the snow line, a_snow. The snow line is located at an approximate distance of (e.g., Kennedy & Kenyon 2008)

$$\begin{eqnarray}&&{a}_{\mathrm{snow}}=2.7{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{2}.\end{eqnarray} \tag{ 12 }$$

Using the estimated mass of GJ 896A, the snow line in this planetary system is located at a radius of a_snow ∼ 0.51 au. Thus the estimated orbit of the planetary companion GJ 896Ab is located outside the snow line. However, given the eccentricity of the orbit, the planet moves around the estimated snow line distance, but most of the time is located outside the snow line. Recent results suggest that stars with a mass of about 0.4 M_⊙, such as GJ 896A, have a 1% probability of having at least one Jovian planet (e.g., Kennedy & Kenyon 2008). Even when the probability of having a Jovian planet is very low, the results presented here show that the main star GJ 896A in the M dwarf binary system GJ 896AB has at least one Jupiter-like planet.

The radio continuum flux density of GJ 896A is clearly variable in time. Figure 10 shows that the source goes through a large variability in short periods of time. The flux variability does not seem to have a clear pattern. The flux density of the source has changed nearly 2 orders of magnitude during the last 16 yr of the radio observations, with variations at scales of months and at scales of a few years. However, further observations will be required to find if the variability has a defined temporal period.

Zoom In Zoom Out Reset image size

Characterizing the full 3D orbital architecture of binary systems containing a planetary companion can aid to investigate the importance of the star–star and star–planet mutual interaction. Combining relative and absolute astrometric data of the binary system, we are able to obtain the 3D orbital architecture of the system, including the planetary companion (see Figure 9). As the full combined astrometric fit (relative plus absolute astrometry) provides the inclination angle and the position angle of the line of nodes of the orbital planes of both the planetary companion GJ 896b and the binary system GJ 896AB, we can directly measure the mutual inclination angle of this system. A first approach would be to estimate the inclination difference (Δi) between the inclination angles of the orbital planes of the planet and the binary system, assuming that their position angle of the line of nodes is equal to 0°. The inclination angles are measured with respect to the plane of the sky (such that i = 0° corresponds to a face-on orbit, and it increases from the east toward the observer). Using the fitted inclination angles we obtain that Δi = 60.9° ± 22.5° (see Table 3 and Figure 9), which is a very large difference. The mutual inclination angle Φ between the two orbital planes can be determined through (e.g., Kopal 1959; Muterspaugh et al. 2006):

$$\begin{eqnarray}&&\cos \,{\rm{\Phi }}={\cos }\,{i}_{{Ab}}\,\cos \,{i}_{{AB}}+\sin \,{i}_{{Ab}}\,\sin \,{i}_{{AB}}\,\cos ({{\rm{\Omega }}}_{{Ab}}-{{\rm{\Omega }}}_{{AB}}),\end{eqnarray} \tag{ 13 }$$

where i_Ab and i_AB are the orbital inclination angles, and Ω_Ab and Ω_AB are the position angles of the line of nodes. The position angle of the line of nodes is measured anticlockwise from the north toward the ascending node. Table 3 contains the inclination angle and the position angle of the line of nodes for the planet GJ 896Ab and the binary system GJ 896AB. There is an ambiguity in the position angles of the line of nodes (Ω or Ω + 180°, where Ω is the fitted angle), which can be disentangled by RV measurements. For the orbital motion of the binary system, recent RV measurements of both stellar components show that GJ 896B is receding (V_obs(B) = +3.34 km s⁻¹; Morin et al. 2008) and GJ 896A is approaching us (V_obs(A) = −0.02 ± 0.31 km s⁻¹; Lindegren et al. 2018); thus the correct position angle of the line of nodes is Ω_AB = 255.09° (Ω_AB + 180°; see Table 3). However, the planetary companion has no measured RV, so the fitted value of the position angle of the line of nodes of the orbital plane of GJ 896Ab could be either Ω_Ab = 456 or 2256; the former is the fitted angle. From these position angles we calculate a mutual inclination angle between the two orbital planes of Φ = 148° for the fitted Ω_Ab , and Φ = 67° for the second possibility (Ω_Ab + 180°). Taking into account the long-term stability of the system (see Section 5.6), we found that the former solution (Φ = 148°) is stable in a very long period of time, while the latter solution (Φ = 67°) is unstable in a short period of time. This result indicates that there is a large mutual inclination angle of Φ = 148° between both orbital planes.

Recent observations suggest that the rotation axis of GJ 896A has an inclination angle of about 60° ± 20° with respect to the line of sight (Morin et al. 2008), and thus, an inclination angle of i_s ∼210° with respect to the plane of the sky (see Figure 9). On the other hand, the inclination angle of the rotation axis of the orbital motion of GJ 896Ab is i_p = 1592 ± 2561 (i_b + 90°). Thus these rotation planes have a difference in their inclination angles of Δ_s−p ∼51° (see Figure 9). This indicates that the orbital motion of the planetary companion and the rotation plane of the star are far from being coplanar. We can also compare the angle of the rotation axis of the host star GJ 896A and the inclination angle, i_bs , of the rotation axis of orbital motion of the binary system GJ 896AB. In this case the inclination angle is i_bs = 220065 ± 0010 (i_AB + 90°). Thus the rotation planes of the star GJ 896A and the binary system have a difference in their inclination angles of Δi_bs−s ∼10°. In addition, the rotation axis of GJ 896B, with respect to the line of sight, is also about 60° ± 20° (Morin et al. 2008), which is basically the same as that of GJ 896A, and thus the difference between the inclination angle of the rotation axis of GJ 896B, with respect of the plane of the sky, compared with the rotation axis of the binary system is the same as that found for the star GJ 896A (∼10°). Thus, this suggests that the rotation planes of both stars are nearly parallel to the orbital plane of the binary system (Δi ∼ 10°), while the orbital planes of the planet and the binary system have a mutual inclination angle of Φ ≃ 148°. However, having similar inclination angles does not necessarily imply alignment between the rotation axis of both stars as the position angles of the line of nodes are unknown. The rotation axis of the stars could be different, and thus they could have larger and different mutual inclination angles with respect to the orbital plane of the binary system.

We applied direct N-body integrations to the full combined orbital solution obtained from the Keplerian orbits of the binary system GJ 896AB and the planetary companion GJ 896Ab (see Table 3). We integrated the orbits for at least 100 Myr using the hybrid integrator in Mercury 6 (Chambers 1999), which uses a mixed-variable symplectic integrator (Wisdom & Holman 1991) with a time step approximately equal to a hundredth ( ≃ 1/100) of the Keplerian orbital period of the planetary companion. During close encounters, Mercury uses a Bulrich–Stoer integrator with an accuracy of 10⁻¹². We identify an unstable system if: (a) the two companions (the planetary companion and low-mass stellar companion) collide; (b) the planet is accreted onto the star (astrocentric distance less than 0.003 au); and (c) the planet is ejected from the system (astrocentric distance exceeds 200 au). The integration time of 100 Myr long exceeds the 10,000 binary periods that is considered as a stability criterion. The simulation using the fitted Keplerian orbital parameters proved to be stable for at least 100 Myr. However, the simulation of the alternative position angle of the line of nodes of the planetary orbit (Ω = Ω_Ab + 180°) turned out to be unstable in very short times, the planet collided with the host star GJ 896A after a few tens of thousand years. In this case, the eccentricity of the planetary orbit increases rapidly due to the interaction between the stellar companion and the planet. After a few tens of thousand years the eccentricity reaches an extreme value of 1, and thus the planet collides against the main star. This indicates that the fitted solution contains the correct angle of the line of nodes of the planetary orbit.

To complement our stability analysis, we also performed N-body long-term integrations using REBOUND (Rein & Liu 2012). We tested the combined orbital solutions using two different integrators: IAS15 (Rein & Spiegel 2015) and Mercurius (Rein et al. 2012). The first one is a 15th-order high-precision nonsymplectic integrator. The second one is a hybrid symplectic integrator. These two integrators allow us to corroborate the results obtained with Mercury 6. For both integrators, we integrate over 10,000 orbits of the binary system GJ 896AB using 20,000 points per orbital period (i.e., one sampling point every 4 days). This allows us to monitor the changes in the orbital parameters of the planet GJ 896Ab with reliable accuracy. Both the best-fit solution reported in Table 3 and the complementary one with Ω = Ω_Ab + 180° were tested. Our results are in agreement with those obtained with Mercury 6. Our best-fit solution is stable over the whole integration time, while the alternative solution becomes unstable in a short period of time. A detailed discussion about the 3D orbital stability of this binary system is out of the scope of this paper, and it will be presented elsewhere.

The present separation between the two stars (∼31.64 au) in this binary system suggests that they were most likely formed in a massive accretion disk by disk fragmentation and not by turbulent fragmentation of the original molecular core. Binary systems formed by turbulent fragmentation are expected to have separations of hundreds or even thousands of au, which better explain the formation of binary systems with wide orbits (e.g., Offner et al. 2016). On the other hand, in the case of disk fragmentation, the lower-mass stellar companion is formed in the outer parts of the disk (probably at a distance of a few hundred au), and then, during the early evolution of the binary system, the stellar companion migrates inward to a closer orbit (e.g., Tobin et al. 2016). The apparently similar spin angles of the two stars (∼210°) and the relatively small, but significant, difference between the spin axis of both stars and the rotation axis of the binary system (≳10°) are consistent with this formation mechanism. As we do not know the position angles of the line of nodes of the rotation plane of both stars, it is possible that the rotation planes are not coplanar, and thus they probably have different mutual inclination angles with respect to the orbital plane of the binary system. Either way, the difference in the inclination angle between the spin axis of the stars and the orbital axis of the binary system (≳10°) suggests that, during the evolution of this binary system, the star–star interaction between both stars may have significantly changed the orientation of the rotation axes of both stars. In addition, the large mutual inclination angle (Φ = 148°) that we find between the orbital planes of the planetary companion GJ 896Ab and the binary system GJ 896AB indicates that something changed the orientation of the orbital plane of the planet, probably from initially being coplanar to be in a retrograde configuration at present time. The most likely origin of this large mutual inclination angle is the star–planet interaction, where the low-mass stellar companion GJ 896B produces some important torque over the orbit of the planetary companion of GJ 896A.

Recent studies have shown that, in the case of planet–planet interaction (planetary systems with at least one external Jovian planet), the mutual inclination angles are ≲10° (see, e.g., Laughlin et al. 2002). In a few cases there is evidence of an inclination between the orbits of the planets of up to ∼40° (see, e.g., McArthur et al. 2010; Dawson & Chiang 2014; Mills & Fabrycky 2017). However, as these studies are based on planetary systems detected with the RV and/or transit techniques, they lack of information about the position angle of the line of nodes, and thus the mutual inclination angles calculated this way are lower limits. In addition, according to the catalog of Exoplanets in Binary Systems (Schwarz et al. 2016), 160 planets have been found associated with 111 binary systems, 79 of which are in S-type orbits (i.e., the planet is orbiting one of the stars; see also Marzari & Thebault 2019). Most of these planets were found using RV and transit techniques, and some of the other planets were found using other techniques, such as imaging and microlensing. In all of them, some of the orbital parameters are missing, in particular, the position angle of the line of nodes of the orbital plane has not been obtained, and in several cases the inclination angle is also unknown. This is the first time that the remarkable full 3D orbital architecture of a binary planetary system has been determined.

The presence of a Jovian planet associated to the main star in a binary system, such as the one we present here, would produce secular perturbations on the lower-mass stellar companion. Such interaction would produce long-term periodic variations in the orbital motion of the secondary (including the eccentricity and the inclination angle), as well as in the orbital parameters of the planetary companion. Given the orbital period of about 229 yr of this binary system, several decades of absolute and relative astrometric observations of the binary system are probably needed to fully constrain the orbital motion of the binary system. As we have mentioned above, we are limited by the time span of the astrometric observations we present here (80 yr for the relative astrometry and 16.6 yr in the case of the absolute astrometry), which is much smaller than the orbital period of the binary system. However, in the analysis we present here, we show that by combining the absolute astrometric observations of the primary and the secondary stars, as well as the relative astrometry of the binary system, in the astrometric fit, the orbital motion of the binary system and the planetary companion are well constrained. We have found that the orbits of the binary system and the planetary companion are somewhat eccentric. A similar astrometric signal due the planetary companion was found using different methods. We obtained a similar astrometric signal using both the periodogram of the astrometric data, and the single-companion astrometric fits of the absolute astrometric data obtained with two different algorithms (least-squares and AGA). This indicates that the astrometric signal is real. The detection of a similar astrometric signal with the full combined astrometric fit (relative plus absolute astrometry) further support the detection of the astrometric signal. In addition, we complemented the analysis of the astrometric observations with the nonlinear least-squares minimization package lmfit (Newville et al. 2020), finding a similar astrometric solution (see Appendix). This strongly indicates that the astrometric signal is real. Furthermore, the fact that the astrometric solution from the different methods is similar indicates that the astrometric signal is consistent with the companion being planetary in origin. A detailed study of the stability of the orbital motions in this system can give important information about the star–star and star–planet interactions. A detailed discussion about the 3D orbital stability of this binary system is out of the scope of this paper, and it will be presented elsewhere.

The astrometric signal that we find is consistent with a planetary companion associated to the M3.5 Dwarf GJ 896A. However, this astrometric signal could be contaminated by the expected astrophysical "jitter" added to the true source position due to stellar activity. It has been estimated that GJ 896A has a radius of ∼0.35 R_⊙ (e.g., Morin et al. 2008; Pearce et al. 2020). Thus, the radius of GJ 896A at a distance of 6.2567 pc is about 0.260 mas. This result indicates that the astrometric signal due to the planetary companion (0.51 mas) is about twice the radius of the host star. In addition, assuming that the radio emission originates within ∼0.194 stellar radius above the photosphere (e.g., Liefke et al. 2008; Crosley & Osten 2018a), the size of the expected "jitter" is about 0.31 mas. However, the analysis that we carried out for the variability of the radio emission of the main star GJ 896A shows that the "flares" seen in several of the observed epochs contribute only with less that 0.12 mas to the expected "jitter" (see Appendix A and panel (d) in Figure 11), compared to the expected 0.31 mas contribution to the "jitter" that we have estimated here. Thus, the expected "jitter" due to the variability of the star is about a factor of 4.3 smaller than the astrometric signal of 0.51 mas observed in GJ 896A (see Table 3). This result supports the detection of the planetary companion GJ 896Ab.

Zoom In Zoom Out Reset image size

To further investigate the validity of the planetary astrometric signal, we have applied several statistical tests comparing the astrometric solutions obtained without and with a planetary companion (see Table 2). The solutions of the best astrometric fits show that when a planetary companion is not included the residuals of the fit have an rms scatter of 0.34 mas and a χ² = 190.91, compared to an rms scatter of 0.20 mas and χ² = 57.41 when the planetary companion is included in the fit. An F-test shows that the probability of the χ² dropping that much (due to the inclusion of the planetary companion) is less than 4% by mere fluctuations of noise. Using the Bayesian information criterion (BIC), we find that the inclusion of the planetary companion is preferred with a Delta BIC = −109.24, this is a significant difference between our best-fit model with the planetary companion and the one without it. Statistically, an absolute value of Delta BIC of more than 10 implies that the model with the planetary companion better reproduces the data without overfitting them by including more free parameters in the model. Similar results were obtained using the Akaike information criterion (AIC) with a Delta AIC = −119.50. We therefore conclude that there is a very high probability that the planetary companion GJ 896Ab is orbiting the main star GJ 896A.

The large proper motions of the binary system may also contribute to the expected "jitter" of the star. The change in position of GJ 896A due to the proper motions of the binary system, Δ_PM (A) in mas, can be estimated by:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{PM}}(A)=\sqrt{{\mu }_{\alpha }^{2}+{\mu }_{\delta }^{2}}\times \left(\delta t/8766\right),\end{eqnarray} \tag{ 14 }$$

where δ t is the time span of each observed epoch in hours (between 3 and 5 hr on source). Using the estimated proper motions of the binary system (μ = 572.71 mas yr⁻¹), we obtain that the maximum expected contribution to the "jitter" by the proper motions of the system is Δ_PM (A) = 0.16 mas in a time span of 2.5 hr (half the observing time), which is about 3 times smaller than the observed astrometric signal. We also obtain that the contribution of the orbital motion of GJ 896A around the center of mass to the expected "jitter" of the star (about 0.009 mas in a time span of 2.5 hr) is smaller than that estimated from the proper motions of the star. Adding in quadrature all these possible contributions, the total expected "jitter" is about 0.2 mas, which is still about a factor of 2.6 smaller than the astrometric signal due to the planetary companion GJ 896Ab. The highest contribution to the expected "jitter" is that due to the proper motions of the star. However, the astrometric signal is observed in both directions (R.A. and decl.), while the proper motion of the star is basically in the east direction. In addition, the contribution of the proper motions of the star to the expected "jitter" probably averages out when we integrate over the time span of the observations. In addition, given the synthesized beam of the images (about 2.7×1.1 mas, in average), the magnitude of the proper motions is in fact too small to even affect the estimated size of the source. This suggests that the contribution of the proper motions to the expected "jitter" is probably much smaller than we have estimated here. Thus, the expected "jitter" that we estimate here is most likely an upper limit. This further indicates that the astrometric signal is real and planetary in origin.