The Full Kepler Phase Curve of the Eclipsing Hot White Dwarf Binary System KOI-964

In our phase-curve analysis of KOI-964, we utilize all available single aperture photometry (SAP) data for the system taken during quarters 0 through 17 (Q0–Q17). Most of the quarters contain ∼90 days of data, with the exception of Q0 (∼10 days), Q1 (∼33 days), and Q17 (∼32 days). The continuous long-cadence photometric observations with 30 minute integrated exposures were interrupted by monthly data downlinks, quarterly field rotations, as well as other occasional cessations in data collection.

The official Kepler data-processing pipeline (e.g., Jenkins et al. 2010, 2017) provides quality flags that indicate when certain spacecraft actions occurred or when nonnominal operation may have yielded unreliable photometry. We discard all exposures with a nonzero quality flag value. We then apply a 16 point-wide (∼8 hr) moving median filter to the photometric series, with transits and secondary eclipses masked, and remove 3σ outliers.

Inspection of the raw SAP light curves reveals clear periodic flux variations attributable to the astrophysical phase-curve signal, as well as long-term low-frequency brightness modulations that indicate the presence of instrumental systematics (correlated noise). We also find occasional flux ramps lasting up to several days immediately following gaps in data collection, such as those corresponding to data downlinks. Detrending these ramp-like features in addition to the other long-term systematic trends requires additional systematics parameters and often incurs correlations between the systematics and astrophysical parameters. We therefore choose to consistently trim away five days worth of data after each gap. An exception is the short first quarter (Q0), because removing five days from the time series would leave just over one orbital period worth of data, with a single primary and secondary eclipse. Without a pair of primary or secondary eclipses to self-consistently constrain the orbital period, this segment would be insufficient to reliably compute the best-fit systematics model in a simultaneous fit with the astrophysical phase-curve model (see below).

The data gaps and trimming divide each quarter's time series into discrete segments. In our initial analysis of the KOI-964 phase curve, we carry out individual fits of each segment separately in order to optimize the systematics model prior to the joint fit. Information on the full list of segments analyzed in this work is given in Table 1.

A total of 11 RV observations of KOI-964 were obtained over nine epochs using the High Resolution Echelle Spectrometer (HIRES) instrument at the Keck Observatory between UT 2014 August 2 and 2014 September 12. The first epoch consists of three consecutive exposures, and the RVs derived from those observations were combined into a single measurement. The HIRES spectra lie across three chips and span the wavelength regimes 364.3–479.5 nm, 497.7–642.1 nm, and 654.3–799 nm, respectively, resulting in 23 total echelle orders, each with a length of 4020 pixels. The first of these chips spans a wavelength range that contains many hydrogen features. The second of these chips contains the spectral features produced by the iodine wavelength calibration cell used to compute the wavelength solution.

The primary in the KOI-964 system is an A-type star, with a mass of ∼2.3 M_⊙ (Carter et al. 2011). Extracting RV measurements for such a target presents two major challenges: first, the spectral features are wide and in most cases take up a large fraction of a single HIRES echelle order, resulting in complications for continuum normalization across those orders; second, because HIRES is not an environmentally stabilized spectrograph, the wavelength solution varies even over the course of a night.

We derive relative RV values from the HIRES echelle spectra using the method of Becker et al. (2015). In summary, we utilize a simultaneous fit to all 23 HIRES orders using a two-dimensional model of the continuum level, which allows featureless orders to act as anchors in the fit. Eight of the orders have discernible spectral features, but the depth and width of these features prevent the determination of the true continuum level. We use the remaining featureless orders to constrain the overall continuum model. When a wavelength solution is not available for a given observation, we extrapolate from an existing solution taken closest in time to the target observation. These solutions can be generated when the iodine cell is used during an observation. The specifics of this extrapolation technique require extrapolating both between chips and between observations, and are described in detail by Becker et al. (2015). We use a Markov Chain Monte Carlo (MCMC) routine to optimize the fit between our template spectrum (i.e., a smoothed version of the first acquired spectrum) and each of the subsequent spectra, simultaneously fitting the blaze function and the relative redshift. Unlike in Becker et al. (2015), we do not fit for the stellar v sin i due to the small number of data epochs and highly broadened spectra, which provide insufficient information content to fit the extra variable. The fit to the full data cube produces relative RVs, where the RV of each observation is measured relative to the smoothed template spectrum. We finally apply a barycentric correction to the extracted RVs. The small number of features on the broadened spectra lead to relatively poor RV precision. The nine epochs of RVs obtained using this method are presented in Table 2.

We model both transits (when the white dwarf passes in front of the larger primary) and secondary eclipses (when the white dwarf is occulted) using the BATMAN package (Kreidberg 2015). Throughout this work, we use the subscripts a and b to refer to the host star and the orbiting white dwarf, respectively. Under the parameterization scheme adopted by ExoTEP, the shape and timing of these eclipse light curves λ(t) are determined by the radius ratio R_b/R_a, the relative brightness of the secondary f_b, the impact parameter b, the scaled orbital semimajor axis a/R_a, a reference midtransit time T₀, the orbital period P, the orbital eccentricity e, and the argument of periastron ω. The midtransit time is determined for the zeroth epoch, which in the case of the joint fit is designated to be the transit event closest to the mean of the combined time series.

When generating the model eclipse light curves λ(t) at the 30 minute cadence of the Kepler time series, we supersample the transit and secondary eclipse light curves in the vicinity of each event at 30 s intervals and average the finely sampled fluxes in order to accurately compute the integrated flux at each exposure. We also account for the nonnegligible relative brightness of the white dwarf and dilute the primary transits accordingly by the factor (1 + f_b).

In the joint fit presented in this work, we allow all of the transit and orbital parameters to vary freely. On the other hand, when initially fitting for the individual segment light curves to optimize the systematics model (Section 3.3), we fix b and a/R_a to the best-fit values from Carter et al. (2011) and e and ω to zero, because the short time baseline of each photometric series does not provide any strong constraints on orbital geometry, and the orbit of the system is consistent with circular (see Section 4).

