NY Vir Eclipse timing observations: Two planet solution

TL;DR

This study updates the circumbinary exoplanet model for NY Virginis using new eclipse timing data, identifying a stable two-planet solution with specific orbital parameters, though some degeneracies remain.

Contribution

It presents the first comprehensive two-planet model for NY Virginis based on combined eclipse timing data and stability analysis.

Findings

A stable two-planet model with specific eccentricities and masses was identified.

The model remains stable for at least 10^8 years under certain eccentricity conditions.

Additional observations are needed to refine planetary parameters due to parameter degeneracies.

Abstract

We report 18 new primary minima timing observations of the short-period eclipsing binary system NY Virginis. We combined these minima with previously published primary minima to update circumbinary exoplanet models in this system based on O-C timing variations. We performed a non-linear least-squares minimization search using a quadratic ephemeris and either one or two exoplanets. The only model with an acceptable fit includes a period derivative $\dot{P} = 2.83\times10^{-12}$ and two planets in eccentric orbits $e = 0.15,0.15$ with minimum masses 2.7 and 5.5 Jovian masses. Analysis of the orbit stability shows that this solution is stable for at least $10^8$ yr, but a small increase in eccentricity ($e\geq0.20$) for either planet renders the orbits unstable in less than $10^6$ years. A number of model parameters are significantly degenerate, so additional observations are required to determine planetary parameters with high statistical confidence.

Figures

Tables

Table 1. Table 1: NY Vir Primary Mid-Eclipse Times

BJD
2458142.913026(9)
2458143.014080(4)
2458143.923178(2)
2458145.943570(3)
2458146.953672(3)
2458147.963801(3)
2458159.984740(4)
2458180.995988(4)
2458190.895673(2)
2458193.825018(2)
2458201.906339(6)
2458202.916509(2)
2458204.936869(4)
2458205.947009(6)
2458206.957149(2)
2458226.756344(1)
2458233.928400(5)
2458236.857910(3)

Table 2. Table 2: Best-fit parameters for the quadratic plus two planet model for NY Vir. 1- σ 𝜎 \sigma uncertainties are given in parentheses.

Parameter	Fitted Values		Unit
Inner Binary
T0	2453174.442647(13)		BJD
P0	0.1010159677(4)		day
dP/dt	$2.83\left(0.25\right)\cdot {10}^{-12}$		s/s
Planets
	LTT 1	LTT 2	units
P	8.64 (0.17)	24.09 (0.65)	yr
$e$	0.15 (0.08)	0.15 (0.01)
K	7.6 (0.7)	31.4 (1.1)	sec
asin(i)	3.55 (0.01)	7.04 (0.25)	AU
T	2453472	2450031	JD
$\omega$	348 (6)	320(4)	deg
Min. Mass	2.66 (0.26)	5.54 (0.2)	MJ
${\chi }_r$	1.12 (p = 0.17)

Equations

$$J{D}_{min}(BJD)=2453174.442647+0.1010159677E$$

$$T(E)=T_0+P_{0}E+\frac{1}{2}\frac{dP}{dt}{E}^2+{\tau }_1(E)+{\tau }_2(E)$$

$$\begin{multlined}\tau_{i}=\frac{K_{i}}{\sqrt{1-\small[e_{i}cos(\omega_{i})\small]^{2}}}\cdot\\ \bigg{[}\frac{1-e_{i}^{2}}{1+e_{i}cos(\nu_{i})}sin(\nu_{i}+\omega_{i})+e_{i}sin(\omega_{i})\bigg{]}\end{multlined}\tau_{i}=\frac{K_{i}}{\sqrt{1-\small[e_{i}cos(\omega_{i})\small]^{2}}}\cdot\\ \bigg{[}\frac{1-e_{i}^{2}}{1+e_{i}cos(\nu_{i})}sin(\nu_{i}+\omega_{i})+e_{i}sin(\omega_{i})\bigg{]}$$

Full text

binaries: close - binaries: eclipsing - stars: individual: NY Virginis - planetary system

1 Introduction

