Eclipse Timing Modeling of Three Post-Common Envelope Binaries: Hybrid Solutions

HW Vir (V = 10.7, P = 2.8 hr) is the prototype of the eponymous HW Vir class of eclipsing subdwarf B stars with low-mass dwarf companions. Kilkenny et al. (1994) first discovered variations in the orbital period and reported a period decrease over nine years. Several authors have carried out subsequent ETV studies of HW Vir (Çakirli & Devlen, 1999; Kiss et al., 2000; Kilkenny et al., 2000; Kilkenny et al., 2003; İbanoǧlu et al., 2004; Lee et al., 2009; Qian et al., 2010; Beuermann et al., 2012b; Skelton, 2012; Lohr et al., 2014; Baran et al., 2018; Esmer et al., 2021; Brown-Sevilla et al., 2021). These papers invoked one or more substellar companions as the primary cause of the observed ETVs. However, as new data became available, these models inevitably needed to be modified.

After ruling out mass transfer, magnetic braking, and apsidal motion, Çakirli & Devlen (1999) concluded that a third body with mass 23 M_Jup and an orbital period of 19 yrs is the most probable cause. İbanoǧlu et al. (2004) also proposed a brown dwarf companion, with a minimum mass of 23 M_Jup  and orbital period of 18.8 yrs. Lee et al. (2009) proposed two companions with periods 15.8 and 9.1 yrs and masses of 19.2 M_Jup and 8.5 M_Jup respectively. Beuermann et al. (2012a) added more ETV data, and showed that they deviated significantly from the solution proposed by Lee et al. (2009). They also showed that the solution proposed by Lee et al. is dynamically unstable. They proposed a new solution consisting of two brown dwarfs (14 M_Jup + 65 M_Jup ) that was stable over more than 10 Myr. However, Baran et al. (2018) showed that the solution of Beuermann et al. does not fit more recent O-C data.

Horner et al. (2012) independently made a stability analysis of the Lee et al. solution, and also found it highly unstable. Indeed, they searched a significant volume of substellar component parameter space, and could not find any stable solutions, concluding "If the [HW Vir ] system does host exoplanets, they must move on orbits differing greatly from those previously proposed". This failure to find a stable solution is reinforced by the recent studies of Esmer et al. (2021) and Brown-Sevilla et al. (2021) using all previously published O-C observations along with new times of minima up to 2019.4. They tried a number of multi-planet models but also failed to find a stable solution.

Drechsel et al. (2001) first confirmed HS0705+6700 (V470 Cam, V= 14.7, P= 2.3 hr) as a post common envelope sdB binary. They fit their eclipsing timing observations with a linear ephemeris. Qian et al. (2009) reported a cyclic change in orbital period and proposed a substellar companion with orbital period of 7.15 years and mass of 40 M_Jup is responsible for the observed ETV. Beuermann et al. (2012b) observed an additional 28 epochs, with a similar solution: A single brown dwarf with mass 31 M_Jup and orbital period of 8.4 years.

In 2013, Qian et al. (2013) fit a single brown dwarf model with a period 8.87 years and mass 32 M_Jup , along with a positive quadratic term $dP/dt=+9.8\times 10^{-9}\textrm{\ days/year}$. They ruled out mass transfer and angular momentum loss due to gravitational radiation and magnetic braking as explanations of this period increase. Pulley et al. (2015) showed that the more recent O-C timing data began to diverge from the Qian model. Later in 2018, they confirmed this departure by adding 65 more times of minima (Pulley et al., 2018).

Bogensberger et al. (2017) proposed a new one planet model with a longer orbital 11.8 year period, eccentricity 0.38, and no quadratic term. Recently, Sale et al. (2020) added 20 additional ETV observations through mid-2018. They found a two component solution that fit all published O-C observations, but unfortunately it was highly unstable over a short timescale (1,000 yr) and hence not physically plausible.

