Prediction of a Red Nova Outburst in KIC 9832227

Stellar evolution for individual stars and multiple star systems consists of both short- and long-lived phases. Transient events are short-lived phases that are sufficiently brief and luminous to be discovered by the time variation of their brightness. Ideally, the study of transient events would include a full characterization of the pre-outburst state. Such data would provide observational tests of outburst mechanisms as well as a fundamental check that the correct pool of progenitors is being connected to that transient type. Unfortunately, it has not been possible, to date, to predict the occurrence of any major transient type. This is not surprising owing to the rarity of events and the low luminosity of precursors. Moreover, theories do not generally predict specific observational signatures of systems in the final years before eruption.

With an ever-increasing observational emphasis on surveys, retrospective identification of precursors has been possible in some fortunate cases, allowing limited but instructive tests of theory. Prediscovery photometry of the Type II supernova SN 1987A in the Large Magellanic Cloud showed a blue supergiant—a surprising result that led to significant reconsideration of the Type II supernova mechanism (cf. Arnett et al. 1989, and references therein). Prediscovery photometry of the Type IIb supernova SN 1993J in M81 showed a similarly surprising color (Aldering et al. 1994). Prediscovery light curves from the Optical Gravitational Lensing Experiment (OGLE) survey of the classical nova V1213 Cen (Nova Centauri 2009) showed a very low quiescent level punctuated with dwarf novae outbursts in marked contrast to the brighter, steadier post-nova behavior (Mróz et al. 2016). This observation directly supports models that suggest the classical nova outburst abruptly alters the mass transfer rate.

Luminous red novae have only recently been identified as a distinct class of stellar transients (Kulkarni et al. 2007). They are characterized by relatively long outbursts, a range of luminosities intermediate between classical novae and supernovae, and red spectra. Prediscovery OGLE light curves of the luminous red nova V1309 Sco (Nova Sco 2008) revealed the light-curve shape of a contact binary system with an orbital period that was decreasing exponentially (Tylenda et al. 2011, hereafter Te11). This was a remarkable confirmation of the hypothesis of Soker & Tylenda (2003) that luminous red novae arise from merging contact binary stars.

Had intensive targeted observations of the precursors been possible in any of these cases, more specific details of the physical state of the precursor could have been determined and used to test theory more fully. For red novae, Te11 suggested the mechanism that triggers outburst is the onset of the Darwin instability. This occurs when the mass ratio is small enough that the companion star can no longer keep the primary star synchronously rotating via tidal interaction. Angular momentum transferred from the binary orbit to the primary's spin changes the orbit more than the spin, leading to a runaway. Specifically, Rasio (1995) found that this occurs at a mass ratio of 0.09 for a main-sequence primary star that is mostly radiative.

The central premise of this paper is that the most valuable use of Te11's discovery may be to consider the exponentially decreasing period as the first predictive signature of an imminent nova outburst. We consider generally how this signature should be sought in the rapidly increasing quantity of survey data that will be available in the near future. We consider specifically the fortuitous identification of a promising candidate red nova precursor, KIC 9832227 (Molnar et al. 2015), laying out the observational evidence in hand and specifying additional observations needed to fill out the story.

Eggleton (2012) describes contact binaries and common envelope evolution as the two great unsolved questions of stellar evolution. There is growing consensus on the mechanisms underlying some of the key stages of contact binary evolution. Initial binary separations are wide (orbital periods of days to months). Dynamic interaction with third star companions combined with tidal friction can reduce the binary period down to a day or two. See, e.g., Fabrycky & Tremaine (2007) for dynamical simulations and Tokovinin et al. (2006) for an observational survey of third star companions of close binaries. Magnetic braking and nuclear evolution may both play roles in bringing the binary from that point into contact. Once in contact, the stellar radii and Roche geometry do not allow for significant evolution of the orbital period. The end result for all orbital periods is a stellar merger. Soker & Tylenda (2003) found that red novae with a wide range of lumonosities can be produced by contact binaries with a range of initial masses.

There is no consensus on the details of mass exchange between the contact binary components or on the specific mechanism that triggers the merger. Eggleton (2012) favors a gradual mass transfer from the less massive (and hence smaller radius) secondary star to the primary driven by structural change of the secondary caused by the energy it receives from the primary. The contact is sustained by magnetic braking or nuclear evolution. This comes to an abrupt end when the mass ratio is extreme enough for the Darwin instability to trigger a merger. Alternatively, Stȩpień (2011) suggests that when a binary system first reaches contact, a brief intense mass transfer event ensues that ends with the originally more massive star becoming the less massive one. A stable contact phase only exists following the reversal. D'Souza et al. (2006) explore dynamic mass transfer without the Darwin instability. Koenigsberger & Moreno (2016) explore merger triggered by a tidal runaway based on a nonequilibrium response to tidal dissipation.

KIC 9832227 was discovered in the Northern Sky Variability Survey (NSVS) (Woźniak et al. 2004) in a catalog for RR Lyrae stars (Kinemuchi et al. 2006). The nature of the variability was unclear in the NSVS data, and the light curve obtained from Kepler further complicated the classification. Originally, KIC 9832227 was listed as a candidate c-type RR Lyrae star, but with a short periodicity, the target could be classified as something else (Kinemuchi 2013). The Kepler light curve revealed a periodicity at 0.23 days (pulsational) or 0.46 days (orbital), but with an unusual behavior of amplitude modulation that appeared to be periodic as well. The Kepler Eclipsing Binary Catalog (Kirk et al. 2016, and references therein) found a morphological parameter of 0.95, meaning that the identification as an eclipsing binary is uncertain.

Molnar et al. (2015, hereafter Me15) presented measurements (color versus orbital phase) that ruled out the pulsation interpretation. Double absorption line spectra presented in the present paper confirm directly the contact binary classification. Furthermore, Me15 found the orbital period determined from Kepler data was shorter than the period determined from prior ground-based observations over a nine-year time baseline—a negative time derivative of the period. The signature behavior of V1309 Sco was a negative period derivative that became more negative over time. Timing measurements of their own from 2013 and 2014 (subsequent to the Kepler data) showed that the period derivative was indeed becoming more extreme. Me15 found that all of the timing data together were well described by the exponentially decreasing period function used for V1309 Sco and identified the system as a red nova precursor candidate on that basis.

However, Me15 concluded by outlining two additional criteria that must be met before the system can be considered a likely precursor for a red nova outburst: (1) the magnitude of the period derivative should exceed the range observed for otherwise similar contact binary systems, and (2) the possibility that the observed period changes are caused by changing light travel time as the binary orbits a distant companion must be ruled out. In this paper we present the photometry presented by Me15 in detail for the first time, along with follow-up photometric and spectroscopic measurements made through 2016 September that address these issues.

To assess the feasibility of searching for red nova precursors, we estimate the number of precursors in our Galaxy that are currently detectable as such. Three values must be determined: the rate of nova outbursts observed, the magnitude of the period derivative required for a precursor to stand out from other systems (the threshold period derivative, ${\dot{P}}_{t}$), and the length of time ${\dot{P}}_{t}$ is exceeded before outburst. The expected number of precursors is simply the product of the outburst rate and the number of years the precursor exceeds ${\dot{P}}_{t}$.

Kochanek et al. (2014) estimated the Galactic rate of stellar mergers as a function of their peak luminosity. They found that events like the V1309 Sco outburst are observed about once per decade.

