11institutetext:  CFisUC, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal 22institutetext: Department of Astronomy and Space
Sciences, Ankara University, 06100 Ankara, Turkey 33institutetext: Institute
of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus
Copernicus University, Grudziadzka 5, 87-100 Torun, Poland. 44institutetext:
Department of Astronomy and Space Science, Chungnam National University, 34134
Daejeon, Republic of Korea 55institutetext: ASD, IMCCE, Observatoire de Paris,
PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France

# Revisiting the analysis of HW Vir eclipse timing data

I. A frequentist data modeling approach and a dynamical stability analysis

Ekrem. M. Esmer 1122 Özgür Baştürk, 22 Tobias C. Hinse Corresponding author:
T. C. Hinse (tchinse@gmail.com)3344 Selim O. Selam 22 Alexandre C. M. Correia
1155

(Received Jun 11, 2020; Accepted Feb 2, 2021)

###### Abstract

Context. HW Vir is a short-period binary that presents eclipse timing
variations. Circumbinary planets have been proposed as a possible explanation,
although the properties of the planets differ in each new study.

Aims. Our aim is to perform robust model selection methods for eclipse timing
variations (ETV) and error calculation techniques based on a frequentist
approach for the case of the HW Vir system.

Methods. We initially performed simultaneous light and radial velocity curve
analysis to derive the masses of the binary. We then analyzed the eclipse
timing variation of the system by fitting multiple models. To select the best
model, we searched the confidence levels for the best model by creating an
$\chi^{2}$ surface grid and bootstrap methods for each pair of parameters. We
searched for stable orbital configurations for our adopted ETV model.

Results. The masses of the binary are found as $0.413\pm 0.008\ M_{\sun}$ and
$0.128\pm 0.004\ M_{\sun}$. Under the assumption of two light time effects
superimposed on a secular change, the minimum masses of the circumbinary
objects are calculated as $25.0_{-2.2}^{+3.5}\ M_{Jup}$ and
$13.9_{-0.45}^{+0.60}\ M_{Jup}$. The projected semi-major axes are found to be
$7.8_{-1.0}^{+1.4}\ au$ and $4.56_{-0.22}^{+0.27}\ au$ in respective order. We
find that this configuration is unstable within a 3$\sigma$ range on the semi-
major axis and eccentricity of the outer circumbinary object.

###### Key Words.:

binaries: eclipsing – stars: individual: HW Vir – planetary systems –
subdwarfs – methods: data analysis

##  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_{\sun}$ and 0.12 $M_{\sun}$ while the
temperatures are 26000 K and 4700 K, respectively. From their single-lined
spectrum, Hilditch et al. (1996) derived the temperatures as 33000 $\pm$ 800 K
and 3700 $\pm$ 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_{2}$ = 275 $\pm$ 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 $\pm$ 0.08 $M_{\sun}$ and $0.15\pm 0.03\,M_{\odot}$ 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 $\mu
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_{\sun}$ for the
primary. However, the mass that they found is considerably smaller than the
canonical helium core ignition mass of $\sim$0.47 $M_{\sun}$ (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_{3}=19.2\ M_{Jup}$ and $M_{4}=8.5\
M_{Jup}$. Orbital periods of these planetary mass companions were given as
$P_{3}$ = 15.84 yr and $P_{4}$ = 9.08 yr with the eccentricities $e_{3}$ =
0.46 and $e_{4}$ = 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_{3}$ = 14.3$\ M_{Jup}$, $M_{4}$ = 65$\ M_{Jup}$). They found
that the orbits should be stable for more than $10^{8}$ yr with the
configuration that they proposed.

In Section 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, and
compared them, as a result, making use of the Durbin-Watson statistics (Durbin
& Watson, 1950), $\chi^{2},$ and $F$-tests. In Section 5, we provide the
details of our error analysis based on an $\chi^{2}$ 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 $BVRI$ filters. The majority of observations were
done through $R$ filters. All of the observations were made under clear to
thin cloud conditions with various moon phases. We used 2x2 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
AstroImageJ111https://www.astro.louisville.edu/software/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.

Telescope | Dates of Observations | Filters | Mean Phot. Errors (mag) | No. of Mid-Eclipse Times
---|---|---|---|---
| 2015 Feb 21, Mar 13, Apr 17, May 15, 26; |  |  |
AUKR - T35 | 2016 Feb 3, 7, 21, 22, Mar 9, Apr 1, 27; | BVRI | 9.2, 8.7, 3.5, 3.3 ($\times 10^{-3}$) | 38 (23 pri; 15 sec)
| 2017 Feb 2, 18, Apr 3, Jun 6, 12 |  |  |
| 2019 Mar 24 |  |  |
TUG - T100 | 2014 Feb 19, May 4; 2015 Apr 24 | R | 0.0020 | 7 (4 pri; 3 sec)
| 2018 May 3 |  |  |
SOAO | 2015 Mar 21, 22; 2016 Mar 21, Apr 1, 13 | R | 0.0054 | 14 (7 pri; 7 sec)
LOAO | 2016 Apr 28, 29, 30 | BVR | 9.2, 9.3, 17 ($\times 10^{-4}$) | 7 (4 pri; 3 sec)
Table 1: Log of photometric observations with ”sec” and ”pri” denoting
secondary and primary eclipses, respectively.

##  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.

Phase | $\textbf{RV}_{1}(\textbf{km\ s}^{-1})$ | $\textbf{RV}_{2}(\textbf{km\ s}^{-1})$
---|---|---
0.05 | -25 $\pm$3 |
0.11 | -52 $\pm$3 |
0.17 | -70 $\pm$3 |
0.24 | -81 $\pm$3 |
0.30 | -78 $\pm$3 | 205 $\pm$10
0.37 | -57 $\pm$3 | 184 $\pm$10
0.43 | -25 $\pm$3 | 126 $\pm$10
0.50 | 7 $\pm$3 |
0.56 | 35 $\pm$3 | -107 $\pm$10
0.63 | 63 $\pm$3 | -174 $\pm$10
0.69 | 84 $\pm$3 | -209 $\pm$10
0.76 | 85 $\pm$3 |
0.82 | 84 $\pm$3 |
0.85 | 72 $\pm$3 |
0.92 | 48 $\pm$3 |
0.98 | 28 $\pm$3 |
Table 2: Radial velocity data from Edelmann (2008). The values in second and
third column are the radial velocity measurements of the primary and the
secondary companion of HW Vir.

The semi-amplitudes of the radial velocities are $82.3\pm 4$ km/s and $83.0\pm
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\sigma$) Radial
velocity of the barycenter of the system ($V\gamma$) was given as $2.9\pm 3.1$
km/s and $-9.1\pm 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$\sigma$) difference in $V\gamma$ 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\gamma$ 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\gamma$, 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.

Figure 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.

