PROPERTIES OF AN ECLIPSING DOUBLE WHITE DWARF BINARY NLTT 11748

We observed NLTT 11748 with ULTRACAM (Dhillon et al. 2007) over 27 eclipses during 2010 and 2012, as summarized in Table 1. ULTRACAM provides simultaneous fast photometry through 3 filters with negligible dead time. During 2010, it was mounted on the 3.5 m New Technology Telescope (NTT) at La Silla Observatory, Chile. We used the u' and g' filters, along with either r' or i'. The integration times were chosen based on the conditions, but were typically 1–2 s for the redder filters and 5–8 s for the u' filter. During 2012, ULTRACAM was mounted on the 4.2 m William Herschel Telescope at the Observatorio del Roque de los Muchachos on the island of La Palma. Here, we used only the u'g'r' filters, although for one observation we discarded the r' data as they were corrupted. Exposure times were 2–3 s for the redder filters and 3–5 s for the u' filter, taking advantage of better conditions, a larger mirror, and a lower airmass toward this northern target. The total observing time for each eclipse was typically less than 40 minutes.

The data were reduced using custom software. We first determined bias and flatfield images appropriate for every observation. Then, we measured aperture photometry for NLTT 11748 and up to six reference sources that were typically somewhat brighter than NLTT 11748 itself. The aperture was sized according to the mean seeing for the observation, but was held fixed for a single observation and a single filter. Finally, a weighted mean of the reference star magnitudes (some of which were removed because of saturation and some of which were removed because of low signal-to-noise, especially in the u' band) was subtracted from the measurements for NLTT 11748. These detrended data were the final relative photometry that we used in all subsequent analysis. Given the short duration of the eclipse (≈3 minutes), any remaining variations in the relative photometry due to transparency or airmass changes could be ignored; we did not attempt to model the out-of-eclipse data (cf. Shporer et al. 2010).

In addition to the ULTRACAM observations, we observed a few eclipses using the Gemini Near-Infrared Imager (NIRI; Hodapp et al. 2003) on the 8 m Gemini-North telescope under program GN-2010B-Q-54. Given the low signal-to-noise and relatively long cadence, the primary eclipse data were not particularly useful in constraining the properties of the system. Instead, we concentrate on the secondary eclipse observations, where the additional wavelength coverage is helpful in constraining the effective temperature of the secondary (Section 4.2). The data were taken on 2010 November 21 and 2010 December 15 using the J-band filter. The total duration of the observations were 37 and 31 minutes around a secondary eclipse, as predicted from our initial ephemeris. Successive exposures happened roughly every 25 s, of which 20 s were actually accumulating data, so our overhead was about 20%.

To reduce the data, we used the nprepare task in IRAF, which adds various meta-data to the FITS files. We then corrected the data for nonlinearities⁹ and applied flatfields computed using the niflat task, which compared dome-flat exposures taken with the lamp on and off to obtain the true flatfield. We used our own routines to perform point spread function photometry with a Moffat (1969) function. This function was held constant for each observation and was fit to bright reference stars. Finally, we subtracted the mean of two bright reference stars to de-trend the data.

As discussed above, our one significant assumption (which was also made in Steinfadt et al. 2010) is that the secondary star follows the mass-radius relation for a CO WD. This seems reasonable, given inferences from observations (Kawka et al. 2010; Kilic et al. 2010) and from evolutionary theory. However, with the high precision of the current dataset, we must examine the choice of mass-radius relation closely. In particular, the zero-temperature model used in Steinfadt et al. (2010) is no longer sufficient. For effective temperatures near 7500 K and masses near 0.7 M_☉, finite-temperature models are larger than the zero-temperature models by roughly 2% (thin envelopes, taken to be 10⁻¹⁰ of the star's mass) to 6% (thick envelopes, taken to be 10⁻⁴ of the star's mass) and moreover the excess depends on mass (Fontaine et al. 2001; Bergeron et al. 2001). This excess is similar to what we observe in a limited series of models computed using Modules for Experiments in Stellar Astrophysics (Paxton et al. 2011, 2013). The envelope thickness can be constrained directly through astroseismology of pulsating WDs (ZZ Ceti stars), with most sources having thick envelope (fractional masses of 10⁻⁶ or above), but extending down such that roughly 10% of the sources have thin envelopes (fractional masses of 10⁻⁷ or below); there is a peak at thicker envelopes, but there is a broad distribution (Romero et al. 2012; consistent with the findings of Tremblay & Bergeron 2008).

We use the parameter r₂ to explore the envelope thickness. Values near 1.02 correspond to thin envelopes (with some slight mass dependence), while values near 1.06 correspond to thick envelopes. Changing r₂ over the range of values discussed above leads to changes in the best-fit physical parameters M₁, M₂, R₁, and R₂. In particular, M₁ was surprisingly sensitive to r₂. Determining a "correct" value for r₂ is beyond the scope of this work, but we can understand how the physical parameters scale with r₂ in a reasonably simple manner. Since we know the period accurately and can also say that sin i ≈ 1 (deeply eclipsing), we know that M₁ + M₂∝a³ (where a is the semi-major axis) from Kepler's third law. We also know K₁, which is the orbital speed of the primary:

$$\begin{equation} K_1 \approx \frac{2\pi a}{P_B} \frac{M_2}{M_1+M_2}, \end{equation} \tag{ 1 }$$

which allows us to constrain aM₂∝M₁ + M₂. Combining these gives M₂∝a². We can further parameterize the mass-radius relation of the secondary:

$$\begin{equation} R_2 \propto r_2 {M_2^\beta }. \end{equation} \tag{ 2 }$$

The duration of the eclipse fixes R₁/a, while the duration of ingress/egress fixes R₂/a (e.g., Winn 2011), so we can also say a∝R₂ or $a\propto r_2 M_2^\beta$. Combining this with M₂∝a² gives $M_2 \propto r_2^{2/(1-2\beta)}$. Near M₂ ≈ 0.7 M_☉, the mass–radius relation has a slope β ≈ −0.78 (compare with the traditional β = −1/3 for lower masses), so $M_2\propto r_2^{0.78}$. Since $a\propto M_2^{1/2}$ and both R₁ and R₂ are ∝a, $R_1\propto r_2^{1/(1-2\beta)}\propto r_2^{0.39}$ and the same for R₂ (at a fixed M₂, R₂∝r₂, as given in Equation(2); however, the result here is for the best-fit value of R₂, which can result in changes of the other parameters including M₂).

To understand how M₁ changes with r₂, we use the Keplerian mass function, which fixes $M_2^3 \propto (M_1+M_2)^2$. If we take the logarithmic derivative of this, we find:

$$\begin{equation} \frac{d\log M_1}{M_1} = \alpha \frac{d\log M_2}{M_2}\protect, \protect \end{equation} \tag{ 3 }$$

with α = 3/2 + M₂/2M₁. Since our mass ratio q is roughly 0.2 (Table 2), we find α ≈ 4.0. From this, we find that $M_1 \propto r_2^{2\alpha /(1-2\beta)}\propto r_2^{3.13}$, which is much steeper than the other dependencies. These relations are borne out by MCMC results (Table 2, Figures 2 and 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The model fit the data well, with a minimum χ² of 22978.8 for 22,566 degrees of freedom (dof).¹³ We show 1D and 2D marginalized confidence contours in Figure 2 and the best-fit light curves in Figure 1. The results are given in Table 2. A linear ephemeris gives a satisfactory fit to the data, although we find a significantly non-zero value for Δt, which we discuss below.

For NIRI, the quality of the data is modest, with typical uncertainties of ±0.03 mag and a cadence of 25 s. Given the quality of the ULTRACAM results, fitting the NIRI data with all parameters free would not add to the results. Instead, we kept the physical parameters fixed at their best-fit values from Table 2 and only fit for the eclipse depth at 1.25 μm.

The results are shown in Figure 4. Despite the modest quality of the data, the fit is good, with χ² = 160.7 for 161 dof. We find a depth d₂(1.25 μm) of 4.2% ± 0.4%, corresponding to a J-band flux ratio d₂(1.25 μm)/(1 − d₂(1.25 μm)) of 4.4%  ±  0.4% (see Section 4.2). This value is slightly off from the predictions based on our fit to the ULTRACAM data (Figure 5), although by less than 2σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To derive the eclipse times in Table 1, we used the results of the full eclipse fit but fit the model (with the shape parameters held fixed) to each observation individually.

Using the measured eclipse times from Table 1, along with those reported by Steinfadt et al. (2010), we computed a linear ephemeris for NLTT 11748 with a constant frequency f_B = 1/P_B. The residuals (Figure 6) are consistent with being flat and with the results from the full eclipse fitting (Table 2), showing no indication of orbital changes. However, we do find a systematic offset between the times of the primary and secondary eclipses, as predicted in Kaplan (2010). The secondary eclipses arrive earlier on average, by Δt = 4.1 ± 0.5 s (after correcting for this, the rms residual for the new data is 1.7 s and the overall χ² for the ephemeris data is 38.9 with 29 dof), which is consistent with Δt = 4.2 ± 0.6 s inferred from the full eclipse fitting (Table 2). The sign is correct for a delay caused by the light-travel delay across the system (the Rømer delay) when the more massive object is smaller, the magnitude of which is (Kaplan 2010)

$$\begin{equation} \Delta t_{\rm LT} = \frac{P_B K_1}{\pi c}\left(1-q\right). \end{equation} \tag{ 4 }$$

If we use Δt = 4.2  ±  0.6 s (from the full eclipse fitting in Table 2), then we infer q_LT = 0.29 ± 0.10. This is fully consistent with our fitted values for q (Table 2; q = 0.192 ± 0.008 for r₂ = 1.00). However, it is also possible that some time delay is caused by a finite eccentricity of the orbit (Kaplan 2010; Winn 2011), with

$$\begin{equation} \Delta t_{e}=\frac{2P_Be}{\pi }\cos \omega, \end{equation} \tag{ 5 }$$

where e is the eccentricity and ω is the argument of periastron.¹⁴ However, the Rømer delay must be present in the eclipse timing with a magnitude $(4.76\pm 0.05)r_2^{2.6}$ s, based on our mass determination. Therefore, instead of using the Rømer delay to constrain the masses, we can use it to constrain the eccentricity. Doing this gives ecos ω = (− 4 ± 5) × 10⁻⁵ (consistent with a circular orbit). This would then be one of the strongest constraints on the eccentricity of any system without a pulsar, as long as the value of ω is not particularly close to π/2 (or 3π/2). This may be testable with long-term monitoring of NLTT 11748, as relativistic apsidal precession (Blandford & Teukolsky 1976) should be $\dot{\omega }\approx 2\deg \,{\rm yr}^{-1}$ (for a nominal e = 10⁻³). As long as any tidal precession is on a longer timescale, the change in ω could separate the Rømer delay from that due to a finite eccentricity. For a system as wide as NLTT 11748  tidal effects are likely to be negligible (Fuller & Lai 2013; Burkart et al. 2013). Further relativistic effects may be harder to disentangle: an orbital period decay $\dot{P}_B$ will be of a magnitude −1 μs yr⁻¹, compared with a period derivative from the Shklovskii (1970) effect¹⁵ of +24 μs yr⁻¹ (current data do not strongly constrain $\dot{P}_B$ because of the short baseline, with $\dot{P}_B=(-1900\pm 900)\,\mu {\rm s}\,{\rm yr}^{-1}$). Therefore, an accurate determination of the distance (whose uncertainties currently dominate the uncertainty in the period derivative) will be necessary before any relativistic $\dot{P}_B$ can be measured.

Zoom In Zoom Out Reset image size

In comparison, the radial velocity constraints on the eccentricity are considerably weaker: Steinfadt et al. (2010) determined e < 0.06 (3σ), while Kilic et al. (2010) and Kawka et al. (2010) both assumed circular orbits. We refit all of the available radial velocity measurements from those three papers. Assuming a circular orbit, we find K₁ = 273.3 ± 0.4 km s⁻¹ (like Kilic et al. 2010, whose data dominate the fit). We also tried the eccentric orbit as parameterized by Damour & Taylor (1991). The fit is consistent with a circular orbit, with ecos ω = 0.001 ± 0.004 and esin ω = −0.004 ± 0.004, which limits e < 0.017 (3σ). The eclipse durations can also weakly constrain the eccentricity, with the ratio of the secondary eclipse to the primary eclipse duration roughly given by 1 + 2esin ω (for e ≪ 1), although as noted by Winn (2011) this is typically less useful than the constraints on ecos ω from eclipse timing. In the future, we can fit for this term directly.

While all data are satisfactorily fit by just a simple linear ephemeris, we can also ask if a third body could be present in the system. Such a body, especially if on an inclined orbit, could significantly speed up the merger of the inner binary and may alter the evolution of the system (Thompson 2011). Any putative tertiary would likely be in a more distant circumbinary orbit, since the interactions necessary to produce the ELM WD would have disrupted closer companions. A tertiary would produce transit timing variations (Holman & Murray 2005; Agol et al. 2005), moving the eclipse times we measure. A full analysis of transit timing variations, including nonlinear orbital interactions, is beyond the scope of this paper. Instead, we did a limited analysis where we considered the system to be sufficiently hierarchical such that the inner binary was unperturbed (consistent with our measurements) and only its center of mass moved due to the presence of the tertiary. We took the ephemeris residuals and fit a variety of periodic models, determining for each trial period what the maximum amplitude could be (marginalizing over phase) such that χ² increased by 1 from the linear ephemeris fit (adding additional terms in general decreases χ², but we wanted to see what the maximum possible amplitude could be). We found that for periods of 1–300 days, the limit on any sinusoidal component was ≲ 1 s (consistent with the rms discussed above) or smaller than the orbit of the inner binary. Therefore, unless it is highly inclined, no tertiary with such a period is possible. As we get to periods that are longer than 300 days, we no longer have sufficient data to constrain a periodic signal, but here the constraints from our polynomial fit also exclude any stellar-mass companion (an amplitude of 1 s at a period of 300 days would require a mass of 0.002 M_☉ if the outer orbit is also edge-on).

As with the NIRI data, we can separate the fitting of the secondary eclipse depth from the rest of the eclipse fitting and derive the eclipse depths as a function of wavelength, d₂(λ) (where we also separate the 2010 and 2012 ULTRACAM observations). Given the eight secondary eclipse depths that we measure, we can determine the ratio of the radius of the secondary to the primary R₂/R₁ as well as the temperature of the secondary, T₂, given measurements of T₁. For that, we use the determination by Kilic et al. (2010)—T₁ = 8690 ± 140 K—largely consistent with the value determined by Kawka et al. (2010) of 8580 ± 50 K. The secondary eclipse depths are related to the wavelength-dependent flux ratios:

$$\begin{equation} f(\lambda) \equiv \frac{F_{2}(\lambda)}{F_1(\lambda)} = \frac{d_2(\lambda)}{1-d_2(\lambda)} = \frac{R_2^2 10^{-m(\lambda,T_2)/2.5}}{R_1^2 10^{-m(\lambda,T_1)/2.5}}, \end{equation} \tag{ 6 }$$

where m(λ, T) is the absolute magnitude of a fiducial WD with a temperature T at wavelength λ.

To determine the secondary's temperature from the eclipse depths, we first sample T₁ from the distribution ${\cal N}(8690,140)$ 1000 times. Then, for each sample, we interpolate the synthetic photometry of Tremblay et al. (2011; using the 0.2 M_☉ model) to determine the photometry for the primary, m(λ, T₁). We then solve for the temperature of the secondary and the radius ratio by minimizing the χ² statistic, comparing our measured f(λ) with the synthetic values (using the 0.7 M_☉ grid for the secondary; the results did not change if we used the 0.6 M_☉ or 0.8 M_☉ grids instead).

The resulting distribution of T₂ and (R₂/R₁) is independent of any assumptions in the global eclipse fitting. We find T₂ = 7643 ± 94 K and R₂/R₁ = 0.255 ± 0.002. However, the fits had an average χ² = 12.6 for 5 dof, mostly coming from a small mismatch between the inferred eclipse depth at g' measured in 2010 versus 2012. If we increase the uncertainties to have a reduced χ² = 1, then we find T₂ = 7643 ± 150 K and R₂/R₁ = 0.255 ± 0.003. Note that the radius ratio here is fully consistent with that inferred from the fit to the rest of the eclipse shape (Section 3 and Table 2).

We show the results in Figure 5 (including the results of the NIRI data analysis), along with the results using blackbodies for the flux distributions instead of the synthetic photometry (as had been done by Steinfadt et al. 2010 and others). It is apparent that the blackbody does not agree well and using the synthetic photometry is vital.

The constraints using the separate wavelength-dependent eclipse depths ended up being slightly less precise than the value from the full eclipse fitting, although the two constraints are entirely consistent. We therefore choose the values from Table 2, where we determined a secondary temperature of T₂ = 7600  ±  120 K. Using the thin (thick) hydrogen atmosphere models for a 0.7 M_☉ CO WD, we find a secondary age τ₂ = 1.70 ± 0.09 Gyr (1.58 ± 0.07 Gyr) by interpolating the cooling curves from Tremblay et al. (2011). So, the envelope uncertainty does not contribute significantly to the uncertainty in the age of the secondary directly. A bigger contribution is through changes in the secondary mass M₂.

Based on the measured J-band photometry (J = 15.84  ±  0.08 from Skrutskie et al. 2006)¹⁶, along with an estimate for the extinction (E(B − V) = 0.10 and R_V = 3.2, from Kilic et al. 2010), we can compare our measured radius with that inferred from a parallax measurement (π = 5.6 ± 0.9 mas; H. Harris 2011, private communication). We again use the Tremblay et al. (2011) synthetic photometry for the bolometric correction at the temperature determined by Kilic et al. (2010) and use A_J = 0.29A_V. Based on these data, we infer R_{1, phot} = 0.049 ± 0.009 R_☉. This is fully consistent (within 0.5σ) with our inferred values from the eclipse fitting. We can use this value along with the radius ratio inferred from the eclipse shape (roughly R₂/R₁ = 0.2567 ± 0.0006) to determine R₂ = 0.013 ± 0.002 R_☉. From here, we can calculate M_{2, thin} = 0.57 ± 0.08 M_☉ and M_{2, thick} = 0.61 ± 0.11 M_☉, which makes use of evolutionary models that were interpolated to the correct effective temperature (Fontaine et al. 2001; Bergeron et al. 2001). These masses are a little lower than our secondary masses calculated by the eclipse fitting, but differ by less than 2σ. Inverting the problem, we infer based on our eclipse fitting for r₂ = 1.02 a distance d = 159 ± 8 pc (π = 6.3 ± 0.3 mas), with the uncertainty dominated by the uncertainty in the photometry.

Our analysis above included the effects of in-binary microlensing. Steinfadt et al. (2010) assumed microlensing would modify the primary eclipse depth, but did not show definitively that it was required for a good fit (cf. Muirhead et al. 2013). We fit the same data with the same procedure as in Section 3, but did not allow for any decrease in the depth of the primary eclipse from lensing. The resulting fit was adequate, although slightly worse than with lensing. With no lensing amplification, to match the depth of the primary eclipse requires a reduced value of R₂/R₁ or a lower inclination. We can do that by increasing R₁ or decreasing R₂, but that is difficult as the eclipse durations fix R₁/a and R₂/a. We end up accomplishing this by increasing the masses, which widens the orbit (increasing a to go along with the increase in R₁ and decreasing R₂ through the mass-radius relation for the WD). For r₂ = 1.00, we find M₁ = 0.16 M_☉ and M₂ = 0.74 M_☉, along with R₁ = 0.044 R_☉ (log (g)₁ = 6.36 ± 0.04). While this combination of parameters does not give as good a fit to the data as the fit with lensing, the difference is not statistically significant, with an increase in χ² of 20 (the reduced χ² increased from 1.018 to 1.019, for a chance of occurrence of about 50%). We would need other independent information (such as a much more precise parallax or a more precise time delay) to break the degeneracies.