Orbital period variation analysis of the HS 0705+6700 post-common envelope binary

Abstract

To detect additional bodies in binary systems, we performed a potent approach of orbital period variation analysis. In this work, we present 90 new mid-eclipse times of a short-period eclipsing binary system. Observations were made using two telescopes from 2014 to 2024, extending the time span of the $O-C$ diagram to 24 yr. The data obtained in the last seven years indicate significant deviations in the $O-C$ diagram from the models obtained in previous studies. We investigated whether this variation could be explained by mechanisms such as the LTT effect or Applegate. To investigate the cyclic behaviour observed in the system with the light travel time effect, we modelled the updated $O-C$ diagram using different models including linear/quadratic terms and additional bodies. The updated $O-C$ diagram is statistically consistent with the most plausible solutions of models that include multiple brown dwarfs close to each other. However, it has been found that the orbit of the system is unstable on short time scales. Using three different theoretical definitions, we have found that the Applegate mechanism cannot explain the variation in the orbital period except for the model containing the fifth body. Therefore, due to the complex nature of the system, further mid-eclipse time is required before any conclusions can be drawn about the existence of additional bodies.

1. Introduction

Celestial systems consisting of binary stars are a popular topic due to their diverse and complex behaviour. The examination of close binary systems, especially those involving subdwarf B stars (sdBs) or white dwarfs, has provided a comprehensive understanding of stellar evolution, mass transfer, and gravitational interactions within such systems. Close binary systems of the HW Vir type consisting of subdwarf B (sdB) stars and low-mass components such as white dwarfs or M dwarfs (dM) typically represent a subset with intriguing properties Maxted et al. 2001; Han et al. 2002, 2003; Morales-Rueda et al. 2003; Kilkenny 2011; Schleicher & Dreizler 2014; Silvotti et al. 2014). The binary systems known as post-common envelope binary (PCEB) provide an important opportunity for understanding various astrophysical phenomena such as the light travel time (LTT) effect that is demonstrated by variations in eclipse timings within binary star systems (Beuermann et al. 2012b; Horner et al. 2012; Almeida, Jablonski, & Rodrigues 2013; Lohr et al. 2014; Marsh 2018). Thus, LTT effect can be used to investigate the possible presence of additional object(s) orbiting close binary systems.

The possible existence of the additional bodies orbiting around a short-period ( $P=2.3$ h) eclipsing PCEB system HS 0705+6700 (V470 Cam, $V=14.7$ mag), consisting of an sdB star and a dM component (Drechsel et al. 2001), has been extensively studied by many researchers. The possible presence of a third body orbiting HS 0705+6700 has been reported (Qian et al. 2009, 2010; Çamurdan, Zengin Çamurdan, & İbanoğlu 2012; Beuermann et al. 2012b; Bogensberger, Clarke, & Lynas-Gray 2017). Qian et al. (2013) and Pulley et al. (2015) detected a positive increase in the orbital period of the system, suggesting that this increase could be due to the presence of a fourth body. Sale et al. (2020) indicated the presence of two circumbinary brown dwarfs orbiting the system with a two-body model containing a quadratic term. However, the proposed model is dynamically unstable over a timescale of $10^3$ yr. Finally, Mai & Mutel (2022) investigated the orbital period variation of the system for both one- and two-body models using different data set rather than those commonly used in literature. Although the two-body model remained stable over a timescale of $10^7$ yr, statistically, the one-body model yielded better results than the two-body model. Pulley et al. (2022) concluded that none of the models in the literature are consistent with the most recent $O-C$ .

This study contributes observationally to new mid-eclipse times to a comprehensive investigation of the mechanisms underlying orbital period variations in HS 0705+6700. This constrains the parameters of potential additional objects and improves our understanding of the structure of the system, which is known to be complex, and its orbital stability.

2. Observations and data analysis