Contact binaries exhibit timing variations for a variety of reasons that are not all understood (see Conroy et al. 2014, for a list). Kubiak et al. (2006) used OGLE data to make a uniform study of 569 contact binaries with data spanning 13 years. They found a limit to the observed magnitude of $\dot{P}$ of $1\times {10}^{-8}$, with typical values less than a tenth of that. Positive and negative values occurred at similar rates, indicating that the dominant period change mechanism was cyclical (such as changes caused by orbital motion around a third star). We adopt $1\times {10}^{-8}$ as our threshold ${\dot{P}}_{t}$.

To determine the time that the period derivative will exceed the threshold, one needs a model of the time evolution of the orbit. Webbink (1976) and Webbink (1977) showed that an unstable phase of evolution can occur in a close binary, resulting in shrinkage of the orbit and merger, but no physically realistic model of the process has yet been constructed. We therefore use the phenomenological expression that Te11 found represented well the time dependence of the orbital period of V1309 Sco. Specifically, they used an exponential formula of the form

$$\begin{eqnarray}&&P(t)={P}_{0}\exp (\tau /[t-{t}_{0}]),\end{eqnarray} \tag{ 1 }$$

where P₀ is the initial orbital period, t₀ is the time at which Equation (1) goes to zero, and τ is the exponential timescale. (Note that the period does not have to go all the way to zero to trigger an outburst—for V1309 Sco the onset of outburst preceded t₀ by 1.8 year.) Hence the period derivative

$$\begin{eqnarray}&&\dot{P}(t)=-{P}_{0}\tau /{(t-{t}_{0})}^{2}\exp (\tau /[t-{t}_{0}]),\end{eqnarray} \tag{ 2 }$$

and $\dot{P}$ exceeds ${\dot{P}}_{t}$ from the time the threshold is reached, t_t, until t₀. Since τ is observed to be much shorter than a year, the exponential factor is nearly unity until the final year and

$$\begin{eqnarray}&&({t}_{0}-{t}_{t})\approx {({P}_{0}\tau /{\dot{P}}_{t})}^{1/2}.\end{eqnarray} \tag{ 3 }$$

For V1309 Sco $({P}_{0},\tau )$ = (1.4456 days, 15.29 days). Hence, when first observed by OGLE in 2001 August, the period derivative, $-2\times {10}^{-6}$, was already more extreme than ${\dot{P}}_{t}$. Using Equation (1), the time spent in excess of the threshold was 129 years. The outburst rate times this lifetime implies there should be 13 such systems in our Galaxy currently exceeding the threshold as they approach merger. As determined in Section 3.3 below, the precursor candidate KIC 9832227 has $({P}_{0},\tau )$ = (0.45796 days, 0.114 days), both lower values than for V1309 Sco, and hence $({t}_{0}-{t}_{t})$ of 6.2 yr. This yields an estimate of 0.6 systems exceeding the threshold on average.

The typical value for $({t}_{0}-{t}_{t})$ probably lies somewhere in between these two cases. The orbital period of V1309 Sco, 1.4 days, is generally considered to be unusually long for contact binaries. Indeed, W UMa systems are sometimes defined to have orbital periods shorter than a day. However, catalogs of ground-based data have a sampling bias against periods near one day. A period histogram of the 261 Kepler binary systems classified as contact binaries (Matijevic et al. morphology between 0.7 and 0.8) shows this period value to be smoothly connected with shorter values. Fully 5% of the systems had longer periods. With shorter main-sequence lifetimes, the more massive stars found in longer-period systems may evolve more quickly, and so may be underrepresented in the histogram relative to their birthrate. Hence, there is no reason to think the period of V1309 Sco is distinctive of systems that merge. The orbital period of KIC 9832227 is in the peak of the period distribution (0.25–0.50 days).

In summary, it is plausible that there are currently observable red nova precursor systems, although the number is statistically small. Furthermore, as $({t}_{0}-{t}_{t})$ is only six years for KIC 9832227, it is desirable in general to identify precursors as soon as they reach ${\dot{P}}_{t}$ to ensure an adequate time to study them before outburst. Hence it is important to use optimal methods of timing analysis.

In Section 2 we reanalyze the V1309 Sco photometry, developing the technique we later apply to KIC 9832227. In particular, we show that emphasizing the measurement of the timing of the orbital phase is more precise than direct estimates of the orbital period. In Section 3 we present the available KIC 9832227 photometry (from ground-based archives, the Kepler spacecraft, and four years of our own ground-based observations). We give special treatment to means of minimizing the effects of time-variable starspots on timing estimates. We find evidence of a small third body, component C, with an orbital period of 1.7 years. We compute an updated fit of timing data to the exponential formula. In Section 4 we present spectroscopic observations from the Apache Point Observatory (APO) and the Wyoming Infrared Observatory (WIRO). In Section 5 we seek a self-consistent model of the light-curve shape and the spectra. In Section 5.1 we estimate temperature and metallicity directly from the spectra and compare with temperatures estimated from Sloan colors (reported in the literature and measured by us). In Section 5.2 we peform a global fit to the spectroscopic data for the stellar rotational and orbital velocities for each of the two data sets. In Section 5.3 we use the light curves and the spectroscopic velocities to determine the physical parameters of the binary system. In Section 5.4 we use the spectra to place an upper limit on the mass of component C. In Section 5.5 we consider and rule out an alternative model of the long-term timing trend involving light travel-time delays due to a massive companion, a hypothetical component D. In Section 6 we briefly discuss three aspects of the merger process for which the case of KIC 9832227 may be instructive: the outburst trigger, the outburst timing relative to the parameter t₀, and changes in light-curve shape prior to outburst. Finally, in Section 7 we summarize our conclusions and outline further work needed.

The Northern Sky Variability Survey (NSVS) monitored KIC 9832227 from 1999 April to 2000 March (Woźniak et al. 2004). It appeared under two object identifications: 5597755 and 5620022, with 238 and 217 unfiltered measurements, respectively, as it straddled the boundary of two program fields. Median measurement uncertainties were 0.034 and 0.044 mag, respectively. We converted the times given in Julian Date (JD) to barycentric Julian Date (BJD) using the Python package astropy and the DE423 solar system ephemeris from JPL.⁴ (In this paper we follow the convention for modified Julian Date used in the Kepler Eclipsing Binary Catalog (Kirk et al. 2016): $\mathrm{MJD}\equiv \mathrm{JD}-2400000.0.$) To ensure a common zero-point for the two data sets, we computed a Fourier series fit to the average light curve folded on the exponential ephemeris of Me15 and computed the average residual for each data set. We subtracted 0.062 mag from the 5620022 data to place them on the same scale as the 5597755 data. We clipped four data points that had residuals relative to this average light curve greater than 0.15 mag. We iterated these calculations so that the fit used in the end reflected both the zero-point correction and the editing.

The All Sky Automated Survey (ASAS) monitored KIC 9832227 from 2006 May to 2008 January (Pigulski et al. 2009). It appeared under the identification 192916+4637.3 for the V-filter data and 192916+4637.4 for the I-filter data, with 83 and 104 measurements, respectively. ASAS computed light curves using several aperture widths. We followed the choices they made for the figures on their web page (number 2 for the V filter and number 0 aperture for the I filter). We converted the times given in Heliocentric Julian Date (HJD) into JD and then into BJD. While the two data sets spanned approximately the same time range, the I-filter data were mostly taken near the end of the time range. As with the NSVS data, we sought highly discrepant measurements and clipped two from the V data set and one from the I data set.

