© ESO 2021

1 Introduction

Studies of variations in mid-eclipse timings can provide invaluable information about various phenomena observed in eclipsing binaries such as mass exchange or loss, apsidal motion, magnetic activity, and the presence of additional bodies. In general, the analysis of eclipse timing variations (ETV) depends heavily on having a sufficient amount of precise photometric observations, which may span over many decades. Circumbinary exoplanet detections are possible due to the effect they induce on the observed eclipsing binary systems. Binaries with extremely short orbital periods are targeted due to (a) their relatively low total mass that enables the detection of a potential, low-amplitude light-time effect (LiTE); and (b) reducing the time spent on the retrieval of a sufficiently large number of observations. Hence, for almost all of the proposed circumbinary planets found with the ETV technique, the host binaries’ orbital periods are as short as a few hours. Some relatively recent ETV studies of short-period binary systems claim detections of gravitationally bound, circumbinary objects with substellar masses. Lee et al. (2009) were among the first to propose circumbinary planets around HW Vir. Sinusoidal trends on eclipse timing variations for some other binaries such as NN Ser (Beuermann et al. 2010), DP Leo (Qian et al. 2010), HU Aqr (Qian et al. 2011), NY Vir (Qian et al. 2012), Kepler-451 (Baran et al. 2015), Kepler-1660 (Getley et al. 2017), and GK Vir (Almeida et al. 2020) were also attributed to circumbinary planets. Some such cases have been re-investigated concerning their orbital stabilities (e.g., HW Vir: Horner et al. 2012, RZ Dra: Hinse et al. 2014, NSVS 14256825: Wittenmyer et al. 2013).

In this paper, we investigate trends that we observed in the ETVs of the short-period binary system HW Vir, based on our own precise photometric follow-up observations with different telescopes, space-borne observations with the Kepler space telescope, Wide Angle Survey for Exoplanets (WASP) telescopes, and timing data compiled from the published literature. We present the technical details relating to our decision-making for the quantitative model selection and perform robust error estimation based on a frequentist approach. In a follow-up paper, we plan to present results using a full Bayesian approach on a very similar dataset.

HW Vir (BD-07 3477) is an Algol-type eclipsing binary with an orbital period of 0.1167 d. The system was first identified as a subdwarf-B star by Berger & Fringant (1980), and later on Menzies & Marang (1986) revealed its binary nature. HW Vir is the prototype of binary systems consisting of a subdwarf-B (sdB) primary and a main-sequence M-type secondary dM. The system has been a topic of research since its discovery for its various aspects. Menzies & Marang (1986) made the first photometric study of the system. They found that the masses of the stars are 0.25 M_⊙ and 0.12 M_⊙ while the temperatures are 26 000 K and 4700 K, respectively. From their single-lined spectrum, Hilditch et al. (1996) derived the temperatures as 33 000 ± 800 K and 3700 ± 700 K for the binary. They did not detect any spectral signature of the secondary star due to the high contrast between the components. Wood & Saffer (1999) were the first to claim the detection of the secondary on the spectra. They derived the semi-amplitude of the radial velocity for the secondary as K₂ = 275 ± 15 km s⁻¹. In order to recover the spectral signatures of the secondary, they corrected their spectra to absolute fluxes by using previous photometric observations and extracted the spectra of the secondary. Edelmann (2008) calculated absolute parameters of the components based on the radial velocity measurements of the much cooler dM-component (secondary),as well as the sdB primary from additional weak absorption lines detected around the secondary eclipse due to reflections on the surface of the less massive dM-star. With these data, the masses were determined as 0.53 ± 0.08 M_⊙ and 0.15 ± 0.03 M_⊙ by Edelmann (2008).

In a more recent study, based on Kepler (K2) data, Baran et al. (2018) showed that the sdB component displayed pulsations with frequencies up to 4600 μHz. They also found the masses of the binary by utilizing the Rømer effect on binaries with different masses (Kaplan 2010). Comparably to the mass derived by Menzies & Marang (1986), they found a mass of 0.26 M_⊙ for the primary. However, the mass that they found is considerably smaller than the canonical helium core ignition mass of ~0.47 M_⊙ (Han et al. 2002; Heber 2009, 2016), which should have taken place for the primary sdB component of HW Vir.

The variations in the mid-eclipse timings of the system was first noticed by Kilkenny et al. (1994). ETVs of HW Vir have been analyzed in many studies to date (Çakirli & Devlen 1999; Kiss et al. 2000; Kilkenny et al. 2003; İbanoǧlu et al. 2004; Qian et al. 2008). Çakirli & Devlen (1999), Kilkenny et al. (2003), and İbanoǧlu et al. (2004) suggested a third body for the cause of ETV within the mass ranges of a substellar object. Lee et al. (2009) ruled out the possible cyclic magnetic activity effect to explain the observed eclipse timings due to the absence of accompanying out-of-transit brightness variations expected by Applegate (1992) and Lanza et al. (1998). Lee et al. (2009) explained the observed timing variations by LiTE induced by two circumbinary planets with masses M₃ = 19.2 M_Jup and M₄ = 8.5 M_Jup. Orbital periods of these planetary mass companions were given as P₃ = 15.84 yr and P₄ = 9.08 yr with the eccentricities e₃ = 0.46 and e₄ = 0.31, respectively, in the same study. Horner et al. (2012) and Beuermann et al. (2012) questioned the validity of these orbital parameters. Horner et al. (2012) performed an orbital stability analysis of the system and reported that the orbits of the circumbinary planets would be dynamically unstable with the parameters given by Lee et al. (2009) on a timescale of a few thousand years. They suggested that the observed variations on mid-eclipse timings of the binary may not have been caused only by the LiTE. Beuermann et al. (2012) also found that the orbits of the planets suggested by Lee et al. (2009) should be unstable. They analyzed the eclipse timings of the binary system using the data set of photometric observations from 1984 to 2012 and gave a new set of parameters for the circumbinary components (M₃ = 14.3 M_Jup, M₄ = 65 M_Jup). They found that the orbits should be stable for more than 10⁸ yr with the configuration that they proposed.

In Sect. 2 of this study, we present our photometric observations of the system carried out by four different observatories. By using our photometric observations and two separate sets of radial velocity data available in the literature, we present the details of a simultaneous light and radial velocity curve analysis in Sect. 3. We formed a combination of two radial velocity data sets by Hilditch et al. (1996) and Wood & Saffer (1999). Since Edelmann (2008) also published the radial velocities for the secondary star, we refer to a separate analysis of his data. After deriving the absolute parameters of the binary from simultaneous light and velocity curve analysis in Sect. 3, we present the results of an ETV analysis based on a long-baseline data set, including the recent observations in the literature as well as our own observations, which we describe in Sect. 4. We attempted fits of the data with various models, extracted the statistical measures of their success, andcompared them, as a result, making use of the Durbin-Watson statistics (Durbin & Watson 1950), χ², and F-tests. In Sect. 5, we provide the details of our error analysis based on an χ² surface search method and bootstrapping to evaluate the uncertainties on each of the model parameters. Finally, in Sect. 6, we describe our search for dynamical stability of the orbital configuration of our adopted ETV model.

2 Observations and data reduction

Photometric observations of HW Vir, starting from 2014 Feb. to 2019 Mar., were carried out with the 0.35 m telescope T35 at Ankara University Kreiken Observatory (AUKR) using Apogee ALTA U47+ CCD, 1 m telescope T100 at TÜBİTAK National Observatory of Turkey (TUG) using the SI 1100 Cryo, UV, AR, BI, 0.6 m telescope at Sobaeksan Optical Astronomy Observatory (SOAO) using a 2k CCD, and the 1 m telescope at Lemmonsan Optical Astronomy Observatory (LOAO) using a 4k CCD. We acquired observations through Bessel BV RI filters. The majority of observationswere done through R filters. All of the observations were made under clear to thin cloud conditions with various moon phases. We used 2 × 2 CCD binning for the observations from TUG, SOAO, and LOAO. All light curves were generated from differential photometry, and comparison stars were checked to be non-variable within nightly observation uncertainties and the limits of our observational setups. For this purpose, we used the ASTROIMAGEJ¹ (Collins et al. 2017) software package, which also allows us to correct our images for instrumental effects (bias-dark-flat corrections), perform differential aperture photometry, and detrend the airmass effect. The log of our photometric observations can be found in Table 1.