We conducted an observation campaign for HS 0705+6700 system between Nov 2014 and Jan 2024 using the 1 m telescope equipped with a 4k $\times$ 4k SI1100 CCD, the $15\times15$ $\unicode{x03BC}$ m pixel size at the TÜBİTAK National Observatory (TUG T100, Antalya, Türkiye) and 50 cm telescope equipped with Apogee Alta U230 2K CCD with $24\times24$ $\unicode{x03BC}$ m pixel size at the Türkiye National Observatories (ATA50, Erzurum, Türkiye). It should be noted that the CCD at the ATA50 telescope was replaced with a CMOS QHY268M Pro I camera (a $3.76\times3.76$ $\unicode{x03BC}$ m pixel size) after April 4, 2022. Our observations were performed in white light to obtain optimal counts, with exposure times ranging from 3 to 25 seconds depending on the seeing. A standard process was used to reduce the CCD frames, i.e. bias subtraction, flat fielding, dark subtraction, and cosmic ray correction. The reduced CCD frames were performed the differential aperture photometry using the same method as in Er, Özdönmez, & Nasiroglu (2021). Thus, we obtained the 90 new primary eclipse light curves of the HS 0705+6700.

The system was observed by TESS in sectors 20, 47 and 60 from Dec 2019 to Jan 2023. For each of the three sectors, the sampling time of an image was between 20 and 120 seconds. To download the photometric images, the Lightkurve packagea (Lightkurve Collaboration et al. 2018) was used, which provides the ability to download TESS data from the public data archive at Barbara A. Mikulski Archive for Space Telescopes (MASTb). From TESS photometric images in all sectors, we obtained 1093 primary eclipse light curves.

To determine the mid-eclipse time, we modelled the primary eclipse light curves obtained from our observations and TESS data with a modified Gaussian profile as in Beuermann et al. (2012b). Modelled eclipse light curves obtained from our observations are shown in Fig. 1. For each of our modeled eclipse light curves, we calculated the root mean square (RMS) from the residuals between observed and modelled light curves. The RMS values range from 0.006 to 0.087 mag with a mean value of 0.017 mag.

Figure 1. The light curves of HS 0705+6700 during the eclipse were observed using the TUG T100 and ATA50 telescopes. A modified Gaussian profile was used to fit the light curves, as explained in Section 2. The date of observation and telescope used are labeled, and the sequence of total eclipses on the same day is indicated with Roman numerals. The calculated RMS values, in units of magnitude, are also provided.

3. Orbital period variations

We converted all mid-eclipse times to barycentric dynamical Julian time (BJD) using the method described in Eastman, Siverd, & Gaudi (2010). Table 1 lists the mid-eclipse times collected from the literature (Drechsel et al. 2001; Niarchos, Gazeas, & Manimanis 2003; Németh, Kiss, & Sarneczky 2005; Kruspe, Schuh, & Traulsen 2007; Qian et al. 2009, 2010; Çamurdan et al. 2012; Beuermann et al. 2012b; Diethelm 2012; Qian et al. 2013; Diethelm 2013; Kubicki 2015; Petropoulou et al. 2015; Bogensberger et al. 2017; Pulley et al. 2018; Faillance et al. 2020; Sale et al. 2020; Mai & Mutel 2022; Pulley et al. 2022), as well as those obtained from our observations and TESS data. Before utilizing this data, we made the subsequent modifications: (i) We excluded the outlier mid-eclipse times which scatter more than three standard deviations from the overall ( $O-C$ ) trend that is calculated every 5 000 cycles. (ii) Since there are no errors in the times and the start times of the exposures were published for the mid-eclipse times in Bogensberger et al. (2017), we used the mid-eclipse times and errors given in Sale et al. (2020) for these mid-eclipse times. (iii) The mid-eclipse times based on data obtained from space-based telescopes are widely scattered in the $O-C$ diagram due to their imprecise nature. Thus, The TESS data was binned into five groups based on their cycles, and only these binned times were used during modelling.

Table 1. The mid-eclipse times of HS 0705+6700, its error, references