For the transits, we employ a standard quadratic limb-darkening law to describe the radial brightness profile of the primary. In the joint fit, both coefficients u₁ and u₂ are allowed to vary, while in the individual segment fits, we fix them to the values from Carter et al. (2011): u₁ = 0.20 and u₂ = 0.2964. For the white dwarf secondary, Carter et al. (2011) fixed the quadratic limb-darkening coefficient to zero and found that the remaining linear coefficient was largely unconstrained and consistent with zero. In our analysis, we fix both coefficients to zero and do not report any significant improvement to the quality of the joint fit when allowing those coefficients to vary.

Following Carter et al. (2011), we model the out-of-occultation brightness variation as a third-order harmonic series in phase (e.g., Carter et al. 2011),

$$\begin{eqnarray}&&\psi (t)=1+\displaystyle \sum _{k=1}^{3}{A}_{k}\sin (k\phi (t))+\displaystyle \sum _{k=1}^{3}{B}_{k}\cos (k\phi (t)),\end{eqnarray} \tag{ 1 }$$

where ϕ(t) = 2π(t − T₀)/P, and we have normalized the flux such that the combined average brightness of both binary components is unity.

To model the instrumental systematics in each segment, we use a generalized polynomial in time:

$$\begin{eqnarray}&&{S}_{n}^{\{i\}}(t)=\displaystyle \sum _{k=0}^{n}{c}_{k}^{\{i\}}{\left(t-{t}_{0}\right)}^{k}.\end{eqnarray} \tag{ 2 }$$

Here, t₀ is the first time stamp in the light curve from segment i, and n is the order of the polynomial model. The full phase-curve model is given by

$$\begin{eqnarray}&&F(t)={S}_{n}^{\{i\}}(t)\times \lambda (t)\times \psi (t).\end{eqnarray} \tag{ 3 }$$

To determine the optimal polynomial order for each segment, we carry out individual fits of each segment to the model in Equation (3) and choose the polynomial order that minimizes the Bayesian information criterion (BIC): $\mathrm{BIC}\equiv k\mathrm{log}N-2\mathrm{log}L$, where k is the number of free parameters in the fit, N is the number of data points, and L is the maximum log likelihood. The optimal polynomial orders determined from our individual segment fits are listed in Table 1. In all cases, when selecting polynomials of similar order, the best-fit astrophysical parameters do not vary by more than 0.3σ. A compilation plot of the systematics-corrected light curves from all segments is provided in Figures 5 and 6 in the Appendix.

The residuals from the light-curve fit for the first segment of Q15 (labeled as 15-0) show severe uncorrected systematics in the form of a quasiperiodic modulation with a characteristic frequency about 10% higher than the orbital frequency. Close inspection of the binned residuals from the individual segment fits reveal several more instances of similar uncorrected systematics. The characteristic frequencies of these residual signals are close to one another, but not identical; the occurrences of these systematics are largely confined to individual segments, with adjacent segments showing clean, featureless residuals. These observations strongly rule out an astrophysical source of these additional photometric modulations, and they are likely to be instrumental systematics. In addition to 15-0, we find significant uncorrected systematics in the segments 3-0, 3-1, 5-0, 7-0, 7-1, 11-0, 11-1, 12-1, 16-0, and 16-1. We trim these segments in the joint light-curve fit presented in this work. However, we find that including them in the photometric series does not affect the measured astrophysical parameters.

When carrying out the joint phase-curve fit of the full Kepler photometric time series, we do not combine the uncorrected light curves and simultaneously fit the systematics model for all segments, because that would include over 200 systematics parameters and incur forbiddingly large computational overheads. Instead, we first remove 4σ outliers from the best-fit model for each segment and divide the photometric series by the best-fit systematics model from the individual analysis. Then, we concatenate the detrended flux arrays to form the combined light curve for our joint analysis. To empirically validate this approach, we have experimented with carrying out joint analyses of select subsets of data (e.g., Q2 and Q3 only) and comparing the results from (a) fits where we compute the full systematics model for all component segments and (b) fits where we use pre-detrended light curves and no systematics modeling. In all the cases we have tested, the astrophysical parameter estimates agree to within 0.2σ. The combined light curve for our joint fit contains 34,293 data points.

In addition to the astrophysical parameters, we fit for a uniform per-point uncertainty σ to ensure that the resultant reduced chi-squared value is near unity and self-consistently derive realistic uncertainties on the astrophysical parameters, given the intrinsic scatter of the light curves. The total number of free parameters in our joint phase-curve analysis is 17.

The combined light curve spans 4 yr and contains (either partially or in full) 229 transits and 237 secondary eclipses. The ExoTEP pipeline utilizes the affine-invariant MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013) to simultaneously compute the posterior distributions of all free parameters. The number of walkers is equal to four times the number of free parameters, and each chain has a length of 30,000 steps. After inspecting the plotted chains by eye, we only use the last 40% of each chain that is extracted to generate the posterior distributions. We also apply the Gelman–Rubin convergence test (Gelman & Rubin 1992) to ensure that the diagnostic value $\hat{R}$ is well below 1.1.

The results of our joint analysis of all 18 Kepler quarters worth of data are listed in Table 3. We report the median and 1σ uncertainties for the astrophysical and noise parameters, as derived from the marginalized one-dimensional posterior distributions. The combined phase-folded light curve is plotted in Figure 1, along with the best-fit full phase-curve model.

Zoom In Zoom Out Reset image size

Table 3.  Results of Joint Phase-curve Analysis

