EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 565 (May 2014) A&A, 565 (2014) A104 Full HTML
Free Access
Issue
A&A
Volume 565, May 2014
Article Number A104
Number of page(s) 8
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/201423799
Published online 19 May 2014

© ESO, 2014

1. Introduction

A hierarchical (Evans 1968) multibody star system is believed to be formed through one or more formation channels. First, van den Berk et al. (2007) considered interaction/capture mechanisms during the formation and dynamical evolution of a globular star cluster. A second mechanism that might explain the existence of such systems is that they might be formed directly from a massive primordial disk involving accretion processes and/or local disk instabilities (Lim & Takakuwa 2006; Duchêne et al. 2007; Marzari et al. 2009). A third mechanism follows a chaotic erosion-process of a nonhierarchical star system by angular momentum and energy exchange via mutual gravitational interactions. In the latter case, and considering an initial triple-system, Reipurth (2000) provided a schematic outline of three stages that might produce a close binary system with a circumbinary disk from redistribution of circumstellar material after chaotic interactions. The formation of the tightly bound central binary is followed by the transport of the third member to a wider orbit as a result of energy conservation. In extreme cases, this can result in the third member being ejected completely from the system, leaving a tightly packed close binary on a quasi-Keplerian orbit.

A particular example of a hierarchical multibody system is a so-called circumbinary system (also known as companions on P-type orbits Schwarz et al. 2011) in which one or more massive objects orbit a binary star system. Such systems have recently been discovered by Kepler and the Planet Hunters community1 (Doyle et al. 2011; Welsh et al. 2012; Orosz et al. 2012a,b). Ofir et al. (2009) presented the results for a search of circumbinary companions based on CoRoT data. The planet orbiting the binary PH-1 (Schwamb et al. 2013; Kostov et al. 2013) is a particularly exotic example of such a system. Here the binary is a member of a quadruple (or quaternary) hierarchical system where two binary pairs form a gravitationally bound star system. Similar in nature, though with no evidence of planetary companions, is the HD 98800 quadruple system (Furlan et al. 2007). Other types of hierarichical star systems reside in so-called S-type (Schwarz et al. 2011) configurations where one body is orbiting one component of a binary pair. Several examples of such systems have been reported in the literature (Neuhäuser et al. 2007; Chauvin et al. 2007).

A well known technique to detect a hierarchical circumbinary systems is to measure and monitor timing variations of the mid-eclipse times of the central binary (also known as times of minimum light). For a detailed description of its application on detecting circumbinary companions of planetary mass we refer to Deeg et al. (2000) and Doyle & Deeg (2004). This technique has recently begun to be applied to the excellent timing data collected by the Kepler mission, resulting in the recent announcement of the first substellar-mass circumbinary companion discovered from that data which orbits KIC002856960 (Lee et al. 2013a). The fundamental principle of the light-travel time effect (LTT2) makes use of the motion of the binary around the total barycenter of the system. Due to the finite speed of light, the eclipses exhibit delays or advances in the timings of minimum light depending on the orbital position of the binary relative to the observer (Irwin 1952, 1959). This method is particularly attractive because it is observationally time-effective, involving only photometric CCD measurements.

In recent times, single and multibody substellar circumbinary companions to known eclipsing binary systems have been proposed using ground-based timing measurements (Lee et al. 2011, 2012, 2013a,b). The same technique was used to detect candidate circumbinary companions of planetary nature: CM Draconis (Deeg et al. 2000, CM Dra, one companion), DP Leonis (Qian et al. 2010a; Beuermann et al. 2011, DP Leo, one companion), HW Virginis (Lee et al. 2009a, HW Vir, two companions), NN Serpentis (Qian et al. 2009; Beuermann et al. 2010, 2013, NN Ser, two companions), UZ Fornazis (Potter et al. 2011, UZ For, two companions) and HU Aquarii (Qian et al. 2011, HU Aqr, two companions) and QS Virginis (Qian et al. 2010b; Almeida et al. 2011, QS Vir, one or two companions). Recently, additional circumbinary companions were proposed from using the LTT effect: RR Caenis (Qian et al. 2012a, RR Cae, one companion), NSVS 14256825 (Almeida et al. 2013, two companions), and NY Virginis (Qian et al. 2012b, NY Vir, one companion).

However, the existence of the proposed multibody systems has been challenged as a result of a number of studies of the dynamical stability of their orbital architectures. The proposed companions around HU Aqr have been studied in detail, and a series of studies have revealed them to be dynamically unfeasible (Horner et al. 2011; Hinse et al. 2012a; Wittenmyer et al. 2012a; Funk et al. 2011, 2012; Goździewski et al. 2012). The same is true of HW Vir (Horner et al. 2012b) and NSVS 14256825 (Wittenmyer et al. 2013; Hinse et al. 2014). Indeed, of the systems studied in this way, the only one to withstand dynamical scrutiny, are the companions around NN Ser (Horner et al. 2012a; Beuermann et al. 2013), although a recent study of the evolution of the central binary in the NN Ser system suggests that it is unlikely that planetary companions on the proposed orbits could have survived the system’s post-main sequence evolution intact (Mustill et al. 2013). Furthermore, Hinse et al. (2012b) showed that the two substellar companions orbiting the SZ Herculis (Lee et al. 2012, SZ Her, two companions) binary also follow highly unstable orbits. The same situation is also seen for the QS Vir system where Horner et al. (2013) was able to show that the proposed two-companion system is highly unstable.

