3.1. Light travel time effect

The $O-C$ diagram (see Figs. 2 and 3) obtained from residuals of the linear fit indicates a cyclic variation, which can be attributed to the LTT effect. To investigate the orbital period variation, we used the models including a quadratic term ( $\unicode{x03B2}=P\dot{P}/2$ ) and/or the LTT term(s) causing from the presence of hypothetical body. The quadratic term is included to the model by adding $\unicode{x03B2} L^2$ . The LTT term is defined by Irwin (Reference Irwin1952), a modified version provided by Goździewski et al. (Reference Goździewski2012) as following;

(2) \begin{equation} \tau_{i}= K_i \left(\sin{\omega_i} (\cos{E_i (t)} - e_i) + \sqrt{1-e_i^2} \cos{\omega_i} \sin{E_i (t)}\right)\end{equation}

This formulation is parameterized by Keplerian orbital elements of the N-body companions orbiting around the mass centre of the system. In Equation (2), $K_i$ is the semi-amplitude of the LTT signal of the i th body, $e_i$ is the eccentricity, $\omega_i$ is the longitude of pericentre, $E_i$ is the eccentric anomaly, $t_{0,i}$ is time of pericentre passage. To prevent weakly constrained values for $e_i$ and $\omega_i$ in quasi-circular and moderately eccentric orbits, we utilized Poincare elements. These elements are represented by $x\equiv e_i cos\omega_i$ and $y\equiv e_i sin\omega_i$ (see Goździewski et al. Reference Goździewski2012, Reference Goździewski2015; Nasiroglu et al. Reference Nasiroglu2017, for more details); Özdönmez, Er, & Nasiroglu Reference Özdönmez, Er and Nasiroglu2023.

Figure 2.

The $O-C$ diagram of HS 0705+6700 is displayed in the upper plot using a linear ephemeris calculated from Equation (1). The TUG T100 is represented by the red-filled circles, while the ATA50 is represented by the orange-filled circles. Data from literature and TESS are shown in different colors and labeled with corresponding abbreviations. The model obtained in our study is represented by the black line, which includes the quadratic term and two additional objects. The grey shaded area represents the $\pm3\unicode{x03C3}$ posterior spread. This is calculated from 1 000 randomly selected parameter samples from the MCMC posterior. The bottom plot shows the residuals of the $O-C$ times for the model obtained in this study. The figure also includes the calculated RMS value of the residuals of the mid-eclipse times for the model.

Figure 3.

The upper plot displays the $O-C$ diagram of HS 0705+6700 using a linear ephemeris calculated from Equation (1). The model including the three LTT term is represented by the black line. Other informations are as in Fig. 2.

This study aims to explain the $O-C$ diagram for each model using various formulations, including those with/without quadratic terms, and with one to three LTT terms. For example, the model containing a quadratic term and two LTT terms are formulated as follows.

(3) \begin{equation} T_{eph}(L) = T_0 + L \times P_{bin} + \unicode{x03B2} L^2 + \tau_{1}(L) + \tau_{2}(L)\end{equation}

Markov Chain Monte Carlo (MCMC) methodology based on a likelihood function ( $\mathcal{L}$ ) was used to express the orbital period variation, following the identical fitting process described in our previous studies (see Nasiroglu et al. Reference Nasiroglu2017; Er et al. Reference Er, Özdönmez and Nasiroglu2021; Özdönmez et al. Reference Özdönmez, Er and Nasiroglu2023) Uniform prior samples have been randomly assigned to all free parameters within the specified ranges $\unicode{x03B2}, K_i, P_i, t_{0,i}, \unicode{x03C3}_f > 0$ days, $x_i,y_i \epsilon$ [–0.75,+0.75], $P_{bin} \epsilon$ [0.08, 0.15] days and $\Delta T_0 \epsilon$ [–10, +10]. The $\mathcal{L}$ function includes the free parameter $\unicode{x03C3}_f$ in units of days to account for systematic uncertainties. This parameter scales the raw uncertainties of eclipsing times ( $\unicode{x03C3}_i$ ) in quadrature. The MCMC method is used to sample the posterior distribution. The sampling process utilized the affine-invariant ensemble sampler implementation from the emcee package, following the approach presented by Goodman & Weare (Reference Goodman and Weare2010), and made available by Foreman-Mackey et al. (Reference Foreman-Mackey, Hogg, Lang and Goodman2013). For MCMC, 512 initial conditions (walkers) were used to observe the dynamics of each distinct variable in models over 30 000–120 000 steps (depending on used models) within chains. The optimal parameter values, along with their corresponding uncertainties, were determined by assessing the 16th, 50th, and 84th percentiles of the marginalized distributions derived from the maximized likelihood (L). The MCMC has been run separately for each models with a variety of formulations. The models, which include only three LTT terms or one quadratic term and two LTT terms, are the most statistically and astrophysically consistent with the $O-C$ diagram. Table 2 presents the most plausible parameters for these models. The best-fitting parameters of the models are shown in the $O-C$ diagram in Figs. 2 and 3. Figs. A1–A2 consists of the 1D and 2D posterior probability distributions of the system parameters sampled by MCMC. It is important to acknowledge that single-body models aren’t consistent with the most recent $O-C$ diagram using all available data. We have discussed the results obtained for all models in Section 4.