Parameter	Value	Error
Fitted Parameters
Rb/Ra	0.080983	0.000077
fb	0.0189365	${}_{-0.0000041}^{+0.0000044}$
T0 (BJDTDB)	2455662.252189	${}_{-0.000022}^{+0.000024}$
P (days)	3.273698741	${}_{-0.000000053}^{+0.000000054}$
b	0.6740	${}_{-0.0015}^{+0.0017}$
a/Ra	7.052	${}_{-0.015}^{+0.014}$
e	0.000271	${}_{-0.000029}^{+0.000153}$
ω (◦)	−19.3	${}_{-36.1}^{+36.2}$
u1	0.176	${}_{-0.031}^{+0.029}$
u2	0.312	${}_{-0.034}^{+0.035}$
A1 (ppm)	99.21	${}_{-0.86}^{+0.95}$
A2 (ppm)	−39.75	${}_{-0.98}^{+0.88}$
A3 (ppm)	0.34	${}_{-0.90}^{+0.87}$
B1 (ppm)	268.4	1.1
B2 (ppm)	−568.2	1.0
B3 (ppm)	−15.0	${}_{-1.1}^{+1.0}$
σ (ppm)	117.74	${}_{-0.47}^{+0.44}$
Derived Parameters
Transit depth (ppm)a	6558	13
i (◦)	84.515	${}_{-0.026}^{+0.023}$
$e\cos \omega$	0.000239	${}_{-0.000012}^{+0.000011}$
$e\sin \omega$	−0.00008	${}_{-0.00026}^{+0.00016}$

Note.

^aCalculated as ${({R}_{b}/{R}_{a})}^{2}$.

Download table as:  ASCIITypeset image

When comparing with the results of Carter et al. (2011), we find that our astrophysical parameter estimates are at least 3.5 times more precise. Our inclination i and a/R_a estimates agree with the previously published values at the 0.4σ level. The Carter et al. (2011) analysis fixed the quadratic coefficient for the primary's limb-darkening law to the values calculated by Sing (2010) for a star with the stellar parameters of the primary in the KOI-964 system—u₂ = 0.2964—while fitting for the linear coefficient and obtaining u₁ = 0.20 ± 0.02. Our fitted values for u₁ and u₂ agree with these estimates at the 0.7σ and 0.5σ levels, respectively. The per-point uncertainty of 118 ppm that we obtain from our joint fit is about 10% smaller than the value reported by Carter et al. (2011), indicating that data in more recent quarters have somewhat less scatter than the Q0 and Q1 photometry.

The period estimate from our joint analysis has an exquisite precision of 5 ms and differs from the value in Carter et al. (2011) by 1.8σ. The radius ratio R_b/R_a is consistent with the previously published value at 0.82σ, while our secondary eclipse depth estimate f_b is smaller by 3.1σ. We have experimented with fitting only the data analyzed by Carter et al. (2011)—Q0 and Q1—and obtain P and f_b values that are consistent with their results at the 0.1σ and 0.6σ levels, respectively. This demonstrates that the addition of 16 more quarters of data has shifted the global estimate of the primary–secondary flux ratio in particular to a significantly different value.

The estimates we obtain for the main parameters of interest—the phase-curve harmonic amplitudes—are consistent with the Carter et al. (2011) values at better than the 1.1σ level. We detect the amplitudes of five harmonic terms at significance levels above 15σ; they are, in order of decreasing amplitude and statistical significance, $\cos 2\phi$, $\cos \phi$, $\sin \phi$, $\sin 2\phi$, and $\cos 3\phi$. The remaining harmonic, $\sin 3\phi$, has an amplitude that is consistent with zero at the 0.4σ level. The astrophysical interpretation of these phase-curve components will be discussed in detail in the following section.

One major discrepancy between our analysis and the results in Carter et al. (2011) is the orbital eccentricity. They report an $e\cos \omega$ value of 0.0029 ± 0.0005, which translates to a significant delay in the midpoint of the secondary eclipse relative to the expectation for a circular orbit: ${\rm{\Delta }}t\equiv 2{Pe}\cos \omega /\pi \,\sim$ 8.7 minutes (note: the value for the time delay reported in their work is erroneously off by a factor of 2). In our joint phase-curve analysis, we obtain $e\cos \omega ={0.000239}_{-0.000012}^{+0.000011}$, more than an order of magnitude smaller than the estimate in Carter et al. (2011). In a separate joint analysis of just Q0 and Q1 data, we find a similar value to the value from our full global fit, leading us to speculate that the reported $e\cos \omega$ value in Carter et al. (2011) may be missing a zero after the decimal point. Our smaller updated constraint translates to a predicted secondary eclipse time delay of 43.0 ± 2.1 s. Assuming the orbital semimajor axis we derive from our stellar isochrone analysis and the fitted a/R_a value (Section 4.4)—0.0620 ± 0.0044 au—the relative light travel time delay between superior and inferior conjunctions is 61.9 ± 4.4 s, which is close (<4σ) to the secondary eclipse time delay predicted by our joint phase-curve analysis. Therefore, we conclude that the orbit of the white dwarf is consistent with circular.

Figure 2 shows the Lomb–Scargle periodogram of the residuals from our best-fit phase-curve model. No significant peaks at integer multiples of the orbital frequency are visible, demonstrating that our phase-curve modeling, which includes Fourier terms up to third order, accounts for the full harmonic content in the astrophysical phase curve phased at the orbital period.

Zoom In Zoom Out Reset image size

There is a cluster of peaks in the power spectrum spanning frequencies somewhat higher than the orbital frequency. These signals are attributable to residual systematics that are not fully removed in the data trimming (see Section 3). When we include all segments in the joint fit, some peaks in this region of frequency space rise to 8σ–9σ significance. With the trimming employed in the fit presented in this work, the noise contribution from these systematics is greatly reduced. Furthermore, the periodicities contributed by the systematics are distinct from the orbital frequency and all associated harmonics, and as such, their presence does not affect the astrophysical phase-curve parameters we measure in the joint fit.

A very strong single peak at a frequency of ∼5.28 period⁻¹ is evident in the periodogram, corresponding to a period of around 0.62 days. We have inspected all available Kepler light curves for stars within ∼1' of KOI-964 and do not find any periodicities at this frequency. Furthermore, the estimated relative contamination within the optimal extraction aperture for KOI-964 provided by the official Kepler data-processing pipeline (given by the CROWDSAP keyword in the header) is 0%. We conclude that this signal is most likely from the KOI-964 system.

To further characterize this signal, we fit a simple sinusoidal function to the unbinned, non-phase-folded residual series:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\beta \cos \left[2\pi \displaystyle \frac{t-{T}_{0}-\tau }{{\rm{\Pi }}}\right].\end{eqnarray} \tag{ 4 }$$

