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/00046361/201423799  
Published online  19 May 2014 
On the RZ Draconis substellar circumbinary companions
Stability study of the proposed substellar circumbinary system
^{1}
Korea Astronomy and Space Science Institute,
776 Daedukdaero, Yuseonggu,
305348
Daejeon,
Republic of Korea
email:
tchinse@gmail.com
^{2}
Computational Engineering and Science Research Centre, University
of Southern Queensland, Toowoomba, 4350
Queensland,
Australia
^{3}
Australian Centre for Astrobiology, University of New South Wales,
Sydney
2052,
Australia
^{4}
School of Physics, University of New South Wales,
Sydney
2052,
Australia
Received:
13
March
2014
Accepted:
12
April
2014
Context. Recently, using the lighttravel time effect, planets and substellar companions have been proposed to orbit around binary star systems (also known as circumbinary companions) as a result of variations in timing of the observed eclipses. For the majority of these systems the proposed orbital architecture features a crossing of the orbital configurations as a result of high eccentricities for one or both of the companions. For such systems, strong mutual gravitational interactions are expected, resulting in catastrophic orbital instabilities, or collisions between the proposed components, on very short timescales.
Aims. We reexamine the primary and secondary eclipse timings of the shortperiod and semidetached binary RZ Draconis (RZ Dra). The proposed companions were reported to have masses of around ≃0.07 and ≃0.18 M_{⊙} with the inner companion on an orbit with moderate eccentricity (0.46), whose apocenter distance crosses the orbit of the outer companion. We show that the companions proposed previously follow highly unstable orbits. In an attempt to find a stable system we searched the underlying χ^{2} parameter space for a bestfit model and carried out an orbitstability study to test possible bestfit models. If the binary period changes are truly due to additional massive companions in a hierarchical configuration, they must follow stable orbits.
Methods. For numerical orbital stability calculations we used well established orbit integration routines. Computations were carried out using a largescale multiCPU computing environment. Our data analysis of times of primary and secondary eclipse is based on the LevenbergMarquardt leastsquares minimization algorithm using the twobody Keplerian lighttravel time effect model.
Results. Despite the wide variety of potential models tested for the RZ Dra system in this work, we found very few models that were stable for even one million years, with the vast majority of systems tested falling apart on timescales of just hundreds of years. It seems therefore likely that the observed timing variations are not solely the result of massive, unseen companions.
Key words: methods: observational / methods: numerical / celestial mechanics / binaries: eclipsing / stars: formation / stars: individual: RZ Draconis
© 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 erosionprocess of a nonhierarchical star system by angular momentum and energy exchange via mutual gravitational interactions. In the latter case, and considering an initial triplesystem, 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 quasiKeplerian orbit.
A particular example of a hierarchical multibody system is a socalled circumbinary system (also known as companions on Ptype 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 community^{1} (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 PH1 (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 socalled Stype (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 mideclipse 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 substellarmass circumbinary companion discovered from that data which orbits KIC002856960 (Lee et al. 2013a). The fundamental principle of the lighttravel time effect (LTT^{2}) 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 timeeffective, involving only photometric CCD measurements.
In recent times, single and multibody substellar circumbinary companions to known eclipsing binary systems have been proposed using groundbased 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 postmain 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 twocompanion system is highly unstable.
Finally, Parsons et al. (2010) presented photometric followup observations of a number of eclipsing postcommonenvelope binaries for which they have been able to rule out previous claims for singleobject circumbinary companions (e.g., Qian et al. 2009, 2010b).
In this work we reexamine the observed timing dataset of RZ Dra as presented by Yang et al. (2010). These authors proposed the existence of two additional lowmass dwarfs from two distinct quasisinusoidal 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 lighttravel time effect and describe the leastsquares methodology in Sect. 4, where we present a new bestfit 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
Orbital and mass parameters of the two substellar companions proposed to orbit the RZ Dra Algoltype 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 Algoltype binary star system RZ Dra might be explained by the presence of two massive, unseen companions moving on eccentric longperiod orbits, as detailed in Table 1. On their nominal bestfit orbits, these two proposed objects can approach one another remarkably closely – with the innermost object (RZ Dra (AB) C)^{3} having a bestfit apastron distance of 18.06 AU, and the outermost (RZ Dra (AB) D) having a nominal bestfit periastron distance of 17.02 AU. In other words, the nominal bestfit models for the two proposed companions have significant overlap – these stellarmass companions can cross one another’s orbits.
Fig. 1 Dynamical stability of the RZ Dra system as a function of the semimajor axis, a, and eccentricity, e of RZ Dra (AB) D, for the solution presented in Yang et al. (2010). The nominal bestfit 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 semimajor 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 bestfit 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 nbody dynamics package mercury (Chambers 1999). We used a constant integration timestep of one day in all our orbit integrations. The error tolerance parameter was set to one part in 10^{12}, which ensures accurate integrations (using the BulirschStoer algorithm) of higheccentricity 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 bestfit values (as detailed in Table 1). We then systematically varied the initial semimajor 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 breakup 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. Fortyone 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.
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 longestlived 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 data^{4} 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 Ephemeris^{5} (Giorgini et al. 1996) and hence were transformed by omitting the lighttravel 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) mideclipse events can be computed from (1)where E denotes the cycle number, T_{0} is the reference epoch of some primary eclipse, and P_{0} measures the eclipsing period of RZ Dra. A linear regression on the 680 recorded eclipse times allows determining P_{0} 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 T_{0} and P_{0} 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 lighttravel 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 mideclipse 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 lighttravel time effect and is given by (2)where T_{O} denotes the measured time of an observed mideclipse, 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 twobody LTT orbits. The quantity τ_{i} is given by the following expression for each companion (Irwin 1952): (3)where K_{b,i} = a_{b,i}sinI_{b,i}/c is the semiamplitude of the lighttime effect (in the O−C diagram) with c measuring the speed of light and I_{b,i} is the lineofsight inclination of the LTT orbit relative to the sky plane, e_{b,i} the orbital eccentricity, f_{b,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 (a_{b,i}sinI_{b,i},e_{b,i},ω_{b,i},T_{b,i},P_{b,i}). The time of pericenter passage T_{b,i} and orbital period P_{b,i} are introduced by expressing the true longitude as a timelike variable via the mean anomaly M = n_{b,i}(T_{O} − T_{b,i}), with n_{b,i} = 2π/P_{b,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. Leastsquares fitting and parameter uncertainties
Fig. 3 Bestfit twoKepler LTT model without secular parameter with reduced . The bestfit 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 bestfit 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 bestfit model follows the methodology as outlined in Hinse et al. (2012a) using a leastsquares minimization algorithm as implemented in MPFIT (Markwardt 2009). We considered several models. First we started to fit a singlecompanion 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 bestfit statistic to be no less than 10^{4} which indicates a onecompanion model that is inadequate to describe the data. We then fitted for a twocompanion model with all twelve parameters allowed to vary freely (680−12 degrees of freedom). A significantly better fit was found with . 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 bestfit 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 onecompanion model with a quadratic trend we refer to Erdem et al. (2011). Our bestfit twocompanion model and resulting bestfit 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 bestfit covariance 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 bestfit model as our initial guess and recorded the resulting parameters. We generated 10^{6} 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 twodimensional jointconfidence parameter scan of the 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 bestfit model
Fig. 4 Twodimensional parameter scans of the statistic in a neighborhood around the bestfit model (marked with an asterisk) considering the twocompanion model. The colorbar is scaled with the value of our bestfit model representing 1. Countour lines show the 1σ (left panel) and 1, 2, 3σ (right panel) confidence level curves around the bestfit 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 bestfit values and created 126 075 test systems in which the initial semimajor axis, eccentricity, mean anomaly, and argument of periastron for RZ Dra (AB) D were systematically varied across the ± 3σ range around their nominal bestfit 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 bestfit 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 nearcontact binaries is not well known. Currently, there is poorly consensus that nearcontact 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 masstransfer is possible via a flare that moves matter from one component to the other. This could cause a sudden periodchange 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 bestfit model, the inner companion apocenter distance is at ≃18 AU and the outer companion pericenter distance is ≃17 AU, implying a crossingorbit architecture. We tested the dynamical feasibility of the proposed bestfit orbits by exploring a large grid of initial conditions using direct nbody 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 breakup 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 anticoplanar 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 bestfit model using the complete set of timing measurements transformed from HJD to BJD time standard. We found a new bestfit 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).
Fig. 5 Dynamical stability of the RZ Dra system as a function of the semimajor axis, a, and eccentricity, e of RZ Dra (AB) D for the alternative solution. The nominal bestfit 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 semimajor 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 bestfit 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 twocompanion 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 quasisinusoidal timing variations could be produced by one or more physical causes: (1) apsidal motion in an elliptical orbit of the binary; (2) a lighttravel 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 10^{51}−10^{52} 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 latetype secondary.
However, the eclipsing pair of RZ Dra is a semidetached binary whose less massive secondary fills its inner Roche lobe, see Yang et al. (2010). In such semidetached 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 + twoLTT 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 + oneLTT 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 noiselevel 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 followup 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.
sometimes also referred to as LITE (Borkovits & Hegedüs 1996).
We adopt the naming convention as proposed by Hessman et al. (2010).
We note that new timing data exist published by Erdem et al. (2011). These data points were not included in the analysis presented here.
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 HighEnd Computing (ICHEC). T.C.H., J.W.L., C.U.L. & J.H.P. acknowledge support from KASI registered under grant number 2013940000/2014140006. Astronomical research at Armagh Observatory is funded by the Department of Culture, Arts and Leisure (DCAL). JPM is supported by Spanish grant AYA 201126202. 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
 Almeida, L. A., & Jablonski, F. 2011, IAU Symp., 276, 495 [NASA ADS] [Google Scholar]
 Almeida, L. A., Jablonski, F., & Rodrigues, C. V. 2013, ApJ, 766, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Beuermann, K., Hessman, F. V., Dreizler, S., et al. 2010, A&A, 521, L60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuermann, K., Buhlmann, J., Diese, J., et al. 2011, A&A, 526, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuermann, K., Dreizler, S., & Hessman, F. V. 2013, A&A, 555, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borkovits, T., & Hegedüs, T., 1996, A&AS, 120, 63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chambers, J. E., 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chauvin, G., Lagrange, A.M., Udry, S., & Mayor, M. 2007, A&A, 475, 723 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deeg, H. J., Doyle, L. R., Kozhevnikov, V. P., et al. 2000, A&A, 358, 5 [Google Scholar]
 Doyle, L. R., & Deeg, H.J. 2004, IAU Symp., 213, 80 [Google Scholar]
 Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Duchêne, G., Bontemps, S., Bouvier, J., et al. 2007, A&A, 476, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935 [NASA ADS] [CrossRef] [Google Scholar]
 Erdem, A., Zola, S., & Winiarski, M. 2011, New Astron., 16, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, D. S., 1968, QJRAS, 9, 388 [NASA ADS] [Google Scholar]
 Funk, B., Eggl, S., Gyergyovits, M., et al. 2011, EPSCDPS, 1725 [Google Scholar]
 Funk, B., Eggl, S., Gyergyovits, M., et al. 2012, IAU Symp., 282, 446 [NASA ADS] [Google Scholar]
 Furlan, E., Sargent, B., Calvet, N., et al. 2007, ApJ, 664, 1176 [NASA ADS] [CrossRef] [Google Scholar]
 Giorgini, J. D., Yeomans, D. K., Chamberlin, A. B., et al., 1996, BAAS, 28, 1158 [NASA ADS] [Google Scholar]
 Goździewski, K., Nasiroglu, I., Słowikowska, A., et al. 2012, MNRAS, 425, 930 [NASA ADS] [CrossRef] [Google Scholar]
 Hessman, F. V., Dhillon, V. S., Winget, D. E., et al. 2010, unpublished [arXiv:1012.0707H] [Google Scholar]
 Hinse, T. C., Lee, J. W., Haghighipour, N., et al. 2012a, MNRAS, 420, 3609 [NASA ADS] [CrossRef] [Google Scholar]
 Hinse, T. C., Goździewski, K., Lee, J. W., et al. 2012b, AJ, 144, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Hinse, T. C., Lee, J. W., Goździewski, K., et al. 2014, MNRAS, 438, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Marshall, J. P., Wittenmyer, R. A., & Tinney, C. G. 2011, MNRAS, 416, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Wittenmyer, R. A., Hinse, T. C., & Tinney, C. G. 2012a, MNRAS, 425, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Hinse, T. C., Wittenmyer, R. A., Marshall, J. P., & Tinney, C. G. 2012b, MNRAS, 427, 2812 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Wittenmyer, R. A., Hinse, T. C., et al. 2013, MNRAS, 435 2033 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Wittenmyer, R. A., Hinse, T. C., & Marshall, J. P. 2014, MNRAS, 439, 1176 [NASA ADS] [CrossRef] [Google Scholar]
 Irwin, J. B. 1952, ApJ, 116, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Irwin, J. B. 1959, AJ, 64, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Kostov, V. B., McCullough, P., Hinse, T. C., et al. 2013, ApJ, 770, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Lanza, A. F. 2006, MNRAS, 369, 1773 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. W., Kim, S.L., Kim, C.H., et al. 2009a, AJ, 137, 3181 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. W., Youn, J.H., Lee, C.U., et al. 2009b, AJ, 138, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. W., Kim, S.L., Lee, C.U., et al. 2011, AJ, 142, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. W., Lee, C.U., Kim, S.L., et al. 2012, AJ, 143, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. W., Kim, S.L., Lee, C.U., et al. 2013a, ApJ, 763, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. W., Hinse, T. C., & Park, J.H. 2013b, AJ, 145, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Lim, J., & Takakuwa, S. 2006, ApJ, 653, 425 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Marshall, J., Horner, J., & Carter, A. 2010, International Journal of Astrobiology, 9, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mustill, A., Marshall, J. P., Villaver, E., et al. 2013, MNRAS, 437, 1404 [Google Scholar]
 Neuhäuser, R., Mugrauer, M., Fukagawa, M., Torres, G., & Schmidt, T. 2007, A&A, 462, 777 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ofir, A., Deeg, H. J., & Lacy, C. H. S. 2009, A&A, 506, 445 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, Science, 337, 1511 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, ApJ, 758, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Parsons, S. G., Marsh, T. R., Copperwheat, C. M., et al. 2010, MNRAS, 407, 2362 [NASA ADS] [CrossRef] [Google Scholar]
 Potter, S. B., RomeoColmenero, E., Ramsay, G., et al. 2011, MNRAS, 416, 2202 [NASA ADS] [CrossRef] [Google Scholar]
 Pribulla, T., Vaňko, M., Ammlervon Eiff, M., et al. 2012, Astron. Nachr., 333, 754 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, S.B., Dai, Z.B., Liao, W.P., et al. 2009, ApJ, 706, L96 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, S.B., Liao, W.P., Zhu, L.Y., & Dai, Z.B. 2010a, ApJ, 708, L66 [Google Scholar]
 Qian, S.B., Liao, W.P., Zhu, L.Y., et al. 2010b, MNRAS, 401, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, S.B., Liu, L., Liao, W.P., et al. 2011, MNRAS, 414, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, S.B., Liu, L., Zhu, L.Y., et al. 2012a, MNRAS, 422, 24. [Google Scholar]
 Qian, S.B., Zhu, L.Y., Dai, Z.B., et al. 2012b, ApJ, 745, 23 [Google Scholar]
 Reipurth, B. 2000, AJ, 120, 3177 [CrossRef] [Google Scholar]
 Robertson, P., Endl, M., Cochran, W. D., et al. 2012a, ApJ, 749, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, P., Horner, J., Wittenmyer, R. A., et al. 2012b, ApJ, 754, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Schwarz, R., Haghighipour, N., Eggl, S., PilatLohinger, E., & Funk, B., 2011, MNRAS, 414, 2763 [NASA ADS] [CrossRef] [Google Scholar]
 Sybilski, P., Konacki, M., & Kozłowski, S. 2010, MNRAS, 405, 657 [NASA ADS] [Google Scholar]
 van den Berk, J., Portegies Zwart, S. F., & McMillan, S. L. W. 2007, MNRAS, 379, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, C. A., & Dhillon, V. S. 2004, MNRAS, 351, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wittenmyer, R. A., Horner, J., Marshall, J. P., Butters, O. W., & Tinney, C. G. 2012a, MNRAS, 419, 3258 [NASA ADS] [CrossRef] [Google Scholar]
 Wittenmyer, R. A., Horner, J., Tuomi, M., et al. 2012b, ApJ, 753, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Wittenmyer, R. A., Horner, J., & Marshall, J. P. 2013, MNRAS, 431, 2150 [NASA ADS] [CrossRef] [Google Scholar]
 Wittenmyer, R. A., Tan, X., Lee, M. H., et al. 2014, ApJ, 780, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, Y.G., Li, H.L., Dai, H.F., & Zhang, L.Y. 2010, AJ, 140, 1687 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Orbital and mass parameters of the two substellar companions proposed to orbit the RZ Dra Algoltype binary system, taken from Yang et al. (2010, their Table 5).
All Figures
Fig. 1 Dynamical stability of the RZ Dra system as a function of the semimajor axis, a, and eccentricity, e of RZ Dra (AB) D, for the solution presented in Yang et al. (2010). The nominal bestfit 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 semimajor 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 bestfit 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 
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 
Fig. 3 Bestfit twoKepler LTT model without secular parameter with reduced . The bestfit 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 bestfit 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 
Fig. 4 Twodimensional parameter scans of the statistic in a neighborhood around the bestfit model (marked with an asterisk) considering the twocompanion model. The colorbar is scaled with the value of our bestfit model representing 1. Countour lines show the 1σ (left panel) and 1, 2, 3σ (right panel) confidence level curves around the bestfit model. (Color version online.) 

In the text 
Fig. 5 Dynamical stability of the RZ Dra system as a function of the semimajor axis, a, and eccentricity, e of RZ Dra (AB) D for the alternative solution. The nominal bestfit 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 semimajor 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 bestfit 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.