The Wide Angle Search for Planets (WASP) survey monitored KIC 9832227 intensively from 2007 May to 2008 August (Butters et al. 2010) and used the identification 1SWASP J192916.01+463719.9. These data were obtained with a broadband filter with a passband from 400 to 700 nm. We converted the times given in JD into BJD. Data were taken in four observing seasons (spanning 2004 May to 2008 August), but the first two seasons had too few measurements to establish an orbital phase. The third season was of generally good quality except for one night (JD 2454234), which had great scatter and which we did not use. Again, we sought highly discrepant measurements, setting the limit at 0.1 mag, which clipped 58 data, leaving 13,298 good measurements in two observing seasons.

The Kepler spacecraft fortuitously observed KIC 9832227 during its primary mission from 2009 May to 2013 January. Originally, this target was classified by the Kepler binary star working group, who classified it as an eclipsing binary star (Prša et al. 2011). Separately, this target was granted continued follow-up Kepler photometry through Guest Observer Director's Discretionary time, under proposal GO program #30024 during quarters 10–13. For all the quarters, this target was regularly observed in long-cadence mode (values averaged in 29.4-minute bins with typical brightness uncertainties of 0.07 mmag) for a total of 53,410 measurements. The Kepler passband spanned 423–897 nm. Details of the Kepler data processing and calibration for this target are provided in Kirk et al. (2016) and references therein. Kepler information on KIC 9832227 can be found in the Kepler Eclipsing Binary Catalog⁵ or at the MAST archive.⁶

KIC 9832227 was observed between 2013 June and 2016 September with Calvin College's two identical (0.4 m OGS Ritchey-Chrétien) telescopes: one operated remotely in Rehoboth, NM, at an elevation of 2024 m, and a second, located on the Calvin campus in Grand Rapids, MI, at an elevation of 242 m. The Rehoboth telescope has an SBIG ST-10XE camera with a plate scale of 1.97 arcsec per pixel (binned 2x2), while the Grand Rapids telescope has an SBIG ST-8XE camera with a plate scale of 1.58 arcsec per pixel. Given better viewing conditions, the New Mexico data have lower noise, but an average anticorrelation in weather conditions between the sites is advantageous for a monitoring program. Altogether, 32,421 images were obtained on 165 distinct nights: 23 in Grand Rapids, and 143 in Rehoboth.

Red filters were used every night to monitor timing, approximating the Kepler passband: Bessel R on 148 nights, and Sloan r on 17 nights, which yielded 28,289 red images in total. For nights during which only the red filter was used, images were obtained every 61 s in New Mexico (every 130 s in Michigan). In addition, Bessel B was used on 9 nights early on in order to establish the color change with phase. These data ruled out the possibility suggested by Kinemuchi (2013) that KIC 9832227 might be a pulsing star with an unusual light curve. Bessel V and I images were also obtained on 12 nights to test the need for altitude-dependent color corrections. We used MaxIm to perform differential photometry of the target relative to the average brightness of four reference stars, chosen for proximity in brightness, position, and color, and hence found no need for color correction. Nightly average brightness over the four years for each reference star brightness were examined to confirm lack of significant variability. The standard deviation of nightly averages for individual reference stars ranged from 3 to 7 mmag. The median random uncertainty on individual measurements of KIC 9832227 was 5 mmag for data taken in New Mexico and 7 mmag for data taken in Grand Rapids. The reference stars used were KIC 9771791, KIC 9832338, KIC 9832411, and GSC 03543:01558. The aperture had a radius of 5 pixels, the gap a width of 2 pixels, and the sky annulus a width of 3 pixels. Finally, on 11 of the nights with Sloan r, we also used Sloan g, i, and z filters in order to determine average colors and a photometric temperature (see Section 5.1).

As a complete light curve cannot be obtained in a single night, we group the observing nights into 38 epochs for purposes of timing analysis. Table 1 lists each epoch name, along with the first and last night included, the number of nights in which observations were made, which observing site(s) was used, how many images were obtained through the R or r filters, which other filters were used, and whether coverage of all orbital phases was complete. Figure 4 shows representative light curves from five epochs spanning 2016 May to July. Successive epoch are offset by 0.1 mag for clarity. Epochs with complete orbital phase coverage are plotted in red, while incomplete ones are shown in gray.

Zoom In Zoom Out Reset image size

The Kepler data set is by far the most precise and continuous of the data sets. Hence it is the natural choice to use for a cross-correlation template. The data set is also valuable for an analysis of short-timescale variations. The results inform the choice of analysis tools and limitations of interpretation.

Figure 5 shows in red a 10th order Fourier fit to all of the Kepler magnitudes along with black dots for the individual measurements (all folded on the exponential ephemeris given in Section 3.3). The variation of the average light curve spans 0.20 mag. The wide spread of the data around the average is due to real changes in the light-curve shape as described by Kinemuchi (2013). Figure 6 illustrates the nature of the variations by showing data from five consecutive orbits beginning on MJD 54983 in red and comparing them to data beginning on MJD 55118 in blue. The asymmetries between the maxima are reversed, as are the asymmetries between the minima. We interpret these asymmetries as due to starspot activity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The strength and arrangements of spots gradually vary with time. At times, a single starspot region dominates and can be seen drifting in longitude. Similar behavior is seen in the Calvin College data. This is illustrated in Figure 7, which shows the residual of the data from Figure 4 relative to the average Kepler light curve (Figure 5). These residuals are characterized by a variation of 0.05 mag with a single maximum that drifts in two months from orbital phase 0.3 in the top light curve to 0.75 in the bottom curve. We defer a detailed analysis of starspots to a later publication, in which we will explore variation of the drift rate over the years. The relevance to this study of orbital timing is that large and time-variable starspots can introduce an additional signal into timing determinations, an effect we wish to eliminate.

Zoom In Zoom Out Reset image size

Conroy et al. (2014) performed an analysis of the eclipse timing of the 1279 close binaries in the Kepler data set. Close binaries were defined as those whose Matijevič et al. (2012) morphology parameters exceeded 0.5. The authors computed O − C (for which they use the term "eclipse timing variations" or ETV) separately for primary and secondary eclipses using a chain of polynomials to identify "knots" (essentially points of ingress and egress). Their (online) plot of KIC 9832227 shows the timing of the two eclipses varying over a range of 35 minutes, largely out of phase with each other, a consequence of drifting and changing starspots. They also plot the timing of the average of the two eclipses, which reduces the range of variation, although there remain abrupt features in which the timing may go up (or down) by five minutes and return over the course of just one or two months. They identified 236 systems with an oscillatory signal in the timing consistent with light travel delays due to a third body. Of these, they listed 35 systems identified as having a likely signature of a third body (which they restrict to third-body periods shorter than 700 days so that two full cycles were covered by the data). They also listed 31 systems that are better fit with a constant period derivative (i.e., O − C plots consistent with a constant nonzero period derivative). Lacking an obvious overall trend in the plot, they did not include KIC 9832227 in either list.