Here, β and Π are the semi-amplitude and period of the periodic signal, the zero point of the time series is set to the median midtransit time from the joint fit T₀ = 2455662.252189, and τ represents the relative phase shift of this modulation. Our MCMC analysis yields Π = 0.620276 ± 0.000011 days, β = 12.1 ± 0.9 ppm, and $\tau ={0.0453}_{-0.0074}^{+0.0079}$ days.

The frequency of this observed periodic signal lies in the range spanned by γ Dor pulsators (e.g., Guzik et al. 2000; Balona et al. 2011). Hundreds of γ Dor pulsators have been discovered in the Kepler era, with typical pulsation frequencies smaller than 4 day⁻¹ and peaked around 1 day⁻¹ (e.g., Tkachenko et al. 2013; Bradley et al. 2015); some of the measured pulsation amplitudes from the Kepler sample are on the order of 10 ppm, consistent with the measured amplitude for the KOI-964 signal. However, γ Dor pulsators tend to be late A-/F-type stars, with effective temperatures in the range 6000–8000 K (e.g., Tkachenko et al. 2013; Bradley et al. 2015), significantly cooler than KOI-964 (T_eff ∼ 10,000 K). Therefore, the origin of the observed pulsation signal on KOI-964 remains uncertain.

We use the RadVel package (Fulton et al. 2018) to model the orbit of the hot white dwarf based on the nine Keck/HIRES RV measurements listed in Table 2. We place Gaussian priors on the orbital period (P') and time of inferior conjunction (i.e., midtransit time; ${T}_{0}^{{\prime} }$) based on the median values and uncertainties from our joint phase-curve fit (Table 3). Here, the prime symbols (') are used to distinguish parameters included in the RV analysis from the analogous parameters calculated in our joint phase-curve fitting. We fit for the RV semi-amplitude K_RV, along with the mean RV offset (γ) and jitter. The linear and second-order acceleration terms ($\dot{\gamma }$, $\ddot{\gamma }$) are fixed to zero.

We run a fit that allows for independent constraints on orbital eccentricity ($e^{\prime}$) and argument of periastron ($\omega ^{\prime}$), as well as a fit with eccentricity fixed to zero. The free-eccentricity RV fit yields $e^{\prime} ={0.082}_{-0.052}^{+0.079}$, which is consistent with zero at the 1.6σ level and in line with the nondetection of significant secondary eclipse timing offset in our joint phase-curve fit. When fitting an RV model with orbital eccentricity fixed to zero, we obtain a similar BIC value to the free-eccentricity fit, while the Akaike information criterion with a small sample size correction (AICc) strongly favors the zero-eccentricity model (ΔAICc = 90.83). Therefore, we present the RV fit results assuming a circular orbit in this work.

The results of our RV fitting are shown in Table 4 and Figure 3. The RV variation, with a fitted semi-amplitude of ${K}_{\mathrm{RV}}={18.5}_{-1.8}^{+2.0}$ km s⁻¹, is detected at the 10.3σ level.

Zoom In Zoom Out Reset image size

Table 4.  Results of RV Analysis

Parameter	Value	Error
$P^{\prime}$ (days)a	3.273698811	${}_{-0.000000048}^{+0.000000049}$
${T}_{0}^{{\prime} }$ (BJDTDB)a	2455665.525884	0.000022
$e^{\prime}$	$\equiv$ 0.0	⋯
${K}_{\mathrm{RV}}$ (km s−1)	18.5	${}_{-1.8}^{+2.0}$
γ (km s−1)b	−1.3	1.3
jitter (km s−1)	3.2	${}_{-1.0}^{+1.5}$

Notes.

^aConstrained by Gaussian priors derived from joint phase-curve fit estimates (Table 3). ^bLinear and second-order RV acceleration terms fixed to zero.

Download table as:  ASCIITypeset image

In order to calculate the mass, radius, and orbital semimajor axis of the transiting white dwarf from the measured RV variation and fitted orbital parameters, we must first obtain a reliable estimate of the host star's mass and radius.

We follow a methodology similar to the one described in Berger et al. (2018) and utilize the stellar classification package isoclassify (Huber et al. 2017). For input, we combine the parallax measurement from Gaia Data Release 2 (0.4257 ± 0.0348 mas; Lindegren et al. 2018) with the 2MASS K-band magnitude (13.056 ± 0.029 mag) and SDSS g-band magnitude (13.000 ± 0.020 mag). We correct the Gaia parallax with a positive offset of 0.03 mas (Berger et al. 2018; Lindegren et al. 2018).

The isoclassify routine combines the distance calculated from the input parallax with the extinction derived from a 3D dust reddening map and interpolated reddening vectors listed in Green et al. (2018) to estimate the absolute magnitude of the host star. Using the "grid method," as opposed to the "direct method" utilized in Berger et al. (2018), then allows for the posteriors of the host star's properties to be self-consistently computed through comparisons with interpolated stellar isochrones from Modules for Experiments in Stellar Astrophysics Isochrones & Stellar Tracks (MIST) grids (Choi et al. 2016).

The presence of the binary companion can bias the retrieved parameters of the host star and lead to an overestimation in the host star's radius and temperature. In the Kepler bandpass, the white dwarf contributes roughly 2% of the total luminosity in the system. As a first-order correction for the contamination, we first use isoclassify on the uncorrected magnitudes to derive the host star's temperature and radius. We then approximate the spectra of the host star and white dwarf companion as blackbodies and utilize the measured flux ratio f_b in the Kepler bandpass and the measured radius ratio R_b/R_a from the joint light-curve fit to compute the white dwarf's effective temperature (Equation (7)). Lastly, while fixing the component radii and the computed white dwarf flux in the Kepler bandpass, we adjust the host star's temperature to match the flux ratio f_b. By calculating the difference in the host star's integrated flux in the g and K bands between the initial and adjusted blackbody spectra, we obtain magnitude corrections of Δg = +0.026 mag and ΔK = + 0.013 mag.