Finally, Parsons et al. (2010) presented photometric follow-up observations of a number of eclipsing post-common-envelope binaries for which they have been able to rule out previous claims for single-object circumbinary companions (e.g., Qian et al. 2009, 2010b).

In this work we re-examine the observed timing dataset of RZ Dra as presented by Yang et al. (2010). These authors proposed the existence of two additional low-mass dwarfs from two distinct quasi-sinusoidal variations in the times of mutual eclipses. Section 2 presents a dynamical stability analysis of the nominal orbital parameters as derived by Yang et al. (2010). We then continue and give an outline of our data analysis based on the light-travel time effect and describe the least-squares methodology in Sect. 4, where we present a new best-fit model of the two proposed companions. A dynamical analysis of our new model is presented in Sect. 5 and we finish with concluding remarks in Sect. 6.

2. Orbital stability of the Yang et al. (2010) model

Table 1

Orbital and mass parameters of the two substellar companions proposed to orbit the RZ Dra Algol-type binary system, taken from Yang et al. (2010, their Table 5).

Yang et al. (2010) proposed that the observed variations in the timings of eclipses between the components of the Algol-type binary star system RZ Dra might be explained by the presence of two massive, unseen companions moving on eccentric long-period orbits, as detailed in Table 1. On their nominal best-fit orbits, these two proposed objects can approach one another remarkably closely – with the innermost object (RZ Dra (AB) C)3 having a best-fit apastron distance of 18.06 AU, and the outermost (RZ Dra (AB) D) having a nominal best-fit periastron distance of 17.02 AU. In other words, the nominal best-fit models for the two proposed companions have significant overlap – these stellar-mass companions can cross one another’s orbits.

thumbnail Fig. 1

Dynamical stability of the RZ Dra system as a function of the semi-major axis, a, and eccentricity, e of RZ Dra (AB) D, for the solution presented in Yang et al. (2010). The nominal best-fit orbit for RZ Dra (AB) D presented in that work is located in the center of the red box, at a = 23.88 AU and e = 0.287. The 1σ uncertainties on the semi-major axis and eccentricity are shown by the red crosshairs radiating from the box. The lifetime at each of the 1681 a − e locations plotted in this figure (41 unique a and e values spanning ± 3σ from the nominal best-fit orbit) is the result of 75 distinct simulations, spanning 5 unique values of mean anomaly, M, and 15 unique values of the argument of periastron, ω. The mean lifetime can be seen to vary between ≃100 and ≃1600 years, demonstrating that the system as proposed in Yang et al. (2010) exhibits extreme dynamical instability. (Color version online.)

As we have shown in previous works (e.g. Horner et al. 2011, 2012b; Wittenmyer et al. 2012a, Wittenmyer et al. 2013), such mutually encountering orbital architectures typically lead to significant dynamical instability, often on timescales of just a few hundred years. Given the age of the systems studied, such instability clearly rules out a planetary origin for the observed variations in these systems, suggesting instead that some other astrophysical effect is the cause of the observed timing variations. It is therefore clearly important to examine the dynamical feasibility of the proposed companions to RZ Dra, to determine whether the “companion hypothesis” of Yang et al. (2010) stands up to close scrutiny.

To study the dynamical stability of the proposed companions (RZ Dra (AB) C and D), we used the HYBRID integration algorithm within the n-body dynamics package mercury (Chambers 1999). We used a constant integration time-step of one day in all our orbit integrations. The error tolerance parameter was set to one part in 1012, which ensures accurate integrations (using the Bulirsch-Stoer algorithm) of high-eccentricity orbits and possible close encounters. As in our earlier work (e.g. Marshall et al. 2010; Horner et al. 2011, 2012b; Robertson et al. 2012a,b; Wittenmyer et al. 2012b), we performed many discrete simulations, each following the dynamical evolution of one potential RZ Dra system. As in those earlier works, we held the initial orbital elements of the better constrained companion (in this case, RZ Dra (AB) C) fixed at their nominal best-fit values (as detailed in Table 1). We then systematically varied the initial semi-major axis, a, eccentricity, e, argument of periastron, ω, and mean anomaly, M, of the orbit of RZ Dra (AB) D across their full ±3σ range, thereby creating a wide variety of test systems, which were then integrated for a period of 100 Myr. In each test system, the dynamical evolution of the two companions was followed until a break-up of the system was detected – either through the ejection of one or the other of the two companions (to a distance of 50 AU – sufficiently distant so that the substellar components could only reach this distance if significant mutual interaction had occurred), through mutual collision between the two companions, or when one of the companions collided with the central bodies. If the test system fell apart in this manner, the time of its demise was recorded.

In total, the dynamical evolution of 126 075 trial systems was considered within our numerical stability analysis. Forty-one distinct values of a were tested, spread evenly across the ±3σ range (i.e. between a = 15.99 and a = 31.77 AU). For each of these a values, 41 unique values of e were tested, again evenly distributed over the ±3σ range (between e = 0.266 and e = 0.308). For each of these 1681 a − e pairs, 15 unique ω values were tested, and at each of these 25 215 a − e − ω trials, 5 unique mean anomalies were tested.

thumbnail Fig. 2

Results from numerical tests considering the orbits of the two companions as shown in Fig. 5. Here we consider an initial condition (IC 1) where the outer companion has been assigned an initial eccentricity of 0.90 (black dot in Fig. 5). All other parameters are as shown in Table 2. For the RADAU algorithm (top panels) we used an initial time step of 0.001 days. For the HYBRID algorithm (bottom panels) we used a constant time step of one day. For both algorithms we used an accuracy parameter of 10-12. Orbital elements are sampled every day and are plotted until the ejection time of RZ Dra(AB)C at around 22 years.