$^{*}$ Full table is available in its entirety in machine-readable form.

We fit the mid-eclipse times by following linear ephemeris;

(1) \begin{equation}\begin{aligned}\ T_{eph}(L) = & T_0 + L \times P_{bin} \\ = & \text{BJD}\: 2451822.758177(54) + L \times 0.0956467117(8)\end{aligned}\end{equation}

In the Equation (1), $T_{eph}$ represents the mid-eclipse time, while $T_0$ , L, and $P_{bin}$ are the initial ephemeris, the cycle and the orbital period of the binary system, respectively.

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 (1952), a modified version provided by Goździewski et al. (2012) 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. 2012, 2015; Nasiroglu et al. 2017, for more details); Özdönmez, Er, & Nasiroglu 2023.

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. 2017; Er et al. 2021; Özdönmez et al. 2023) 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 (2010), and made available by Foreman-Mackey et al. (2013). 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. (2001) 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 1992). 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 2009), Finiteshell two-zone model (Völschow et al. 2016), and Spin–orbit coupling model model (Lanza 2020). 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. (2016) 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. 2009), 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 c (Rein & Liu, 2012), which includes a Mean Exponential Growth factor of Nearby Orbits (MEGNO, (Cincotta & Simó, 2000)) indicator and a Wisdom-Holman symplectic integrator (WHFAST, (Rein & Tamayo, 2015)). 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. 2009; Völschow et al. 2016; Lanza 2020).

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.

4. Discussion and conclusions

We present 90 new primary mid-eclipse times for HS 0705+6700 from 2014 to 2024. By combining our mid-eclipse times with those obtained from TESS data in this study and from the literature, we analysed the detected orbital period variation in the derived $O-C$ diagram. Our data covers a time span of 10 yr, which extends the time span of the $O-C$ diagram by about 2 yr, over a total of 24 yr.

All possible models were used to test the new $O-C$ diagram. These include models with quadratic/nonquadratic terms and with one to three LTT terms. Studies in the literature explaining the orbital period variation of the HS 0705+6700 with models containing only a single body are based on $O-C$ diagrams prior to 2017. However, Mai & Mutel (2022) reported a single body model consistent with their $O-C$ diagram, without using the eclipse times around 2004 (cycle 12 000) obtained by Németh et al. (2005). We could not find a valid reason, such as large uncertainty, to exclude all of these times. We could not find any plausible model that includes only one LTT term explaining the orbital variability of the most recent $O-C$ diagram using all available data. The statistical coherence of the $O-C$ diagram is maximized when using a model including more than one LTT terms. Exceptionally, the model including only two LTT signals resulted in very high semi-amplitudes, implying M-type stars with $M_3= 138\: M_{\text{Jup}}$ and $M_4= 151\: M_{\text{Jup}}$ . It is astrophysically uncommon for a stable quadruple system to contain stars so close together ( $a<5$ au). Although the RMS for this model is 21.95 s, the posterior probability distributions show bimodality with two solutions of the parameters, so the parameters have high uncertainties. Thus, this model is statistically and astrophysically less likely to explain the current $O-C$ diagram.

The system parameters for the two most plausible models are listed in Table 2. For the first model with a quadratic and two LTT terms, the RMS value is calculated as $20.13$ s. For the inner and outer bodies of the model, the semi-amplitudes of the LTT signals were determined to be $K_3=85.17$ and $K_4=76.69$ s, while the orbital periods are $P_3=8.08$ yr and $P_4=9.60$ yr, and the semi-major axis is $a_{3}\sin{i_3}=3.48$ and $a_{4}\sin{i_4}=3.96$ au. These parameters yielded minimum masses of 33.30 and 27.45 $M\text{jup}$ , implying brown dwarfs. The quadratic term with a positive coefficient ( $\unicode{x03B2}=1.16\times10^{-12}$ ) obtained for this model can be associated with a long-term perturbation caused by an additional body. Thus, we also investigated the $O-C$ diagram with three LTT terms without the quadratic term. The final model provides the best RMS value of $\sim$ 20 s of all the models. This model includes third brown dwarf with a minimum mass of 28.34 $M_{\text{Jup}}$ and a semi-major axis of 10.03 au. In the system configuration of this model, the other brown dwarfs have minimum masses of 28.57 $M_{\text{Jup}}$ and 20.96 $M_{\text{Jup}}$ (see Table 2).