For the LOAO observations, the technical staff encountered a software issue in recording the correct time stamp. Technical follow-up tests have subsequently identified and resolved the problem. An offset of 15 min was added to the recorded time stamp in order to obtain the correct HJD_UTC time stamp.

The Kepler Space Telescope observed HW Vir in K2 Campaign 10 in mid-2016. We collected only short cadence photometric data for ETV modeling, since the sampling rate of the long cadence data is not sufficient to detect even the secondary eclipse at all. We calculated 1011 mid-eclipse times with XTREMA software (Bahar et al. 2015), which makes use of Kwee-van Woerden method (Kwee & van Woerden 1956). These data consist of 499 primary and 512 secondary mid-eclipse times.

3 Simultaneous analysis of light and radial velocity curves

In order to derive the absolute physical parameters of the binary, we performed a simultaneous analysis of light and radial velocity curves of HW Vir. We selected our photometric observations with T35 on 2017 Feb. 2 (filter B), 2016 Mar. 9 (filter V), 2016 Feb. 7 (filter R), and 2016 Apr. 27 (filter I). The orbital phase coverage is complete and data is relatively precise (see Table 1 for mean photometric errors) compared to other observations.

The comparison and check stars for these specific observations selected for modeling were TYC 5528-596-1 and TYC 5528-655-1, which are also found to be non-variable within the observational limits. Since the observations of the four nights mentioned were selected for simultaneous light and velocity curve analysis, computations of the orbital phases and normalization were performed individually. Each light curve was normalized with respect to the mean brightness around the brightest phases of the system, which corresponds to the orbital phase between 0.43 and 0.57.

In their work, Lee et al. (2009) informed us of the absence of any long-term change in the light curve of HW Vir. We confirm their claim based on a comparison of the light curves spaced by two years apart from each other. Therefore, we decided to use the light curves from different observation nights in the same light curve analysis. The light curves display two eclipses with significantly different depths. This is due to the large difference between the luminosities of the stars. Other than the eclipse features, there is the reflection effect highly prominent around the second eclipse where the hot primary eclipses the cooler secondary, that is, phase 0.5.

In order to derive system parameters through a simultaneous analysis, we collected radial velocity observations from the literature. The radial velocity data from Hilditch et al. (1996, hereafter H96) and Wood & Saffer (1999, hereafter WS99) were selected for our simultaneous analysis. Both radial velocity data sets are single lined due to the high contrast in the luminosities of the binary. On the other hand, Edelmann (2008) published a few radial velocity measurements for the secondary component, around the orbital phases of 0.5, by making use of the weak absorption lines that are wavelength shifted compared to the sdB’s strong absorption lines. They suggested that these weak lines have formed due to the heating by reflection from the hemisphere of the secondary star facing the hot primary component. Since it was the only study so far including radial velocity data for the secondary, we decided to use the data from Edelmann (2008, hereafter ED08, see Table 2) for a separate analysis. For all the simultaneous light and velocity curve analyses, we initially worked on data from WS99+H96 and adopted the results from this analysis. Then, we analyzed data from ED08 to investigate the consequences of determining the radial velocity as in Edelmann (2008) on the physical parameters of the binary. We gathered the data from H. Edelmann by private communication.

The semi-amplitudes of the radial velocities are 82.3 ± 4 and 83.0 ± 1.2 km s⁻¹ in WS99 and H96 datasets, respectively, which indicates that the amplitude does not change over time (agreement level of 0.17σ) Radial velocity of the barycenter of the system (V γ) was given as2.9 ± 3.1 and − 9.1 ± 0.9 km s⁻¹ by Wood & Saffer (1999) and Hilditch et al. (1996), respectively. Since there are about a few years of difference between the two observations of the two study, we think that such a significant (3.72σ) difference in V γ can be the result of the uncertainties in the parameter and the lack of velocity observations of the secondary, not an actual change in the systemic velocity. The V γ parameter does not have any impact on the binary masses. Therefore, we shifted the data of Hilditch et al. (1996) by the difference of V γ, to remove the velocity shift between the two datasets. We phased and combined them and used this combined set for the first simultaneous light and velocity curve analysis. The same procedure was followed for the Edelmann (2008) radial velocity data forming the basis of a separate analysis. Both datasets can be seen in Fig. 1.