The results of our simulations can be seen in Fig. 1. It is immediately apparent that the system as proposed in Yang et al. (2010) is extremely unstable – with mean lifetimes ranging between ≃100 years and ≃1600 years. Indeed, the longest-lived of the 126 075 systems tested fell apart after just 173 000 years. A remarkable 21% of the test architectures were unstable on timescales shorter than 100 years (25 805 of the 126 075 systems tested). More than 88% of the systems fell apart within 1000 years (111 225 of the systems tested). While it is clear from the figure that the stability of the proposed system does increase somewhat as the separation of the companions is increased, this effect is clearly insufficient to allow the proposed companions to survive on sufficiently long timescales for their existence to be reasonable. The likelihood of by chance observing the RZ Dra system within the last thousand years (or even the last hundred thousand years) before the destruction or ejection of the companions seems low.

We tested the accuracy of our numerical computations by comparing the results from several HYBRID integrations with results obtained from accurate RADAU integrations. In Fig. 2 we show the results of integrating an orbit (black dot in Fig. 5) of the outer companion with initially high eccentricty. For the chosen integration time step and accuracy parameters we conclude that the HYBRID integrations are reliable since the two orbits seem to follow the same time evolution until the inner companion is ejected after some 20 years. Orbits integrated with the RADAU algorithm are generally considered as producing reliable results. We carried out similar spot tests for other initial conditions with similar outcomes.

As was the case for a number of other proposed circumbinary systems (Horner et al. 2011; Wittenmyer et al. 2012a, HU Aqr); HW Vir (Horner et al. 2012b); and NSVS 14 256 825 Wittenmyer et al. 2013), the RZ Dra system proposed by Yang et al. (2010) does not pass dynamical scrutiny. If massive companions do exist within that system, they must clearly move on orbits far different from those proposed in the discovery work. The likelihood of observing the RZ Dra system in the last 1000 years (or even 100 000 years) before the destruction or ejection of its companions by chance is vanishingly small because the host system’s lifetime extends several hundred million years.

3. Data analysis and LTT model

We used the same timing data4 set as Yang et al. (2010) along with their adopted timing precisions. To remove systematic offsets in the measured times of minimum light we transformed the heliocentric julian dates (HJD) time stamps (all assumed to be in the UTC time standard) to barycentric julian dates (BJD) times in the barycentric dynamical time (TDB) time standard using the transformation routines in Eastman et al. (2010). Timing measurements earlier than JD 2 433 266.0 are limited by the DE405 JPL Planetary Ephemeris5 (Giorgini et al. 1996) and hence were transformed by omitting the light-travel time corrections from the Einstein and Shapiro effects (J. Eastman, priv. comm.). Since these effects contribute with timing precision to much less than one second, we judged them to be negligible in this work and thus omitted them.

For an idealized, unperturbed and isolated binary system, the linear ephemeris of future/past (primary) mid-eclipse events can be computed from TC(E)=T0+P0E,\begin{equation} T_{\rm C}(E) = T_{0} + P_{0}E, \label{ephemeriseq} \end{equation}(1)where E denotes the cycle number, T0 is the reference epoch of some primary eclipse, and P0 measures the eclipsing period of RZ Dra. A linear regression on the 680 recorded eclipse times allows determining P0 to a high precision. We chose to place the reference epoch at the same date as in Yang et al. (2010). This is relatively close to the middle of the overall data set and avoids parameter correlation between T0 and P0 during the fitting process. In the following we briefly outline the LTT model we used in this work.

3.1. Analytic LTT model

We used a model similar to the formulation of a single light-travel time orbit introduced by Irwin (1952). In this model, the two binary components are assumed to represent one single object with a total mass equal to the sum of the masses of the two stars. If a circumbinary companion exists, then the combined binary mass follows an orbit around the system barycenter. The eclipses are then given by Eq. (1). This defines the LTT orbit of the binary. The underlying reference system has its origin at the center of the LTT orbit.

Following Irwin (1952), if the observed mid-eclipse times exhibit a sinusoidal variation (due to an unseen companion), then the quantity O−C measures deviation from the linear ephemeris possibly attributable to the light-travel time effect and is given by (O−C)(E)=TO(E)−TC(E)=∑i=12τi,\begin{equation} ({\rm O}-{\rm C} )(E) = T_{\rm O}(E) - T_{\rm C}(E) = \sum_{i=1}^2\tau_{i}, \label{omc} \end{equation}(2)where TO denotes the measured time of an observed mid-eclipse, summed over i circumbinary companions. When plotting the quantity “O − C” as a function of time, we will use the nation “eclipse timing diagram”.

We note that τ1 + τ2 is the combined LTT effect from two separate two-body LTT orbits. The quantity τi is given by the following expression for each companion (Irwin 1952): τi=Kb,i[1−eb,i21+eb,icosfb,isin(fb,i+ωb,i)+eb,isinωb,i],\begin{equation} \tau_{i} = K_{b,i}\left[\frac{1-e_{b,i}^2}{1+e_{b,i}\cos f_{b,i}} \sin (f_{b,i}+\omega_{b,i}) + e_{b,i}\sin \omega_{b,i}\right], \label{taueq} \end{equation}(3)where Kb,i = ab,isinIb,i/c is the semi-amplitude of the light-time effect (in the O−C diagram) with c measuring the speed of light and Ib,i is the line-of-sight inclination of the LTT orbit relative to the sky plane, eb,i the orbital eccentricity, fb,i the true anomaly and ωb,i the argument of pericenter. The five model parameters for a single LTT orbit are given by the set (ab,isinIb,i,eb,i,ωb,i,Tb,i,Pb,i). The time of pericenter passage Tb,i and orbital period Pb,i are introduced by expressing the true longitude as a time-like variable via the mean anomaly M = nb,i(TO − Tb,i), with nb,i = 2π/Pb,i denoting the mean motion of the combined binary in its LTT orbit. Computing the true anomaly as a function of time (or cycle number) requires the solution of Kepler’s equation. We refer to Hinse et al. (2012a) for more details.