We note in passing that we checked whether there are contact systems in the Kepler catalog that exceed the threshold ${\dot{P}}_{t}$ that we set using OGLE data. The Kepler data set has a much shorter time baseline than the OGLE data, so the constant period derivatives found may not represent sustained values as they may turn out to be portions of third-body oscillations with longer periods. We downloaded the Kepler data for the 31 systems that Conroy et al. (2014) listed as having a constant period derivative and computed the values of the derivatives. Of the systems that they also classified as contact binaries, the only system with a period derivative of either sign with a magnitude exceeding ${\dot{P}}_{t}$ was KIC 9840412, for which $\dot{P}$ was $-3.3\times {10}^{-8}$. We observed this system for two nights in 2016 July (using the equipment described in Section 3.1) and found the observed orbital phase 19.3 ± 0.5 minutes behind the phase expected, indicating that the magnitude of $\dot{P}$ for this system has subsequently decreased. Hence, the Kepler data set contains no exceptions that would require adjustment of ${\dot{P}}_{t}$ to a more extreme value.

For our first pass evaluating the timing of the Kepler data of KIC 9832227, we made a template like that of Figure 5, but folded on the linear ephemeris of Conroy et al. (2014). We then performed a set of cross correlations to determine the O − C for segments of the data spanning five consecutive orbits. While optimal timing of wide binaries is done by analysis of eclipses, cross correlation makes sense for contact binaries as the light is continuously variable. Furthermore, similar to the averaging done by Conroy et al. of primary and secondary eclipse times, cross correlation averages over the entire cycle and removes the first-order timing distortions of starspots. This is plotted with blue circles in Figure 8(a) and looks very similar to the average O − C figure of Conroy et al. We connect adjoining measurements with lines to emphasize that observational noise is negligible in this plot. The few-months-long excursions of up to five minutes are resolved in time, with a series of independent measurements that proceed smoothly. These excursions are residual effects of changing starspot patterns. We also plot a quadratic fit to the cross correlations. This shows a long-term trend with a negative period derivative, $-5.4\times {10}^{-9}$, which is among the most extreme contact binary period derivatives found by Conroy et al.

Zoom In Zoom Out Reset image size

We wish to remove the excursions as much as possible in order to see the underlying trend more clearly. We explore the potential for doing this with Fourier components by plotting the amplitude of the first six components versus time for five orbit data segments (Figure 9(a)). The modulation due to orbital geometry is principally in the even-numbered orders, plotted with solid lines. These are stronger and less variable than the odd-numbered orders (plotted with dashed lines), which are largely due to the time-variable spots. We also plot the amplitudes of the Fourier components of the residual of the data to the folded average light curve (Figure 9(b)), thereby removing the contribution from orbital geometry. Since light modulation of starspots is a consequence of rotation (and therefore gradual), we expect most of the power in the first order, with progressively less power for increasing orders, as Figure 9(b) shows.

Zoom In Zoom Out Reset image size

Hence use of even-numbered Fourier orders for timing should suppress the starspot contribution. The best order to use will be a tradeoff between the improved separation of signals at higher orders and the increased observational noise at higher orders. We plot the second-order determination in red with squares in Figure 8 along with the cross correlations. The similarity of shape indicates that this order is no better at separation than the cross correlation. The vertical offset indicates a 3-minute error in the time of zero phase reported by Conroy et al.

In Figure 8(b) we plot the second, fourth, and sixth orders with five-minute offsets inserted between them so that the orders may be seen separately. The short-timescale excursions seen in the second order are only half as long in the fourth order and are essentially gone in the sixth order. The observational scatter that is invisible in the second order becomes significant by the sixth order. Nonetheless, the sixth order is the cleanest indication of long-term trends. In addition to the general curvature seen in the cross correlation, there is also a broad dip around 2011.

For our second pass, we plot the sixth order again (Figure 10), but this time relative to the exponential ephemeris of Section 3.3. For the Kepler epoch, this essentially means removing a period derivative of $-3.4\times {10}^{-9}$. The parameters of that ephemeris are determined by comparison of the Kepler timing to earlier and later measurements. Hence the relatively flat result is an indication that the period and period derivative of the exponential ephemeris apply within the Kepler data.

Zoom In Zoom Out Reset image size

This figure also shows a small oscillation of the O − C around the average, which we fit with a simple sine wave of semiamplitude 1.69 ± 0.03 minutes and a period of 590 ± 8 days and plot as a blue line in Figure 10. Conroy et al. (2014) interpreted O − C diagrams like this as caused by light travel-time variation arising from the orbit of the binary around a third body. Although the time signature of a third body may be more complex than a sine wave (for nonzero eccentricity), the data are not of adequate quality to justify a more complex fit. We designate this third-body component C. Assuming a binary mass of 1.71 ${M}_{\odot }$ (as determined in Section 5.3), then ${m}_{C}\sin {i}_{C}$ is 0.11 ${M}_{\odot }$, where m_C is the mass of component C and i_C is the inclination of its orbit to our line of sight. Finally, our best estimate of the time of zero orbital phase in the Kepler data set comes from setting the average of the sine wave to zero: BJD₀ = MJD 55688.49913(2).

Me15 fit the KIC 9832227 timing data taken through 2014 to the exponential formula of Equation (1) and found the best-fit parameters to be $({P}_{0},\tau ,{t}_{0})$ = (0.4579603 days, 0.089 days, 2020.8(15)). Updating this fit using data through 2016 September, we find $({P}_{0},\tau ,{t}_{0})$ = (0.4579615(9) days, 0.114(14) days, 2022.2(7)). The numbers in parentheses indicate the uncertainty on the final digit(s). The new values are consistent with the previous ones to within their uncertainties. The most relevant change is the reduction in the uncertainty on t₀, which is sensitive to data taken near t₀. In this section we present the details of the fitting and plot the results in several ways to give deeper insight.

Figure 11 shows the residuals of all of the orbital timing data relative to our new exponential ephemeris and the BJD₀ from the preceding section. The three pre-Kepler data sets are relatively sparse, so a single cross correlation was done for each observing season. The horizontal error bars indicate the time range covered for each timing measurement. An advantage of averaging over data that are randomly spread over a long time range is that it should also average out residual effects of variable starspots. The Kepler timing estimates come from the sixth-order phases shown in Figure 10, but averaged in groups of seven to reduce the crowding. The Calvin estimates use second-order phases rather than sixth order as there is less random noise. As a consequence, these values are precisely measured, but may contain starspot-induced excursions like those seen in the second-order Kepler data (Figure 8(b)). The 25 epochs with complete phase coverage were used in the fit. The 13 epochs without complete phase coverage show similar values but are omitted as systematic effects of starspots are greater in them. The effect of component C exceeds the measurement uncertainties in both the Kepler and Calvin data, and the residual effect of starspots exceeds the measurement uncertainties in the Calvin data. These two effects should have no systematic influence on the exponential fit parameters as they are small in amplitude compared to the long-term drift in the timing and oscillate on short timescales. However, they do increase the effective measurement noise. Hence the uncertainties used in the fit for the Kepler and Calvin data sets were based on their respective residual scatter.

Zoom In Zoom Out Reset image size