NY Vir (= PG1336-018, GSC 04966-00491; V=13) is an eclipsing binary comprising a hot sub-dwarf B star and a M5 dwarf with an orbital period of 2.42 hours (Kilkenny et al., 1998; Vucković et al., 2007). It is a member of the Post-Common Envelope Binary (PCEB) class, in which the primary star, initially in a wide ($\sim$1 AU) binary, has reached the giant branch and the secondary has been engulfed in the giant star’s envelope. This results in both energy and angular momentum transfer onto the outer shell of the giant. The secondary then spirals inward, so that the binary separation is approximately one solar radius (Paczynski, 1976; Heber, 2009). Long-term eclipse timing observations of these systems are very useful to probe their evolution and dynamics, such as angular momentum transfer and detection of additional components, including exoplanets. This approach for detecting exoplanets has been applied to several short-period eclipsing binaries, including HW Vir (Lee et al., 2009), NN Ser (Beuermann et al., 2010; Caceres et al., 2014), and NSVS 14256825 (Almeida et al., 2013).

Kilkenny et al. (1998) initially fit primary eclipse timing observations of NY Vir with a linear (constant period) ephemeris. However, later observations deviated from this linear prediction significantly. Both Kilkenny (2011) and Çamurdan et al. (2012) found that the O-C residuals could be fit with a negative quadratic function i.e., a slow decrease in orbital period ($\dot{P}\approx-4\cdot 10^{-8}$ day/yr) which Çamurdan et al. (2012) ascribed to angular momentum loss from the binary system. This loss cannot be caused by mass transfer between the components, since the eclipsing pair is a detached system and both stars are nearly spherical (Lee et al., 2014). Rather, the momentum loss is more likely caused by magnetic braking in the cool secondary (Rappaport et al., 1983).

As more timing observations became available, the observed O-C timing residuals began to diverge from a quadratic ephemeris. Qian et al. (2012) added a sinusoidal term to the quadratic model, implying the existence of an orbiting circumbinary third body that would cause periodic motion of the binary center of mass with a light travel time modulation of 6.3s. Lee et al. (2014) observed NY Vir over an additional 39 epochs. They fit two models: A single planet solution with a quadratic term, and a two-planet solution with a constant period. The best fit was a two-planet model, with masses 2.8 and 4.5 Jovian masses and orbital periods 8 and 27 years respectively.

Unfortunately, these models have recently been invalidated by more recent observations. Pulley et al. (2018) reported new eclipse timings for seven sdB binaries, including 15 new times of primary minima for NY Vir. This extended the time baseline by three years to 2017.7. They found that that neither the one-planet model of Qian et al. (2012) nor the two-planet model of Lee et al. (2014) can fit the full O-C residual time history. They did not attempt a new planetary model fit, arguing instead that more timing observations are required to determine any planetary parameters with high confidence.

Other mechanisms that could account for eclipse timing variations have also been discussed. Both Qian et al. (2012) and Lee et al. (2014) considered the possibility that gravitational radiation from the binary may be responsible for the decreasing period, but concluded that the magnitude is two dex lower than the observed effect. Another possibility for the quadratic term is gravitational coupling of the orbit to changes of oblateness of a magnetically active star (Applegate, 1992). Völschow et al. (2016) considered the Applegate mechanism for a group of sixteen close binaries, including NY Vir. They found that although it could marginally account for observed period variations in several of the eleven listed PCEB binaries, it was insufficient to explain the period variations in the majority of these systems, including NY Vir, where the energy discrepancy was more than a factor of 100 too small.

Finally, not all planetary models that can account for the observed O-C residuals have stable orbits. (Tinney et al., 2011) determined that the two-planet solution for the eclipsing polar HU Aqr proposed by Zhao et al. (2011) is dynamically unstable on a timescale of 5000 yr. Similarly, (Tinney et al., 2012) found that the planetary solutions for NN Ser proposed by Beuermann et al. (2010) are also unstable, although this claim has been disputed (Caceres et al., 2014).

In this paper, we report new eclipse timing observations of NY Vir that extend the time base to epoch 2018.5. We compare the O-C residuals to several models, including one, two, and multiple planet systems, and evaluate the uniqueness of the derived model parameters. We also present the results of a dynamical model to determine the stability of the derived planetary orbits.

2 Observations

We observed NY Vir at primary eclipse on 18 nights between 24 Jan 2018 and 18 April 2018. The observations were made with the 0.5m telescope at the Iowa Robotic observatory111http://astro.physics.uiowa.edu/iro, located at Winer Observatory222winer.org in southeastern Arizona. We used a SBIG 6303e CCD camera (2048x3072 x 12 micron pixels) and a Sloan r*′* filter. At each epoch, we took a sequence (typically 20) of 15 sec exposures centered on the predicted time of primary eclipse.