4. Least-squares fitting and parameter uncertainties

Table 2

Best-fit parameters for the LTT orbits of RZ Dra corresponding to Fig. 3.

thumbnail Fig. 3

Best-fit two-Kepler LTT model without secular parameter with reduced χr2=12.7\hbox{$\chi_{\rm r}^2=12.7$}. The best-fit model parameters are shown in Table 2. In the top panel the LTT signal with smaller/larger amplitude corresponds to the inner/outer companion. The bottom panel shows the residuals between the data and best-fit model. We note that, judging by eye, there is a slight asymmetry in the residuals with O−C > 0 for E> 0.

The overall procedure of finding a best-fit model follows the methodology as outlined in Hinse et al. (2012a) using a least-squares minimization algorithm as implemented in MPFIT (Markwardt 2009). We considered several models. First we started to fit a single-companion LTT model to the data with seven freely varying parameters (680−7 degrees of freedom). Considering a few ten thousand initial random guesses (all parameters randomized as in Hinse et al. (2012a)) we found the reduced best-fit statistic χr2\hbox{$\chi^2_{\rm r}$} to be no less than 104 which indicates a one-companion model that is inadequate to describe the data. We then fitted for a two-companion model with all twelve parameters allowed to vary freely (680−12 degrees of freedom). A significantly better fit was found with χr2≃13\hbox{$\chi^2_{\rm r}\simeq 13$}. To investigate the possibility of a secular trend in the timing data set (representing a linear period decrease by mass transfer) we conducted an independent best-fit model search by adding an additional term to Eq. (2), that is quadratic in time (or cycle number). However, upon inspecting the residual plots we were not convinced about the necessity of adding an additional free parameter. We judged that no obvious secular parabolic trend was present in the dataset. For a one-companion model with a quadratic trend we refer to Erdem et al. (2011). Our best-fit two-companion model and resulting best-fit parameters are shown in Fig. 3 and Table 2.

During this study (in conjunction with previous studies, e.g. Horner et al. (2012b)) we found that the formal parameter errors as obtained from the best-fit co-variance matrix are too optimistic. We therefore determined parameter errors from Monte Carlo bootstrap experiments (Goździewski et al. 2012). For each bootstrap dataset we used the best-fit model as our initial guess and recorded the resulting parameters. We generated 106 bootstrap simulations. From the resulting distributions we determined the 1σ uncertainty in each parameter to be the 68.2% percentile interval. A graphical representation of our confidence intervals is given in Fig. 4 showing a two-dimensional joint-confidence parameter scan of the χr2\hbox{$\chi^2_{\rm r}$} statistic. Here we consider the dependency between the period and eccentricity parameters of the outer companion as well as the reference epoch and the binary period. The 1-σ confidence level agrees well with the outcome of our bootstrap experiment.

5. Dynamical stability of our new best-fit model

thumbnail Fig. 4

Two-dimensional parameter scans of the χr2\hbox{$\chi^2_{\rm r}$} statistic in a neighborhood around the best-fit model (marked with an asterisk) considering the two-companion model. The color-bar is scaled with the value of our best-fit model representing 1. Countour lines show the 1-σ (left panel) and 1-, 2-, 3-σ (right panel) confidence level curves around the best-fit model. (Color version online.)

After deriving a new Keplerian model for the observed timing variations in the RZ Dra system, we once again carried out a detailed dynamical study of the proposed orbits. Following the same procedure as described in Sect. 2, we held the initial orbital elements of RZ Dra (AB) C fixed at their nominal best-fit values and created 126 075 test systems in which the initial semi-major axis, eccentricity, mean anomaly, and argument of periastron for RZ Dra (AB) D were systematically varied across the ± 3σ range around their nominal best-fit values. Once again, we considered 41 discrete values of a and e, 15 unique values of ω, and 5 unique values of M. Owing to the high values of a and e allowed by the ± 3σ range, the initial apastron distance for RZ Dra (AB) D can range as high as 88.06 AU – and so we set the ejection distance as 150 AU in this case, so that objects are only considered to be ejected if they have experienced significant orbital perturbations.

The results of our simulations are shown in Fig. 5. As was the case for the companions proposed by Yang et al. (2010), the new model parameters lead to extremely unstable orbits. The entire region spanned by the ± 1σ uncertainties features mean lifetimes measured in thousands of years, and the only (very small) regions where the stability is measured in hundreds of thousands of years, or even just over a million years, are located at a large distance from the nominal best-fit model.

On the basis of these results, we conclude that the observed variations in the timing of the eclipses in the RZ Dra system are not the result of the gravitational influence of two unseen companions – such companions would be so massive and move on such extreme orbits that they would destabilize one another on timescales of hundreds of years – far shorter than the typical lifetime of the host binary system. We recall that the nature of near-contact binaries is not well known. Currently, there is poorly consensus that near-contact binaries are intermediate objects observed between two contact phases. In contact binaries mass and energy exchange between the two components can be significant, which might and probably does affect the interior structure of the stars. The orbital distance between the two stellar surfaces is just a few percent of the stellar radii. We therefore conjecture that stellar mass-transfer is possible via a flare that moves matter from one component to the other. This could cause a sudden period-change of the eclipsing binary, resulting in the observation of ETVs.