We used PHOEBE v0.31222http://phoebe-project.org/1.0 (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_{1}$) to 28500 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_{1}$) to unity, and that of the
secondary ($g_{2}$) to 0.32. We left the remaining parameters (semi-major axis
($a$), mass ratio ($q$), systemic velocity ($V_{\gamma}$), orbital inclination
($i$), surface temperature of the secondary ($T_{2}$), surface potentials of
both of the components ($\Omega_{1,2}$), luminosity ($L_{1}$), and limb
darkening of the primary ($ld_{1}$) as free parameters. We assumed synchronous
rotation concerning the high projected rotational velocity of the sdB
component calculated as 74 $\pm$ 2 km s-1 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 $VRI$ 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 ($BVRI$)
separately, but this time fixing the physical parameters ($T_{1,2}$, $q$, $a$,
$i$, $\Omega_{1,2}$, $V_{\gamma}$) that we derived from $VRI$ 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_{2}$) 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_{2}$ to 0.32. Then, we realized that the albedo
$A_{2}$ and the decoupled luminosity of the secondary $L_{2}$ 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_{2}=1.0$), the other with a fixed
luminosity ($L_{2}$), to the same values derived in the $VRI$ 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{dX^{2}+dY^{2}}$. Table 4
clearly shows that the output absolute parameters are significantly different
from each other in terms of standard deviations ($\sigma$), 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 $\sim$0.3 and phase $\sim$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.

Figure 2: Observed light curves with best-fit models and their residuals in various pass-bands. Parameter | WS99+H96 | ED08 | Agreement
---|---|---|---
|  |  | Level ($\sigma$)
$i$ (∘) | 81.00(7) | 80.9(1) | 0.82
$T_{1}$ (K) | 28500* | 28500* | -
$T_{2}$ (K) | 3903(7) | 3909(9) | 0.53
$\Omega_{1}$ | 5.02(5) | 5.04(6) | 0.26
$\Omega_{2}$ | 2.89(2) | 3.09(2) | 7.07
$A_{1}$ | 1* | 1* | -
$A_{2}$ | 1* | 1* | -
$g_{1}$ | 1* | 1* | -
$g_{2}$ | 0.32* | 0.32* | -
$q=m_{2}/m_{1}$ | 0.310(4) | 0.357(4) | 8.31
Luminosities
$L_{1}$ / $L_{T}$ [V] | 0.998106(3) | 0.998107(4) | 0.2
$L_{1}$ / $L_{T}$ [R] | 0.996305(5) | 0.996307(7) | 0.23
$L_{1}$ / $L_{T}$ [I] | 0.98920(2) | 0.9892(2) | 0
$L_{2}$ / $L_{T}$ [V] | 0.001894 | 0.001893 | -
$L_{2}$ / $L_{T}$ [R] | 0.003695 | 0.003693 | -
$L_{2}$ / $L_{T}$ [I] | 0.01080 | 0.0108 | -
Table 3: Results of the simultaneous light curve and velocity curve analysis of HW Vir. Numbers in parentheses denotes the uncertainty on the last digit. The asterisk (*) denotes fixed parameter during the fitting process. $L_{2}$ / $L_{T}$ are not derived parameters ($L_{2}$ / $L_{T}$ = 1 - $L_{1}$ / $L_{T}$). See text for details. Parameter | WS99+H96 | ED08 | Agreement
---|---|---|---
|  |  | Level ($\sigma$)
$a$ ($R_{\sun}$) | 0.818(3) | 0.750(3) | 16.03
$M_{1}$ ($M_{\sun}$) | 0.413(8) | 0.307(6) | 10.6
$M_{2}$ ($M_{\sun}$) | 0.128(4) | 0.110(3) | 3.6
$R_{1}$ ($R_{\sun}$) | 0.175(2) | 0.161(3) | 3.9
$R_{2}$ ($R_{\sun}$) | 0.166(3) | 0.153(3) | 3.06
$L_{1}$ ($L_{\sun}$) | 18.0(5) | 15.3(5) | 3.82
$L_{2}$ ($L_{\sun}$) | 0.006(2) | 0.005(2) | 0.35
$M_{bol,1}$ (mag) | 1.61(3) | 1.79(4) | 3.6
$M_{bol,2}$ (mag) | 10.35(5) | 10.52(5) | 2.40
$\log g_{1}$ (cgs) | 5.570(3) | 5.512(6) | 8.65
$\log g_{2}$ (cgs) | 5.105(2) | 5.108(2) | 1.06
Table 4: Absolute parameters of HW Vir for two different radial velocity
datasets. Numbers in parentheses denotes the uncertainty on the last digit.
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 10.

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; Wood et al., 1993; Kilkenny et al., 1994; Gurol & Selam, 1994;
Kilkenny et al., 2000; Ogloza et al., 2000; Selam et al., 1999; Wood & Saffer,
1999; Çakirli & Devlen, 1999; Agerer & Hubscher, 2000; Agerer et al., 1999;
Kilkenny et al., 2003; Gurol et al., 2003; Lee et al., 2009; Blaettler &
Diethelm, 2000; Kiss et al., 2000; Agerer & Hubscher, 2002; İbanoǧlu et al.,
2004; Hubscher, 2005; Agerer & Hubscher, 2003; Kotkova & Wolf, 2006; Dvorak,
2005; Zejda et al., 2006; Nagai, 2006; Hubscher et al., 2005; Diethelm, 2005;
Dvorak, 2006; Qian et al., 2008; Dvorak, 2008; Nagai, 2007; Brát et al., 2007;
Dvorak, 2009; Beuermann et al., 2012; Nagai, 2009; Nelson, 2009; Brat et al.,
2009; Parimucha et al., 2009; Nagai, 2010; Brat et al., 2011; Diethelm, 2011;
Hoňková et al., 2013; Nagai, 2013, 2014; Kubicki, 2015; Basturk et al., 2014;
Hubscher, 2017; Kubicki, 2017; 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 $\sim$12000 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.

$Min\leavevmode\nobreak\ I=BJD_{TDB}\leavevmode\nobreak\
2,450,280.28472(4)+0.116719504(1)\leavevmode\nobreak\ E,$ (1)

where, $MinI$ 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 ($\approx$ 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 $\sigma_{i}=6.48\times 10^{-5}\ d=5.6\ s$ and assigned
this value as uncertainties.

###  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 Astropy333http://www.astropy.org package (Astropy
Collaboration et al., 2013; Price-Whelan et al., 2018) functions to detect the
periodicities. We found the highest peak in the period spectrum at
$f_{1}=1.20\times 10^{-4}\,\rm d^{-1}$ with an amplitude of 126.4 s and false
alarm probability (FAP) of $2.4\times 10^{-75}$. There are some other peaks
around $2.7\times 10^{-3}\,\rm d^{-1}$, 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_{2}=1.89\times 10^{-4}\,\rm d^{-1}$
with an amplitude of 93.3 s and FAP of $1.6\times 10^{-94}$ on the residual
spectrum (see Fig. 3). Since we only removed $f_{1}$ 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_{1}$ and $f_{2}$ corresponding to periods of
8320 d ($\sim$22.8 yr) and 5298 d ($\sim$14.5 yr), respectively. Considering
the total time span of the data ($\sim$13200 d), it should cover more than one
full cycle of the variation with $f_{1}$ and more than two cycles with
$f_{2}$. As a result, we decided to use these first two peaks in a first
attempt to model the observed ETV.

Figure 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.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, $\beta$, to be added if needed. Therefore ETV models
($T_{C}$) with $\beta$ would have a form of

$T_{C}=(T_{0}+\Delta T_{0})+(P_{0}+\Delta P_{0})E+\beta E^{2}+\tau(E).$ (2)