Table 5 lists the estimates of the host star's properties we obtain from isoclassify, using the corrected magnitudes. All of these values are consistent to within 1σ with the corresponding values derived from using the uncorrected magnitudes. Our value for R_a is in good agreement with the values computed by Berger et al. (2018) using the "direct method," while the stellar mass M_a is consistent with the value in Mathur et al. (2017)—${2.760}_{-0.674}^{+0.289}$ M_⊙—which was derived without including the Gaia parallax.

Table 5.  KOI-964 System Properties

Parameter	Value	Error
(1) A-type host stara
Ra (R⊙)	1.89	${}_{-0.13}^{+0.14}$
Ma (M⊙)	2.23	0.12
${\rho }_{a}$ (g/cc)b	0.5597	${}_{-0.0069}^{+0.0071}$
Teff,a (K)	9940	${}_{-230}^{+260}$
$\mathrm{log}g$	4.226	${}_{-0.054}^{+0.046}$
$[\mathrm{Fe}/{\rm{H}}]$	−0.09	${}_{-0.14}^{+0.15}$
Luminosity, L1 (L⊙)	31.7	${}_{-4.5}^{+5.8}$
Distance, d (pc)	2220	${}_{-140}^{+170}$
Age (Gyr)	0.21	${}_{-0.08}^{+0.11}$
(2) Transiting white dwarf companionb
Rb (R⊙)	0.153	0.011
q ≡ Mb/Ma	0.106	0.012
Mb (M⊙)	0.236	${}_{-0.027}^{+0.028}$
ρb (g/cc)	93	23
Teff,b (K)	15080	400
a (au)	0.0620	0.0044

Notes.

^aComputed with isoclassify using Gaia parallax, 2MASS K-band photometry, and the SDSS g-band magnitude. ^bDerived from R_a, M_a, and the results of our RV and joint phase-curve analyses.

Download table as:  ASCIITypeset image

Having estimated the host star's radius R_a, we use the best-fit values for R_b/R_a and a/R_a from the joint phase-curve analysis (Table 3) to calculate the radius and orbital semimajor axis of the white dwarf companion: R_b = 0.153 ± 0.011 R_⊙ and a = 0.0620 ± 0.0044 au.

We separately calculate the mass ratio q ≡ M_b/M_a and the white dwarf's mass M_b using the measured RV semi-amplitude K_RV and the primary's mass derived above. In the case of a circular orbit, the quantities q and M_b are related to K_RV by the following expression:

$$\begin{eqnarray}\begin{array}{rcl}{K}_{\mathrm{RV}} & = & q\sin i{\left[\displaystyle \frac{2\pi G}{P}\displaystyle \frac{{M}_{a}}{{\left(1+q\right)}^{2}}\right]}^{1/3}\\ & = & \ {M}_{b}\sin i{\left[\displaystyle \frac{2\pi G}{P}\displaystyle \frac{1}{{\left({M}_{a}+{M}_{b}\right)}^{2}}\right]}^{1/3}.\end{array}\end{eqnarray} \tag{ 5 }$$

To construct the posteriors of q and M_b, we randomly sample combinations of (P, i, M_a, K_RV) from the respective posteriors and numerically solve for q and M_b. The resultant estimates are q = 0.106 ± 0.012 and ${M}_{b}={0.236}_{-0.027}^{+0.028}$ M_⊙. The deduced properties of the white dwarf companion are summarized in Table 5.

The density of the host star ρ_a can be calculated in two ways. The first method is direct, ${\rho }_{a}={M}_{a}/(4\pi {R}_{a}^{3}/3)$, while the second method utilizes Kepler's third law and the mass ratio q to infer ρ_a from the orbital period and semimajor axis:

$$\begin{eqnarray}&&{\rho }_{a}=\displaystyle \frac{3\pi }{{{GP}}^{2}}\displaystyle \frac{1}{1+q}{\left(\displaystyle \frac{a}{{R}_{a}}\right)}^{3}.\end{eqnarray} \tag{ 6 }$$

Using these two methods, we obtain the statistically consistent estimates ${0.450}_{-0.079}^{+0.082}$ g/cc and ${0.5597}_{-0.0069}^{+0.0071}$ g/cc. The second method relies on quantities with significantly smaller relative uncertainties and therefore yields the more precise density estimate, which we take and list in Table 5. For the white dwarf's density, we can only utilize the direct method.

Lastly, we estimate the temperature of the white dwarf companion, T_eff,b. The secondary eclipse depth gives the relative fluxes of the two binary components and is related to their emission spectra F_ν,a and F_ν,b by

$$\begin{eqnarray}&&{f}_{b}={\left(\displaystyle \frac{{R}_{b}}{{R}_{a}}\right)}^{2}\displaystyle \frac{\langle {F}_{\nu ,b}\rangle }{\langle {F}_{\nu ,a}\rangle },\end{eqnarray} \tag{ 7 }$$

where the spectra are averaged over the Kepler bandpass, weighted at each frequency by the corresponding value in the Kepler response function. We model the host star's emission using PHOENIX stellar spectra (Husser et al. 2013), while the white dwarf's flux is represented as a simple blackbody.

To facilitate estimating T_eff,b and reliably propagating the uncertainties in stellar properties to the uncertainty in T_eff,b, we derive empirical analytic expressions for $\langle {F}_{\nu ,a}\rangle$ and $\langle {F}_{\nu ,b}\rangle$ as functions of (T_eff,a, [Fe/H], $\mathrm{log}g$) and T_eff,b, respectively. Using all available PHOENIX model stellar spectra spanning the ranges T_eff,a = [8000, 12,000] K, [Fe/H] = [−1.0, 0.5], and $\mathrm{log}g=[3.5,5.0]$, we compute $\langle {F}_{\nu ,a}\rangle$ at each point in the three-dimensional grid and fit a generalized linear polynomial in the dependent variables. Similarly, we fit a cubic polynomial in T_eff,b to the array of $\langle {F}_{\nu ,b}\rangle$ values calculated for blackbody spectra across the temperature range T_eff,b = [10,000, 20,000] K at 50 K intervals. We obtain the following relationships (units of $\langle {F}_{\nu }\rangle$ are ${10}^{15}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}$):