6. Discussion and conclusion

Based on measurements of primary and secondary eclipse times, Yang et al. (2010) proposed an interpretation of the eclipse timing variation due to two circumbinary companions. The two bodies they proposed are very low mass stars with minimum mass 0.07 M⊙ for the inner and 0.18 M⊙ for the outer companion. According to their best-fit model, the inner companion apocenter distance is at ≃18 AU and the outer companion pericenter distance is ≃17 AU, implying a crossing-orbit architecture. We tested the dynamical feasibility of the proposed best-fit orbits by exploring a large grid of initial conditions using direct n-body integrations. The results of our simulations exploring the (a,e)-space are shown in Fig. 1. All of the tested orbits were highly unstable and resulted in either the break-up of the system or collisions between the two bodies. We found that the mean lifetimes range between 100 and 1600 years. We did not consider mutually inclined orbits between the two companions. We refer to Horner et al. (2011); Hinse et al. (2012b); Horner et al. (2013, 2014); Wittenmyer et al. (2014) who investigated the effect of mutually inclined orbits which resulted in little improvement of the overall dynamical stability. However, the only exception to this finding were the scenarios for which the two planets were placed on anti-coplanar orbits – in other words, where they moved in the same plane, but with a mutual orbital inclination of 180 degrees. This setup predictably led to extremely stable systems whenever the two planets were not placed on orbits that crossed one another (which led to extreme instability). Such an orbital architecture, while of theoretical interest, seems highly physically implausible, and so needs not to be considered further at this time.

We then searched for and determined a new best-fit model using the complete set of timing measurements transformed from HJD to BJD time standard. We found a new best-fit model with parameters shown in Table 2. Compared with the orbital parameters given in Yang et al. (2010), our new model assigned higher eccentricities to the orbits of the companions and increased the mass of the outer companion to 0.4 M⊙. We then tested the our new model parameters for dynamical stability and found lifetimes of only a few hundred years (Fig. 5).

thumbnail Fig. 5

Dynamical stability of the RZ Dra system as a function of the semi-major axis, a, and eccentricity, e of RZ Dra (AB) D for the alternative solution. The nominal best-fit orbit for RZ Dra (AB) D is located in the center of the red box, at a = 26.6 AU and e = 0.62. The 1σ uncertainties on the semi-major axis and eccentricity are shown by a crosshair. Almost all the tested models, particularly those within the ± 1σ uncertainties, are dynamically unstable on remarkably short timescales. A few regions are stable on timescales of around a million years, but these are all located at the extreme edges of the plot, almost ± 3σ away from the nominal best-fit orbit. The black dot is the location of a test orbit as discussed in Sect. 2. This figure is available in color in electronic form.

On the basis of our dynamical simulations we conclude that the two-companion hypothesis around RZ Dra does not pass scrutiny. At this point we would like to emphasize the robustness of our stability analysis. First we carried out tests that demonstrate that the results obtained from HYBRID integrations are reliable. Second, our dynamical setup replaces the two binary components as a single massive object. In that sense, our setup considers a system that favors dynamical stability by ignoring the gravitational perturbations from an additional body (significant in mass). If no stable orbits are found in this simplified system, it is generally hard to conceive how stable orbits are ensured by adding an additional pertubing force to the system. From an intuitive point of view, adding more perturbing forces will increase the effect and possibility of chaos and hence favors orbital instabilities. Another aspect of replacing the two binary components by a single body is a question of being consistent with the LTT model, which assumes the binary to be a single massive object.

It is clearly important, therefore, to consider what other mechanisms might account for the observed variations. In Fig. 3 we have some indication from an asymmetric distribution of the residuals of additional effects that might be causing a period change. The orbital periods of many binary systems have varied because of some combination of secular and/or cyclical variations. Generally, the quasi-sinusoidal timing variations could be produced by one or more physical causes: (1) apsidal motion in an elliptical orbit of the binary; (2) a light-travel time (LTT) effect (or several) due to additional companion(s); or (3) cyclical changes of magnetic activity of the component stars with deep convective envelope. The secular variations can be interpreted as being caused by either mass transfer between the two component stars or by angular momentum loss (AML).

Because the binary has a circular orbit, the cyclical variations of RZ Dra cannot be described by apsidal motion. Furthermore, Yang et al. (2010) ruled out the magnetic activity cycles because the variations of the gravitational quadrupole moment (ΔQ) are two orders of magnitude smaller than typical values of 1051−1052 for close binaries. This finding is also supported by a recent study (Lanza 2006) that indicated that the magnetic mechanism (Applegate model) is not sufficiently effective to explain the period modulation of close binaries with a late-type secondary.

However, the eclipsing pair of RZ Dra is a semi-detached binary whose less massive secondary fills its inner Roche lobe, see Yang et al. (2010). In such semi-detached binaries, a secular variation (quadratic term) might be produced through mass transfer from the secondary to the primary star, AML due to a magnetic stellar wind, or the combination of these two mechanisms. The mass transfer causes a period increase (upward parabola), while the AML causes a period decrease (downward parabola). We considered a quadratic + two-LTT model (see Sect. 4), but we were unable to convincingly detect either one of these trends in the residuals. However, a secular variation may be hidden in the timing data set, and this system may be in a weak phase of mass transfer. We refer to Erdem et al. (2011), who considered a quadratic + one-LTT model.