After the normal CCD image calibration (dark subtraction, flat fielding), we created light curves using airmass-corrected differential photometry. We found each epoch’s time of minimum and associated uncertainty by fitting a Gaussian profile to the primary eclipse light curve using a nonlinear least-squares (Levenberg-Marquardt ) algorithm. We then determine the observed mid-eclipse times. Table 1 lists the eclipse times and uncertainties.

3 Additional Data and Analysis

We combined the present observations with previously published times of primary minima listed in Kilkenny et al. (1998),Kilkenny et al. (2000), Kilkenny (2011), Çamurdan et al. (2012), Qian et al. (2012), Lee et al. (2014), Pulley et al. (2016), Pulley et al. (2018) and Baştürk & Esmer (2018). In addition, we added several minima found in the AAVSO (aavso.org) and WASP (Butters et al., 2010) public archives. We made a few modifications to these data before use: (a) Kilkenny et al. (2000) did not publish exact uncertainties of the minima. They reported that the uncertainties are smaller than 0.00005 days. We used 0.00005 days for all their uncertainties. (b) The individual WASP times of minima have large ($\sigma>$10s) uncertainties. These data are concentrated in three seasons, so we took the median O-C value for each season, with an uncertainty given by the standard deviation of all times for that season. (c) Lee et al. (2014) have several outliers i.e., closely-spaced times of minima that differ from by more than three standard deviations from the overall trend. We removed these outliers.

The aggregate dataset consisted of 104 times of primary minima spanning 23 years, from 1996.4 to 2018.5. We calculated a new best-fit linear ephemeris using all of the available time of minima,

[TABLE]

to calculate the O-C residuals and cycle number for all observed primary minima times. The O-C diagram with all of the available data is shown in Figure 1.

It is clear that the O-C plot is complex, and excludes both a simple quadratic ephemerides model or strictly periodic variations, as might be expected from e.g., pure apsidal motion or a modulation by a third body in a circular orbit. Hence, we used a family of models that included both a quadratic term and either one or two orbiting bodies in elliptical orbits. The predicted times of minima are

[TABLE]

where the light travel time delay caused by the reflex motion of the binary by each perturbing body is given by (Irwin, 1952, 1959),

[TABLE]

where the orbit number dependence $E$ is computed by solving Kepler’s equation for the eccentric anomaly, followed by the true anomaly $\nu$.

We used a non-linear least-squares algorithm (Python library LMFIT) to fit the O-C data starting with a simple quadratic model, then adding one or two planets. For the simple quadratic and quadratic plus one planet models, we used a brute-force grid search, followed by a downhill-simplex (Nelder-Mead) chi-square minimization parameter search starting at the brute-force minimum grid point. For the two- and three-planet searches, the number of free parameters was large enough that a brute force grid search was computationally prohibitive, so we used a Monte Carlo randomized selection of initial parameter start values and a downhill-simplex minimization search for each set of initial parameters.

4 Results

Figure 1 shows the O-C timing residuals versus epoch computed using the linear ephemeris in equation (1). Also shown are the quadratic plus one-planet model of Qian et al. (2012) (red dotted line), the two-planet model of Lee et al. (2014) (green dotted line), and our best-fit model, a quadratic plus two planet solution (black solid line), summarized in Table 2. The reduced chi-square of the fit is $\chi^{2}=1.12$, p =0.17. Both the quadratic only and quadratic plus one planet models, for which a brute-force minimization search was used, resulted in model fits with reduced chi-square greater than 1.9, which can be rejected with high confidence ($p<10^{-6}$).

Our two-planet solution consists of a quadratic (period derivative) term $dP/dt=(2.83\pm 0.28)\cdot 10^{-12}$ s/s and two planets with minimum masses 2.66$\pm$0.36 MJ and 5.54$\pm$0.28 MJ, and orbital periods 8.64$\pm$0.17 yr and 24.09$\pm$0.65 yr respectively. Note that to convert light travel times to planetary masses, we have assumed a combined stellar binary mass of 0.59 Msun (Vucković et al., 2009).