Figure 12 shows the timing data relative to a constant reference period, the standard O − C plot. As with V1309 Sco (Figure 2), the reference period was chosen to be the actual period near the middle of the data set. The sharper curvature in later years than in earlier years is again the mark of an increasingly negative period derivative with time. The exponential fit of Me15 is marked by a black dashed line, while our revised fit is given as a solid purple line. The 18 new Calvin data points from 2015–2016 (filled purple circles) follow the prediction of the exponential line. As forecast by Me15, the magnitude of the period derivative exceeded ${\dot{P}}_{t}$ in 2015. Unlike the case of KIC 9840412 (described in the preceding section), there is no evidence to date of a reversal in the timing trend. After presentation of our spectroscopic data in Section 4, we discuss in detail the limits that can be placed on alternative interpretations of the long-term trend in timing as light travel-time delay (Section 5.5).

Zoom In Zoom Out Reset image size

Although the exponential model parameters were fit directly to the O − C timing data, it is instructive to plot period versus time as well to see directly the amount of change in the period. Figure 13 shows the exponential formula period with time along with three period determinations made by separately fitting the slope of the O − C plot for pre-Kepler, Kepler, and Calvin data, respectively. As before, uncertainties in the Kepler and Calvin data are determined using the actual scatter in the timing data. The points are plotted at the times of the midpoint of the time ranges of the respective data sets. The periods thus derived, listed in Table 2 and plotted in Figure 13, show very significant changes relative to the uncertainties. It is likewise instructive to plot the period derivative with time, Figure 14, to see directly the change in this quantity. Two estimates of the period derivative are plotted that are based on the differences in the periods in the preceding plot. These are marked at the times midway between the times of the period estimates. The apparent disagreement with the model curve is an artifact of the curvature of the model over the time range used to make the estimate. The increasingly negative period derivative with time is the essential signature of the merging star model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1.1. Photometric Estimates

5.1.2. Spectroscopic Determination

Spectra were analyzed for ${T}_{\mathrm{eff}}$ and [Fe/H] determination using the iSpec Python tool developed by Blanco-Cuaresma et al. (2014b). This tool provides a number of different line databases and atmospheric models to allow the user to determine atmospheric parameters of the input stellar spectrum by fixing input parameters that are well known while fitting for unknown parameters. The tool also allows for determination of radial velocity by cross-correlating the data against any of a number of template spectra. Procedures were generally followed as described by Blanco-Cuaresma et al. (2014a), the iSpec user manual, and S. Blanco-Cuaresma (2015, private communication).

The normalization of the APO spectra was checked over the reduced wavelength range and adjusted as needed. Normalizing the extracted orders of echelle spectra is notoriously challenging, particularly in low-S/N spectra, and in some instances, small adjustments to the IRAF normalization were needed. The renormalization was done with iSpec using several spline and low-order polynomial fitting models. These were repeated to identify systematic effects on the final parameter fit values.

APO spectra from orbital phase ϕ = 0.0 were coadded into a master ϕ = 0.0 APO spectrum with maximum S/N. The same was done for all phase ϕ = 0.5 APO spectra, as well as WIRO spectra from these phases. These four master spectra were fit using the MARCS model atmosphere of Gustafsson et al. (2008), solar abundances from Grevesse et al. (2007), and the most recent atomic line list including hyperfine structure transitions (Heiter et al. 2015) from the Gaia ESO survey (Gilmore et al. 2012). Owing to the high rotational speed of this contact system, spectral features were Doppler-broadened ($v\sin i\sim 150$ km s⁻¹) to the extent where the data were not sensitive to macroturbulence and microturbulence, so default values of 0.0 km s⁻¹ and 1.0 km s⁻¹, respectively, were adopted. Small increases in these parameters did not affect the results. A precise surface gravity could not be determined either, and so $\mathrm{log}\,g$ was fixed at 4.0. Results were consistent with tests using 4.5 as well.

Within iSpec the initial value of ${T}_{\mathrm{eff}}$ was varied systematically between 5800 K and 5900 K (spanning the KIC value) to check for consistency in the output fit values, while initial values for [Fe/H] were varied between −0.05 and +0.05. The best-fit spectra for phases 0.0 and 0.5 are shown in Figures 15 and 16, where the surface temperature and metallicity for the system were determined to be ${T}_{\mathrm{eff}}$ = 5828 ± 42 K (random) ± 40 K (systematic) and [Fe/H] = −0.04 ± 0.03 (random) ± 0.04 (systematic). The largest source of systematic uncertainty was the renormalization of the echelle spectra, while the random error is a combination of fitting error reported by iSpec and the scatter in the output fit values reported while the parameters were being incrementally varied.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Throughout this work, we use linear limb-darkening coefficients from van Hamme (1993),⁷ using $\mathrm{log}g=4.2$ (as found for the primary star in Section 5.3), the temperature and metallicity given above, and the relevant filter or wavelength range. For the APO and WIRO spectra, the coefficients were 0.55 and 0.511, respectively.

An effective way to digest the velocity information in a spectrum is to compute the broadening function, a least-squares deconvolution of the data that represents spectral power as a function of radial velocity. Rucinski (2004) reviewed the advantages of this approach over the simpler cross-correlation function. Like the cross correlation, this technique requires a template spectrum for comparison. Fortunately, in the case of a contact binary system, both components have very similar atmospheric temperature and surface gravity, so one template works equally well for both stars. Coelho et al. (2005) published a grid of model atmospheres, from which we selected the one with solar abundances and ($T,\mathrm{log}g$) = (5750 K, 4.0). To compute broadening functions, we selected spectral regions with relatively good S/N ratios and that were free from strong terrestrial lines. For the APO data we used the wavelength range 373.0–686.6 nm, smoothed with a 3.5 km s⁻¹ Gaussian with 3.5 km s⁻¹ bins. For the WIRO data we used the wavelength range 539.5–669.4 nm, with 7.0 km s⁻¹ bins. We made a sequence of three analysis passes.

5.2.1. First Pass: Fitting Gaussians to Individual Spectra

The first pass follows the typical application of the broadening function (e.g., Rucinski 2002). A pair of Gaussians are fit individually to each spectrum. Then the Gaussian center velocities as a function of orbital phase are fit to find the radial velocity semiamplitudes of the two stars, K_A and K_B, and the heliocentric radial velocity of the system barycenter, γ, where the subscripts A and B refer to the primary and secondary stars, respectively (components A and B). Figure 17 shows the radial velocity curves for each component. The sinusoidal fits shown yield (K_A, K_B, γ) = (48.9, 226.6, −22.2 ± 0.8) km s⁻¹. This implies a stellar mass ratio ${m}_{B}/{m}_{A}={K}_{A}/{K}_{B}=0.216$.

Zoom In Zoom Out Reset image size

5.2.2. Second Pass: Fitting Gaussians to the Combined Data Set

5.2.3. Final Pass: Fitting Rotationally Broadened Spheres to the Combined Data Set