Here, $\Delta T_{0}$ and $\Delta P_{0}$ are the correction terms on the
reference mid-eclipse time ($T_{0}$) and the orbital period of the binary
($P_{0}$), respectively, and $\tau(E)$ is a LiTE model(s). The models without
the quadratic term will have a form similar to the one above, but without the
$\beta E^{2}$ 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\approx 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^{2}$, $\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_{bin}+M_{3})^{2}}=\Bigg{[}\frac{173.15A}{\sqrt[]{1-e_{3}^{2}\cos^{2}\omega_{3}}}\Bigg{]}^{3}\frac{1}{P_{3}^{2}}.$
(3)

Here, $M_{3}$ is the mass, $i_{3}$ is the orbital inclination, $e_{3}$ is the
eccentricity, $P_{3}$ is the orbital period, $\omega_{3}$ 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_{3}$ 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.

Model Component | Analytic Expression
---|---
Linear (Lin) | $\Delta PE+\Delta T_{ref}$
Quadratic (Quad) | $\beta E^{2}$
Light-time effect1 (LiTE) | $\frac{A}{\sqrt[]{1-e^{2}\cos^{2}\omega}}$
| $\times\leavevmode\nobreak\ \big{[}\frac{1-e^{2}}{1+e\cos
v}\sin(v+\omega)+e\sin\omega\big{]}$
Model (No. of param.) | Short Hand Notation
Lin (2) | Model 1
Lin+Quad (3) | Model 2
Lin+LiTE (7) | Model 3
Lin+Quad+LiTE (8) | Model 4
Lin+ 2$\times$LiTE (12) | Model 5
Lin+Quad+ 2$\times$LiTE (13) | Model 6
Table 5: Analytical expressions of the model components and model list. We
note that short-hand notations are ordered by the increased number of
parameters. (1 Irwin 1959)

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 $\nu$ 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^{6}d\approx 2750yr$) and extremely
eccentric orbit, which coincidentally should have passed from its superior
conjunction relative to Earth within the time frame of the dataset ($\sim$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 10499 d (29 yr) for binned dataset with almost the same
eccentricity ($e\sim 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
$\sim$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, $\omega$, 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.

Figure 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 terms of fit statistics and best-fit parameter values, one of the
satisfactory ETV models has been achieved for the case of Model 6 (Fig. 4 and
Fig. 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 $\omega$ 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 $\nu$ in order
to have an $\chi_{\nu}^{2}$ of unity. However, for seasonal binned data, the
ratio $N$ / $\nu$ $\approx$ 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 $\sim 0.2kms^{-1}$ with the assumption of $i=90^{\circ}$.
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\times 10^{4}$ d-1, which corresponds to a period of 1690.7 d and
a small FAP value of $6\times 10^{-4}$. However, the amplitude of the peak is
3.97 s, which is less than the mean standard deviation ($\sim$5.6 s) of the
timing data and less than the RMS for Model 6 fit to the original data
($\sim$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.

Figure 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. | Parameter | Unit | Original Dataset | Seasonal Binned Dataset
---|---|---|---|---
Model 1 | $\chi_{\nu}^{2}$ |  | $89.28$ | $125.43$
RMS | s | $52.764$ | $60.674$
$T_{0}$ | $BJD_{TDB}$ | $2,450,280.28471\pm 4\times\>10^{-5}$ | $2,450,280.2850\pm 2\times 10^{-4}$
$P_{0}$ | $d$ | $0.116719504\pm 1\times\>10^{-9}$ | $0.116719502\pm 4\times 10^{-9}$
Model 2 | $\chi_{\nu}^{2}$ |  | $79.52$ | $118.26$
RMS | s | $49.747$ | $57.924$
$T_{0}$ | $BJD_{TDB}$ | $2,450,280.28489\pm 5\times 10^{-5}$ | $2,450,280.28513\pm 2\times\>10^{-4}$
$P_{0}$ | $d$ | $0.116719513\pm 1\times\>10^{-9}$ | $0.116719511\pm 6\times\>10^{-9}$
$\beta$ | $d\ cycle^{-2}$ | $-2.71\times 10^{-13}\pm 3\times 10^{-14}$ | $-2.3\times 10^{-13}\pm 1\times 10^{-13}$
Model 3 | $\chi_{\nu}^{2}$ |  | $5.58$ | $9.34$
RMS | s | $13.126$ | $15.114$
$T_{0}$ | $BJD_{TDB}$ | $2,450,280.3\pm 0.1$ | $2,450,280.3\pm 0.1$
$P_{0}$ | $d$ | $0.11671961\pm 3\times 10^{-8}$ | $0.11671959\pm 1\times 10^{-8}$
$T_{0,LiTE}$ | $BJD_{TDB}$ | $2,449,352\pm 104$ | $3,595,768\pm 8\times 10^{6}$
$P_{LiTE}$ | $d$ | $1.0\times 10^{6}\pm 7\times 10^{6}$ | $1.2\times 10^{6}\pm 8\times 10^{6}$
$e$ | - | $0.98\pm 0.1$ | $0.98\pm 0.1$
$A_{LiTE}$ | $s$ | $2530\pm 1.2\times 10^{4}$ | $2527\pm 1\times 10^{4}$
$\omega$ | ∘ | $133\pm 7$ | $129\pm 11$
$f(m)$ | $M_{Jup}$ | $0.04\pm 0.9$ | $0.03\pm 0.6$
$M\sin\\!i$ | $M_{Jup}$ | $25\pm 170$ | $21\pm 140$
$a\sin\\!i$ | $au$ | $156\pm 1300$ | $170\pm 1400$
Model 4 | $\chi_{\nu}^{2}$ |  | $5.26$ | $6.69$
RMS | s | $12.726$ | $12.533$
$T_{0}$ | $BJD_{TDB}$ | $2,450,280.28554\pm 2\times 10^{-5}$ | $2,450,280.28546\pm 8\times 10^{-5}$
$P_{0}$ | $d$ | $0.1167195268\pm 8\times 10^{-10}$ | $0.11671954\pm 1\times 10^{-8}$
$\beta$ | $d\ cycle^{-2}$ | $-7.21\times 10^{-13}\pm 1\times 10^{-14}$ | $-7.76\times 10^{-13}\pm 7\times 10^{-14}$
$T_{0,LiTE}$ | $BJD_{TDB}$ | $2,449,127\pm 51$ | $2,449,662\pm 281$
$P_{LiTE}$ | $d$ | $9399\pm 103$ | $10499\pm 1708$
$e$ | - | $0.52\pm 0.03$ | $0.51\pm 0.08$
$A_{LiTE}$ | $s$ | $101\pm 2$ | $107\pm 14$
$\omega$ | ∘ | $115\pm 2$ | $145\pm 15$
$f(m)$ | $M_{Jup}$ | $0.0140\pm 0.0008$ | $0.017\pm 0.09$
$M\sin\\!i$ | $M_{Jup}$ | $16.8\pm 0.4$ | $17.8\pm 3.2$
$a\sin\\!i$ | $au$ | $7.0\pm 0.2$ | $7.5\pm 1.7$
Model 5 | $\chi_{\nu}^{2}$ |  | $1.04$ | $0.68$
RMS | s | $5.625$ | $3.640$
$T_{0}$ | $BJD_{TDB}$ | $2,450,280.28456\pm 4\times 10^{-5}$ | $2,450,280.28458\pm 6\times 10^{-5}$
$P_{0}$ | $d$ | $0.1167195065\pm 4\times 10^{-10}$ | $0.1167195059\pm 8\times 10^{-10}$
$T_{0,LiTE}$ | $BJD_{TDB}$ | $2,446,408\pm 448$ | $2,453,385\pm 670$ | $2,446,482\pm 1885$ | $2,453,621\pm 3199$
$P_{LiTE}$ | $d$ | $7800\pm 253$ | $7448\pm 369$ | $7645\pm 1200$ | $7442\pm 1368$
$e$ | - | $0.19\pm 0.04$ | $0.18\pm 0.02$ | $0.18\pm 0.2$ | $0.17\pm 0.02$
$A_{LiTE}$ | $s$ | $480\pm 801$ | $467\pm 801$ | $786\pm 9761$ | $773\pm 9761$
$\omega$ | ∘ | $242\pm 36$ | $20\pm 34$ | $239\pm 158$ | $33\pm 166$
$f(m)$ | $M_{Jup}$ | $2\pm 10$ | $2\pm 11$ | $10\pm 350$ | $10\pm 400$
$M\sin\\!i$ | $M_{Jup}$ | $97\pm 180$ | $98\pm 190$ | $173\pm 250$ | $175\pm 260$
$a\sin\\!i$ | $au$ | $6\pm 14$ | $5\pm 14$ | $5\pm 100$ | $5\pm 100$
Model 6 | $\chi_{\nu}^{2}$ |  | $1.07$ | $0.68$
RMS | s | $5.699$ | $3.556$
$T_{0}$ | $BJD_{TDB}$ | $2,450,280.28569\pm 8\times 10^{-5}$ | $2,450,280.2855\pm 1\times 10^{-4}$
$P_{0}$ | $d$ | $0.116719556\pm 3\times 10^{-9}$ | $0.116719549\pm 5\times 10^{-9}$
$\beta$ | $d\ cycle^{-2}$ | $-1.40\>\times\>10^{-12}\ \pm\ 9\ \times\ 10^{-14}$ | $-1.17\times 10^{-12}\pm 1\times 10^{-13}$
$T_{0,LiTE}$ | $BJD_{TDB}$ | $2,444,369\pm 305$ | $2,452,958\pm 68$ | $2,444,744\pm 490$ | $2,453,024\pm 166$
$P_{LiTE}$ | $d$ | $11391\pm\ 299$ | $4958\pm\ 30$ | $11006\pm 478$ | $5040\pm 58$
$e$ | - | $0.45\pm 0.01$ | $0.27\pm 0.02$ | $0.50\pm 0.03$ | $0.25\pm 0.05$
$A_{LiTE}$ | $s$ | $152\pm 9$ | $54\pm 1$ | $133\pm 13$ | $55\pm 2$
$\omega$ | ∘ | $0\pm\ 2$ | $13\pm 5$ | $358\pm 4$ | $13\pm 13$
$f(m)$ | $M_{Jup}$ | $0.043\pm 0.008$ | $0.0079\pm 0.0005$ | $0.034\pm 0.010$ | $0.0082\pm 0.0011$
$M\sin\\!i$ | $M_{Jup}$ | $24.6\pm 1.5$ | $13.89\pm 0.33$ | $22.8\pm 2.4$ | $14.0\pm 0.6$
$a\sin\\!i$ | $au$ | $7.8\pm 0.7$ | $4.56\pm 0.16$ | $7.7\pm 1.1$ | $4.61\pm 0.30$
Table 6: Parameter values and formal uncertainties for the ETV models of HW
Vir.

###  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
$\sim$$N(0,\sigma^{2}$), 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,

$\textit{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$ $\approx$ 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 $\beta$ 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.

Figure 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: $\chi^{2}$ test

We compared the $p$-value, $p(\chi^{2},\nu)$, of each model calculated from an
$\chi^{2}$ distribution for a given number of degrees of freedom, with a
critical probability indicated by the significance level ($\alpha$) to test
the null-hypothesis stating that the model under investigation fits the data
well ($p<\alpha$). In this case, the probability of finding the observed
$\chi^{2}$ is too low at the given significance level and cannot be explained
by a random process. When $\alpha<p\leq 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 $\chi^{2}$ is overestimated, resulting in an unrealistically
small value of $\chi^{2}$ (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 $\chi^{2}$ is close to the mean of the
($\chi^{2},\nu$) distribution: $\chi^{2}\approx\nu$ with probability
$p(\chi^{2},\nu)\approx 0.5$. In the case of considerably large $\chi^{2}$
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 $\chi^{2}$ to a
critical probability, $\alpha$. If the $\chi^{2}$ is considerably less than
the mean value of the $\chi^{2}$ distribution, the probability of the
$\chi^{2}$ 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
$\alpha$ 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(\chi^{2})<0.01,$ we judge very strong evidence
against $H_{0}$; if $0.01<p<0.05,$ we judge strong evidence against $H_{0}$;
if $0.05<p<0.10,$ we judge some weak evidence against $H_{0}$; and if
$p>0.10,$ we judge little or no evidence against $H_{0}$.

Table 7 shows the calculated $\chi^{2}$ for each model that we considered in
this study. To calculate the probabilities $p(\chi^{2},\nu)$ for each model,
we made use of our own PYTHON code. For Models 1 - 4, the corresponding
$p$-values are smaller than $10^{-200}$, 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 ($\alpha$ = 0.1), we see no evidence
against Models 5 and 6, since $p_{model}>\alpha$. In a relative sense, we
therefore interpret this result as a non-rejection of the null-hypothesis for
Models 5 and 6.

| $\nu$ | $\chi^{2}$ | $\chi_{\nu}^{2}$ |  $p$-value
---|---|---|---|---
Model 1 | $507$ | $45263.8$ | $89.28$ | $<10^{-300}$
Model 2 | $506$ | $40236.1$ | $79.52$ | $<10^{-300}$
Model 3 | $502$ | $2801.1$ | $5.58$ | $<10^{-300}$
Model 4 | $501$ | $2633.0$ | $5.26$ | $2.2\times 10^{-285}$
Model 5 | $497$ | $514.5$ | $1.04$ | $0.285$
Model 6 | $496$ | $528.0$ | $1.07$ | $0.155$
Table 7: Results of the $\chi^{2}$-test ($\alpha$ = 0.1; see text for more
details).

###  4.5 Model selection: $F$-Test

In the previous section, we determined the most likely model by comparing the
minimum $\chi^{2}$ of each model with the $\chi^{2}$ probability density
distribution for $\nu$ degrees of freedom under the assumption of the null-
hypothesis. We find an $\chi^{2}$ value for the Model 5, indicating the
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 $\chi^{2}$ 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 $\chi^{2}$ 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

$\textit{F}=\frac{(\chi^{2}_{1}-\chi^{2}_{2})/(\nu_{2}-\nu_{1})}{\chi^{2}_{2}/\nu_{2}},$
(5)

where $\nu_{i}$ are the corresponding degrees of freedom, and $\chi^{2}_{i}$
are the $\chi^{2}$ 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 $\chi^{2}$ and $\nu$ values from the best-fit results from Sect. 4.2
(see Table 7), for each model pair ordered by the corresponding $\nu$. We then
calculated the probability of F-statistics in an $F$-distribution defined by
the two degrees of freedom, $P(F;v_{2}-v_{1},v_{2})$. Similar to the
$\chi^{2}$ test, we chose the critical probability, $\alpha$, as 0.1 (10%),
and the condition of rejecting the null-hypothesis is in the cases of
$p<\alpha$, while if $p>\alpha$, 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 $\alpha$ of 0.1, thus the 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 the fit, 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^{\circ}$, 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.

$F$-test between |  $F$-value |  $p$-value
---|---|---
Model 1 - Model 2 | $63.2$ | $1.2\times 10^{-14}$
Model 2 - Model 3 | $1677.3$ | $7.8\times 10^{-289}$
Model 3 - Model 4 | $32.0$ | $2.6\times 10^{-8}$
Model 4 - Model 5 | $511.7$ | $1.2\times 10^{-174}$
Model 5 - Model 6 | $-12.8$ | $1$
Table 8: Results of the $F$-test (see text for more details).

##  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 independent algorithmic methods as part of a comparative
analysis study. The first method explores the projected $\chi^{2}$ 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 $\chi^{2}$ surface search

Parameter uncertainties can be determined by exploring the projected
$\chi^{2}$ surface around $\chi^{2}_{min}$. 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^{2}_{min}$ \+ 1 is met. This method is
then applied to two-parameter projections of $\chi^{2}$ 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 $\chi^{2}$ value while letting the
remaining parameters to vary freely. The upper and lower standard deviations
were then determined from $\chi^{2}_{min}$ \+ 1. Any parameter correlation is
propagated and accounted for in the calculated confidence intervals. For
complicated correlations, the $\chi^{2}$ 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 $\chi^{2}$ surface
was $75\leavevmode\nobreak\ \times\leavevmode\nobreak\ 75,$ and we calculated
a total of 438750 fits. For each fit, we calculated and stored the resulting
$\chi^{2}$ corresponding to a given parameter pair. To map out the $\chi^{2}$
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 $\chi^{2}$ 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^{2}_{min}$. For each grid point, we stored the resulting $\chi^{2}$
value.

We converted the resulting $\chi^{2}$ surfaces to a grid of probability
density as $P(x,y)\propto exp(-\chi^{2}/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 $\int 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(\alpha,\beta,...)$, where ($\alpha,\beta,...$)
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(\alpha+\sigma_{\alpha},\beta,...)-f(\alpha,\beta,...)$. This procedure and
calculation of differences were also determined for the other parameters
$f(\alpha,\beta+\sigma_{\beta},...)-f(\alpha,\beta,...)$. The final
uncertainty for X is then obtained from the quartiles as

$\sigma_{X}^{2}=[f(\alpha\pm\sigma_{\alpha},\beta,...)-f(\alpha,\beta,...)]^{2}\\\
+[f(\alpha,\beta\pm\sigma_{\beta},...)-f(\alpha,\beta,...)]^{2}+...$

Since the uncertainties calculated by the $\chi^{2}$ 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 of error 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 $\chi^{2}$ 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. 10 in the appendix, and the
values that correspond to the 68.3% confidence levels are shown in Table 9.
$\chi^{2}$ 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 derived uncertainties 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\times 10^{4}$, $2.5\times
10^{4}$, $1\times 10^{5},$ and $2\times 10^{5}$ random samples and plotted the
resulting distributions. We found no difference between $1\times 10^{5}$ and
$2\times 10^{5}$ bootstrap samples. Clearly, choosing $1\times 10^{4}$ 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. 11 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.

| Parameter | Unit |  $\chi^{2}$ Surface Search | Bootstrap | Agreement Level ($\sigma$)
---|---|---|---|---|---
Lin | $T_{0}$ | $BJD_{TDB}$ | $2,450,280.2857\ _{-1.2\times 10^{-4}}^{+1.5\times 10^{-4}}$ | $2,450,280.2857\ _{-1.3\times 10^{-4}}^{+2.2\times 10^{-4}}$ | $0.068$
$P_{0}$ | $d$ | $0.116719556\ _{-7.4\times 10^{-9}}^{+7.3\times 10^{-9}}$ | $0.116719556\ _{-4.45\times 10^{-9}}^{+4.73\times 10^{-9}}$ | $0.001$
Quad | $\beta$ | $d/cycle^{2}$ |  $-1.42\times 10^{-12}$$\ {}_{-1.9\times 10^{-13}}^{+1.6\times 10^{-13}}$ |  $-1.46\times 10^{-12}$$\ {}_{-2.02\times 10^{-13}}^{+1.35\times 10^{-13}}$ | $0.171$
$LiTE_{1}$ | $T_{0,1}$ | $BJD_{TDB}$ | $2,444,243\ _{-776}^{+407}$ | $2,444,230\ _{-798}^{+266}$ | $0.012$
$P_{1}$ | $d$ | $11494\ _{-483}^{+691}$ | $11310\ _{-356}^{+801}$ | $0.174$
$e_{1}$ | - | $0.450\ _{-0.011}^{+0.011}$ | $0.448\ _{-0.020}^{+0.021}$ | $0.091$
$A_{1}$ | $s$ | $155\ _{-14}^{+22}$ | $154\ _{-13}^{+20}$ | $0.064$
$\omega_{1}$ | ∘ | $359\ _{-2}^{+2}$ | $0\ _{-5}^{+5}$ | $0.070$
| $f(m)_{1}$ | $M_{jup}$ | $0.045\ _{-0.018}^{+0.022}$ | $0.045\ _{-0.011}^{+0.021}$ | $0.012$
| $M_{1}\sin\\!i$ | $M_{jup}$ | $25.1\ _{-2.5}^{+3.7}$ | $25.0\ _{-2.2}^{+3.5}$ | $0.013$
| $a_{1}\sin\\!i$ | $au$ | $7.9\ _{-1.1}^{+1.5}$ | $7.8\ _{-1.0}^{+1.4}$ | $0.041$
$LiTE_{2}$ | $T_{0,2}$ | $BJD_{TDB}$ | $2,452,958\ _{-80}^{+78}$ | $2,452,939\ _{-107}^{+101}$ | $0.137$
$P_{2}$ | $d$ | $4963\ _{-31}^{+32}$ | $4949\ _{-68}^{+74}$ | $0.167$
$e_{2}$ | - | $0.271\ _{-0.022}^{+0.022}$ | $0.279\ _{-0.025}^{+0.026}$ | $0.212$
$A_{2}$ | $s$ | $54\ _{-1}^{+1}$ | $54\ _{-1}^{+2}$ | $0.066$
$\omega_{2}$ | ∘ | $13\ _{-6}^{+6}$ | $12\ _{-8}^{+8}$ | $0.139$
| $f(m)_{2}$ | $M_{jup}$ | $0.00798\ _{-0.00048}^{+0.00051}$ | $0.00802\ _{-0.00070}^{+0.00103}$ | $0.029$
| $M_{2}\sin\\!i$ | $M_{jup}$ | $13.9\ _{-0.3}^{+0.3}$ | $13.9\ _{-0.45}^{+0.60}$ | $0.029$
| $a_{2}\sin\\!i$ | $au$ | $4.57\ _{-0.16}^{+0.16}$ | $4.56\ _{-0.22}^{+0.27}$ | $0.026$
Table 9: Results of best-fit parameters and associated uncertainties for Model
6 as determined from the $\chi^{2}$ surface search and bootstrap method.
Further details are given in the text. The agreement level is calculated as
$ABS(X-Y)$ / $\sqrt{dX^{2}+dY^{2}}$. Results of an $\chi^{2}$ surface search
and bootstrap search with their agreements in units of $\sigma$.

##  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 $\chi^{2}$ 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

$\textit{$D$}=\frac{\lvert n_{1}-n_{2}\rvert}{n_{1}},$ (6)

where $n_{1}$ is the frequency of the mean motion calculated from the first
half of each integration, and $n_{2}$ is the same for the second half. In the
case of regular motion, the normalized stability index should be $D<10^{-6}$
(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
$\times 10^{6}$ d. We calculated $n_{1}$ and $n_{2}$ 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\times 10^{4}$ 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^{\circ}$. 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^{-6}$) around $a\sim 10.5\ au$, $e\sim
0.12$ and another one around $a>12.2\ au$, $e<0.1$ (Figs. 7, 8, and 9). None
of the orbital configurations of Model 6 are in stable regions within a
3$\sigma$ 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 Fig. 7, 8 and 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.

Figure 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 $\sigma$ uncertainties on the best-fit values (see text for more
details). Figure 8: Stability analysis of Model 6 from $\chi^{2}$ surface
search in Sect. 5.1 of Model 6. See Fig. 7 for the description of the figure
(see text for more details). Figure 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).

##  7 Conclusions

From the simultaneous light and velocity curve analysis in Sect. 3, mass of
the primary is found to be 0.413 $M_{\sun}$ and 0.307 $M_{\sun}$ for the
WS99+H96 and ED08 datasets, respectively. Both of the masses are below the
canonical helium ignition mass of $\sim$0.47 $M_{\sun}$. 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_{\sun}$. 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_{\sun}$
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_{\sun}$. 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_{\sun}$ and 0.110 $M_{\sun}$ in
our analysis based on the WS99+H96 and ED08 datasets, respectively, for the
mass of the secondary. Literature values for the secondary mass change between
0.1 $M_{\sun}$ from Baran et al. (2018) and 0.15 $M_{\sun}$ 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 $\sim$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\approx 0.8\ R_{\sun}$) observed under rather
high inclination ($i\approx 81^{o}$).

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 $\beta$ coefficient of Model 6 in the ETV analysis corresponds to a
constant period change with a linear rate of $-9.16\times 10^{-9}d\ yr^{-1}$.
The angular momentum change can be derived (e.g., Brinkworth et al., 2006) as
$dJ/dt=-3.25\times 10^{35}erg$. The 3$\sigma$ confidence interval of the
angular momentum change should be within the range of $-4.60\times 10^{35}erg$
and $-2.36\times 10^{35}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\times
10^{32}erg$. The angular momentum loss calculated for the $\beta$ 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\times
10^{36}$ and $-4.16\times 10^{33}$ (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 by simultaneously 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\approx 39\ kG$. The
luminosity variation of the secondary with $L_{2}=0.006\ L_{\odot}$ should be
$\Delta L=0.011\ L_{\odot}$. 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 not carefully 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.

###### 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.

## References

  * Agerer et al. (1999) Agerer, F., Dahm, M., & Hubscher, J. 1999, Information Bulletin on Variable Stars, 4712, 1
  * Agerer & Hubscher (2000) Agerer, F. & Hubscher, J. 2000, Information Bulletin on Variable Stars, 4912, 1
  * Agerer & Hubscher (2002) Agerer, F. & Hubscher, J. 2002, Information Bulletin on Variable Stars, 5296, 1
  * Agerer & Hubscher (2003) Agerer, F. & Hubscher, J. 2003, Information Bulletin on Variable Stars, 5484, 1
  * Almeida et al. (2020) Almeida, L. A., Pereira, E. S., Borges, G. M., et al. 2020, MNRAS, 497, 4022
  * Applegate (1992) Applegate, J. H. 1992, ApJ, 385, 621
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bahar et al. (2015) Bahar, E., Şenavcı, H. V., & Baştürk, Ö. 2015, in Astronomical Society of the Pacific Conference Series, Vol. 496, Living Together: Planets, Host Stars and Binaries, ed. S. M. Rucinski, G. Torres, & M. Zejda, 288
  * Baran et al. (2018) Baran, A. S., Østensen, R. H., Telting, J. H., et al. 2018, MNRAS, 481, 2721
  * Baran et al. (2015) Baran, A. S., Zola, S., Blokesz, A., Østensen, R. H., & Silvotti, R. 2015, A&A, 577, A146
  * Basturk et al. (2014) Basturk, O., Bahar, E., Senavci, H. V., et al. 2014, Information Bulletin on Variable Stars, 6125, 1
  * Berger & Fringant (1980) Berger, J. & Fringant, A.-M. 1980, A&A, 85, 367
  * Beuermann et al. (2012) Beuermann, K., Dreizler, S., Hessman, F. V., & Deller, J. 2012, A&A, 543, A138
  * Beuermann et al. (2010) Beuermann, K., Hessman, F. V., Dreizler, S., et al. 2010, A&A, 521, L60
  * Bevington & Robinson (2003) Bevington, P. R. & Robinson, D. K. 2003, Data reduction and error analysis for the physical sciences; 3rd ed. (New York, NY: McGraw-Hill)
  * Blaettler & Diethelm (2000) Blaettler, E. & Diethelm, R. 2000, Bulletin der Bedeckungsveraenderlichen-Beobachter der Schweizerischen Astronomischen Gesellschaft, 122, 1
  * Brat et al. (2009) Brat, L., Trnka, J., Lehky, M., et al. 2009, Open European Journal on Variable Stars, 107, 1
  * Brat et al. (2011) Brat, L., Trnka, J., Smelcer, L., et al. 2011, Open European Journal on Variable Stars, 137, 1
  * Brát et al. (2007) Brát, L., Zejda, M., & Svoboda, P. 2007, Open European Journal on Variable Stars, 0074, 1
  * Brinkworth et al. (2006) Brinkworth, C. S., Marsh, T. R., Dhillon, V. S., & Knigge, C. 2006, MNRAS, 365, 287
  * Çakirli & Devlen (1999) Çakirli , Ö. & Devlen, A. 1999, Astronomy and Astrophysics Supplement Series, 136, 27
  * Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793
  * Collins et al. (2017) Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77
  * Correia et al. (2010) Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, A&A, 511, A21
  * Correia et al. (2005) Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751
  * Couetdic et al. (2010) Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10
  * Diethelm (2005) Diethelm, R. 2005, Information Bulletin on Variable Stars, 5653, 1
  * Diethelm (2011) Diethelm, R. 2011, Information Bulletin on Variable Stars, 5992, 1
  * Durbin & Watson (1950) Durbin, J. & Watson, G. S. 1950, Biometrika, 37, 409
  * Dvorak (2005) Dvorak, S. W. 2005, Information Bulletin on Variable Stars, 5603, 1
  * Dvorak (2006) Dvorak, S. W. 2006, Information Bulletin on Variable Stars, 5677, 1
  * Dvorak (2008) Dvorak, S. W. 2008, Information Bulletin on Variable Stars, 5814, 1
  * Dvorak (2009) Dvorak, S. W. 2009, Information Bulletin on Variable Stars, 5870, 1
  * Eastman et al. (2010) Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935
  * Edelmann (2008) Edelmann, H. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 392, Hot Subdwarf Stars and Related Objects, ed. U. Heber, C. S. Jeffery, & R. Napiwotzki, 187
  * Efron & Tibshirani (1986) Efron, B. & Tibshirani, R. 1986, Statistical Science, 1
  * Gastineau & Laskar (2011) Gastineau, M. & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194
  * Getley et al. (2017) Getley, A. K., Carter, B., King, R., & O’Toole, S. 2017, MNRAS, 468, 2932
  * Gurol et al. (2003) Gurol, B., Gurdemir, L., Caglar, A., et al. 2003, Information Bulletin on Variable Stars, 5443, 1
  * Gurol & Selam (1994) Gurol, B. & Selam, S. 1994, Information Bulletin on Variable Stars, 4109, 1
  * Han et al. (2002) Han, Z., Podsiadlowski, P., Maxted, P. F. L., Marsh, T. R., & Ivanova, N. 2002, MNRAS, 336, 449
  * Heber (2009) Heber, U. 2009, ARA&A, 47, 211
  * Heber (2016) Heber, U. 2016, PASP, 128, 082001
  * Hilditch et al. (1996) Hilditch, R. W., Harries, T. J., & Hill, G. 1996, MNRAS, 279, 1380
  * Hinse et al. (2014) Hinse, T. C., Horner, J., Lee, J. W., et al. 2014, A&A, 565, A104
  * Horner et al. (2012) Horner, J., Hinse, T. C., Wittenmyer, R. A., Marshall, J. P., & Tinney, C. G. 2012, MNRAS, 427, 2812
  * Hoňková et al. (2013) Hoňková, K., Juryšek, J., Lehký, M., et al. 2013, Open European Journal on Variable Stars, 160, 1
  * Hubscher (2005) Hubscher, J. 2005, Information Bulletin on Variable Stars, 5643, 1
  * Hubscher (2017) Hubscher, J. 2017, Information Bulletin on Variable Stars, 6196, 1
  * Hubscher et al. (2005) Hubscher, J., Paschke, A., & Walter, F. 2005, Information Bulletin on Variable Stars, 5657, 1
  * Hughes & Hase (2010) Hughes, I. G. & Hase, T. P. A. 2010, Measurements and their Uncertainties: A Practical Guide to Modern Error Analysis (Oxford University Press)
  * İbanoǧlu et al. (2004) İbanoǧlu, C., Çakırlı, Ö., Taș, G., & Evren, S. 2004, A&A, 414, 1043
  * Irwin (1959) Irwin, J. B. 1959, AJ, 64, 149
  * Juryšek et al. (2017) Juryšek, J., Hoňková, K., Šmelcer, L., et al. 2017, Open European Journal on Variable Stars, 179, 1
  * Kaplan (2010) Kaplan, D. L. 2010, ApJ, 717, L108
  * Kilkenny et al. (1991) Kilkenny, D., Harrop-Allin, M., & Marang, F. 1991, Information Bulletin on Variable Stars, 3569, 1
  * Kilkenny et al. (2000) Kilkenny, D., Keuris, S., Marang, F., et al. 2000, The Observatory, 120, 48
  * Kilkenny et al. (1994) Kilkenny, D., Marang, F., & Menzies, J. W. 1994, MNRAS, 267, 535
  * Kilkenny et al. (2003) Kilkenny, D., van Wyk, F., & Marang, F. 2003, The Observatory, 123, 31
  * Kiss et al. (2000) Kiss, L. L., Csák, B., Szatmáry, K., Furész, G., & Sziládi, K. 2000, A&A, 364, 199
  * Kotkova & Wolf (2006) Kotkova, L. & Wolf, M. 2006, Information Bulletin on Variable Stars, 5676, 1
  * Kubicki (2015) Kubicki, D. 2015, Information Bulletin on Variable Stars, 6133, 1
  * Kubicki (2017) Kubicki, D. 2017, Information Bulletin on Variable Stars, 6232, 1
  * Kwee & van Woerden (1956) Kwee, K. K. & van Woerden, H. 1956, Bull. Astron. Inst. Netherlands, 12, 327
  * Lanza et al. (1998) Lanza, A. F., Rodono, M., & Rosner, R. 1998, MNRAS, 296, 893
  * Laskar (1990) Laskar, J. 1990, Icarus, 88, 266
  * Laskar (1993) Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257
  * Lee et al. (2009) Lee, J. W., Kim, S.-L., Kim, C.-H., et al. 2009, AJ, 137, 3181
  * Levenberg (1944) Levenberg, K. 1944, Quarterly Journal of Applied Mathmatics, II, 164
  * Lohr et al. (2014) Lohr, M. E., Norton, A. J., Anderson, D. R., et al. 2014, A&A, 566, A128
  * Mandel (1964) Mandel, J. 1964, The Statistical Analysis of Experimental Data (Interscience Publishers)
  * Marang & Kilkenny (1989) Marang, F. & Kilkenny, D. 1989, Information Bulletin on Variable Stars, 3390, 1
  * Marquardt (1963) Marquardt, D. W. 1963, SIAM Journal on Applied Mathematics, 11, 431
  * Menzies & Marang (1986) Menzies, J. W. & Marang, F. 1986, in IAU Symposium, Vol. 118, Instrumentation and Research Programmes for Small Telescopes, ed. J. B. Hearnshaw & P. L. Cottrell, 305
  * Mikulášek et al. (2014) Mikulášek, Z., Chrastina, M., Liška, J., et al. 2014, Contributions of the Astronomical Observatory Skalnate Pleso, 43, 382
  * Nagai (2006) Nagai, K. 2006, VSOLJ Variable Star Bulletin, 44, 1
  * Nagai (2007) Nagai, K. 2007, VSOLJ Variable Star Bulletin, 45, 1
  * Nagai (2009) Nagai, K. 2009, VSOLJ Variable Star Bulletin, 48, 1
  * Nagai (2010) Nagai, K. 2010, VSOLJ Variable Star Bulletin, 50, 1
  * Nagai (2013) Nagai, K. 2013, VSOLJ Variable Star Bulletin, 55, 1
  * Nagai (2014) Nagai, K. 2014, VSOLJ Variable Star Bulletin, 56, 1
  * Nelson (2009) Nelson, R. H. 2009, Information Bulletin on Variable Stars, 5875, 1
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python¶
  * Ogloza et al. (2000) Ogloza, W., Drozdz, M., & Zola, S. 2000, Information Bulletin on Variable Stars, 4877, 1
  * Paczyński (1967) Paczyński, B. 1967, Acta Astron., 17, 287
  * Parimucha et al. (2009) Parimucha, S., Dubovsky, P., Baludansky, D., et al. 2009, Information Bulletin on Variable Stars, 5898, 1
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in C. The art of scientific computing
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123
  * Prsa et al. (2011) Prsa, A., Matijevic, G., Latkovic, O., Vilardell, F., & Wils, P. 2011, PHOEBE: PHysics Of Eclipsing BinariEs, Astrophysics Source Code Library
  * Qian et al. (2008) Qian, S. B., Dai, Z. B., Zhu, L. Y., et al. 2008, ApJ, 689, L49
  * Qian et al. (2010) Qian, S.-B., Liao, W.-P., Zhu, L.-Y., & Dai, Z.-B. 2010, ApJ, 708, L66
  * Qian et al. (2011) Qian, S.-B., Liu, L., Liao, W.-P., et al. 2011, MNRAS, 414, L16
  * Qian et al. (2012) Qian, S. B., Zhu, L. Y., Dai, Z. B., et al. 2012, ApJ, 745, L23
  * Rappaport et al. (1983) Rappaport, S., Verbunt, F., & Joss, P. C. 1983, ApJ, 275, 713
  * Selam et al. (1999) Selam, S. O., Gurol, B., & Muyesseroglu, Z. 1999, Information Bulletin on Variable Stars, 4670, 1
  * von Essen et al. (2016) von Essen, C., Cellone, S., Mallonn, M., Tingley, B., & Marcussen, M. 2016, arXiv e-prints, arXiv:1607.03680
  * Wittenmyer et al. (2013) Wittenmyer, R. A., Horner, J., & Marshall, J. P. 2013, MNRAS, 431, 2150
  * Wolf et al. (2018) Wolf, M., Kučáková, H., Zasche, P., et al. 2018, A&A, 620, A72
  * Wood & Saffer (1999) Wood, J. H. & Saffer, R. 1999, MNRAS, 305, 820
  * Wood et al. (1993) Wood, J. H., Zhang, E.-H., & Robinson, E. L. 1993, MNRAS, 261, 103
  * Zejda et al. (2006) Zejda, M., Mikulasek, Z., & Wolf, M. 2006, Information Bulletin on Variable Stars, 5741, 1

##  Appendix A Mid-Eclipse Timings

space

Table 10: List of new primary (p) and secondary (s) eclipse timings as presented in this work (see text in Sect. 4 for details). Mid-Eclipse Times | Error | Type | Filter | Tel.
---|---|---|---|---
$(BJD_{TDB})$ | (d) | p/s |  |
2456708.494970 |  $\pm$0.000024 | p | R | TUG
2456708.553380 |  $\pm$0.000072 | s | R | TUG
2456708.611720 |  $\pm$0.00002 | p | R | TUG
2456813.367760 |  $\pm$0.000095 | s | R | TUG
2457075.461120 |  $\pm$0.000063 | p | R | AUKR
2457095.478600 |  $\pm$0.00013 | s | R | AUKR
2457095.536860 |  $\pm$0.000081 | p | R | AUKR
2457103.123850 |  $\pm$0.000092 | p | R | SOAO
2457103.182120 |  $\pm$0.00014 | s | R | SOAO
2457103.240600 |  $\pm$0.000054 | p | R | SOAO
2457103.298927 |  $\pm$0.00036 | s | R | SOAO
2457104.115700 |  $\pm$0.000042 | s | R | SOAO
2457104.174010 |  $\pm$0.000049 | p | R | SOAO
2457104.232810 |  $\pm$0.00012 | s | R | SOAO
2457104.290880 |  $\pm$0.000036 | p | R | SOAO
2457130.377900 |  $\pm$0.000054 | s | R | AUKR
2457130.435820 |  $\pm$0.000057 | p | R | AUKR
2457137.322470 |  $\pm$0.000022 | p | R | TUG
2457158.448580 |  $\pm$0.000064 | p | R | AUKR
2457169.303600 |  $\pm$0.000044 | p | R | AUKR
2457422.468150 |  $\pm$0.000095 | p | R | AUKR
2457422.584870 |  $\pm$0.000078 | p | R | AUKR
2457426.553410 |  $\pm$0.000064 | p | R | AUKR
2457426.611750 |  $\pm$0.000094 | s | R | AUKR
2457440.442900 |  $\pm$0.00016 | p | R | AUKR
2457440.501550 |  $\pm$0.00014 | s | R | AUKR
2457440.559700 |  $\pm$0.00008 | p | R | AUKR
2457441.493442 |  $\pm$0.000022 | p | R | AUKR
2457441.551859 |  $\pm$0.000069 | s | R | AUKR
2457441.610196 |  $\pm$0.000022 | p | R | AUKR
2457457.484029 |  $\pm$0.000034 | p | V | AUKR
2457457.542536 |  $\pm$0.000075 | s | V | AUKR
2457457.600768 |  $\pm$0.000045 | p | V | AUKR
2457469.097705 |  $\pm$0.00023 | s | R | SOAO
2457469.155938 |  $\pm$0.00019 | p | R | SOAO
2457469.214385 |  $\pm$0.00012 | s | R | SOAO
2457469.272720 |  $\pm$0.000048 | p | R | SOAO
2457469.331065 |  $\pm$0.00023 | s | R | SOAO
2457480.069250 |  $\pm$0.000088 | s | R | SOAO
2457480.185980 |  $\pm$0.00017 | s | R | SOAO
2457480.244330 |  $\pm$0.000094 | p | R | SOAO
2457480.302610 |  $\pm$0.00010 | s | R | SOAO
2457480.361057 |  $\pm$0.000077 | p | V | AUKR
2457480.419463 |  $\pm$0.00010 | s | I | AUKR
2457492.091380 |  $\pm$0.00002 | s | R | SOAO
2457492.149730 |  $\pm$0.000019 | p | R | SOAO
2457506.331194 |  $\pm$0.000083 | s | I | AUKR
2457506.389490 |  $\pm$0.000011 | p | I | AUKR
2457506.447964 |  $\pm$0.000036 | s | V | AUKR
2457506.506248 |  $\pm$0.000013 | p | V | AUKR
2457506.856370 |  $\pm$0.000005 | p | B | LOAO
2457507.673350 |  $\pm$0.000011 | p | V | LOAO
2457507.731700 |  $\pm$0.000038 | s | V | LOAO
2457507.790090 |  $\pm$0.00001 | p | V | LOAO
2457507.848350 |  $\pm$0.000031 | s | V | LOAO
2457508.782240 |  $\pm$0.000021 | s | R | LOAO
2457508.840580 |  $\pm$0.000009 | p | R | LOAO
2457787.566715 |  $\pm$0.000035 | p | B | AUKR
2457787.625322 |  $\pm$0.00022 | s | B | AUKR
2457803.381877 |  $\pm$0.00018 | s | B | AUKR
2457803.440585 |  $\pm$0.000037 | p | B | AUKR
2457803.498890 |  $\pm$0.00017 | s | B | AUKR
2457803.557359 |  $\pm$0.00004 | p | B | AUKR
2457803.615831 |  $\pm$0.000095 | s | B | AUKR
2457847.502204 |  $\pm$0.000069 | s | R | AUKR
2457847.560560 |  $\pm$0.00005 | p | R | AUKR
2457911.347827 |  $\pm$0.00015 | s | R | AUKR
2457917.358775 |  $\pm$0.00005 | p | R | AUKR
2458242.305835 |  $\pm$0.00001 | p | R | TUG
2458242.364198 |  $\pm$0.000032 | p | R | TUG
2458567.369609 |  $\pm$0.000026 | p | R | AUKR

##  Appendix B Parameter uncertainties for Model 6

Figure 10: Marginalized probability density distributions ($P(x,y)\propto
exp(-\chi^{2}/2)$) for the parameters of Model 6 as obtained from the
$\chi^{2}$-grid method. The 68.3% confidence interval and median values are
shown with vertical lines (see text for more details). Figure 11: Histograms
of parameters from the bootstrap method. $\pm 68.3\%$ confidence limits and
median values are shown with vertical lines.

##  Appendix C $\chi^{2}$ surface grids for Model 6

Figure 12: $\chi^{2}$ surface plots for the pairs of parameters of Model 6.
The vertical and horizontal dashed lines represent $\pm 1\sigma$ values.
Figure 13: (continued) $\chi^{2}$ surface plots for the pairs of parameters of
Model 6.

Figure 14: (continued) $\chi^{2}$ surface plots for the pairs of parameters of
Model 6. Figure 15: (continued) $\chi^{2}$ surface plots for the pairs of
parameters of Model 6.