5 Discussion

5.1 Comparison with previous planetary solutions

The planetary model described in this paper most closely resembles the two-planet solution of Lee et al. (2014). However, Lee et al. (2014)’s model had no period derivative term, and does not fit the observed O-C time differences since epoch 2015.0 (cf. Fig. 1). Likewise, the one-planet plus period-derivative model of Qian et al. (2012) has a single-planet mass (2.3 MJ) and period (7.9 yr) similar to our inner planet solution, but has a significantly larger negative period derivative. As with Lee et al. (2014), the Qian et al. (2012) model does not match recent O-C data.

Baştürk & Esmer (2018) recently published a quadratic plus one-planet model that claims to fit all published O-C data through epoch 2017.4 with a remarkably low reduced chi-square value ($\chi^{2}$= 0.69). However, we could not confirm their results. We were unable to reproduce their published O-C data (Fig. 1) using their linear ephemeris. This is in part because the published linear period (P = 0.101016 $\pm$ 0.000001 day) has insufficient precision: A difference of 0.000001 days over 60,000 periods changes the O-C plot by 0.06 days, over 100x larger than the O-C departures shown in their Fig. 1. We also tried comparing their one-planet model with our O-C data, generated with our linear ephemeris (Eqn. 1). This resulted in a very poor fit ($\chi_{r}^{2}$ = 3.1, $p\ll 10^{-9}$). We also performed a brute-force chi-squared minimization grid search using their model as a starting point, but were unable to find any one-planet solution with a reduced chi-square less than 1.9 ($p<10^{-6}$).

Finally, we note that the planetary masses are similar to those derived for NY Vir by (Schleicher & Dreizler, 2014) using a second generation common-envelope ejecta model for planet formation (6.2 MJ , 2.5 MJ ), with the larger mass being the outer planet in both cases.

5.2 How well-determined are the fitted model parameters?

Although the two-planet model described in this paper fits all available O-C timing observations with an acceptable goodness of fit, it would be naive to accept the fitted model parameters with their associated formal uncertainties as a unique solution, owing to several well-known difficulties associated with non-linear fitting algorithms (e.g., Transtrum et al., 2010). These include the often poorly-constrained choice of initial parameter space coordinates, parameters non-orthogonality (correlation), and overfitting data by adding additional free parameters (e.g., additional planets) that are unjustified by the dataset. We now discuss each of these briefly in the context of the present model.

As described in section 3, we used a Monte Carlo randomized selection of initial parameter values as a starting point for the simplex chi-squared minimization search in the multi-dimensional parameter space. For each of the 10 independent parameters (4 per planet, plus period and period derivative), we chose a plausible range, then divided it into 10 uniformly sampled values. This resulted in $10^{10}$ possible starting points, which is clearly computationally impossible to evaluate. Hence we randomly chose a single coordinate as a starting point, then used the Neader-Mead downhill-simplex algorithm to find the nearest local minimum based on reduced chi-square. This procedure was repeated for 50 trials, and the smallest chi-square solution was used as the fitted parameter set. Although more than 30% of the trials found the same minimum, this does not guarantee that the minimum is indeed a global minimum. However, since the reduced chi-square of the chosen fit (1.12, p= 0.17) was already in good agreement, alternate solutions with even lower chi-square values would be only marginally statistically better.

To evaluate possible correlations between parameters, we used the LMFIT Python library implementation of the Nelder-Mead algorithm (method leastsq), which reports the solution cross-correlation matrix. There were about 20 high ($\rho>0.5$) correlation coefficients, especially between planetary parameters describing the outer planet. For example, the correlation coefficients between the the outer planet eccentricity and the period (P2), inner planet light-travel time (K1), and outer planet argument of periastron ($\omega_{2}$) all exceeded 0.8. This degeneracy can in principle be broken by additional data, as shown in Fig. 2, where we use four different fixed values of the outer planet’s orbital eccentricity, fitting for the other parameters and using the resulting models to predict O-C values 20 years into the future. It is clear that for this example, the degeneracy will not be broken for at least 10 years.