Furthermore, the presence of systematic residuals shown in Fig. 3 might be caused by an apparent phase shift of the real conjunctions due to asymmetrical eclipse minima originating from starspot activity (Lee et al. 2009b). The effects of starspots on timing measurements of eclipsing binaries were also studied by Watson & Dhillon (2004).

At present we cannot rule out that most of the timing measurements have been underestimated and hence the plotted errorbars in Fig. 3 might be much larger than stated in the literature. To obtain an idea of the timing uncertainty we scrutinized Table 4 in Yang et al. (2010) and noticed the visual recording of the same secondary eclipse (presumably by two observers) on the night of HJD 2 442 984 (July 24, 1976; we refer to the VO ASCII data file available online). The first observer measured the secondary eclipse to be at HJD 2 442 984.637, the second observer measured the same eclipse event to occur at HJD 2 442 984.641. These observations suggests a larger uncertainty since the times of minimum light from the visual observations differ by more than 5 min (over 300 s). This assumes, of course, that the entry in the corresponding VO file is neither a duplicated entry or typing error. If times of minimum light truly have a precision of 0.003 days, then most of the variation seen in Fig. 3 would be within the noise-level of timing estimates. In general, Eastman et al. (2010) recommended that uncertainties of at least 60 s should be used for the timing precision if the time standard has not been specified explicitly.

We therefore encourage future follow-up observations of eclipsing binaries to obtain as precise timing measurements as possible for RZ Dra and other systems that are mainly characterized via visual measurements of the minima. Several monitoring programs are currently running or in planning (Sybilski et al. 2010; Pribulla et al. 2012) within the context of searching for circumbinary companions of planetary, substellar, and stellar nature.

2

sometimes also referred to as LITE (Borkovits & Hegedüs 1996).

3

We adopt the naming convention as proposed by Hessman et al. (2010).

4

We note that new timing data exist published by Erdem et al. (2011). These data points were not included in the analysis presented here.

5

http://ssd.jpl.nasa.gov/?horizons

Acknowledgments

We would like to thank the anonymous referee for improving this research paper. Research by TCH is carried out at the Korea Astronomy and Space Science Institute (KASI) under the KRCF (Korea Research Council of Fundamental Science and Technology) Young Scientist Research Fellowship Program. Numerical simulations were carried out on the “Pluto” computing cluster at KASI and the SFI/HEA Irish Centre for High-End Computing (ICHEC). T.C.H., J.W.L., C.U.L. & J.H.P. acknowledge support from KASI registered under grant number 2013-9-400-00/2014-1-400-06. Astronomical research at Armagh Observatory is funded by the Department of Culture, Arts and Leisure (DCAL). JPM is supported by Spanish grant AYA 2011-26202. This work was supported by iVEC through the use of advanced computing resources located at the Murdoch University, in Western Australia. This research has made use of NASA’s Astrophysics Data System (ADS).