Table 2.

System Parameters of HS 0705+6700 for two models.

The minimum masses of additional bodies can be determined from the following mass function,

(4) \begin{equation}f(M_i)=\frac{(M_i\sin{i_i})^3}{(M_{i}+M_{bin})^2}=\frac{4\pi^2(a_{12}\sin{i_i})^3}{G P_i^2}\end{equation}

where G is the gravitational constant, $M_{bin}$ is the total mass of binary, $i_i$ is the inclination of the ith body’s orbit, $a_{12}\sin{i_i}$ is the projected semi-major axis of the binary system around the barycentre of the system, $P_i$ is the orbital period, and $M_i$ is the mass of the ith body’s. We used of $M_{bin}=\:\sim0.617\: M_\odot$ reported by Drechsel et al. (Reference Drechsel2001) for the stellar binary mass.

3.2. Applegate mechanism

The orbital period variation of binary star systems can also be attributed to the magnetic cycle of the low-mass active stars in the system. This is known as the Applegate mechanism (Applegate Reference Applegate1992). Changes in the shape of a magnetically active component can contribute to the orbital period variation of the system. To test the magnetic mechanism on the $O-C$ signals, we calculated the energy ratios ( $\Delta E/E_{sec}$ ) through three different Applegate approaches, as follows: Thin-shell model (Tian, Xiang, & Tao Reference Tian, Xiang and Tao2009), Finiteshell two-zone model (Völschow et al. Reference Völschow, Schleicher, Perdelwitz and Banerjee2016), and Spin–orbit coupling model model (Lanza Reference Lanza2020). For all Applegate models, the magnetic activity could occur under the condition that $\Delta E/E_{sec}$ is less than 1 (i.e. $\Delta E \ll E_{sec}$ ). For the calculations of the energy ratios, the parameters obtained in Table 2 and determined by Völschow et al. (Reference Völschow, Schleicher, Perdelwitz and Banerjee2016) were used. Those calculated for the smaller LTT signal are listed in Table 3. For the LTT signal of the fifth body in the model with three LTT terms, the $\Delta E/E_{sec}$ is calculated to be 0.66 using the thin shell model (Tian et al. Reference Tian, Xiang and Tao2009), depending on solar-like magnetic cycles in the secondary star. The other calculated energy ratios are much higher than the threshold. Thus, only the LTT signal of the fifth body can be attributed to the magnetic cycle in the case of the thin-shell model.

3.3. Orbit stability analysis

To investigate the orbital stability of HS 0705+6700 in our models, we used the N-body orbital integration package of the REBOUND Footnote c (Rein & Liu, Reference Rein and Liu2012), which includes a Mean Exponential Growth factor of Nearby Orbits (MEGNO, (Cincotta & Simó, Reference Cincotta and Simó2000)) indicator and a Wisdom-Holman symplectic integrator (WHFAST, (Rein & Tamayo, Reference Rein and Tamayo2015)). Using N-body integration, REBOUND simulates the motion of celestial objects and provides two significant insights: First, the MEGNO chaotic parameter surface is mapped, yielding an indicator $\lt Y \gt$ that assesses the chaotic behaviour of the system over a range of semi-major axis and eccentricity values over a given period of time. A stable system is indicated by $\lt Y\gt \leq2$ , while values greater than 2 indicate chaotic (unstable) orbital configurations. A value of 10 is assigned to $\lt Y\gt$ when a particle is ejected or collides. Second, the orbital stability timeline integrates the orbits for a given time and shows the variations in parameters such as semi-major axis and eccentricity as a function of time. This is useful for understanding planetary interactions, predicting system escape or collision, and determining the stability period of orbits.

Table 3.

The energy ratios ( $\Delta E$ / $E_{sec}$ ) for the formulation of corresponding Applegate mechanisms (Tian et al. Reference Tian, Xiang and Tao2009; Völschow et al. Reference Völschow, Schleicher, Perdelwitz and Banerjee2016; Lanza Reference Lanza2020).

In both simulation scenarios, the central binary star was treated as a singular mass, and all orbital trajectories were confined to a co-planar configuration. We set the optimal timestep for WHFast to be roughly 0.1% of the shortest orbital period of the additional bodies. It was also assumed that the limit distance for escaping from the system is 20 AU. Dynamic stability simulations were performed using the model parameters to obtain both the MEGNO value and the orbital stability timeline. It has been found that all system configurations constructed from the system parameters of the models in Table 2 are unstable even <2 000 yr. In addition, the stability tests were performed under the assumption that the detected signal of the fifth body was raised from the magnetic cycle, yet the stable system configuration on longer time scales can not be constructed.