The sinusoidal variation in $O-C$ is attributed to the magnetic cycle, as is the case for the LTT effect. The studies in the literature searched for the magnetic cycle through Applegate mechanism (see Völschow et al. 2016; Pulley et al. 2022), but it was reported that the magnetic cycle is not a possible explanation for the orbital period variation of the HS 0705+6700. Our investigation of the magnetic cycle for orbital period variation includes the three different modified Applegate models using the parameters for the LTT term with the smallest amplitude (see Table 3). In the case of the fifth body in the model including three LTT terms, the energy ratios are calculated close to the required energy limit only for the thin-shell magnetic mechanism. It suggests that the magnetic cycle is potentially responsible for the periodic signal in the $O-C$ caused by the fifth body in the system. The orbital period variation in the other models cannot be explained by the Applegate mechanisms alone due to much higher energy ratios than the threshold limit.

Although it is possible to obtain a statistical model explaining the $O-C$ diagram, it is important to ensure that the orbits in the system remain stable for at least a few thousand years. We investigated the stability of the orbital configurations constructed for the model parameters in Table 2. The orbital configurations of all models remain unstable and disrupt the system configuration within 2 000 yr. This agrees with those reported by Sale et al. (2020). The orbital structure of the HS 0705+6700 system appears to be highly complex, according to these results.

The variation in light travel time resulting from the reflex motion of the centre of mass of a HW Vir binary system can be attributed to the presence of one or more orbiting sub-star objects (Beuermann et al. 2012a; Heber, 2016; Baran et al. 2018; Esmer et al. 2021; Brown-Sevilla et al. 2021). For instance, the existence of additional objects orbiting HW Vir has been postulated on the basis of analysis of eclipse timing variations (Beuermann et al. 2012a; Esmer et al. 2021). It is not sufficient to identify the existence of an additional body in motion within the system with the use of the LTT alone. Consequently, Baycroft employed the catalogue of Hipparcos and Gaia proper motion anomalies to demonstrate the existence of a slight indication of a circumbinary companion orbiting HW Vir (Baycroft, Triaud, & Kervella 2023). It has been reported by Baycroft et al. (2023) that the eventual publication of the complete Gaia epoch astronomy will be an important method to confirm the existence of possible additional components around HW Vir and similar systems. Furthermore, to understand the evolution of binary systems, it is important to investigate additional bodies in evolved star systems and explore potential formation scenarios. There are studies indicating the existence of brown dwarfs in common post-envelope binaries (PCEBs) (Perets, 2011; Zorotovic & Schreiber, 2013; Schaffenroth et al. 2015). It was reported that additional bodies can form before the common envelope (CE) phase and orbits evolve due to changes in gravitational potential for a pure first-generation scenario, while additional bodies can be formed from material ejected during the CE in the second-generation scenario. Additionally, a hybrid scenario consists of a combination of the first and second-generation formation scenarios (Schleicher et al. 2015). The formation of the brown dwarf(s) as an additional body(ies) is more likely with this hybrid scenario. It is possible that the brown dwarfs in the HS 0705+6700 system formed before CE and evolved during CE. Furthermore, additional planets may have formed from the second-generation disk. For the latest $O-C$ diagram, the brown dwarfs within the complex and chaotic orbital configuration of HS 0705+6700 are found to be responsible for the observed period variations through the LTT effect. Thus, further observations of this system are needed to ultimately understand the orbital configuration, formation, and evolution of the system.