A better approach still is to build on the second pass approach by also using a physically motivated model: two rotating spheres. For a synchronously rotating contact binary, rotational broadening is comparable to the orbital velocity and larger than other sources of line broadening. Since the wing of a rotationally broadened line cuts off more precipitously than that of a Gaussian broadened line, forcing a Gaussian shape will systematically shift the measured radial velocity of the smaller star. The fitting parameters we use now are K_A, K_B, γ, ${v}_{A}\sin i$, ${v}_{B}\sin i$, R, and A, where the new parameters are the projected stellar rotational velocities, the ratio of the integrated signal from the second star to that of the first, and the sum of the integrated signal from both stars. In addition to this approach providing a better fit, three of the new parameters (the rotational velocities and the signal ratio) are physically meaningful values that will be used as a consistency check when modeling the system geometry. A final element of the model is convolution with a fixed Gaussian (to account for instrumental effects). We determined the appropriate Gaussian widths by fitting Gaussians to our standard star spectra: 5.5 km s⁻¹ and 31.7 km s⁻¹ for APO and WIRO, respectively. These values are consistent with the combined effects of instrumental resolution and smoothing. The model spheres included limb darkening, with coefficients 0.629 and 0.529 for APO and WIRO, respectively.

To exemplify the results, Figure 18 shows model fits to one APO and one WIRO broadening function. The orbital phases are labeled in the upper right corner.

Zoom In Zoom Out Reset image size

The final APO fit yields (K_A, K_B, γ,${v}_{A}\sin i$, ${v}_{B}\sin i$, R) = (49.1, 215.5, −22.9, 152.0, 88.9, 0.304), while the WIRO fit was (49.8, 216.0, −22.2, 150.3, 96.7, 0.311). (All units are km s⁻¹ except for the dimensionless ratio.) The residuals to this model were 12% smaller than the second pass residuals for the APO data (10% for the WIRO data), indicating a significantly better model shape. Likewise, the agreement between the two data sets was better than for the Gaussian model, consistent with elimination of a systematic problem. The lone exception is the rotational velocity of the secondary star. For this we conclude that the WIRO data (which have a lower spectral resolution) are not sensitive to this parameter. Otherwise, the differences between the two fits are consistent with random uncertainties on fit parameters. In Section 5.3 we show that the rotational velocities and the ratio parameter are also very consistent with the geometric model of all the data.

The most important result of the velocity study is the mass ratio, which is 0.228/0.231 for the APO/WIRO data. These values are ∼6% greater than our first pass value.

In Table 4 we summarize the physical parameters of KIC 9832227 as inferred from the time-averaged Kepler light curve (Figure 5) and the spectroscopic temperature and radial velocities. First, we used the APO mass ratio in a fit of the light curve made with Binary Maker 3.0 (Bradstreet & Steelman 2002) to find the system inclination (i), fillout parameter (f), temperature ratio (${T}_{B}/{T}_{A}$), and the parameters of one starspot. Fixed input parameters include 0.549 for the limb-darkening coefficient, 0.32 for the gravity-brightening coefficient (appropriate for convective atmospheres Lucy 1967), and 0.5 for the reflection coefficient (again appropriate for convective atmospheres Ruciński 1969). As we show in Section 5.4, component C is so small that no spectral lines can be detected from it. Accordingly, we do not expect its continuum emission to significantly affect this light-curve analysis.

Table 4.  KIC 9832227 Binary Physical Parameters

Stellar Parameters
	Primary	Secondary	Total
$m/{M}_{\odot }$	1.395(11)	0.318(5)	1.714(12)
$L/{L}_{\odot }$	2.609(55)	0.789(17)	3.398(72)
${R}_{V}/{R}_{\odot }$	1.581(9)	0.830(5)
T (K)	5800(59)	5920(63)	5828(58)
$v\sin i$ (km s−1)	149.7(8)	84.7(5)
$\mathrm{log}g$	4.19(1)	4.10(1)
System Parameters
i (deg)	53.19(10)
f	0.430(15)
TB/TA	1.021(6)
mass ratio	0.228(3)
area ratio	0.279(6)
light ratio (R)	0.307(10)
distance (pc)	565(16)
E(B–V) (mag)	0.03(2)
a (${R}_{\odot }$)	2.992(15)
P0 (days)	0.4579615(9)
τ (days)	0.114(14)
t0 (year)	2022.2(7)
BJD0	MJD 55688.49913(2)
Coordinates
R.A., Decl. (J2000)	19h29m15954	+46°37'1988
Galactic (l, b)	788007	133557

Download table as:  ASCIITypeset image

The inclination is largely determined by the depth of the orbital modulation. The fillout parameter is a convenient equivalent of the surface equipotential scaled such that it is zero for a system precisely filling its Roche lobe and unity for a system reaching the outer critical Roche potential. It is determined by how rounded the eclipses are. The value, 0.43, indicates an intermediate depth of contact between the two stars. The temperature ratio is determined by the relative depths of the eclipses. Temperature ratios exceeding unity (1.021 in this case) are typical of contact binaries with shorter orbital periods. Note that we have chosen orbital phase zero to be the point at which the secondary star eclipses the primary one, although due to the temperature ratio, this is not the deeper of the two eclipses. The need for a starspot is evident from the asymmetric maxima. It is not generally possible to determine a unique starspot configuration from the light curve. We fit a single spot on the primary star with radius of 96 and colatitude 60°. We solved for the spot longitude (276°) and fractional temperature (${T}_{\mathrm{spot}}/{T}_{A}=0.825$). These are physically reasonable values, and they yielded an excellent fit to the light-curve shape.

The difference between the time-averaged data and the best-fit model implies a standard error on the data of only 0.0015 mag. The uncertainties given for the Binary Maker parameters (i, f, and ${T}_{B}/{T}_{A}$) are based on the assumption of random errors of this size. As no correction was made for the small smoothing effect of the 29.4-minute binning of the Kepler data, there are also systematic effects that are probably comparable to the random uncertainties.