$$\begin{eqnarray}\begin{array}{rcl}\langle {F}_{\nu ,a}\rangle & = & 4.709+10.525{\tau }_{a}\\ & & \ +\ 0.284[\mathrm{Fe}/{\rm{H}}]-0.0081(\mathrm{log}g-4.0),\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\langle {F}_{\nu ,b}\rangle =4.856+13.934{\tau }_{b}+6.261{\tau }_{b}^{2}-1.864{\tau }_{b}^{3},\end{eqnarray} \tag{ 9 }$$

where ${\tau }_{i}\equiv \left(\tfrac{{T}_{\mathrm{eff},i}}{{\rm{10,000}}\,{\rm{K}}}-1\right)$. Following a similar Monte Carlo sampling method used previously to compute the white dwarf's mass, we randomly sample from the posteriors of (f_b, R_b/R_a, T_eff,a, [Fe/H], $\mathrm{log}g$) and numerically compute the corresponding values of T_eff,b to obtain T_eff,b = 15,080 ± 400 K.

The large radius and high temperature of the white dwarf companion indicate that the object must be young and still cooling. When compared to the radius of a degenerate He star of the same mass, our measured value is roughly eight times larger. As mentioned in Carter et al. (2011), cooling models of He white dwarfs show that lower-mass objects with relatively H-rich atmospheres can remain bloated and hot for longer than their more massive, H-poor counterparts (Hansen & Phinney 1998; Nelson et al. 2004). Our measured mass ${M}_{b}={0.236}_{-0.027}^{+0.028}$ M_⊙ is consistent with the smaller of the two discrepant white dwarf masses that Carter et al. (2011) derived from their phase-curve analysis. They indicated that a H-rich 0.21 M_⊙ white dwarf can remain large and hot for ∼0.15 Gyr. The young host star age we deduce from the isochrone analysis (${0.21}_{-0.08}^{+0.11}$ Gyr) is therefore consistent with the predicted evolutionary track of the hot white dwarf companion.

As the two components of the binary system orbit around their center of mass, the apparent spectral intensity of both objects at a given wavelength modulates due to the periodic red- and blueshifting of the spectra, as well as variations in the photon emission rate and light aberration (e.g., Shakura & Postnov 1987; Loeb & Gaudi 2003; Zucker et al. 2007; Shporer et al. 2010). The composite effect is commonly referred to as Doppler boosting. This temporal modulation is driven by the projected RVs of the two components, v_r,a and v_r,b, and as such, the harmonic contribution of Doppler boosting to the phase curve is carried by the sine function of the orbital phase. The amplitude ${A}_{1}^{\mathrm{DB}}$ of the Doppler boosting is related to the system properties through the following expression (e.g., Loeb & Gaudi 2003; Shporer et al. 2010; Carter et al. 2011):

$$\begin{eqnarray}\begin{array}{rcl}{A}_{1}^{\mathrm{DB}} & = & {\alpha }_{a}\left(\displaystyle \frac{{v}_{r,a}}{c}\right)+{\alpha }_{b}{f}_{b}\left(\displaystyle \frac{{v}_{r,b}}{c}\right)\\ & = & \ \left(\displaystyle \frac{2\pi a}{{Pc}}\right)\left[{\alpha }_{a}{\left(1+\displaystyle \frac{1}{q}\right)}^{-1}-{\alpha }_{b}{\left(1+q\right)}^{-1}{f}_{b}\right].\end{array}\end{eqnarray} \tag{ 10 }$$

Here, q is the mass ratio as defined previously, f_b is the secondary eclipse depth, and c is the speed of light.

The prefactors α_i, in short, reflect the relative change in the objects' integrated fluxes through the observed bandpass due to Doppler shifting and depend on the shape of the emission spectra:

$$\begin{eqnarray}&&{\alpha }_{i}=3-\left\langle \displaystyle \frac{d\mathrm{log}{F}_{\nu ,i}}{d\mathrm{log}\nu }\right\rangle .\end{eqnarray} \tag{ 11 }$$

F_ν is the object's emission spectrum as a function of frequency ν, and the derivative is averaged over the Kepler bandpass.

We empirically model the dependence of α_a and α_b on stellar temperature, metallicity, and surface gravity using the same method described in Section 4.4. In the case of α_a, we find that including quadratic and cubic terms in temperature as well as a correlation term between the temperature and surface gravity greatly improves the fit:

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{a} & = & 2.077-1.543{\tau }_{a}+4.478{\tau }_{a}^{2}-8.790{\tau }_{a}^{3}-0.021[\mathrm{Fe}/{\rm{H}}]\\ & & \ +\ 0.043(\mathrm{log}g-4.0)-0.311{\tau }_{a}(\mathrm{log}g-4.0),\end{array}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\alpha }_{b}=2.635-1.852{\tau }_{b}+1.433{\tau }_{b}^{2}-0.502{\tau }_{b}^{3},\end{eqnarray} \tag{ 13 }$$

where, as before, ${\tau }_{i}\equiv \left(\tfrac{{T}_{\mathrm{eff},i}}{{\rm{10,000}}\,{\rm{K}}}-1\right)$.

We calculate the predicted harmonic amplitude due to Doppler boosting using Equations (10), (12), and (13) by sampling the posteriors for the dependent variables—P, f_b, a, q, T_eff,a, [Fe/H], log g, and T_eff,b (Tables 3 and 5). The resultant estimate is ${A}_{1}^{\mathrm{DB}}={114}_{-16}^{+18}$ ppm, which is consistent with our best-fit phase-curve amplitude (${A}_{1}={99.21}_{-0.86}^{+0.95}$ ppm) at the 0.9σ level.

Among the three main processes that contribute to the out-of-eclipse phase-curve modulation, Doppler boosting is solely responsible for the fundamental harmonic of the sine term. The consistency between the predicted value of A₁ computed using the RV-derived mass ratio q and the observed photometric amplitude indicates that the formalism described above adequately captures the physical mechanisms that drive the Doppler-boosting signal. We can then use the very precise measurement of A₁ from our joint phase-curve fit and Equation (12) to derive an independent estimate of the mass ratio, q*. Following a sampling method similar to that before, we obtain ${q}^{* }={0.0932}_{-0.0058}^{+0.0068}$, which is consistent with the value derived from our RV analysis (Table 5) at the 0.9σ level and twice as precise.