We used PHOEBE v0.31² (Prsa et al. 2011) for the simultaneous light and radial velocity curve analysis for the WS99+H96 dataset. We selected the modeling option for an unconstrained binary system. We fixed the surface temperature of the primary component (T₁) to 28 500 K (Wood & Saffer 1999; İbanoǧlu et al. 2004), albedos for both components (A_1,2) to unity, because of the intense reflection on the secondary, gravity brightening of the primary (g₁) to unity, and that of the secondary (g₂) to 0.32. We left the remaining parameters (semi-major axis (a), mass ratio (q), systemic velocity (V_γ), orbital inclination (i), surface temperature of the secondary (T₂), surface potentials of both of the components (Ω_1,2), luminosity (L₁), and limb darkening of the primary (ld₁) as free parameters. We assumed synchronous rotation concerning the high projected rotational velocity of the sdB component calculated as 74 ± 2 km s⁻¹ by Edelmann (2008) as consistent with the tidal locking of the close system.

Our first trial for the simultaneous analysis yielded a considerable difference for the parameters concerning the luminosities between the model curve and data especially for the B filter. This difference was probably caused by the radiative properties of the hot sdB companion and the intense reflected light from the secondary companion. Therefore, we decided to derive the physical parameters of the binary using only V RI filters and followed the same procedure, which resulted in an acceptable model fit and parameter values. However, there were deviations from the model in certain orbital phases. For the V filter, while the model expected higher luminosities around the secondary minimum than the observed, the deviation is in the opposite direction in the R band for the same phase. The model for the I filter was much better than the other filters, which is most probably due to the strong reflection effect on the secondary component. To better understand the extent of the reflection effect, we analyzed each filter (BV RI) separately, but this time fixing the physical parameters (T_1,2, q, a, i, Ω_1,2, V_γ) that we derived from V RI solution, and adjusting the parameters concerning the luminosities (L_1,2, ld_1,2, A_1,2, g_1,2). The luminosity of the secondary (L₂) was decoupled from the temperature in this step, due to the additional luminosity coming from the reflection. The first problem we encountered in this step was the convergence of the gravity brightening of the secondary to non-physical values. Thus, we fixed g₂ to 0.32. Then, we realized that the albedo A₂ and the decoupled luminosity of the secondary L₂ were affecting the model in the same manner. Therefore, we split the analysis for each filter into two groups, one with a fixed albedo (A₂ = 1.0), the other with a fixed luminosity (L₂), to the same values derived in the V RI solution. The albedo values converged to nonphysical values when we fixed the luminosities. Therefore, we decided to reject these results and adopt only the results for the analyses with the fixed albedos. The resultant luminosities for the individual filters can be found in Table 3, while the absolute parameters can be found in Table 4. Synthetic light and radial velocity curves can be seen in Fig. 2 superimposed on the observational data.

We followed the same modeling approach for the ED08 dataset. We compared the models based on the level of agreement between their results as well as the absolute parameters derived from the analysis making use of the WS99+H96 and ED08 datasets by computing ABS(X − Y) / $\sqrt{\text{d}{X}^2+\text{d}{Y}^2}$. Table 4 clearly shows that the output absolute parameters are significantly different from each other in terms of standard deviations (σ), thus the two models strongly disagree. Radial velocity measurements of the primary for both datasets are consistent with each other. However, models for the radial velocity of the primary differ in a way that the amplitude for the model based on the ED08 dataset is slightly larger than that of WS99+H96.

The method that Edelmann (2008) followed to calculate the radial velocities of the secondary is based on the detection of the weak absorption lines from the reflected hemisphere on the secondary. The reflected hemisphere is mostly visible to the observer around phase 0.5, and the projected surface with respect to the observer decreases when phase differs. Therefore, the signal from the weak absorption lines are at maximum around 0.5. It is plausible to expect that the radial velocities of the secondary are relatively more precise and accurate around this phase, while this precision and accuracy will decrease toward phases 0.25 and 0.75. Hence,it is expected that the velocity measurements of the secondary will deviate more from the true velocities when phase diverges from 0.5. The lower panel of Fig. 1 demonstrates that the secondary velocities at phase ~0.3 and phase ~0.7 are below the model curve, while the velocities that are closest to phase 0.5 are above. We suspect the velocity model for the secondary has smaller amplitude than it should due to the decreased accuracy of the measurements at orbital phases, which diverge from 0.5. This explains the relatively large mass ratio found for EB08 dataset compared to WS99+H96. Thus, we interpret the strong disagreement between the parameter values as the effect of the secondary’s radial velocities (only available in the ED08 dataset) on the mass ratio, which is then propagated to the remaining parameters that depend on this mass ratio.

The total mass of the binary provides key information for the minimum mass of a potential circumbinary object(s) suggested in earlier works. We decided to adopt only one of the solutions of two datasets due to the strong disagreement between the results explained above. While the measurement technique of the radial velocities of the secondary by Edelmann (2008) is plausible and an interesting attempt, we suspect that the results may contain some systematic tendency making the mass ratio closer to unity. Therefore, we decided to adopt the masses of the binary components found in the analysis based on the WS99+H96 radial velocity dataset in our models for the eclipse timing variations, which we present in the next section.

Fig. 1 Radial velocity data along with synthetic curves from simultaneous light and radial velocity analysis. Combined radial velocity data of Wood & Saffer (1999) (white) and shifted Hilditch et al. (1996) (black) is above. Radial velocities from Edelmann (2008) is below. The white triangle marker represents the measurements for the secondary companion. See text for details.

4 Timing data modeling

We derived 66 mid-eclipse times in total, which consist of 38 primary and 28 secondary eclipses, from 28 nights of observations. In order to determine the mid-eclipse times from the light curves of the system, we used the Kwee-van Woerden method (Kwee & van Woerden 1956). All mid-eclipse times were then converted from HJD_UTC to the timescale of barycentric dynamical time BJD_TDB (Eastman et al. 2010), which we list in Table A.1.

In addition to our own observations, we collected 418 mid-eclipse times of HW Vir from the literature based on individual ground-based CCD and photoelectric observations (Menzies & Marang 1986; Marang & Kilkenny 1989; Kilkenny et al. 1991, 1994, 2000, 2003; Wood et al. 1993; Gurol & Selam 1994; Ogloza et al. 2000; Selam et al. 1999; Wood & Saffer 1999; Çakirli & Devlen 1999; Agerer & Hubscher 2000; Gurol et al. 2003; Agerer et al. 1999; Lee et al. 2009; Blaettler & Diethelm 2000; Kiss et al. 2000; Agerer & Hubscher 2002, 2003; İbanoǧlu et al. 2004; Hubscher 2005, 2017; Kotkova & Wolf 2006; Dvorak 2005, 2006, 2008, 2009; Zejda et al. 2006; Nagai 2006, 2007, 2009, 2010, 2013, 2014; Hubscher et al. 2005; Diethelm 2005, 2011; Qian et al. 2008; Brát et al. 2007; Beuermann et al. 2012; Nelson 2009; Brat et al. 2009, 2011; Parimucha et al. 2009; Hoňková et al. 2013; Kubicki 2015, 2017; Basturk et al. 2014; Juryšek et al. 2017; Baran et al. 2018). In total, 336 primary and 82 secondary minima timings have been reported within these studies. We also converted them to the BJD_TDB time standard. Timing data derived from visual observations were not included in this study due to high uncertainties.

There are two more sources in which one can find long-term photometric observations of HW Vir, SuperWASP (SWASP) survey, and Kepler Space Telescope’s K2 Campaign 10. Lohr et al. (2014) published the mid-eclipse times of HW Vir along with some other post-common envelope binaries using SWASP data. These data include two types of timings: good times of minima (179 mid-eclipse times) and extra times (94 mid-eclipse times), as the authors name them. Good minima times have errors comparable with other modern photometric observations and almost one order of magnitude more precise compared to the extra timings. Kepler K2 data has a very short time coverage of 75 d, compared to the whole ETV dataset of ~12 000 d (see Fig. 4). Therefore we used weighted average to create one point resembling all of the K2 data.

In order to draw an ETV diagram, we first corrected the linear ephemeris calculated by Horner et al. (2012) based on a linear fit to the dataset, which we provide in Eq. (1). $$\text{Min} I={\text{BJD}}_{\text{TDB}} 2,450,280.28472(4)+0.116719504(1) E,$$(1)

where, Min I refers to the time of primary eclipse and E refers to the epoch. We found the root mean square (RMS) of residuals of linear fit as 52.76 s and the reduced chi-squared ${\chi }_{\nu }^2$ as 89.28. The RMS of the linear fit is found to be an order of magnitude larger than the average standard errors of the observations (≈5.6 s), which is a direct indication of an ETV variation. We formed the ETV diagram in the usual manner by subtracting the observed times of minima from that computed from the linear ephemeris that we corrected (see Fig. 4). The resultant ETV diagram has clear signs of potentially multiple cyclic trends ans is also easily distinguishable by eye. The amplitude of the trend appears to be changing due to the nature of the responsible mechanism(s), or there is a combination of more than one cyclic variation. We used only primary eclipse timings in the fitting procedure due to their smaller scatter over the general trend, and a potential apsidal motion is ruled out because secondary mid-eclipse timings are following the same trend with the primaries in ETV diagram. This is exactly what we expect for a circular orbit. We discarded a few of the timing data from the literature, where the published uncertainties are high, thus unreliable. A small portion of the early literature data, which follow the general trend, have mid-eclipse timing measurements without quoted uncertainties. In order to find an uncertainty estimate of these data, we calculated the mean standard errors of the data as σ_i = 6.48 × 10⁻⁵ d = 5.6 s and assigned this value as uncertainties.

Fig. 2 Observed light curves with best-fit models and their residuals in various pass-bands.

Fig. 3 Lomb-Scargle spectrum of the ETV data. Upper panel is for the whole dataset, and the lower panel is for the residuals from the highest frequency of the whole dataset. We plot the 1% FAP line for reference.

4.1 Analysis of the orbital period variation

Since the mid-eclipse timings data form an unevenly distributed dataset, we made use of a Lomb-Scargle (LS) periodogram which we computed with our own PYTHON code based on the ASTROPY³ package (Astropy Collaboration 2013, 2018) functions to detect the periodicities. We found the highest peak in the period spectrum at f₁ = 1.20 × 10⁻⁴ d⁻¹ with an amplitude of 126.4 s and false alarm probability (FAP) of 2.4 × 10⁻⁷⁵. There are some other peaks around 2.7 × 10⁻³ d⁻¹, which are due to the annual sampling of the ETV data. We then removed the highest frequency from the original dataset and found a second frequency at f₂ = 1.89 × 10⁻⁴ d⁻¹ with an amplitude of 93.3 s andFAP of 1.6 × 10⁻⁹⁴ on the residual spectrum (see Fig. 3). Since we only removed f₁ from the ETV data, the frequencies arose from the annual sampling remained in the second spectrum as well.

We also removed f2 from the ETV dataset to check any remaining periodicity. The highest peak on the residual spectrum corresponds to a period more than 100 times longer than the total time span of the dataset. This is consistent with a period decay, a downward parabola in the ETV, or a very long period sinusoidal. The amplitude of this frequency is almost 2 d, which is unlikely to be a part of a long period sinusoidal, since the orbital period of the binary is only 2.8 h. Therefore, we interpret this signal as a secular change.

In order to avoid overfitting the data and for simplicity, we decided to stick to a maximum of two cyclic models in modeling the ETV. Therefore, we investigated the potential of the first two peaks in our frequency analysis with LS periodograms, found at f₁ and f₂ corresponding to periods of 8320 d (~22.8 yr) and 5298 d (~14.5 yr), respectively. Considering the total time span of the data (~13 200 d), it should cover more than one full cycle of the variation with f₁ and more than two cycles with f₂. As a result, we decided to use these first two peaks in a first attempt to model the observed ETV.

4.2 ETV modeling

In order to account for the accumulation of the uncertainties on the ephemeris parameters, we added a linear expression to each of the models that we think have the potential to represent the observed ETV trend. The lowest significant frequency found in the LS periodogram could be interpreted as a secular change of the eclipse timings. Therefore, we decided to select the models to allow for a quadratic term, β, to be added if needed. Therefore ETV models (T_C) with β would have a form of $$T_{\text{C}}=(T_0+\Delta T_0)+(P_0+\Delta P_0)E+\beta E^2+\tau (E).$$(2)

Here, ΔT₀ and ΔP₀ are the correction terms on the reference mid-eclipse time (T₀) and the orbital period of the binary (P₀), respectively,and τ(E) is a LiTE model(s). The models without the quadratic term will have a form similar to the one above, but without the βE² term. By using two different groups of models with and without a quadratic term, it can be tested if a secular change will improve the fits based on the statistical significance of the results.

From a visual inspection of the phased light curve (Fig. 2) and the ETV diagram (Fig. 4), we do not see any evidence for apsidal motion, so we discarded it. The signals found in LS analysis may be due to LiTE caused by additional bodies or due to magnetic activity. LiTE would lead to a periodic variation with constant amplitude, while magnetic activity causes cyclic changes with variable amplitudes both in the short term (modulation due to stellar spots with rotation) and in the long run (magnetic activity cycle). We fit a sinusoidal model without the eccentricity term with the assumption of magnetic activity modulation of a secondary star as described in Applegate (1992). The corresponding period is 77 yr, and the magnetic field should be B ≈ 39 kG, which is an unrealistic level of activity expected from an M-type dwarf. Therefore, we decided not to assume magnetic-activity-induced, cyclic orbital period variations in our models. For the sake of simplicity, we used the analytical expression of LiTE (Irwin 1959) for an initial model of the cyclic trends and discuss the choice between the mechanisms in Sect. 7.

The models that we used to fit the ETV data are, linear (Model 1), linear + quadratic (Model 2), linear + LiTE (Model 3), linear + quadratic + LiTE (Model 4), linear + two LiTE (Model 5) and linear + quadratic + two LiTE (Model 6). The analytical expressions for each component of the models and the short-hand notations are shown in Table 5. These models are fit separately, and corresponding goodness-of-fit statistics (RMS, χ², ${\chi }_{\nu }^2$) are calculated for each of them. For all the models that include LiTE, the mass function for a circumbinary object can be calculated as described in, for example, Wolf et al. (2018): $$f(m_3)=\frac{M_3^3{\sin }^3{i}_3}{{(M_{\text{bin}}+M_3)}^2}={[\frac{173.15A}{\sqrt[]{1-e_3^2{\cos }^2{\omega }_3}}]}^3\frac{1}{P_3^2}.$$(3)

Here, M₃ is the mass, i₃ is the orbital inclination, e₃ is the eccentricity, P₃ is the orbital period, ω₃ is the argument of periastron of the circumbinary object, and A is the semi-amplitude of LiTE, while M_bin is the total mass of the binary stars. P₃ is in units of years, masses are in units of solar masses, and the constant 173.15 arises from the conversion of units. Assuming the stellar binary as a single object with the total mass of the binary, the minimum mass of the circumbinary object, M sin i, regarding each LiTE can be found using an iterative method if the mass of the binary is already known. For the case of HW Vir, we used the masses that we derived from simultaneous light and radial velocity curve analysis in Sect. 3.

We performed the fitting procedure on both the original dataset and the mid-eclipse times binned to remove the seasonal variations. Since the cyclic trends in the ETV of HW Vir that we found in frequency analysis are of the order of tens of years, seasonal averaging of one year should not diminish the periodic signals and is therefore plausible. By using seasonal averaged data, we minimized the effect of the outliers on the parameter values and their uncertainties, which made it possible to compare with the results of the fitting to the original dataset.

We used our own PYTHON code that makes use of the Levenberg-Marquardt algorithm for a nonlinear least-squares fitting (Levenberg 1944; Marquardt 1963) based on the LMFIT package functions by Newville et al. (2014) for the fitting procedure. Using initial parameter values, the Levenberg-Marquardt algorithm can find local minima in the parameter space. Since the Levenberg-Marquardt algorithm requires initial parameter guesses, we used the amplitude and period values derived from a Lomb-Scargle periodogram to derive initial guess values for the models including LiTE. Initial values for the remaining parameters were found from several trial test-fits and a visual inspection as implemented in a spreadsheet program. The results containing the best-fit parameter values, RMS, and ${\chi }_{\nu }^2$ for every model fit to the whole range of data, as well as seasonal binned data, can be found in Table 6. Best-fit curves for each model fit can be seen in Fig. 4 overplotted on the full dataset and the averaged data. In the following, we report the results of model fits to the original dataset and to the seasonal binned dataset.

Model 1 is the linear model that we used to update the ephemeris information as we mentioned earlier. RMS and ${\chi }_{\nu }^2$ for Models 1 and 2, for fitting original and seasonal binned datasets are large enough to indicate the need for further improvement to the models. The difference between the goodness-of-fit statistics for two datasets arises from (i) the difference of the model parameters, and (ii) the difference of ν and standard errors of data points between two datasets.

With the addition of LiTE to Model 3, fit statistics decreased dramatically compared to Models 1 and 2. However, the best-fit results look far from convincing. For the original dataset, the model fit indicates an additional body with a very wide, long period (10⁶ d ≈ 2750 yr) and extremely eccentric orbit, which coincidentally should have passed from its superior conjunction relative to Earth within the time frame of the dataset (~35 yr). The results for a seasonal binned dataset also indicate a similar, highly unlikely configuration with a relatively longer period.

The addition of a quadratic term to Model 4 resulted in almost the same fit statistics. The orbital period was found to be 9399 d (26 yr) for the original dataset and 10 499 d (29 yr) for binned dataset with almost the same eccentricity (e ~ 0.5). While the best-fit parameters are in an acceptable range, the fit statistics imply that Model 4 is not a plausible explanation to the ETV.

We introduced a second LiTE with Model 5, for which the fit statistics were found to be within acceptable limits, and in fact the best of all of the models. ${\chi }_{\nu }^2$ for the original data is around unity, and RMS is ~5.6 s. However, periods, ETV amplitudes (A_LiTE), eccentricities and corresponding minimum masses (M sin i), and semi-major axes (a sin i) calculated for the original dataset are almost identical for each of the two LiTE models included in Model 5. This infers that there should be two objects with almost identical orbital and physical properties according to the best-fit solution of Model 5. The only considerable difference is the argument of periastron, ω, and time of periastron passage, T_0,LiTE. M sin i values indicate that the masses of the additional bodies are at the limit of stellar masses, and the uncertainties of M sin i and a sin i are enormous. The minimum masses calculated from the best-fit parameters values for binned dataset are even higher than the mass of the secondary companion of the binary. We highly doubt this solution represents the case for HW Vir.

In terms of fit statistics and best-fit parameter values, one of the satisfactory ETV models has been achieved forthe case of Model 6 (Figs. 4 and 5). Both the fits for original and seasonal binned dataset gives almost the same solution. LiTE results of Model 6 correspond to M sin i of 24.6 and 13.89 M_Jup, respectively. While the periods of Model 6 differ from the results of LS analysis, they still fall in a time span shorter than that of the dataset. The orbits are eccentric and seem to be co-oriented in space as ω values are almost the same.

The best-fit values for the binned dataset are almost the same as original dataset. The orbits are slightly closer to each other, and the outer orbit is somewhat more eccentric compared to the best-fit values of the original dataset. ${\chi }_{\nu }^2$ statistics of the binned dataset indicate overfitting or overestimated uncertainties. As Hughes & Hase (2010, p. 107) mentioned, the number of data points (N) has to be similar to ν in order to have an ${\chi }_{\nu }^2$ of unity. However, for seasonal binned data, the ratio N/ν ≈ 1.68, while the same ratio for the original dataset is 1.03. Therefore, the ${\chi }_{\nu }^2$ statistic for the binned dataset may not be a suitable goodness-of-fit indicator.

The radial velocity semi-amplitudes corresponding the two LiTEs of the Model 6 solution are both ~ 0.2 km s⁻¹ with the assumption of i = 90°. Therefore, we do not expect a relatively small semi-amplitude to be detectable from the radial velocity datasets that we used in Sect. 3 in terms of their observational uncertainties, the number of observations, and that the system is a single lined eclipsing binary.

Finally, we performed a frequency analysis on the residuals of Model 6. The LS periodogram of the residuals from the original dataset has its highest peak at f_res = 5.9 × 10⁴ d⁻¹, which corresponds to a period of 1690.7 d and a small FAP value of 6 × 10⁻⁴. However, the amplitude of the peak is 3.97 s, which is less than the mean standard deviation (~5.6 s) of the timing data and less than the RMS for Model 6 fit to the original data (~5.7). We did not find any similar frequency with a FAP less than 0.1 on the LS periodogram of the residuals of the binned dataset. Therefore, we interpret f_res to resemble the sampling frequency of the original dataset or its multiples.

Fig. 4 ETV diagram of HW Vir with best-fit curves of all models. The ETV diagram is corrected by the results of the linear fit. Color-filled circles represents the primary mid-eclipse timings, while white-filled ones represent the secondaries. The black markers are for data from the literature, blue ones are for SWASP, green for Kepler K2, red for our own observations, and yellow markers represent seasonal binned data.

4.3 Durbin-Watson test

Often within the framework of a least-squares minimization, the residuals are qualitatively assessed to judge any trends or structures (auto-correlations) that might be present in the data. The nature of a structure is either astrophysical, and/or a non-Gaussian measurement process is involved. If the measurement errors are assumed to follow a Gaussian distribution with ~ N(0, σ²), the presence of auto-correlation (change of variance with time) could be interpreted as an additional signature of astrophysical origin. We applied the Durbin-Watson (DW) statistic (Hughes & Hase 2010) to quantitatively assess the quality of a given model fit. In addition, we drew so-called lag-plots to visualize the underlying concept of the DW statistic for each residual set. A lag plot is constructed by plotting the residuals of a given model against a lagged residual by selecting a lag interval k = 1. Lag plots exhibiting no correlation suggest that the residuals do follow a random normal distribution. Lag plots with correlated (positive or negative) residuals indicate some degree of auto-correlation. The DW statistic is given as, $$\text{DW}=\frac{{\sum }_{i=2}^N{[R_i-R_{i-1}]}^2}{{\sum }_{i=2}^N{[R_i]}^2},$$(4)

where R_i is the residuals in the original order and R_i−1 is the k = 1 lagged residual. The range of DW is 0 to 4 with the two extrema corresponding to anti-correlation and auto-correlation, respectively, and a DW = 2 indicates no auto-correlation in the investigated residuals.

With the increased model complexity in our analysis, the DW statistic converges to 2 (DW ≈ 2) and the auto-correlation seems to diminish (Fig. 6). The lag plots of Models 1–4 show clear signs of deviation from normal distribution around the origin. The DW statistics of Models 1–4 are, 0.055, 0.057, 0.326, and 0.340, in respective order, all of which are far from DW = 2. Adding a β parameter has no significant effect on either the DW statistic or the auto-correlation of the residuals, when comparing Models 1 and 2 as well as Models 3 and 4.

The models with two LiTEs have almost the same DW statistics, while Model 5 has a DW statistic slightly closer to 2. However, the difference is so small that any conclusion from DW statistics cannot be achieved. The auto-correlation of the residuals in lag plots for Models 5 and 6 looks almost the same, and these are both close to being normally distributed, as can be seen in the figure. With DW values deviating from 2, Models 5 and 6 can be interpreted as follows: either (i) an additional trend of likely astrophysical origin is present if the timing errors distribute normally, or (ii) the deviation from DW = 2 is due to measurement errors and/or the presence of additional astrophysical signals in the data if timing errors are not normally distributed.

The ETV dataset of HW Vir consists of observations performed in different atmospheric conditions with various instruments having different characteristics. Mid-eclipse times are also measured following different methods from data reduced making use of different packages. Moreover, variable atmospheric conditions can introduce time-correlated, so-called red noise, even on the observations from the same instrument. All of these factors, and perhaps others, may be responsible for timing errors not being normally distributed (for a detailed discussion on potential error sources, see von Essen et al. 2016). Even if the observational errors of ETV data are not normally distributed, we cannot discard the possibility of assuming wrong ETV models and/or additional astrophysical signals over assumed models on the ETV.

Fig. 5 ETV diagram of HW Vir with best-fit curve of Model 6 and its residuals. Color-filled circles represent the primary mid-eclipse timings, while white-filled ones represents the secondaries. The black markers represent data from the literature, blue ones are for SWASP, green for Kepler K2, red for our own observations, and yellow markers represent seasonal binned data.

Fig. 6 Lag plots corresponding the six models considered in this study. Lag plots (a) for Models 1–4 (red, green, blue, and yellow, in respective order), (b) for Model 5, and (c) for Model 6. Plots (b) and (c) show the Durbin-Watson (DW) statistic. We refer the reader to the text for more details.

4.4 Model comparison: χ² test

We compared the p-value, p(χ², ν), of each model calculated from an χ² distribution for a given number of degrees of freedom, with a critical probability indicated by the significance level (α) to test the null-hypothesis stating that the model under investigation fits the data well (p < α). In this case, the probability of finding the observed χ² is too low at the given significance level and cannot be explained by a random process. When α <p ≤ 0.5, the discrepancies between the model and the data are random; there is no evidence to reject the null-hypothesis. Finally, if 0.5 < p < 1, the standard error used in the determination of χ² is overestimated, resulting in an unrealistically small value of χ² (Hughes & Hase 2010).

At the heart of a statistical significance test is the assumption that data are independent and normally (random) distributed. Assuming the null-hypothesis, the expectation value of χ² is close to the mean of the (χ², ν) distribution: χ² ≈ ν with probability p(χ², ν) ≈ 0.5. In the case of considerablylarge χ² values, the model function that we assumed to fit the data well in the null-hypothesis is unlikely to explain the observations, which can be quantitatively decided by comparing the probability of the χ² to a critical probability, α. If the χ² is considerably less than the mean value of the χ² distribution, the probability of the χ² becomes close to unity, which eventually means the standard errors of the data are overestimated; in other words, we are fitting the noise of the data.

The critical probability, alpha, sets the significance level (1-alpha) and is somewhat arbitrary depending on the nature of the underlying problem. Often α is chosen to be 0.001 (0.1%), 0.01 (1%), or 0.05 (5%). The significance levels of rejection and non-rejection are 99.9, 99, and 95%, respectively. The p-value itself provides a measure of the strength of the evidence against the null-hypothesis. As a rough guideline, we adopt the following criteria: if p(χ²) < 0.01, we judge very strong evidence against H₀; if 0.01 < p < 0.05, we judge strong evidence against H₀; if 0.05 < p < 0.10, we judge some weak evidence against H₀; and if p > 0.10, we judge little or no evidence against H₀.

Table 7 shows the calculated χ² for each model that we considered in this study. To calculate the probabilities p(χ², ν) for each model, we made use of our own PYTHON code. For Models 1–4, the corresponding p-values are smaller than 10⁻²⁰⁰, thus with a high level of significance, and under the mentioned assumptions, we can reject the null-hypothesis for these models.

For Models 5 and 6, we find the p-value to be 0.285 and 0.155, respectively. Considering the significance level (α = 0.1), we see no evidence against Models 5 and 6, since p_model > α. In a relative sense, we therefore interpret this result as a non-rejection of the null-hypothesis for Models 5 and 6.

4.5 Model selection: F-Test

In the previous section, we determined the most likely model by comparing the minimum χ² of each model with the χ² probability density distribution for ν degrees of freedom under the assumption of the null-hypothesis. We find an χ² value for the Model 5, indicatingthe greatest probability amongst other candidates, and hence deem it to be very significant not to reject the null-hypothesis based on the p-value. While evaluating a low χ² is a necessary requirement for assessing the goodness of fit, this method also has limitations favoring models overfitting the data. The more the adjustable parameters are included, the better the overall goodness of fit will be (lower variance). Therefore, one drawback of the χ² test is that it does not have a mechanism to distinguish models with increasing numbers of adjustable model parameters (for the same number of data points).

In this section, we describe the application of the F-test to our problem and report the results. This test aims to quantitatively select one model in favor of another model. Intuitively, a model (2) with more adjustable parameters should describe a given dataset at least as well as a model (1) with fewer adjustable parameters. Hence, model 2 will provide lower residual errors. The question is then: how significant is the improvement of model (2) over model (1)? The test statistic quantifying this problem is given as $$F=\frac{({\chi }_1^2-{\chi }_2^2)/({\nu }_2-{\nu }_1)}{{\chi }_2^2/{\nu }_2},$$(5)

where ν_i are the corresponding degrees of freedom, and ${\chi }_i^2$ are the χ² statistics for any two models that are being compared.

With the null-hypothesis stating that the model (2) does not describe the data significantly better than model (1), the F-statistic follows the F probability density distribution. To assess the rejection probability, we again used the calculated p-value for the observed F-statistic as outlined in the previous section.

By using our own PYTHON code, we calculated F-statistics as given in Eq. (5) by adopting χ² and ν values from the best-fit results from Sect. 4.2 (see Table 7), for each model pair ordered by the corresponding ν. We then calculated the probability of F-statistics in an F-distribution defined by the two degrees of freedom, P(F;v₂ − v₁, v₂). Similar to the χ² test, we chose the critical probability, α, as 0.1 (10%), and the condition of rejecting the null-hypothesis is in the cases of p < α, while if p > α, there is no reason to state that the simpler model does not fit the data statistically well compared to the more complex one. The model pairs and their corresponding F-statistics and p-values can be seen in Table 8.

Other than for Models 5 and 6, p-values of all of the remaining F-tests were closer to zero than being comparable to α of 0.1, thusthe null-hypotheses for all of these cases were rejected. F-value being negative for the test between Models 5 and 6 means that not only did the additional parameter not improve thefit, but also the more complex model has poorer statistics than the simpler one. Therefore, the p-value is unity and there is no reason for rejecting the null-hypothesis.

In all of the model selection and comparison techniques performed in this study, Models 5 and 6 are significantly better compared to the other simpler models. The statistical significance of Model 5 over Model 6 could be misleading. Such orbital architecture is also expected to be unstable even within a few orbital revolutions. Nevertheless, we checked the stability of Model 5 in Sect. 6. From a physical point of view, we expect that Model 5 is not feasible. On the other hand, the additional objects of Model 6 have much smaller minimum masses with a relatively larger separation between the two. The observational signature of these two bodies may expected to be undetectable, if their orbital inclinations are relatively close to 90°, thus their true masses are close to the minimum masses. The eccentricities are larger compared to Model 5 and should be subject to the dynamical analysis. However, we do not see a direct reason to reject Model 6. Therefore, we adopt Model 6 as a possible explanation to the ETV of HW Vir.

5 Error analysis

Preliminary uncertainties for each parameter of Model 6 have already been assigned as shown in Table 6. However, these are formal errors because they have been derived from the covariance matrix as a result of the LMFIT least-squares minimization process. In general, formal uncertainties are not reliable, since parameter correlations are not taken into account.

In order to compute realistic parameter uncertainties for Model 6, we made use of two different and independentalgorithmic methods as part of a comparative analysis study. The first method explores the projected χ² surface for any two parameters. We believe this method is somewhat overlooked in the astronomical data analysis work, and we provide an in-depth description of the underlying machinery. The second method is based on the bootstrap Monte Carlo resampling method, which is widely practiced in data analysis in astronomy.

5.1 χ² surface search

Parameter uncertainties can be determined by exploring the projected χ² surface around ${\chi }_{\text{min}}^2$. In the following, we chose to adopt the technique outlined in Hughes & Hase (2010, p. 74). The 68.3% confidence level for a given parameter is calculated numerically by successive iterations until the condition ${\chi }_{\text{min}}^2$ + 1 is met.This method is then applied to two-parameter projections of χ² in order to account for the parameter correlations that are otherwise neglected when quoting formal uncertainties. For a systematic change of any two parameters around their best-fit values we monitored the χ² value while letting the remaining parameters to vary freely. The upper and lower standard deviations were then determined from ${\chi }_{\text{min}}^2$ + 1. Any parameter correlation is propagated and accounted for in the calculated confidence intervals. For complicated correlations, the χ² topology leads to asymmetric error bars around the best-fit values (Bevington & Robinson 2003).

There are 13 parameters in the analytic expression of Model 6. This corresponds to a total of 78 parameter pairs. We selected each pair within these sets and fixed them to 75 different values around the best-fit model considering appropriate ranges. Thus, the resolution of the χ² surface was 75  × 75, and we calculated a total of 438 750 fits. For each fit, we calculated and stored the resulting χ² corresponding to a given parameter pair. To map out the χ² topology in some more detail, we chose to adopt an exponential parameter gradient around the best-fit model. This approach provides a higher resolution of the χ² landscape for values close to the best-fit model and a lower resolution far from it. The reason for this choice is to reduce the computation time and assign weights to the topology change in the vicinity of ${\chi }_{\text{min}}^2$. For each grid point, we stored the resulting χ² value.

We converted the resulting χ² surfaces to a grid of probability density as P(x, y) ∝exp(−χ²∕2), where the latter quantity is the maximum likelihood and x and y represent the two parameters. We normalized the probabilities to impose the constraint ∫ P(x, y)dxdy = 1. We created 12 marginalized probability density distributions for each parameter by summing the probability values that correspond to the same parameter value on the grid of each pair. We then summed these 12 probability density distributions to create a final normalized probability distribution for 13 parameters. Finally, we calculated the 68.3% confidence intervals as well as median values from the final probability distributions.

In addition to the model parameters, we calculated derived parameters of the mass function (f(m)), minimum mass (M sin i), and projected semi-major axis (a sin i). We calculated their confidence intervals by following the functional approach of error propagation for multi-variable functions (see Hughes & Hase 2010). The procedure is relatively straightforward, but there is merit in some detailed explanation for the purpose of clarity. The derived parameter is defined as X = f(α, β, …), where (α, β, …) are the median values of parameters determined from the marginalization process. To obtain a conservative estimate we then added the upper limit in uncertainty to the median of a given parameter and calculated the difference f(α + σ_α, β, …) − f(α, β, …). This procedure and calculation of differences were also determined for the other parameters f(α, β + σ_β, …) − f(α, β, …). The final uncertainty for X is then obtained from the quartiles as $$\begin{array}{lll}{\sigma }_X^2\hfill & =\hfill & {[f(\alpha \pm {\sigma }_{\alpha },\beta ,\mathrm{...})-f(\alpha ,\beta ,\mathrm{...})]}^2\hfill \\ \hfill & \hfill & +{[f(\alpha ,\beta \pm {\sigma }_{\beta },\mathrm{...})-f(\alpha ,\beta ,\mathrm{...})]}^2+\mathrm{...}\hfill \end{array}$$

Since the uncertainties calculated by the χ² surface search are not necessarily normally distributed and their distribution can be asymmetric, we needed to calculate the uncertainties both for the positive and negative sides separately by using the formula above. According to Hughes & Hase (2010, p. 41), this procedure oferror propagation assumes the variables that are needed to calculate the derived parameters to be uncorrelated and independent. To determine the derived parameters, we adopted the uncertainties derived from 68.3% confidence intervals of χ² surface. Therefore, the propagated uncertainties are not formal and the effects of correlations on the error surface were already taken into account. The normalized probability distributions for each parameter are shown in Fig. B.1 in the appendix, and the values that correspond to the 68.3% confidence levels are shown in Table 9. χ² surface plots are shown in Appendix C.

5.2 Bootstrap method

Another method to determine parameter uncertainties is known as the re-sampling bootstrap method (Efron & Tibshirani 1986). Again, we aim to determine the parameter uncertainties in Model 6.

The boostrap method is also referred to as the ‘quick and dirty Monte Carlo method’ in Press et al. (1992), who found wide-spread acceptance in frequentist astrostatistics to make inferences of the parent distribution for each parameter. The deriveduncertainties can then be compared to the uncertainties obtained from the marginalization process as described in the previous section. The underlying algorithm is based on random sampling from the original dataset with a replacement. A given simulated (bootstrapped) dataset has the same number of data points as the original dataset. The replacement condition introduces the possibility that data points can be duplicated in a given bootstrap set. Since the dimension of the dataset is equal, this implies that some data is missing among simulated datasets.

There is, however, one question remains to be addressed concerning the number of random samples. To answer this question, we performed a systematic numerical experiment where we calculated parameter distributions for a range of the number of random samples. We considered 1 × 10⁴, 2.5 × 10⁴, 1 × 10⁵, and 2 × 10⁵ random samples and plotted the resulting distributions. We found no difference between 1 × 10⁵ and 2 × 10⁵ bootstrap samples. Clearly, choosing 1 × 10⁴ is too small a number to warrant a reliable large-number statistic of the final result.

As described in the previous section, we again determined the 15.9, 50, and 84.2% quartiles from each parameter distribution in order to obtain the median and 68.2% confidence interval. In Fig. B.2 in the appendix, we give the histograms obtained from the bootstrap method. Table 9 summarizes the uncertainties for Model 6 parameters. We would like to emphasize the relative strength of the bootstrap method, which is proven to work on any parent distribution, even without knowing the nature of it at all. Hence, we recommend taking the usage of uncertainties from bootstrap into account in any related study in order to be cautious in further interpretation rather than stating only the formal errors.

6 Stability of Model 6

The orbital parameters of Model 6 derived from the ETV analysis correspond to two planetary-mass companions in eccentric orbits that overlap. The additional companions thus undergo strong mutual gravitational interactions, and system stability is only possible if a resonant mechanism is present to prevent close encounters. Therefore, a dynamical stability search of the orbits derived from Model 6 is needed.

We performed a stability search for the results from the LM fit, a bootstrap search, and an χ² surface search. We used frequency map analysis (FMA; Laskar 1990, 1993), similarly to the studies of Correia et al. (2010) and Couetdic et al. (2010). We calculated the normalized stability index D as $$D=\frac{\text{|}{n}_1-n_2\text{|}}{n_1},$$(6)

where n₁ is the frequency of the mean motion calculated from the first half of each integration, and n₂ is the same for the second half. In the case of regular motion, the normalized stability index should be D < 10⁻⁶ (Correia et al. 2005).

For the integrations, we used the MERCURY6 code (Chambers 1999) with the RADAU algorithm. We set the timestep as 100 d and total integration time as 3 × 10⁶ d. We calculated n₁ and n₂ frequencies using a TRIP code (Gastineau & Laskar 2011). We integrated the system by varying the semi-major axis and eccentricity of the outer companion while keeping the orbital parameters of the inner companion fixed since the uncertainties on the inner companion are smaller compared to those in outer orbit. The orbital period of the inner circumbinary orbit is more than 4 × 10⁴ longer than the orbital period of the central binary. In this case, interactions between an individual star and a circumbinary object can be neglected. Therefore, we assumed that the central binary as a single object with a mass of the total binary mass. We also assumed that the orbits of the two circumbinary objects are coplanar and inclinations are 90°. Therefore, we fixed the masses and the other orbital parameters to the values in Table 9. We set the step size for a as 0.04 au and 0.005 for e.

Almost all of the orbits within the searched a–e range are unstable other than a small stable region (D < 10⁻⁶) around a ~ 10.5 au, e ~ 0.12 and another one around a > 12.2 au, e < 0.1 (Figs. 7–9). None of the orbital configurations of Model 6 are in stable regions within a 3σ uncertainty range. In the presence of the inner companion with the parameters derived from Model 6, the outer companion should have stable orbits for a > 10 au and e < 0.2 (black regions in Figs. 7–9). We tried to fit the ETV by using the parameters that correspond to the stable regions mentioned above. However, none of these trials corresponded to satisfactory fits. Therefore, we conclude that Model 6 is not plausible in terms of explaining the ETV of the HW Vir system.

We also checked if the circumbinary stellar mass objects of Model 5 have any stable orbits. By assuming the central binary as a single object and keeping the two orbits coplanar, we simulated the best-fit orbital configuration of Model 5. As expected, such massive objects are interacting with each other very strongly, and the outer object is being ejected in less than one orbital revolution. Therefore, Model 5 can be ruled out as well.

Fig. 7 Stability analysis of Model 6 in Sect. 4.2 (Table 6). The color scale represents the normalized stability index in log scale. The black circle represents the best-fit values of the outer orbit. The errorbars represent 1, 2, and 3σ uncertainties on the best-fit values (see text for more details).

Fig. 8 Stability analysis of Model 6 from χ2 surface search in Sect. 5.1 of Model 6. See Fig. 7 for the description of the figure (see text for more details).

7 Conclusions

From the simultaneous light and velocity curve analysis in Sect. 3, mass of the primary is found to be 0.413 M_⊙ and 0.307 M_⊙ for the WS99+H96 and ED08 datasets, respectively. Both of the masses are below the canonical helium ignition mass of ~0.47 M_⊙. However, similar values below the canonical mass have been given in the literature previously. Menzies & Marang (1986) found an even smaller mass for the primary companion as 0.25 M_⊙. A recent study of Baran et al. (2018) used the method of Kaplan (2010), which makes use of LiTE, and derived a mass of 0.26 M_⊙ for the primary of HW Vir. On the other hand, Edelmann (2008) found a larger value for the mass of the primary as 0.53 M_⊙. This wide range of the masses of the sdB component in the HW Vir system is still an open question requiring further investigation, which is beyond the scope of this work. HW Vir is a challenging case for light curve analysis due to its compact orbital configuration and high contrast in the luminosities of the individual components. We do not intend to speculate about whether there is a core helium ignition process in the primary, since our intention is to investigate eclipse timing variation of the binary.

Due to its very low luminosity compared to the sdB primary, little is known about the secondary star. We found 0.128 M_⊙ and 0.110 M_⊙ in our analysis based on theWS99+H96 and ED08 datasets, respectively, for the mass of the secondary. Literature values for the secondary mass change between 0.1 M_⊙ from Baran et al. (2018) and 0.15 M_⊙ from Edelmann (2008). Almost all these masses are within the range of a main-sequence M star. However, the temperature values that we found are somewhat higher for an M-type main-sequence star, probably due to the reflection effect. For both WS99+H96 and ED08 datasets, we found surface temperatures of ~3900 K for the secondary. The properties of the secondary companion in HW Vir, as well as many other secondaries in sdB+dM binaries, will be tentative until precise spectrographic observations can be obtained to distinguish it from those of the hot primary. The sizes of the primary and the secondary stars are found to be almost the same for both of the radial velocity datasets, which is consistent with the partial primary and secondary eclipses, in spite of the tight orbital configuration (a ≈ 0.8 R_⊙) observed under rather high inclination (i ≈ 81°).

Our ETV analysis for the seasonally binned data gave almost the same results as the original unbinned data. Since the binning size is considerably shorter than the cyclic trends on ETVs, we did not encounter any loss of information; for example, LiTE amplitude, period etc. Moreover, taking weighted averages of the data minimized the effects of the outliers, which can impact the results of the least-squares minimization methods (Mandel 1964). From our results, we can conclude that using seasonally binned data is permitted for the case of HW Vir.

The β coefficient of Model 6 in the ETV analysis corresponds to a constant period change with a linear rate of − 9.16 × 10⁻⁹ d yr⁻¹. The angular momentum change can be derived (e.g., Brinkworth et al. 2006) as d J∕dt = −3.25 × 10³⁵ erg. The 3σ confidence interval of the angular momentum change should be within the range of − 4.60 × 10³⁵ and − 2.36 × 10³⁵ erg. The possible mechanisms of angular momentum loss may be gravitational radiation or magnetic stellar wind breaking. By using the absolute parameters derived in Sect. 3 and the formulation from Paczyński (1967), the gravitational radiation of the binary should be − 6.57 × 10³² erg. The angular momentum loss calculated for the β coefficient is three orders of magnitude larger as it is caused by gravitational radiation. On the other hand, magnetic stellar wind breaking due to the secondary star of the HW Vir system will be within the range of − 5.48 × 10³⁶ and − 4.16 × 10³³ (Rappaport et al. 1983), which makes it possible to explain the rate of constant period change in ETV.

The circumbinary objects from Model 6, in a wide range of a and e parameters are found to be unstable, unless the outer object has a > 10 au, for which no ETV solution is found. The stable configuration that Beuermann et al. (2012) found was shown to diverge from the more recent timing data in Baran et al. (2018). Therefore, circumbinary objects alone cannot explain the ETV periodicities. One should also keep in mind that the analytical expression of LiTE (Irwin 1959) does not consider the mutual interactions in a system that consists of more than two bodies. A better solution would be to analyze the ETV bysimultaneously solving n-body interactions between the bodies in the case of multiple circumbinary objects.

Another explanation for the apparent ETV trend might be the effect of the magnetic activity modulation of the cooler companion. Our initial sinusoidal fit to the ETV resulted in a period of 77 yr. By using the formulation of Applegate (1992), we calculate the magnetic field as B ≈ 39 kG. The luminosity variation of the secondary with L₂ = 0.006 L_⊙ should be ΔL = 0.011 L_⊙. Other than the amount of change being larger than the luminosity level of the secondary, such a change was not reported by any observer. Therefore, magnetic activity of the secondary can be ruled out at the level of period change calculated in this study. However, the possibility of post-common envelope evolution and/or a close orbital configuration of the binary affecting the magnetic activity of the cooler companion cannot be ruled out. We would like to encourage interested researchers to investigate the cooler companion and its effect on ETVs in HW Vir-like binaries.

The fit statistics and parameter uncertainties heavily depend on timing uncertainties. Mikulášek et al. (2014) questioned the reliability of the uncertainties from the widely used Kwee-van Woerden method (Kwee & van Woerden 1956) in determining the eclipse timings. From their comparison with the least-squares method, they concluded that the uncertainties from the Kwee-van Woerden method are systematically underestimated. Analyses based on underestimated uncertainties may result in the underestimation or overestimation of the parameter uncertainties depending on the model, if they are notcarefully examined. Therefore, we suggest carefully investigating error estimates on the mid-eclipse timings.

As Baran et al. (2018) showed, the primary star in the HW Vir system displays pulsations that are not detectable from the ground. Since many proposed circumbinary planets are in systems with components expected to show pulsations with different timescales and amplitudes, the effects of pulsations on the measurement of the eclipse timings should also be investigated, even when they only contribute to the noise budget, especially in ground-based observations.

Fig. 9 Stability analysis of Model 6 from bootstrap method in Sect. 5.2 of Model 6. See Fig. 7 for the description of the figure (see text for more details).

Acknowledgements

We thank TÜBİTAK National Observatory for a partial support in using T100 telescope with the project number 17BT100-1208 and 17BT100-378. We also thank Ankara University Kreiken Observatory (AUKR) for the observation time and all the student observers who helped in the observations of the target. E.M.E. acknowledges support from TÜBİTAK (2214-A, No. 1059B141800521). O.B. thanks The Scientific and Technological Council of Turkey (TÜBİTAK) for their support through the research grand 118F042. T.C.H. acknowledges staff at SOAO and LOAO observatories for assistance with the observations and fruitful discussions. T.C.H. acknowledges financial support from the National Research Foundation (NRF; No. 2019R1I1A1A01059609). A.C. acknowledges support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal.

Appendix A Mid-eclipse timings

Appendix B Parameter uncertainties for Model 6

Fig. B.1 Marginalized probability density distributions (P(x, y) ∝exp(−χ2∕2)) for the parameters of Model 6 as obtained from the χ2-grid method. The 68.3% confidence interval and median values are shown with vertical lines (see text for more details).

Fig. B.2 Histograms of parameters from the bootstrap method. ± 68.3% confidence limits and median values are shown with vertical lines.

Appendix C χ² surface grids for Model 6

Fig. C.1 χ2 surface plots for the pairs of parameters of Model 6. The vertical and horizontal dashed lines represent ± 1σ values.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

References

All Tables

All Figures

	Fig. 1 Radial velocity data along with synthetic curves from simultaneous light and radial velocity analysis. Combined radial velocity data of Wood & Saffer (1999) (white) and shifted Hilditch et al. (1996) (black) is above. Radial velocities from Edelmann (2008) is below. The white triangle marker represents the measurements for the secondary companion. See text for details.
In the text

	Fig. 2 Observed light curves with best-fit models and their residuals in various pass-bands.
In the text

	Fig. 3 Lomb-Scargle spectrum of the ETV data. Upper panel is for the whole dataset, and the lower panel is for the residuals from the highest frequency of the whole dataset. We plot the 1% FAP line for reference.
In the text

Fig. 4 ETV diagram of HW Vir with best-fit curves of all models. The ETV diagram is corrected by the results of the linear fit. Color-filled circles represents the primary mid-eclipse timings, while white-filled ones represent the secondaries. The black markers are for data from the literature, blue ones are for SWASP, green for Kepler K2, red for our own observations, and yellow markers represent seasonal binned data.

In the text

	Fig. 5 ETV diagram of HW Vir with best-fit curve of Model 6 and its residuals. Color-filled circles represent the primary mid-eclipse timings, while white-filled ones represents the secondaries. The black markers represent data from the literature, blue ones are for SWASP, green for Kepler K2, red for our own observations, and yellow markers represent seasonal binned data.
In the text

	Fig. 6 Lag plots corresponding the six models considered in this study. Lag plots (a) for Models 1–4 (red, green, blue, and yellow, in respective order), (b) for Model 5, and (c) for Model 6. Plots (b) and (c) show the Durbin-Watson (DW) statistic. We refer the reader to the text for more details.
In the text

	Fig. 7 Stability analysis of Model 6 in Sect. 4.2 (Table 6). The color scale represents the normalized stability index in log scale. The black circle represents the best-fit values of the outer orbit. The errorbars represent 1, 2, and 3σ uncertainties on the best-fit values (see text for more details).
In the text

	Fig. 8 Stability analysis of Model 6 from χ2 surface search in Sect. 5.1 of Model 6. See Fig. 7 for the description of the figure (see text for more details).
In the text

	Fig. 9 Stability analysis of Model 6 from bootstrap method in Sect. 5.2 of Model 6. See Fig. 7 for the description of the figure (see text for more details).
In the text

	Fig. B.1 Marginalized probability density distributions (P(x, y) ∝exp(−χ2∕2)) for the parameters of Model 6 as obtained from the χ2-grid method. The 68.3% confidence interval and median values are shown with vertical lines (see text for more details).
In the text

	Fig. B.2 Histograms of parameters from the bootstrap method. ± 68.3% confidence limits and median values are shown with vertical lines.
In the text

	Fig. C.1 χ2 surface plots for the pairs of parameters of Model 6. The vertical and horizontal dashed lines represent ± 1σ values.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text