Once the inclination is known, the radial velocities and the orbital period can be used to find the total mass, the semimajor axis (a), and the individual stellar masses. The semimajor axis in turn sets the scale for the stellar size. Since the stars are not spherical, several distinct equivalent radii may be of use. The radius (R_V) in the table is the radius of a sphere with the same volume as the star. The surface gravity (g) makes use of this radius. The rotational velocity ($v\sin i$) in the table makes use of the back radius (the distance from the star's center of mass to the equatorial point farthest from its companion. The APO model fit rotational velocities lie between these values and the values one would obtain measuring radius from the center of mass to the inner Lagrange point, a confirmation of the self-consistency of the physical parameters. The area ratio is the ratio of the surface areas of the two stars. The individual stellar temperatures are determined in an iterative process such that the brightness-weighted average of the temperatures is the average temperature determined from spectroscopy. The light ratio is the ratio of stellar brightnesses taking into account both the different surface areas and temperatures. This is within 1% of the independently determined ratio from the APO spectroscopy, another confirmation of the self-consistency of the physical parameters.

The distance was estimated starting from the absolute magnitude of the Sun in each of the four Sloan filters (Hill et al. 2010), scaling for the Planck intensities at those wavelengths and for the surface areas, and comparing them to our measured maximum brightness in those filters. This approach requires a correction for interstellar extinction. Maps of Galactic extinction provide an upper limit to the extinction as our target is within the Galaxy. We estimated interstellar extinction by varying E(B–V) until the dereddened colors corresponded with our spectroscopic temperature using the color-temperature relations of Boyajian et al. (2013) and the extinction curve of Yuan et al. (2013).

In the table the uncertainties in the temperatures include the random and systematic uncertainties added in quadrature. In general, uncertainties on all table quantities were computed by propagating the uncertainties of the two radial velocities and the three Binary Maker parameters as appropriate for each quantity.

For completeness, the table also includes the time of zero orbital phase from Section 3.2 and the parameters of the exponential timing fit from Section 3.3. We conclude the table with ICRS and equivalent galactic coordinates from the UCAC4 astrometric catalog, in which the star is designated 684–067807 (Zacharias et al. 2013).

In Section 3.2 we found evidence in the timing data for a component C with an orbital period P_C = 590 days and mass ${m}_{C}\sin {i}_{C}$ is 0.11 ${M}_{\odot }$. If the orbital plane of component C is more nearly face-on, the mass would be much greater than 0.11 ${M}_{\odot }$. We can place an upper limit on the m_C by looking for its signal in the broadening function.

As noted by Rucinski (2002), it is often easier to see the signature of a third body in a broadening function than that of a close binary as the tidally enforced corotation of the binary components makes their signal broad and short. For our considerations, we use the empirical main-sequence relations of Demircan & Kahraman (1991) to estimate radius, temperature, and luminosity from hypothesized mass. A star hotter than the binary can be immediately ruled out as it would dominate the photometric colors. Likewise, a star similar in temperature to the binary would show a strong narrow spike in the broadening functions in Figure 18. On the other hand, there is no hope of seeing the spectral signature of an $0.11\,{M}_{\odot }$ star, which would have a luminosity just 0.04% of the binary luminosity. We show below that stars of 0.5 or 0.6 ${M}_{\odot }$, with 1.9 or 3.8% of the binary luminosity, respectively, define the limits of what could be detected in the broadening function.

We take three steps to optimize the narrow nonvariable signal expected: (1) use the higher resolution APO spectra,(2) average together all of the spectra from the night with the greatest number of images (2015 June 21), and (3) deconvolve with a template the temperature of the hypothesized star (4000 K and 4250 K for 0.5 or 0.6 ${M}_{\odot }$, respectively). The signal from the binary will be smeared in velocity from the time averaging and diminished in amplitude because their spectra differ from the cooler template. Figure 19 shows the broadening functions optimized for 0.5 ${M}_{\odot }$ (a) and 0.6 ${M}_{\odot }$ (b) as black lines. The binary stars produce the very broad peaks at the center of each plot as well as the spikes at the right and left edges, which are artifacts of the temperature mismatch.

Zoom In Zoom Out Reset image size

If component C were 0.5 or 0.6 ${M}_{\odot }$, we would expect to see a spike near the system velocity of the binary in Figure 19. The question is how the size of the expected spike compares to noise fluctuations in that portion of the broadening function. We address this question with a simulation. We add the template to the June 21 data, binning to match the wavelength resolution and scaling the amplitude for the correct relative strength (for the given temperature and stellar surface area) at each wavelength. No additional noise is added as the Poisson noise of the small additional component is much lower than the random noise that is already present in the data. The broadening functions of these simulated companions are plotted in Figure 19 with gray lines. The black and gray lines differ only near the velocity of the hypothesized companions. In both cases, the gray spike far exceeds the random variation of the observed function. Hence, we conclude that m_C is lower than 0.5 ${M}_{\odot }$, and i_C is greater than 13°.

While the long-term trend in the timing data of KIC 9832227 remains consistent with an exponential formula, we cannot consider it a likely interpretation until we can reject the more mundane alternative of light travel delay due to a massive companion star orbiting the KIC 9832227 triple system with a period of 25+ years. (We designate this hypothetical companion component D.) This is especially true since the $\dot{P}$ threshold has only just been reached. There is nothing to exclude such a body a priori, although it is not required here to explain the formation of the close binary as component C is adequate for that.

Unlike the exponential formula, which can only have one shape in the timing plot with two scale factors, the various orbital parameters allow companion orbit models to exhibit a wide range of shapes. For our calculations, we assume that the companion orbital plane is inclined by 90° (i.e., includes our line of sight, hence maximizing its effect), the mass of the triple is 1.82 ${M}_{\odot }$ (using the minimum mass for component C), and the reference period used in the timing plot is the intrinsic binary period. The latter assumption is for specificity and has negligible affect on our conclusions. For each assumed companion mass (m_D), we fit the orbital period (P), eccentricity (e) angle between the line of sight and the line of apsides (ω), and time of periastron passage (t_p). These models account for the recent epoch of large period derivative by having an eccentric orbit with the time of periastron near the current date. They predict an abrupt turnaround in the timing plot in the near future and ultimately periodic behavior in the long run on the orbital period. Best-fit models for m_D = 0.8, and 1.0 ${M}_{\odot }$ are shown as red and green lines in Figure 12. The model parameters, (m_D, P, e, t_p, ω) are (0.8 ${M}_{\odot }$, 25.1 yr, 0.609, 43°, 2016.5) and (1.0 ${M}_{\odot }$,30.7 yr, 0.569, 53°, 2017.2). The 0.8 ${M}_{\odot }$ model is a marginal fit, differing by 5 minutes for the most recent data (which could be accounted for by fortuitous starspots). The 1.0 ${M}_{\odot }$ fits the existing data well.

We conclude that the minimum mass for a hypothesized component D is 0.8 ${M}_{\odot }$. As discussed in the preceding section, any nondegenerate star of this mass would easily be detected in the broadband colors or the spectroscopic broadening function. We are unaware of any known degenerate companions to triple systems, so this hypothesis itself is unlikely.

Nonetheless, this alternative model can be tested with additional data. The timing prediction rapidly diverges from the exponential formula in the next year. Additionally, the barycenter velocity of the triple system is predicted to rise sharply in the next few years and then level off as the companion moves away from periastron. Figure 20 shows the radial velocity of the triple system as it orbits a $1.0\,{M}_{\odot }$ component D (solid black line). The radial velocity of component D is shown as a dashed red line.

Zoom In Zoom Out Reset image size

While a number of papers have modeled the physics of the observed outburst of V1309 Sco, none cover the full decade preceding the outburst. Hence there is no physical basis for the exponential formula (Equation (1)) used to describe the early period changes, nor any insight into what sets the timescale, τ. Te11 suggest the relevance of both the onset of the Darwin instability and dynamically unstable mass loss. Whether the tidal viscosity timescale is adequate for the Darwin instability in the case of a contact system has been discussed back and forth by Pejcha (2014), Nandez et al. (2014), and Koenigsberger & Moreno (2016). We note here only that the observed mass ratio for KIC 9832227, 0.23, is far from the expected value for the onset of the Darwin instability, 0.1. As models are developed covering the decade approaching outburst, a wider range of mass ratios should be considered.

As noted in the introduction, the beginning of the V1309 Sco outburst preceded t₀ by 1.8 years. What might this time difference be in other systems? If the trigger is a certain rate of energy deposition, the rate of loss of orbital energy scales with $\dot{P}$. For fixed $\dot{P}$, the time offset scales with ${({P}_{0}\tau )}^{0.5}$ (Equation (3)). Alternatively, if the trigger is a certain fraction of orbital energy deposited, it may occur after the same change in $P/{P}_{0}$, so the time offset would scale with τ (Equation (1)). At the time of the observed outburst of V1309 Sco, $\dot{P}\approx -5\times {10}^{-5}$ and P was just 2.2% less than P₀. Given the much smaller τ found for KIC 9832227, these two milestones are not reached until a month or so before t₀. Hence for KIC 9832227, we expect t₀ to be a good estimate of the actual time of outburst.

We note for context that even the currently observed value of $\dot{P}\sim -1\times {10}^{-8}$ implies a high energy dissipation rate. As an illustrative calculation, conservative mass transfer from the primary to the secondary at a rate of $1.1\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ would produce the observed $\dot{P}$ and infer a loss of orbital energy at a rate of 19 ${L}_{\odot }$, more than five times the observed system luminosity. As exemplified by the case of V1309 Sco, only a small fraction of the lost energy is emitted from the surface as light. The primary star might also store the lost energy as increased thermal or gravitational potential energy.

In the final observing season before merger, the shape of the light curve of V1309 Sco changed significantly, with one of the two maxima gradually disappearing and the average system brightness decreasing. Te11 suggested that as the system approached merger, the primary rotated subsynchronously. Energy dissipation due to friction occurring near the L1 point would then heat the portion of the primary that would then rotate into view near orbital phase 0.25, increasing the strength of that maximum. We suggest alternatively that mass lost through the L2 point would accumulate in a disk that would be initially denser on the side viewed at orbital phase 0.25. This material could cool and obscure the system at this phase, decreasing the strength of that maximum. Without pre-outburst spectroscopy of V1309 Sco, we do not know for sure which maximum was at orbital phase 0.25, and so the observations do not decide between these two suggestions. Should a similar transition occur in KIC 9832227, however, we will know how to interpret it. Shu et al. (1979) show that mass lost through the L2 point escapes from the system for binaries with mass ratios in the range 0.064–0.78 (a range including the value of 0.23 we find for KIC 9832227), resulting in disk-like outflows. Pejcha et al. (2016) model radiation from such outflows, showing that results similar to red nova outbursts can be obtained.

The central premise of this paper is to use the exponentially decreasing orbital period found by Te11 to precede the red nova outburst of V1309 Sco as a signature to identify the precursor to the next outburst. The distinctive shape of the exponential formula (fixed except for scale factors) reduces the chance of false positives. We identify a threshold period derivative (${\dot{P}}_{t}\approx -1\times {10}^{-8}$) that should be exceeded to qualify as a serious precursor candidate. We estimate that our galaxy should have roughly 1–10 observable precursors of events at least as large as the V1309 Sco event (which brightened by 10 mag).

We develop the observational case for the precursor candidate suggested by Me15, KIC 9832227. The development exemplifies analytical tools that are needed to identify a precursor in a timely fashion. Determination of orbital phase versus time (the O − C plot) makes the best use of timing information. Computations relative to a tentative exponential ephemeris avoid smearing a signal in the presence of a rapidly changing period. Fourier filtering can help separate the competing signal of time-variable starspots (common in contact binaries with shorter periods) from that of true orbital changes. Spectroscopic analysis of radial velocities is more precise when a global solution is determined and avoids significant systematic errors when a model with realistic rotational broadening is used.

Me15 found that timing data for KIC 9832227 spanning 1999–2014 were consistent with an exponential fit with parameters. Photometric and spectroscopic observations presented in this paper make two strong tests of the precursor hypothesis in addition to establishing numerous physical parameters of the system. The most important test passed is that timing data through 2016 September are consistent with the extrapolation of the earlier fit. Along with this, we find that ${\dot{P}}_{t}$ has been exceeded as of late 2015. As a byproduct of the timing analysis, we discovered evidence for a small third star (component C with mass 0.11–0.5 ${M}_{\odot }$ and P_C = 590 days). A companion of this sort is expected as its interaction with the binary may have drawn the binary together in the first place. The physical parameters of the binary system are listed in Table 4 as deduced from a fit to the average light curve and the orbital velocities. Of particular note is the mass ratio, 0.23, which is inconsistent with the theoretical expectation of 0.1 based on the assumption that the onset of the Darwin instability is the trigger for merger.

The second test passed is that we rule out an alternative interpretation of the long-term trend in the timing data in which a luminous companion (a component D) orbits the triple system and inducing light travel-time delays. Specifically, we find that a component D with a minimum mass of 0.8 ${M}_{\odot }$ is required to match the timing data, but that the absence of an observed spectral signature rules out a main-sequence star of 0.5 ${M}_{\odot }$ or greater. The possibility of a degenerate component D remains, although such a system would in itself be a unique discovery.

Passing these two tests greatly strengthens the case for this candidate. However fortuitous the identification of this precursor candidate seems, the merger hypothesis has had predictive power and we currently have no alternative explanation for its timing behavior. The prediction of a naked-eye nova (second magnitude if it brightens as much a V1309 Sco did) in the year 2022.2 ± 0.6 makes it essential to use the available time to fullest advantage. A full characterization of the system properties (and perhaps their evolution over the next few years) will provide much more specific constraints of physical models of the outburst process.

Additional work is needed to develop the physical theory of mergers, observe the properties of the KIC 9832227 system, and survey other contact binary systems.

While much progress has been made on the physics of a red nova outburst, no detailed consideration has yet been given to the decade preceding outburst. This is difficult because of the long timescale, but might be tractable if treated as a series of quasi-steady-state stages. It is important to develop a physical basis for the exponential formula, and in particular to estimate the timescale as a function of orbital period and mass ratio. In the present analysis the timescale was fit to the data as the only significant free parameter. A theoretical constraint on that parameter could be a third significant test of precursor status. It would also be instructive to develop the theory of the geometric distribution of energy dissipation as merger approaches. Such a theory could be tested by observation of long-term trends in light-curve shape and system luminosity.

A variety of additional observations of KIC 9832227 are needed. First priority is continued monitoring of the orbital timing. With each year closer to t₀, the prediction of the merger hypothesis becomes more distinctive and hence a more stringent test of the hypothesis. Along with this, based on the example of V1309 Sco, we may expect changes in the average brightness and the light-curve shape when very close to t₀.

A series of high-precision spectroscopic observations of the system velocity of the binary can be used to confirm the existence of component C. Based on the timing variation, we expect sinusoidal variation with a semiamplitude of 1.9 km s⁻¹ on the 590-day period. The same observations can also test the alternative hypothesis of a degenerate component D. This hypothesis predicts a gradual increase in system velocity of 10 km s⁻¹ over the next several years, followed by a long period of no change.

A detailed study of the time variation of starspots would be a probe of the convective layer. As dissipation of orbital energy increases with the approach of merger, one might expect a disruption in the present pattern of convection. High-quality spectroscopic measurements would be a valuable supplement to the photometric light curves as they can resolve the ambiguity of which star has a spot. In addition, radio and X-ray emission might be expected from the corona surrounding the contact binary, and this could provide an independent probe of magnetic activity.

Reanalysis of the photometric and spectroscopic data with a more realistic physical model than used in this paper would improve estimates of physical properties and their uncertainties. The soon-to-be-released second version of the Phoebe analysis package⁸ (as described by Prsa et al. 2016) will be well suited for this.

Finally, more observations of the orbital timing of other contact binaries are needed. For example, if new timing measurement of the systems in the Kepler eclipsing binary catalog were to show that the observation of an increasingly negative period derivative in KIC 9832227 is unique, we might conclude that some new physical behavior is being exhibited, independent of the merging star hypothesis. Moreover, other precursor candidates should be sought in new large surveys now coming online. The key is to identify contact binaries promptly and follow up candidates with large period derivatives with intensive monitoring to ensure that precursors are identified before they merge.