In a binary system of two self-luminous objects, the radiation emitted by one component is incident on the other and is subsequently scattered or absorbed and re-emitted. The regions near the substellar points on the mutually facing hemispheres are therefore expected to be brighter than the other regions. Over the course of an orbit, the viewing phase of the secondary and the position of the illuminated region on the primary both change, imparting a periodic brightness variation to the phase curve. The maximum illumination of the primary is observed midtransit, while the primary-facing hemisphere of the secondary is fully oriented toward the observer mid-eclipse. Both of these variations contribute primarily to the fundamental cosine mode in the phase-curve harmonic series, albeit with opposite signs.

Under the assumption of radiative equilibrium, i.e., when all incident radiation is re-emitted at the effective temperature of the illuminated object, the amplitudes of the mutual illumination modulation are given by Kopal (1959)

$$\begin{eqnarray}\begin{array}{rcl}{B}_{1}^{\mathrm{ILL}} & = & \displaystyle \frac{17}{16}{\left(\displaystyle \frac{{R}_{a}}{a}\right)}^{2}\left[\displaystyle \frac{1}{4}\displaystyle \frac{{R}_{a}}{a}+\displaystyle \frac{1}{3}\right]{\left(\displaystyle \frac{{T}_{\mathrm{eff},b}}{{T}_{\mathrm{eff},a}}\right)}^{4}{\left(\displaystyle \frac{{R}_{b}}{{R}_{a}}\right)}^{2}\\ & & -\ \displaystyle \frac{17}{16}{\left(\displaystyle \frac{{R}_{b}}{a}\right)}^{2}\left[\displaystyle \frac{1}{4}\displaystyle \frac{{R}_{b}}{a}+\displaystyle \frac{1}{3}\right]\displaystyle \frac{{\beta }_{b}}{{\beta }_{a}},\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}\begin{array}{rcl}{B}_{2}^{\mathrm{ILL}} & = & \displaystyle \frac{17}{16}{\left(\displaystyle \frac{{R}_{a}}{a}\right)}^{2}\left[\displaystyle \frac{3}{16}\displaystyle \frac{{R}_{a}}{a}+\displaystyle \frac{16}{27{\pi }^{2}}\right]{\left(\displaystyle \frac{{T}_{\mathrm{eff},b}}{{T}_{\mathrm{eff},a}}\right)}^{4}{\left(\displaystyle \frac{{R}_{b}}{{R}_{a}}\right)}^{2}\\ & & -\ \displaystyle \frac{17}{16}{\left(\displaystyle \frac{{R}_{b}}{a}\right)}^{2}\left[\displaystyle \frac{3}{16}\displaystyle \frac{{R}_{b}}{a}+\displaystyle \frac{16}{27{\pi }^{2}}\right]\displaystyle \frac{{\beta }_{b}}{{\beta }_{a}},\end{array}\end{eqnarray} \tag{ 15 }$$

where the bolometric correction to optical wavelengths β_i is approximated by Allen (1964),

$$\begin{eqnarray}&&\mathrm{log}{\beta }_{i}\simeq 17.0-4\mathrm{log}{T}_{\mathrm{eff},i}-\displaystyle \frac{{\rm{11,600}}\,{\rm{K}}}{{T}_{\mathrm{eff},i}}.\end{eqnarray} \tag{ 16 }$$

With the measured and derived values for a/R_a, R_b/R_a, T_eff,a, and T_eff,b as input, Equations (14)–(16) predict the following amplitudes for the mutual illumination variation: ${B}_{1}^{\mathrm{ILL}}={261}_{-39}^{+44}$ ppm and ${B}_{2}^{\mathrm{ILL}}={62.0}_{-9.0}^{+10.3}$ ppm.

The gravitational interaction between the primary and the secondary produces prolate deviations from sphericity on both components, with the long axis of the resultant ellipsoidal shape lying very nearly along the line connecting the two components. The ellipsoidal distortion incurs variations in the sky-projected areas of both components as a function of orbital phase, yielding modulations in the apparent flux with maxima occurring during quadrature, i.e., a variation at the first harmonic of the cosine, $\cos 2\phi$.

The detailed physical formalism of the ellipsoidal modulation was derived in Kopal (1959) and can be described as a series of cosines, with the three leading-order terms having the following amplitudes:

$$\begin{eqnarray}&&{B}_{1}^{\mathrm{ELP}}=-{Z}^{(3)}q{\left(\displaystyle \frac{{R}_{a}}{a}\right)}^{4}(15{\sin }^{3}i-12\sin i),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{B}_{2}^{\mathrm{ELP}}=-{Z}^{(2)}q{\left(\displaystyle \frac{{R}_{a}}{a}\right)}^{3}{\sin }^{2}i,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{B}_{3}^{\mathrm{ELP}}=-5{Z}^{(3)}q{\left(\displaystyle \frac{{R}_{a}}{a}\right)}^{4}{\sin }^{3}i,\end{eqnarray} \tag{ 19 }$$

where, in line with the notation of Morris (1985) and Morris & Naftilan (1993),