HS2231+2441 (V = 14.2, P = 2.7 hr) is an eclipsing white dwarf + brown dwarf binary (Almeida et al., 2017) first reported by Østensen et al. (2007). They fit two years of primary eclipse timing observations with a linear (constant period) ephemeris. Qian et al. (2010), using more extensive timing data, reported that the orbital period shows both a cyclic change and a period derivative. They proposed a 14 M_Jup component to explain the cyclic variation and magnetic braking to explain the period decrease. However, as more timing observations became available, Lohr et al. (2014) found that a linear function might provide a better fit to the O-C curve. Almeida et al. (2017) reported eight more eclipse times and updated the linear ephemeris. Most recently, Pulley et al. (2018) extended the timing database with 26 more timing minima. They compared the quadratic ephemeris by Qian et al. (2010) and linear ephemeris proposed by Lohr et al. (2014) concluded there is no statistically significant difference between the two.

Our eclipse timing dataset compromised 90 new time of primary eclipse minima (Table 2) together with previously published times from the literature. Before use, we made a few modifications.

1. 1.

   For HW Vir and HS2231+2441, we excluded a portion of archival SuperWASP times of minima that were rejected by Lohr et al. (2014)’s outlier methodology and by the size of their uncertainties. In particular, we did not include Lohr et al. (2014)’s listed extra times of minima due to their precision, and discarded timing data for which uncertainties are between 8 sec $\sim$40 sec. We also omitted the individual superWASP times of minima for HS2231+2441, since they had large ($\sigma>20s$) uncertainties.

2. 2.

   We also removed O-C eclipse timings obtained from the Mercator observation in Baran et al. (2018) for HW Vir since they deviate from the overall O-C trend line by 50 seconds. We also excluded outliers that had O-C values deviating by more than 400 seconds from the overall O-C trend.

For a binary system with orbiting substellar components, the predicted primary mid-eclipse time as a function of cycle number $E$ can be written as

where $T_{0}$ is the reference epoch, $P_{0}$ is the linear eclipse period, where $\tau_{i}$ are the time-dependent contributions from each perturbing body. and

is a quadratic term can would result from non-planetary effects, such as apsidal motion, mass exchange, magnetic braking, or gravitational radiation. We used the usual formulation of Irwin (1952, 1959) to calculate the time signature of perturbing bodies on eclipse timing as a function of each body’s orbital properties.

A physically credible model must not only fit all the current and previously published ETV datasets, but it is also subject to the prior constraint that the substellar components must be stable to orbital instabilities over timescales at least comparable to the lifetime of the PCEB phase ($\sim$ 100 Myr) (Han et al., 2002; Han et al., 2003; Zorotovic et al., 2011). We used the Hill stability criterion (Gladman, 1993; Chambers et al., 1996) to predict whether any two components in orbit about a massive central body will avoid mutual close encounters over many orbits.

For circular orbits, stability occurs if the semi-major axis difference between pairs of components exceeds a critical distance $\Delta$ given by (Chambers et al., 1996),

where $R_{H}$ is the Hill radius,

where $m_{1},m_{2}$ and $a_{1},a_{2}$ are the masses and semi-major axes of the components respectively, and $M$ is the mass of the central body.

For non-circular orbits, a more general stability criterion is (Hasegawa & Nakazawa, 1990)

where

and $e$ is the larger eccentricity. For low eccentricity orbits and $K\ll 1$, we can expand the second term to obtain

The second term in parentheses is much less than 1 for the components considered here ($K\ll 1$) so the simpler expression equation 3 was used in our calculations.

We can define a dimensionless parameter $\beta$ that characterizes each planet pair spacing in terms of the normalized mutual Hill radii as (cf. Smullen et al., 2016)

Stable component pairs correspond to $\beta>1$ while unstable pairs have $\beta<1$.

For systems with more than two substellar components, instability can occur even if all components lie outside the Hill radius of their nearest companion, although the timescale for close encounters increases as the initial component separations become larger (Chambers et al., 1996). Hence, we used the Hill parameter as an optimization prior but ran numerical orbit simulations to test orbit stability after optimization, as described below.

The model fitting scheme comprised four steps (Fig. 3):

1. 1.

   Create an initial model of between one and four sub-stellar components, with multi-component models spaced at least one Hill radius apart (i.e., $\beta>1$ for all pairs),

2. 2.

   Search for the minimum of an objective function proportional to the product of the summed $\chi^{2}$ of the model vs. O-C time series and an ansatz ’penalty’ function that sharply increases when the smallest Hill parameter $\beta$ becomes less than 1.0 for any component pair,

3. 3.

   Use Bayesian inference to solve for the final fitted parameters using a Markov chain Monte Carlo (MCMC) sampling, a weighted least-squares objective function, and an orbit stability prior ($\beta>1$),