References

  1. Almeida, L. A., & Jablonski, F. 2011, IAU Symp., 276, 495 [NASA ADS] [Google Scholar]
  2. Almeida, L. A., Jablonski, F., & Rodrigues, C. V. 2013, ApJ, 766, 11 [NASA ADS] [CrossRef] [Google Scholar]
  3. Beuermann, K., Hessman, F. V., Dreizler, S., et al. 2010, A&A, 521, L60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Beuermann, K., Buhlmann, J., Diese, J., et al. 2011, A&A, 526, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Beuermann, K., Dreizler, S., & Hessman, F. V. 2013, A&A, 555, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Borkovits, T., & Hegedüs, T., 1996, A&AS, 120, 63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Chambers, J. E., 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  8. Chauvin, G., Lagrange, A.-M., Udry, S., & Mayor, M. 2007, A&A, 475, 723 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Deeg, H. J., Doyle, L. R., Kozhevnikov, V. P., et al. 2000, A&A, 358, 5 [Google Scholar]
  10. Doyle, L. R., & Deeg, H.-J. 2004, IAU Symp., 213, 80 [Google Scholar]
  11. Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
  12. Duchêne, G., Bontemps, S., Bouvier, J., et al. 2007, A&A, 476, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935 [NASA ADS] [CrossRef] [Google Scholar]
  14. Erdem, A., Zola, S., & Winiarski, M. 2011, New Astron., 16, 6 [NASA ADS] [CrossRef] [Google Scholar]
  15. Evans, D. S., 1968, QJRAS, 9, 388 [NASA ADS] [Google Scholar]
  16. Funk, B., Eggl, S., Gyergyovits, M., et al. 2011, EPSC-DPS, 1725 [Google Scholar]
  17. Funk, B., Eggl, S., Gyergyovits, M., et al. 2012, IAU Symp., 282, 446 [NASA ADS] [Google Scholar]
  18. Furlan, E., Sargent, B., Calvet, N., et al. 2007, ApJ, 664, 1176 [NASA ADS] [CrossRef] [Google Scholar]
  19. Giorgini, J. D., Yeomans, D. K., Chamberlin, A. B., et al., 1996, BAAS, 28, 1158 [NASA ADS] [Google Scholar]
  20. Goździewski, K., Nasiroglu, I., Słowikowska, A., et al. 2012, MNRAS, 425, 930 [NASA ADS] [CrossRef] [Google Scholar]
  21. Hessman, F. V., Dhillon, V. S., Winget, D. E., et al. 2010, unpublished [arXiv:1012.0707H] [Google Scholar]
  22. Hinse, T. C., Lee, J. W., Haghighipour, N., et al. 2012a, MNRAS, 420, 3609 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hinse, T. C., Goździewski, K., Lee, J. W., et al. 2012b, AJ, 144, 34 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hinse, T. C., Lee, J. W., Goździewski, K., et al. 2014, MNRAS, 438, 307 [NASA ADS] [CrossRef] [Google Scholar]
  25. Horner, J., Marshall, J. P., Wittenmyer, R. A., & Tinney, C. G. 2011, MNRAS, 416, L11 [NASA ADS] [CrossRef] [Google Scholar]
  26. Horner, J., Wittenmyer, R. A., Hinse, T. C., & Tinney, C. G. 2012a, MNRAS, 425, 749 [NASA ADS] [CrossRef] [Google Scholar]
  27. Horner, J., Hinse, T. C., Wittenmyer, R. A., Marshall, J. P., & Tinney, C. G. 2012b, MNRAS, 427, 2812 [NASA ADS] [CrossRef] [Google Scholar]
  28. Horner, J., Wittenmyer, R. A., Hinse, T. C., et al. 2013, MNRAS, 435 2033 [NASA ADS] [CrossRef] [Google Scholar]
  29. Horner, J., Wittenmyer, R. A., Hinse, T. C., & Marshall, J. P. 2014, MNRAS, 439, 1176 [NASA ADS] [CrossRef] [Google Scholar]
  30. Irwin, J. B. 1952, ApJ, 116, 211 [NASA ADS] [CrossRef] [Google Scholar]
  31. Irwin, J. B. 1959, AJ, 64, 149 [NASA ADS] [CrossRef] [Google Scholar]
  32. Kostov, V. B., McCullough, P., Hinse, T. C., et al. 2013, ApJ, 770, 52 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lanza, A. F. 2006, MNRAS, 369, 1773 [NASA ADS] [CrossRef] [Google Scholar]
  34. Lee, J. W., Kim, S.-L., Kim, C.-H., et al. 2009a, AJ, 137, 3181 [NASA ADS] [CrossRef] [Google Scholar]
  35. Lee, J. W., Youn, J.-H., Lee, C.-U., et al. 2009b, AJ, 138, 478 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lee, J. W., Kim, S.-L., Lee, C.-U., et al. 2011, AJ, 142, 12 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lee, J. W., Lee, C.-U., Kim, S.-L., et al. 2012, AJ, 143, 34 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lee, J. W., Kim, S.-L., Lee, C.-U., et al. 2013a, ApJ, 763, 74 [NASA ADS] [CrossRef] [Google Scholar]
  39. Lee, J. W., Hinse, T. C., & Park, J.-H. 2013b, AJ, 145, 100 [NASA ADS] [CrossRef] [Google Scholar]
  40. Lim, J., & Takakuwa, S. 2006, ApJ, 653, 425 [NASA ADS] [CrossRef] [Google Scholar]
  41. Markwardt, C. B. 2009, Astronomical Data Analysis Software and Systems XVIII, eds. D. A., Bohlender, D., Durand, & P. Dowler, Conf. Ser., 411, 251 [Google Scholar]
  42. Marshall, J., Horner, J., & Carter, A. 2010, International Journal of Astrobiology, 9, 259 [NASA ADS] [CrossRef] [Google Scholar]
  43. Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Mustill, A., Marshall, J. P., Villaver, E., et al. 2013, MNRAS, 437, 1404 [Google Scholar]
  45. Neuhäuser, R., Mugrauer, M., Fukagawa, M., Torres, G., & Schmidt, T. 2007, A&A, 462, 777 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Ofir, A., Deeg, H. J., & Lacy, C. H. S. 2009, A&A, 506, 445 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, Science, 337, 1511 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  48. Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, ApJ, 758, 87 [NASA ADS] [CrossRef] [Google Scholar]
  49. Parsons, S. G., Marsh, T. R., Copperwheat, C. M., et al. 2010, MNRAS, 407, 2362 [NASA ADS] [CrossRef] [Google Scholar]
  50. Potter, S. B., Romeo-Colmenero, E., Ramsay, G., et al. 2011, MNRAS, 416, 2202 [NASA ADS] [CrossRef] [Google Scholar]
  51. Pribulla, T., Vaňko, M., Ammler-von Eiff, M., et al. 2012, Astron. Nachr., 333, 754 [NASA ADS] [CrossRef] [Google Scholar]
  52. Qian, S.-B., Dai, Z.-B., Liao, W.-P., et al. 2009, ApJ, 706, L96 [NASA ADS] [CrossRef] [Google Scholar]
  53. Qian, S.-B., Liao, W.-P., Zhu, L.-Y., & Dai, Z.-B. 2010a, ApJ, 708, L66 [Google Scholar]
  54. Qian, S.-B., Liao, W.-P., Zhu, L.-Y., et al. 2010b, MNRAS, 401, L34 [NASA ADS] [CrossRef] [Google Scholar]
  55. Qian, S.-B., Liu, L., Liao, W.-P., et al. 2011, MNRAS, 414, L16 [NASA ADS] [CrossRef] [Google Scholar]
  56. Qian, S.-B., Liu, L., Zhu, L.-Y., et al. 2012a, MNRAS, 422, 24. [Google Scholar]
  57. Qian, S.-B., Zhu, L.-Y., Dai, Z.-B., et al. 2012b, ApJ, 745, 23 [Google Scholar]
  58. Reipurth, B. 2000, AJ, 120, 3177 [CrossRef] [Google Scholar]
  59. Robertson, P., Endl, M., Cochran, W. D., et al. 2012a, ApJ, 749, 39 [NASA ADS] [CrossRef] [Google Scholar]
  60. Robertson, P., Horner, J., Wittenmyer, R. A., et al. 2012b, ApJ, 754, 50 [NASA ADS] [CrossRef] [Google Scholar]
  61. Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127 [NASA ADS] [CrossRef] [Google Scholar]
  62. Schwarz, R., Haghighipour, N., Eggl, S., Pilat-Lohinger, E., & Funk, B., 2011, MNRAS, 414, 2763 [NASA ADS] [CrossRef] [Google Scholar]
  63. Sybilski, P., Konacki, M., & Kozłowski, S. 2010, MNRAS, 405, 657 [NASA ADS] [Google Scholar]
  64. van den Berk, J., Portegies Zwart, S. F., & McMillan, S. L. W. 2007, MNRAS, 379, 111 [NASA ADS] [CrossRef] [Google Scholar]
  65. Watson, C. A., & Dhillon, V. S. 2004, MNRAS, 351, 110 [NASA ADS] [CrossRef] [Google Scholar]
  66. Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  67. Wittenmyer, R. A., Horner, J., Marshall, J. P., Butters, O. W., & Tinney, C. G. 2012a, MNRAS, 419, 3258 [NASA ADS] [CrossRef] [Google Scholar]
  68. Wittenmyer, R. A., Horner, J., Tuomi, M., et al. 2012b, ApJ, 753, 169 [NASA ADS] [CrossRef] [Google Scholar]
  69. Wittenmyer, R. A., Horner, J., & Marshall, J. P. 2013, MNRAS, 431, 2150 [NASA ADS] [CrossRef] [Google Scholar]
  70. Wittenmyer, R. A., Tan, X., Lee, M. H., et al. 2014, ApJ, 780, 140 [NASA ADS] [CrossRef] [Google Scholar]
  71. Yang, Y.-G., Li, H.-L., Dai, H.-F., & Zhang, L.-Y. 2010, AJ, 140, 1687 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Orbital and mass parameters of the two substellar companions proposed to orbit the RZ Dra Algol-type binary system, taken from Yang et al. (2010, their Table 5).