$$\begin{eqnarray}\begin{array}{rcl}{Z}^{(2)} & = & \displaystyle \frac{3}{4(3-{u}_{1})}\left[\displaystyle \frac{1}{5}(15+{u}_{1})(1+{\gamma }_{1})\right.\\ & & \ +\ \left.\displaystyle \frac{5}{16}(1-{u}_{1})(3+{\gamma }_{1})(6-7{\sin }^{2}i){\left(\displaystyle \frac{{R}_{a}}{a}\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{Z}^{(3)}=\displaystyle \frac{5{u}_{1}}{32(3-{u}_{1})}(2+{\gamma }_{1}).\end{eqnarray} \tag{ 21 }$$

Here, u₁ and γ₁ are the limb-darkening and gravity-darkening coefficients, respectively, assuming a linear law.

There is a straightforward relationship between γ₁ and T_eff,a according to von Zeipel's law for early-type stars (Morris 1985):

$$\begin{eqnarray}&&{\gamma }_{1}=0.25\left\langle \displaystyle \frac{C/\lambda {T}_{\mathrm{eff},a}}{1-\exp (-C/\lambda {T}_{\mathrm{eff},a})}\right\rangle ,\end{eqnarray} \tag{ 22 }$$

where C = 1.43879 cm K and λ is wavelength. For the linear limb-darkening coefficient, we fit an analytic function to the values computed by Sing (2010) in the Kepler bandpass for the same grid of host star parameter values we used previously in Section 5.1. Likewise, we calculate γ₁ for a range of host star temperatures and fit a polynomial through the resultant array of values. We obtain the following empirical relationships:

$$\begin{eqnarray}&&{\gamma }_{1}=0.660-0.501{\tau }_{a}+0.590{\tau }_{a}^{2}-0.634{\tau }_{a}^{3},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\begin{array}{rcl}{u}_{1} & = & 0.445-0.325{\tau }_{a}+0.407{\tau }_{a}^{2}\\ & & \ +\ 0.0014[\mathrm{Fe}/{\rm{H}}]-0.0050(\mathrm{log}g-4.0).\end{array}\end{eqnarray} \tag{ 24 }$$

Using Equations (17)–(21) and (23)–(24) and sampling the posteriors of (T_eff,a, [Fe/H], $\mathrm{log}g$, q, a/R_a, i) from Tables 3 and 5, we calculate the predicted ellipsoidal distortion amplitudes: ${B}_{1}^{\mathrm{ELP}}=-8.9\pm 1.0\,\mathrm{ppm}$, ${B}_{2}^{\mathrm{ELP}}=-451\pm 51\,\mathrm{ppm}$, and ${B}_{3}^{\mathrm{ELP}}=-15.3\pm 1.8\,\mathrm{ppm}$.

We can instead use the more precise mass ratio ${q}^{* }\,={0.0932}_{-0.0058}^{+0.0068}$ derived from the Doppler-boosting term A₁ in the phase curve (see Section 5.1) to produce a separate set of predicted ellipsoidal distortion amplitudes. An analogous calculation yields ${B}_{1}^{\mathrm{ELP},* }=-{7.80}_{-0.58}^{+0.56}$ ppm, ${B}_{2}^{\mathrm{ELP},* }=-{396}_{-28}^{+27}$ ppm, and ${B}_{3}^{\mathrm{ELP},* }$ = $-13.5\pm 1.0\,\mathrm{ppm}$.

We summarize the measured values and theoretical predictions for the phase-curve harmonic amplitudes in Table 6 and Figure 4. We compute the total predicted values of all five harmonic terms with observed amplitudes significantly different from zero—A₁, A₂, B₁, B₂, and B₃. For the ellipsoidal distortion amplitudes, the table lists values assuming both the RV-derived mass ratio q and the mass ratio q* derived from the Doppler-boosting term; in Figure 4 only the prediction calculated from the RV-derived mass ratio is shown.

Zoom In Zoom Out Reset image size

Table 6.  Predicted and Observed Phase-curve Amplitudes (in ppm)

	Doppler	Mutual	Ellipsoidal	Ellipsoidal	Predicted	Predicted	Observed
	Boosting	Illumination	Distortion	Distortion (with q*)a		(with q*)a
A1	${114}_{-16}^{+18}$	⋯	⋯	⋯	${114}_{-16}^{+18}$	$\equiv {99.21}_{-0.86}^{+0.95}$	${99.21}_{-0.86}^{+0.95}$
A2	⋯	⋯	⋯	⋯	0	0	$-{39.75}_{-0.98}^{+0.88}$
B1	⋯	${261}_{-39}^{+44}$	−8.9 ± 1.0	$-{7.80}_{-0.58}^{+0.56}$	252 ± 42	253 ± 42	$268.4\pm 1.1$
B2	⋯	${62.0}_{-9.0}^{+10.3}$	−451 ± 51	$-{396}_{-28}^{+27}$	−389 ± 52	−334 ± 29	$-568.2\pm 1.0$
B3	⋯	⋯	−15.3 ± 1.8	−13.5 ± 1.0	−15.3 ± 1.8	−13.5 ± 1.0	$-{15.0}_{-1.1}^{+1.0}$

Note.

^aThese values are calculated using the mass ratio estimate ${q}^{* }$ inferred from the measured Doppler-boosting phase-curve amplitude ${A}_{1}={99.21}_{-0.86}^{+0.95}$ (see Section 5.1).

Download table as:  ASCIITypeset image

The observed amplitudes A₁, B₁, and B₃ all lie within 0.9σ of the predicted values. Doppler boosting is the only physical process that produces variation at the fundamental harmonic of the sine (A₁), while mutual illumination is the dominant contributor to the B₁ term. The consistency between the observed and predicted amplitudes indicates that the phase variations due to Doppler boosting and mutual illumination in the KOI-964 system are both well described by the theoretical formalism presented in Sections 5.1 and 5.2.

The first harmonic of the cosine at the orbital phase (B₂ term) is a combination of contributions from mutual illumination and ellipsoidal distortion, with the latter being the predominant source of the variation. From Table 6, we see that the total predicted amplitude, with the contribution from ellipsoidal distortion calculated using the mass ratio measured from RVs, differs from the observed value by 3.4σ. This discrepancy was first noted in Carter et al. (2011) based on the inconsistent mass ratios derived from Bayesian retrievals of the phase-curve amplitudes including Doppler boosting or ellipsoidal distortion. We can perform an analogous consistency test by calculating the predicted B₂ value assuming the more precise mass ratio q* derived from the observed Doppler-boosting variation. As shown in Table 6, the inconsistency between this predicted value and the observed amplitude is more severe—a 7.0σ discrepancy.

5.4.1. Dynamical Tide of the Host Star

5.4.2. Effects of Rapid Stellar Rotation