An alternate approach to model fitting, the random-walk Markov Chain Monte Carlo (MCMC) algorithm, is less sensitive to initial parameter guesses, provides a more nuanced description of parameter uncertainty, and has a much higher likelihood of finding a global minimum in chi-square space. MCMC has recently been used to evaluate possible exoplanet parameters derived from light-time variations in several close binary systems e.g., NSVS 14256825 (Almeida et al., 2013; Nasiroglu et al., 2017), NN Ser (Caceres et al., 2014), and HU Aqr (Goździewski et al., 2015). However, given that the frequentist (least-squares) technique found an acceptable fit at a [possibly local] chi-square minimum and the high degeneracy of the fitted parameters, finding a slightly better fit would not significantly better constrain the exoplanet parameters.

Finally, we consider overfitting, i.e., whether adding additional parameters (e.g., a third planet) is justified. One such test was developed by Schwarz (Schwarz, 1978). In this test, one calculates a dimensionless parameter, the Bayesian information criterion (BIC), for two models with a differing number of parameters. If the difference in BIC is less than 2 (Kass & Raftery, 1995), the more complex model is not justified. In our case, adding a third planet reduced the BIC from 52.8 (two-planet solution) to 50.4 (three-planet solution), which barely exceeds the unjustified threshold. We conclude that although a three planet model fits the O-C data slightly better, additional observations will be needed before one could confidently distinguish between these alternatives.

5.3 Are the planetary orbits stable?

An obvious question for any multiple planet solution is whether the orbits are dynamically stable over the timescale associated with the system. This in turn depends on whether the planets were ‘first generation’ i.e., formed coeval with the parent stars or ‘second generation’ i.e., from the ejecta of the common envelope of the close binary (Bear & Soker, 2014; Nixon et al., 2018). Schleicher & Dreizler (2014) considered models for secondary planet formation from disc ejecta for twelve post-common-envelope binaries with evidence of planets from O-C variations, including NY Vir . They found that some planetary systems were more likely first generation, whereas others, including NY Vir , may have formed from common envelope disk material. We note that for NY Vir, their two-planet mass estimates 6.2 MJ and 2.4 MJ are in good agreement with our results.

Bear & Soker (2014) also analyzed the same twelve systems, focusing on angular momentum evolution. They evaluated the ratio of planetary system angular momentum to the initial binary system angular momentum. For NY Vir , the ratio was 0.05, which is well within the 0.20 limit they adopted for plausible second generation planet formation. If the planets are indeed second-generation, the timescale for planetary orbit stability is likely to be $t\ll 10^{8}$ yr, the timescale of post-common envelope evolution (e.g., Hu et al., 2007), rather than the age of the constituent stars.

Lee et al. (2014) analyzed orbit stability of a two-planet solution for NY Vir using a numerical orbit integrator over a timescale of $10^{6}$ yr. They found that although the system exhibited large-scale instabilities, there were islands of stability in parameter space. In particular, for the outer planet orbit, there was a stable island near semi-major axis $a\sim 7$ AU and eccentricity $e\sim 0.15$, which is within one uncertainty of our solution for the outer planet.

We analyzed orbital stability with Mercury6 program. We used RADAU algorithm with an integration period of 1 day. The results are shown in Figure 3. The solutions given in Table 2 is stable for at least $10^{8}$ years, which meets the requirement of SG planets. We did more analysis with different sets of parameters. Since the mass and semi-major axis of the planets of different solutions of NY Vir are similar, the eccentricities of the two planets determine that stability of the system. More specifically, when eccentricities of both planets reach 0.2, or when eccentricity of one of the planets reaches 0.3, this system becomes unstable.

6 Summary

We report new eclipse timing observations of the short period sdB binary NY Vir . The new times of primary minimum do not agree with predictions from previously published one and two planet models calculated from O-C deviations from a linear ephemeris. We found a new stable two-planet plus a period derivative solution that is in agreement with all historical O-C timing data. However, we find that most model parameters (planetary masses, orbital parameters, and period derivative) are highly degenerate, with large cross-correlations between parameters. The degeneracy can be broken by additional timing data, but the required timescale for breaking the degeneracies is at least a decade.

7 Acknowledgments

This research made use of the LMFIT non-linear least-square minimization Python library, and Astropy333http://www.astropy.org a community-developed core Python package for Astronomy. We also made used of the Mercury orbit stability analysis program Chambers (1999).The Iowa Robotic Observatory is funded by the College of Liberal Arts at the University of Iowa. The project was also funded in part by NSF award 1517412.