In the text
Table 2

Best-fit parameters for the LTT orbits of RZ Dra corresponding to Fig. 3.

In the text

All Figures

thumbnail Fig. 1

Dynamical stability of the RZ Dra system as a function of the semi-major axis, a, and eccentricity, e of RZ Dra (AB) D, for the solution presented in Yang et al. (2010). The nominal best-fit orbit for RZ Dra (AB) D presented in that work is located in the center of the red box, at a = 23.88 AU and e = 0.287. The 1σ uncertainties on the semi-major axis and eccentricity are shown by the red crosshairs radiating from the box. The lifetime at each of the 1681 a − e locations plotted in this figure (41 unique a and e values spanning ± 3σ from the nominal best-fit orbit) is the result of 75 distinct simulations, spanning 5 unique values of mean anomaly, M, and 15 unique values of the argument of periastron, ω. The mean lifetime can be seen to vary between ≃100 and ≃1600 years, demonstrating that the system as proposed in Yang et al. (2010) exhibits extreme dynamical instability. (Color version online.)
In the text
thumbnail Fig. 2

Results from numerical tests considering the orbits of the two companions as shown in Fig. 5. Here we consider an initial condition (IC 1) where the outer companion has been assigned an initial eccentricity of 0.90 (black dot in Fig. 5). All other parameters are as shown in Table 2. For the RADAU algorithm (top panels) we used an initial time step of 0.001 days. For the HYBRID algorithm (bottom panels) we used a constant time step of one day. For both algorithms we used an accuracy parameter of 10-12. Orbital elements are sampled every day and are plotted until the ejection time of RZ Dra(AB)C at around 22 years.
In the text
thumbnail Fig. 3

Best-fit two-Kepler LTT model without secular parameter with reduced χr2=12.7\hbox{$\chi_{\rm r}^2=12.7$}. The best-fit model parameters are shown in Table 2. In the top panel the LTT signal with smaller/larger amplitude corresponds to the inner/outer companion. The bottom panel shows the residuals between the data and best-fit model. We note that, judging by eye, there is a slight asymmetry in the residuals with O−C > 0 for E> 0.
In the text
thumbnail Fig. 4

Two-dimensional parameter scans of the χr2\hbox{$\chi^2_{\rm r}$} statistic in a neighborhood around the best-fit model (marked with an asterisk) considering the two-companion model. The color-bar is scaled with the value of our best-fit model representing 1. Countour lines show the 1-σ (left panel) and 1-, 2-, 3-σ (right panel) confidence level curves around the best-fit model. (Color version online.)
In the text
thumbnail Fig. 5

Dynamical stability of the RZ Dra system as a function of the semi-major axis, a, and eccentricity, e of RZ Dra (AB) D for the alternative solution. The nominal best-fit orbit for RZ Dra (AB) D is located in the center of the red box, at a = 26.6 AU and e = 0.62. The 1σ uncertainties on the semi-major axis and eccentricity are shown by a crosshair. Almost all the tested models, particularly those within the ± 1σ uncertainties, are dynamically unstable on remarkably short timescales. A few regions are stable on timescales of around a million years, but these are all located at the extreme edges of the plot, almost ± 3σ away from the nominal best-fit orbit. The black dot is the location of a test orbit as discussed in Sect. 2. This figure is available in color in electronic form.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.