4. 4.

   Run orbit stability tests using both the N-body orbital integration package REBOUND (Rein & Liu, 2012; Rein & Spiegel, 2015) and the chaotic indicator program MEGNO (Cincotta & Simó, 2000; Hinse et al., 2010).

We first chose the number of substellar companions to model. Examination of the O-C plot for both HS0705+6700 and HW Virginis indicates that a large number of components, although certainly possible, is not justified by the data, since the model parameters are already highly degenerate for few component models, as we will see. Hence we restrict our model selection to one through four component models. Each component is characterized by five parameters: mass, semi-major axis, eccentricity, longitude of periapsis $\tilde{\omega}$, and time of periastron passage . We restricted the model to co-planar orbits for the sake of computational simplicity. We also added a period derivative parameter to account for the possible contribution of non-orbiting components, such as apsidal motion or gravitational radiation.

We started model fitting using a downhill simplex (Nelder-Mead) minimization (Python library LMFIT³³3https://lmfit.github.io/lmfit-py to fit the O-C curve. However, instead of minimizing the sum-squared differences between the O-C data and the model, we incorporate the orbit stability constraint by minimizing the product,

Where we define an empirical Hill stability function

where the Hill parameter $\beta_{min}$ is the smallest normalized Hill parameter (equation 8) of all component pairs. This stability function asymptotes to unity for $\beta>1$ but increases sharply for $\beta<1$ (see Fig.2).

The resulting model was used as a starting point for Bayesian inference modeling. We used the EMCEE (Foreman-Mackey et al., 2013) Python implementation⁴⁴4https://github.com/dfm/emcee with 4096 walkers and 4000 iterations. The likelihood function was summed-squared difference between the model and O-C data weighted by the data uncertainties. We used three priors: (i) stability parameter $\beta>1.05$, (ii) all masses and semi-major axes positive-definite, and (iii) eccentricities restricted to the range $0<e<1$.

Using a model fitting scheme with a Bayesian methodology as the final step has two important advantages over the more commonly-used frequentist least-squares minimization approach. First, the orbit stability constraint (Hill parameter) can be incorporated naturally as a Bayesian prior without the ad hoc ansatz function we used in the LMFIT minimization. Second, the Bayesian posterior distribution provides a more comprehensive description of the parameter uncertainties and correlations.

For each model, we tested dynamical stability by numerically integrating the component orbits using WHFAST, a Wisdom-Holman symplectic integrator (Rein & Tamayo, 2015). We used a time step of 36.5 days and integrated for 10 Myr for each model. In addition, we tested for chaotic behavior using the MEGNO algorithm (Mean Exponential Growth factor of Nearby Orbits, Cincotta & Simó (2000)) also with a time interval of 10 Myr. Both WHFAST and MEGNO were called from the Python package REBOUND (Rein & Liu, 2012) which also provided convenient tools for plotting orbits and time histories of orbital elements.

After adding 26 additional primary eclipse timing observations (listed in Table 2) to previously published observations, we fit the linear ephemeris of Almeida et al. (2017) with a slightly changed orbital period,

The O-C diagram with all of the available data is shown in Figure 4.

We find that there is no statistically significant period variation for HS22315+2441 spanning 16 years, from 2004 to 2020, but we cannot rule out the variations smaller than 15 seconds. A linear fit without a orbiting planet is consistent with the full data set.

Figure 5 Shows the O-C timing residuals versus epoch computed using linear ephemeris,

along with our one component model (solid blue line) and the one-component model of Bogensberger et al. 2017 (solid green line). Table 4 lists the model parameters of our model.

We also found a best-fit two component model with masses 21.8 M_Jup  and 1 M_Jup  and semi-major axes 4.9 AU and 37.3 AU, respectively. The two component model was stable for at least $10^{7}$ years as verified by REBOUND and MEGNO simulations, but the fit to the timing data was no better than the one component model. Therefore, we do not discuss this model further.

It is clear that the previously published model of Bogensberger et al. deviates strongly from the observed timing data starting shortly after the publication of the model. As previously noted, this discouraging trend has been seen for several other eclipsing sdB binary exoplanet detection. Also, we did not plot the recent two component solution proposed by Sale et al. (2020) since we could not determine their ephemeris or interpret orbital parameters of their model from their paper.

An expanded view of the O-C fits for our model is shown in Fig. 6, along with the $\pm 1\sigma$ posterior spread from the EMCEE solutions. Fig. 7 shows the two-dimensional projections of the sample covariances for our one component model. They demonstrate that our model have well-defined parameters with small fractional uncertainties.

Figure 8 Shows the O-C timing residuals versus epoch using the linear ephemeris

On the right of each panel are plots of eccentricity and semi-major axes of all model components vs. epoch as computed by numerical integration using REBOUND with time steps of 0.1 yr. The plots also show two model fits:

• •

  The upper panel shows a three component model computed using the minimization scheme described in section 3, with a Bayesian prior $\beta>1.05$ for all component pairs, where $\beta$ is the normalized Hill constraint. This solution is stable for at least 10 Myr, but is a very poor fit to the timing data ($\chi^{2}\approx 93$).

• •

  The lower panel shows a 4-component model, also computed using the minimization scheme in section 3, but without the Hill stability prior. This provides a much better fit ($\chi^{2}\approx 19$) but is highly unstable over short time internals (one hundred years).

We do not list the orbital elements for either model since they are implausible: They are either very poor fits to the timing data (3-component stable model) or are highly unstable (four-component model). This unfortunate scenario (acceptable O-C model fits are invariably unstable) was also found by (Brown-Sevilla et al., 2021) who analyzed a series of best-fit six multi-component models for HW Vir. In the next section we consider possible alternatives to pure substellar component models.

As we have seen, the best-fit stable orbit model fits for both HS0705 and HW Vir are not fully in agreement with the observed O-C data, with residual average differences of order 10 sec and 20 sec over 10-20 year timescales for HS0705 and HW  Vir respectively (Figs. 5,8). We now consider the question of whether the ETV variations may result from a combination of substellar component perturbation and Applegate-Lanza effects.

The change in fractional angular velocity $\Omega$ of the secondary is related to the observed fractional period changes by (Lanza, 2020, eqn. 56)

Where $m$ in the reduced mass of the binary, $a$ is the semi-major axis of the binary, $I_{s}$ is the moment of inertia about the active (secondary) star spin axis, and $\Delta P/P=\Delta t/T$ is the fractional change in the period, or equivalently the O-C deviations $t$ over the timescale $T$ of the O-C data . The consequent change in rotational energy is

Combining these equations, the dependence on $I_{p}$ cancels, giving,

For both HS0705+6700 and HW Vir, the average timing residuals of the substellar companion model fits to the observed O-C data are both of order 10 sec over 10 yr ($\Delta t/3T\sim 10^{-8}$), while the bracketed term evaluates to $3\cdot 10^{40}$ J, so the change in rotational energy is

We can compare this with the maximum energy available from luminosity changes over the timescale for O-C variability in the secondary star. From Table 3,

Remarkably, this is the same as the energy loss required to account for residual changes in the O-C model fits, but would be only $\sim 10\%$ of the entire O-C deviations from a constant period were due to an Applegate-Lanza mechanism. This suggests that the observed O-C variations could be explained by a combination of substellar companions accounting for the bulk of the ETV’s, with a Lanza-Applegate mechanism contributing about 10% of the total.

Another relevant dynamical timescale is the tidal synchronization time i.e., the timescale in which tides on the secondary dissipate the kinetic energy of asynchronous rotation (Lanza, 2020),

where Q is a dimensionless tidal quality factor that characterizes the efficiency of tidal energy dissipation in the active (secondary) component, $M_{T}$ is the total mass of the binary system, $m_{p}$ is the primary mass, $a$ is the binary semi-major axis, and $R$ is the active star radius. Following Lanza (2020), we assume $Q\sim 10^{5}-10^{6}$ corresponding to strong tidal coupling due to nearly synchronous rotation. For the moment of inertia $I_{s}$, we use the numerical expression of Criss & Hofmeister (2015) for a fully convective secondary star with a polytropic index n = 1.5,

Using stellar parameters list in Table 3, we find tidal synchronization timescales $t_{s}$ = 4 – 40 yr and 10 – 100 yr for HS0705 and HW Vir respectively. These times are comparable to the observed O-C variations, so tidal synchronization may dictate the timescale for spin-orbit coupling rather than magnetic activity timescales.