KOI 1224, A Fourth Bloated Hot White Dwarf Companion Found With Kepler
We present an analysis and interpretation of the Kepler binary system KOI 1224. This is the fourth binary found with Kepler that consists of a thermally bloated, hot white dwarf in a close orbit with a more or less normal star of spectral class A or F. As we show, KOI 1224 contains a white dwarf with K, mass , and radius , and an F-star companion of mass that is somewhat beyond its terminal-age main sequence. The orbital period is quite short at 2.69802 days. The ingredients that are used in the analysis are the Kepler binary light curve, including the detection of the Doppler boosting effect; the NUV and FUV fluxes from the Galex images of this object; an estimate of the spectral type of the F-star companion; and evolutionary models of the companion designed to match its effective temperature and mean density. The light curve is modelled with a new code named Icarus which we describe in detail. Its features include the full treatment of orbital phase-resolved spectroscopy, Doppler boosting, irradiation effects and transits/eclipses, which are particularly suited to irradiated eclipsing binaries. We interpret the KOI 1224 system in terms of its likely evolutionary history. We infer that this type of system, containing a bloated hot white dwarf, is the direct descendant of an Algol-type binary. In spite of this basic understanding of the origin of KOI 1224, we discuss a number of problems associated with producing this type of system with this short of an short orbital period.
Subject headings:binaries: eclipsing — stars: evolution — stars: individual (KOI~1224 (catalog KOI 1224)) — techniques: photometric — white dwarfs
Close binary systems containing a single compact star (i.e., white dwarf, neutron star, or black hole) have the potential to enhance our knowledge of stellar evolution, perhaps even more so than for double degenerate systems whose evolution is yet more complicated. In almost all cases, the compact star originated from the remnant core of the original primary star in the binary system. The primary’s envelope was lost via some combination of Roche-lobe overflow mass transfer, stellar wind, common envelope, or explosive process (e.g., a supernova).
In cases where the original secondary star is now transferring mass back to the compact star, the systems may be highly detectable via the release of gravitational potential energy in the form of optical, UV, or X-ray radiation. However, systems in which the compact star has not yet started to accrete may be much more difficult to discover. This is especially true in the case of white dwarfs in orbit with intrinsically much brighter stars that outshine them (e.g., Regulus; Gies et al., 2008; Rappaport et al., 2009). Even in the case where the observer is in a fortuitous orientation to possibly see eclipses, such events may be of very small amplitude, e.g., at the level if the white dwarf has cooled to near the temperature of the parent star. Due to the extraordinary photometric precision of the Kepler mission, however, eclipses and transits of this depth are readily detectable.
Since the Kepler mission has been launched, there have in fact been a number of such binaries containing a white dwarf that have been discovered (Rowe et al., 2010; van Kerkwijk et al., 2010; Carter et al., 2011). In those systems the white dwarf is still sufficiently hot (i.e., K) and thermally bloated (i.e., up to 10 times their degenerate radii) that the eclipses and transits are in the range of up to a couple of percent of the system light. We report on the study of a fourth such system with a number of similarities to the others, KOI 1224 (KIC 6606653; Borucki et al. 2011; Prša et al. 2011); however, it has the shortest among the orbital periods at 2.69802 days111We note that a new binary system, 1SWASP J024743.37251549.2, recently discovered with the WASP survey (Maxted et al., 2011) contains a bloated hot white dwarf in an even shorter period binary with days., and the primary star in this system has the lowest effective temperature and is the most evolved away from the main sequence of any of the previously detected hot white dwarf systems.
We report here an analysis of the Kepler Q1 and Q2 public data for KOI 1224, as well as supplementary observations from the GALEX satellite. We make use of stellar evolution models to better estimate the mass and evolutionary state of the current primary F star in the system.
We analyze the Kepler light curve using Icarus, a newly developed state-of-the-art binary light curve synthesis code which is aimed at the study of detached and semi-detached binaries. This code has the advantage of handling spectroscopic data as well as Doppler boosting, irradiation effects and transits/eclipses, hence making it particularly attractive for irradiated binaries. An eclipsing binary such as KOI 1224 represents an ideal testbed to benchmark the synthesis code because its results can be easily compared with those of a semi-analytic analysis. The synthesis code is described in some detail in this work.
The organization of the paper is as follows. In §2 we introduce the new binary light curve synthesis code. Some detailed specific aspects of the code are given in Appendix A. After presenting KOI 1224 and the available data in §3, we proceed with the light curve modelling in §4, first semi-analytically and then with our synthesis code. We discuss the probable evolution of KOI 1224, its implications for understanding mass transfer, and its relation to Algol systems in §5. We summarize our findings and draw some general conclusions in §6.
2. Binary Light Curve Synthesis Code
We have implemented our own light curve synthesis code with the purpose of modelling the light curves and spectra of detached and semi-detached binaries. Our code, called Icarus222Icarus is freely available at http://www.astro.utoronto.ca/$∼$breton/icarus., is most similar to the ELC code by Orosz & Hauschildt (2000). Essentially, the code constructs a finite-element surface grid for a star having some pre-determined physical and binary properties by solving the gravitational equipotential equation. Each surface element is characterized by a set of physical properties (i.e., temperature and an effective gravity), and the observed flux is obtained by integrating the specific intensity emerging from the surface visible by an observer located in a given direction. In the event that both stars contribute significantly to the total observed light, we model the second star of the system by inverting the mass ratio and adding a half orbital shift to it.
The core implementation is made in Python and relies heavily on the Scipy and Numpy libraries333Freely available at http://scipy.org/. In order to bolster the execution speed, critical components are written in C and included directly in the Python code with the Scipy Weave module444See http://www.scipy.org/Weave. The modular, object-oriented approach that we have adopted should facilitate adding features to the code. Hence, the code currently includes various effects such as eclipses/transits, irradiation from a companion, ellipsoidal light variations, and Doppler boosting. As opposed to other synthesis codes, such as the Wilson-Devinney code (Wilson & Devinney, 1971), or its modern incarnation in PHOEBE (Prša & Zwitter, 2005), which bundle together the binary synthesis and the parameter optimization in a single piece of software, ours resembles more a suite of routines that can be glued together in order to perform binary modelling.
The main module of our code is the physics engine, which is responsible for synthesizing the star. A separate module deals with the atmosphere model and provides utility functions that return the intensity given a set of properties – typically the temperature, surface gravity, emission angle and, sometimes, the velocity. This layer works independently of the binary light curve physics and serves only as a frontend to retrieve the specific intensities from a synthetic atmosphere model lookup table or an analytic model (e.g., blackbody). These atmosphere models can be integrated over a passband or used with their full spectral content. In the same spirit, the physics engine does not tamper with the returned values and leaves it to the user layer to alter them if necessary (e.g., to resample the spectrum or apply reddening). The user layer exists to connect the physics engine, the atmosphere models module, and the experimental data. It typically contains (1) a data reader function; (2) a utility function that takes user input parameters, recasts them into the physics engine input parameter format, and post-processes the model data to match the observed ones; and (3) a function that returns the fit residuals or goodness-of-fit. Finally, on top of the user layer, a fitting layer can be added. The function that returns the goodness-of-fit can be fed within the framework of a minimization procedure, a brute-force search, a Bayesian optimization, or any other algorithm suited to the needs of the problem.
The main innovation of our light curve synthesis code is that it can handle phase-resolved spectroscopic light curves and hence it allows one to perform modelling without the need for fitting radial velocities separately using a template and having to introduce a correction factor for the displacement of the light-center with respect to the barycenter. This capability will be demonstrated in a forthcoming paper. Thorough details about our code are provided in Appendix A.
3. Data and Data Reduction
We have obtained the KOI 1224 data from the Kepler public archive555http://archive.stsci.edu/kepler/publiclightcurves.html. The original raw light curve, presented in Figure 1, includes the first and second quarters (Q1 and Q2) data, which were collected using a 30-min integration time. We removed the offset between the two quarters and corrected for three jumps in the flux. Next we cleaned the raw light curve from outliers using the following method. We removed the systematic trend using a 7th-degree polynomial and folded the light curve at the orbital period. We calculated the running median and running standard deviation statistics using a 20-point window in order to identify outliers deviating by more than 1.5 standard deviations of the data scatter around the median, which we flagged and discarded for the rest of the analysis. This threshold was chosen in order to remove the obvious visible outliers as some data points were clearly beyond the regular data scatter. In total, the outliers removal discarded 699 out of the 5723 data points.
3.1. Ultraviolet Galex Data and Other Photometry
In order to constrain our modelling, we started with information from the Kepler Input Catalog (Brown et al., 2011), which lists optical and infrared 2MASS photometry as well as estimates of the temperature and reddening (see Table 1). The reddening is inferred using a simple exponential dust model, which, for this source, is consistent with the total dust column of expected along this line of sight from the COBE/DIRBE and IRAS based dust maps of Schlegel et al. (1998). Comparing the and colours (dereddened using the coefficients of Schlegel et al. 1998) to those inferred from spectrophotometry of MK standards (Covey et al., 2007), we infer a spectral type F6V, with an uncertainty of 1 subclass, corresponding to a temperature of (Cox, 2000). This is consistent with that of the Kepler Input Catalog, within the expected uncertainty of K (Brown et al., 2011). From a spectrum taken at the 1.6-m Observatoire du Mont-Megantic (L. Nelson, private communication) we infer a spectral type F5-6 III-IV (more likely luminosity class IV), which implies a temperature similar to that quoted above (Gray et al., 2001). Below, we will use .
We checked other archives for additional information, and found that KOI 1224 was detected by GALEX (see Table 1), with and . The estimated intrinsic color is bluer than expected for the primary (; Bianchi et al. 2007; Vennes et al. 2011), but we cannot exclude that it could significantly contribute to the NUV flux given the uncertainties in the temperature, the reddening and the extinction curve. However, no contribution from the primary is expected to the FUV flux, and indeed, with ,666For a standard extinction curve, (Rey et al., 2007). Given the small difference, the intrinsic color is secure. the source is bluer in the UV than any of the brighter stellar sources near it. Thus, the FUV emission is almost certainly from the hot white dwarf secondary.
In order to estimate a temperature for the white dwarf, we can use the fact that when the white dwarf is eclipsed, the flux in the Kepler bandpass decreases by 1.4%. Here, the 1.4% drop includes a 10% correction for a “third light” (e.g., a fainter blended star within the Kepler postage stamp), and hence it represents the white dwarf’s contribution to the system’s luminosity. From color terms of Brown et al. (2011), the -band flux should decrease about 10% less. We thus estimate and infer and . We can compare this with predicted GALEX and optical colors for white dwarfs from Vennes et al. (2011). For their lowest mass, , they list for to 16000 K in steps of 1000 K, with the color depending somewhat on gravity below 13000 K (here, we converted the V-band magnitudes of Vennes et al. (2011) to -band using , valid for temperatures in this range; Fukugita et al. 1996). Thus, we infer a temperature of K.
For the above temperature, the models also predict , which is much bluer than the observed color. Apparently, therefore, either our estimate is inaccurate and the true temperature is , or the primary contributes of the NUV flux, i.e., is about a magnitude brighter than expected (perhaps because of stellar activity associated with rapid rotation). From the ratio of the transit to eclipse depths (see §4.1), it seems clear the latter is the correct interpretation, and the white dwarf effective temperature is K. This is further confirmed by our formal fits to the light curve (see §4.2).
4. Modelling the Light Curve of KOI 1224
4.1. Semi-Analytic Modelling
We first performed a semi-analytical analysis of the KOI 1224 light curve in order to obtain rough estimates of the system parameters. We fitted the raw data (see Figure 1) for three out-of-eclipse orbital harmonics as well as a 7th-degree polynomial and three “pulsations” (at 3.49, 1.75 and 1.17 d). The polynomial function aims to remove the long-term instrumental variations, while the pulsations were detected as significant signals in a power spectrum. The 3.49-day pulsation is strongest and might reflect the rotation period of the primary, as seemed to be the case for KOI 74, where a possible pulsation period of matched the rotation period inferred from its rotational velocity. If so, it would imply the star is rotating sub-synchronously, perhaps because the star is expanding during the course of its evolution. We then removed the systematic trend from the data and folded them at the orbital period.
We measured various fiducial points on the white dwarf eclipse profile using a detailed plot of the light curve (Figure 2) in order to estimate the radius of the primary star in units of the semi-major axis, i.e., . We found the phase of the first contact, the start of the full eclipse, the end of the full eclipse, and the last contact to be , , and , respectively. From these values, measured in orbital cycles from the primary’s inferior conjunction, we can infer and (implicitly assuming , since we have no other way to infer it from simple geometry).
The first three out-of-eclipse orbital harmonics capture information about the system properties that can be derived from ellipsoidal light variations, irradiation effects, and Doppler boosting (see the folded out-of-eclipse light curve in Figure 3). From our fit we found the fractional amplitudes , and along with the corresponding phases , and . The phases are in orbital cycles measured from the ascending node of the white dwarf for and from the white dwarf eclipse for . Because the ingress and egress time scales are comparable to that of the Kepler integration time, it would be futile to make any inference about without a proper numerical treatment, as in §4.2.
In the simple, semi-analytic analysis that follows, we correct amplitudes and up by a factor of 10% to account for the third light in the system, and up by an additional 6% to account for the much smaller expected Doppler boosting effect of the white dwarf which decreases the total amplitude of (see, e.g., Carter et al., 2011, equation 11). We label these “corrected” amplitudes and .
From the color information of KOI 1224 presented in §3.1, we estimate at 6350 K and . Using the above expression we find that the harmonic amplitude implies a projected velocity amplitude km/s, where the cited uncertainty includes 3% each from the amplitude and . From we find the mass function:
The second orbital harmonic is due to the ellipsoidal light variations, and has an expected amplitude:
where is a pre-factor of order unity that depends on the limb and gravity darkening:
(Morris, 1985)777Here we dropped the factor from the original equation because it is negligible. Using the corrected amplitude, and the geometrically measured value of , and considering that , we find a mass ratio . We also find that if we take the linear limb darkening for Kepler’s bandpass at K from Sing (2010), and by integrating Equation 10 from Morris (1985) over the Kepler bandpass. Hence, we estimate that . One also sees that the ellipsoidal term is phased as expected with the ascending node.
Combining the mass function found from the amplitude and the mass ratio from the term, we find
Second order terms from ellipsoidal light variations arise because the nose of the star (i.e., near the L point) is darker than its back side. This results in the third harmonic, with a relative amplitude . This values differs slightly with that of 0.06 calculated with the equations of Morris (1985), assuming the same coefficient as above.
Lastly, there should also be a contribution to the fundamental from the ellipsoidal light variation effect, with an amplitude of , with a maximum at the white dwarf’s eclipse. Instead, the fundamental modulation is a bit shifted towards the white dwarf transit, with amplitude . This indicates that the shift of the fundamental modulation is mainly due to irradiation, rather than second-order ellipsoidal light variations. Assuming the second-order contribution due to the ellipsoidal light variation is correct, we infer an irradiation term of . According to harmonic decomposition of the reflection effects in close binaries by Kopal (1959), we can write the fractional amplitude of the first harmonic term due to reflection as:
where is the primary’s albedo, is the ratio of the secondary’s luminosity to the total luminosity. In this estimate we have ignored the leading-order term (due to the irradiation of the secondary by the primary) because it is nearly two orders of magnitude smaller than the term – as we find in §4.2 (see Table 2). From this we infer , which is quite consistent with the drop in the light curve during the white dwarf eclipse.
In the next section we discuss our numerical modelling from which we obtain more accurate results for the system masses. The main reason for this improved accuracy is that we utilze measured properties of the primary star to directly determine its mass. In that case, the use of in determining the mass ratio is less critical, and, more importantly, then depends principally on the cube root of the mass function, and so is only linearly dependent on (rather than on its cube).
4.2. Numerical Modelling
In order to fully describe the light curves, our code includes the following list of free parameters (with indices 1 and 2 referring to the F-star primary and the white dwarf secondary, respectively), which are also summarized in Table 2: the mass ratio, ; the filling factor888See §A.2 for definition., ; the (polar) temperature ratio, ; the primary’s (polar) temperature, ; the irradiation temperature ratio, (see §A.4 for the definition of ); the orbital inclination, ; and a third light contribution, .
|Orbital period, (d)||2.69802(2)|
|Epoch of primary’s inferior conjunction, (MJD)||55030.038(1)|
|Filling factor, (fraction of aafootnotemark: )||0.389(6)|
|Filling factor, (fraction of aafootnotemark: )||0.034(2)|
|Orbital inclination, (degree)||85(1)|
|Third-light, (fraction of total light)||0.10(5)|
|Extra noise, (fraction of average flux errors)||3.69(4)|
|Mass primary, ()||1.59(7)|
|Mass secondary, ()||0.20(2)|
|Semi-major axis, ()||9.9(2)|
|Projected velocity amplitude, (km s)||21.0(1.7)|
|Relative volume-averaged radius,||0.270(5)|
|Relative volume-averaged radius,||0.0103(4)|
|Volume-average density, (g cm)||0.118(6)|
|Volume-average density, (g cm)||270(30)|
|Polar surface gravity, ()||3.79(2)|
|Polar surface gravity, ()||5.72(5)|
One more parameter is required in order to fully specify the system dynamics and solve for the masses. This is required in order to determine the projected velocity amplitude of the primary, , which, in turn, allows us to calculate the Doppler boosting amplitude. Instead of introducing this quantity as an additional free parameter to the fit, we use the primary mass, , itself; this is possible because the mean stellar density of the primary, , and its effective temperature 999For this purpose, we approximated the effective temperature by , the polar temperature. determine a nearly unique mass and age, provided that the star has not evolved far up the giant branch. The average density of the primary star, in turn, can be found using Kepler’s third law, which depends only on known quantities such as the volume averaged stellar radius 101010 is found after solving for the equipotential surface, which does not require knowing the masses in the system.:
We have used the newly developed MESA stellar evolution code (Paxton et al., 2011) to evolve a sequence of stars between 1.1 and 2.0 in steps of 0.1 . These tracks are shown in Figure 4 in the plane. Given that there is an uncertainty in both and , and the fact that the evolution tracks cross each other near the TAMS, we have converted this discrete set of tracks to a nearly continuous mass distribution in the plane (see Figure 4). To generate this finely gridded () mass array we generated a probability weighting function associated with each point on the evolution tracks from a 2-dimensional Gaussian distribution in and – as described in the caption to Figure 4. If more than one probability was assigned to a given point in the plane for any given track, we took the highest of the probabilities. Finally, we used these probabilities to compute a weighted “mean” mass at each location in the plane.
We fixed the orbital ephemerides at the values reported in Table 2 and, because the rotational periods of the stars are unknown, we have assumed that they are rotating co-synchronously with the orbit for the calculation of the equipotential surfaces. We also used a gravity darkening coefficient for both stars based on empirical values from Che et al. (2011)111111This is appropriate for the primary star of K, and the hotter white dwarf is sufficiently small that its gravity darkening has a negligible effect.. As we pointed out in §4.1, irradiation effects on the white dwarf due to the primary are negligible (2% of the primary’s irradiation by the secondary) and hence can be safely neglected. From exploratory fitting, we also realized that the flux uncertainties were insufficient to account for the root-mean-square of the fit residuals (i.e., the was large even though the fits visually looked good). Hence, we introduced an extra noise parameter, , which was added in quadrature to the flux uncertainties.
The process of calculating the goodness of fit of our 9-parameter model works as follows. We evaluate the synthetic light curve using a variable sampling scheme (i.e., more densely sampled near the transit/eclipse and less densely in the smoother ellipsoidal light variations region). This allows us to boost the computing speed while still capturing the critical details. The 30-min integration time of Kepler produces a non-negligible smearing that needs to be accounted for. Hence, we resample the synthetic light curve at a higher resolution using a second-order spline in order to perform a 30-min boxcar filtering. The synthetic light curve is then interpolated at the phase of each observed data point. Next, the systematic trend in the data is removed after fitting a 7th-order polynomial and three sinusoidal pulsations having unknown phases, amplitudes and harmonically related periods , and , where the fundamental is d. Finally, the residuals are calculated using the detrended light curve.
4.2.1 MCMC Algorithm
We performed the parameter optimization using a Markov Chain Monte Carlo (MCMC) algorithm since it is well suited to making parameter inferences in high-dimensional problems. Parameter correlations are sometimes complicated and hence we used a Gibbs sampler in our MCMC because it allows for an easier tuning of the proposal distribution’s acceptance rate (see e.g. Walsh, 2004; Gregory, 2005, for more information about MCMC and Gibbs sampling).
We sampled the orbital inclination using flat priors in , such that it eliminates projection effects. Based on the Kepler Input Catalog estimation of the third light (see §3.1), which is listed as , we used a logistic prior for this parameter:
In such way, the probability of values in the range is close to unity before dropping to 50% at 0.15, and then becomes negligible beyond 0.2. For the extra noise parameter, we used modified Jeffrey’s priors , which provides equal weighting per decade. This choice is motivated by the scaling nature of the extra noise parameter. The constant of 0.1 linearizes the prior when becomes smaller than 10% of the average flux measurement errors. In a Bayesian parameter estimation framework, the extra noise parameter has an effect similar to normalizing the reduced to unity in standard fitting. Since Kepler uses a single broadband filter and the data contain no color information, we also made use of the photometry constraints (see §3.1) to add a Gaussian prior of K on the primary’s temperature. Flat priors were assumed for the remaining parameters.
We ran a total of 6 independent MCMC chains, each consisting of 400,000 steps. A burn-in period of 20,000 steps was removed from each chain and we kept only every 360 iterations – a process called thinning – in order to reduce auto-correlation effects. We compared the summary statistics of each MCMC to the others in order to ensure that convergence has been reached. For all model parameters, the Gelman-Rubins convergence diagnostic (Gelman & Rubin, 1992) yielded , which indicates that the chains have reached stationary distributions. The final results presented in Table 2 and Figure 5 were obtained by merging all the chains together (after burn-in and thinning), while Figures 2 and 3 display an example of the best-fit light curve.
We find constituent masses of and for the primary and white dwarf, respectively. The radius of the primary is , indicating that it has evolved substantially off the zero-age main sequence, consistent with its low mean density of 0.12 g cm. The white dwarf’s radius, , is 5 times its degenerate radius, and indicates that it is still actively cooling (possibly punctuated with episodic shell flashes). In turn, this indicates that its cooling age is of order 120 – 600 Myr for masses in the range of (see, e.g., Driebe et al., 1999; Nelson et al., 2004). We discuss below what these properties may imply for the formation of KOI 1224.
5. The Origin of KOI 1224
We have presented results for a system found with the Kepler mission that contains a thermally bloated hot white dwarf: KOI 1224 (KIC 6606653). The system consists of a normal F6 star that has evolved to near the TAMS in a 2.698-day binary with a white dwarf of mass and radius . We discuss in this section how such a binary may have evolved.
There are two plausible evolutionary paths to the formation of the current KOI 1224 binary system. The transfer of mass from the white dwarf progenitor to the primordial secondary was either dynamically unstable, leading to a common envelope phase, or it was stable, even if not completely conservative. The former case leads to a dramatic shrinkage of the orbital separation while the latter results in only modest changes in the orbit. We consider each of these, in turn, and in more detail. For more discussion of the formation and evolution of these systems see Farmer & Agol (2003), Di Stefano (2011), van Kerkwijk et al. (2010) and Carter et al. (2011).
5.1. Common Envelope Scenario
In the common-envelope scenario, the primary star develops a degenerate He core of mass , starts to transfer mass to the secondary, and initiates a dynamically unstable process leading to a common-envelope phase. The secondary spirals into the common envelope of the primary, ejects the envelope, and this results in a very close orbit of the secondary with the He core of the primary (i.e., the current white dwarf of the system). The final orbital separation can be related to the initial (pre-CE) separation by
(see, e.g. de Kool, 1990; Podsiadlowski et al., 2003) where , , and are the core, envelope, and total mass of the primordial primary star, respectively, is the mass of the secondary, is the Roche-lobe radius of the primary in units of the orbital separation, and the product encapsulates the energy ejection efficiency and binding energy of the common envelope. The latter CE parameter is often taken to be 1, but is more likely to for an evolved star (see, e.g. Dewi & Tauris, 2000; Tauris & Dewi, 2001; Podsiadlowski et al., 2003). In this latter case, the orbit must shrink by a factor of 100. For a final orbital period of 2.7 days, the initial period would have to have been years. Such a wide initial orbit for the progenitor binary would necessarily imply the development of a much more massive He (or likely CO) white dwarf of . Thus, given the fact that the current white dwarf in KOI 1224 is much lower in mass, this seems to essentially rule out a common-envelope scenario.
5.2. Stable Mass-Transfer Scenario
Stable, but not necessarily conservative, mass transfer from the progenitor of the white dwarf to the secondary seems the much more likely route to the production of KOI 1224, though, as we discuss, this formation scenario is not without its own difficulties. In this scenario, when the primary fills its Roche lobe and commences mass transfer to the secondary, that mass transfer will proceed on the thermal timescale of the envelope of the donor star (see, e.g., Podsiadlowski et al., 2002; Lin et al., 2011) until well after the masses of the two stars have become equal. Thermal timescale mass transfer results since the donor star is more massive than the accretor and it is not too evolved. The remainder of the mass transfer is driven by the nuclear evolution of the primary and is generally much slower than the thermal timescale transfer (Lin et al., 2011).
Depending on the relative masses of the primary and secondary at the time mass transfer commences, the outer envelope of the accreting star may swell up to the point where further accretion is inhibited and a substantial fraction of the transferred mass is ejected from the system (see, e.g. Kippenhahn & Meyer-Hofmeister, 1977; van Rensbergen et al., 2010, 2011). If we define the mass retention fraction by the accretor to be , the value of is likely to be low at the start of the mass transfer process and higher toward the end.
As we have shown, the age of the white dwarf in KOI 1224 is likely to be Myr because it is still quite hot (see, e.g. Driebe et al., 1999; Nelson et al., 2004). The fact that the current primary star in KOI 1224 has a mass of 1.6 and has evolved to at least the TAMS (and probably somewhat beyond) implies an age of Gyr for a single, isolated star of that mass. These two facts taken together imply that the current primary must have been already substantially evolved at the end of the mass-transfer phase when the white dwarf was unveiled. From all of this we can infer that the primordial secondary (now the current primary) was not substantially less massive than the primordial primary.
Therefore, while there likely was a brief phase of thermal timescale mass transfer, the bulk of the mass transfer should have been at a slower rate and therefore could have been mostly conservative. For purposes of simplicity and for presenting illustrative examples, we therefore parameterize the evolution of KOI 1224 as having a constant mass retention fraction, . We further take the specific angular momentum of any ejected material to be in units of that of the binary system.
Figure 6 shows how the period of such a binary evolves as mass is transferred from the primary, for a range of plausible values of . In all cases, the orbit shrinks as mass transfer gets under way. However, for the orbital period eventually dramatically increases above the original orbital period. Whereas, for smaller values of () the orbit either does not grow or continues to shrink with mass transfer. Regardless of the evolutionary path, the orbit cannot get close enough so that the accreting secondary overflows its Roche-lobe, otherwise the two stars would be in danger of merging.
The range of possible progenitor masses is explored in Figure 7. The shaded region indicates the most plausible pairs of progenitor masses, and . The constraints leading to this region are that (i) the ratio of nuclear and thermal timescales of the primordial pair lies between 0.5 and 0.9; (ii) the initial orbital period was shorter than 10 days; and (iii) the primary expanded by at least 30% in radius by the time that Roche lobe overflow commenced. The red curves are contours of constant orbital period at the start of mass transfer, while the blue curves are contours of constant mass retention fraction, . We again assumed a single value for the angular momentum loss parameter of .
The rationale for constraint (i) above is that the secondary should have a nuclear evolution time that is not too much longer than that of the primary in order for it to have evolved off the main sequence only after the mass transfer was completed, but before the white dwarf had a chance to cool below 14,000 K. The requirement (ii) for an initially short orbital period ( days) is to avoid developing a white dwarf mass that is too large (i.e., remaining below , our upper limit on ) to allow its progenitor to fit within its orbit at the time mass transfer commenced (see the discussion below). Finally, constraint (iii) allows room in the initial binary for the primary star to develop a He core before thermal timescale mass transfer reduces the primary’s mass below 1 and impedes its further evolution.
Finally, in regard to the evolutionary history of KOI 1224, we note that there is a theoretical relationship between the mass of the white dwarf and the orbital period in systems where the progenitor of the white dwarf was a low-mass star (i.e., ; see, e.g., Rappaport et al., 1995; Lin et al., 2011) and where the mass transfer was stable. This is conveniently summarized in Figure 5 of Lin et al. (2011) and their equation (1). In their Figure 5 we can see that stars with initial mass (the red and green points in that plot) produce a steep functional relation between orbital period and the final core mass in the evolution of binary systems involving neutron star accretors. This same relation also somewhat holds for donor stars of initial mass up to 4 that commence mass transfer in so-called “late case A” or “case AB” mass transfer, i.e., when the donor is near the end of its main-sequence phase. The fact that the accretors in the current problem are main sequence stars does not alter the core mass–radius relation on which this effect is based.
If we utilize the expression relating orbital period to core mass:
(Lin et al., 2011) we can infer that a final orbital period of 2.7 days should correspond to a white dwarf mass of with a spread in the theoretical values of . This is highly consistent with the mass of the white dwarf we find for KOI 1224. One caveat, however, is that the relation is based on the core mass–radius relation for giants (see, e.g. Han et al., 1994; Rappaport et al., 1995). In order for the latter to hold, the mass transfer must be effected in at least a semi-controlled way (see §6).
In this paper we presented KOI 1224 (KIC 6606653), a new Kepler eclipsing system consisting of a bloated white dwarf having a normal, slightly evolved companion star in a compact -d orbit. Archival analysis of the optical and near-infrared photometry, as well as inspection of a preliminary classification spectrum (L. Nelson, private communication), reveal that the primary is likely an F6 spectral class star, while the GALEX UV flux allows us to infer a white dwarf temperature of , a value that is comparable with that found by our numerical modelling when using the primary’s temperature as a prior.
We modelled the light curve of this binary with Icarus, a new light curve synthesis code, which uses proper atmosphere models and a physical description of the binary parameters that allows us to reproduce various effects visible in the light curve such as irradiation, ellipsoidal light variations, and Doppler boosting. The derived binary parameters agree relatively well with those inferred from a harmonic analysis of the Kepler data using semi-analytic approximations of the above effects.
From our light curve fitting, we found masses of and for the F-star primary and the white dwarf secondary, respectively, in a 2.698-day orbit. These values pose a modest challenge in terms of binary evolution. Indeed, the low mass of the white dwarf likely excludes a common envelope scenario, which would require a more evolved progenitor – and therefore a more massive white dwarf – in order to produce a compact system with the observed parameters. At the same time, the orbital separation of a system like KOI 1224 is expected to increase significantly in a stable Roche-lobe overflow evolution unless the fraction of mass accreted by the current primary remains relatively low, . In any case, it appears that the system had to start its mass transfer evolution at an orbital period that is not longer than 10 days (as discussed above).
KOI 1224 is the 4th such system found with Kepler. The orbital periods are: 2.7 d (KOI 1224), 3.3 d (KHWD3), 5.2 d (KOI 74), and 24 d (KOI 81) (current work; Carter et al. 2011; van Kerkwijk et al. 2010). The corresponding white dwarf masses are 0.21, 0.26, 0.21, and 0.3 , respectively, with uncertainties. Unfortunately, the theoretical relationship between orbital period and white dwarf mass (discussed above) is too steep, i.e., , where , to test effectively on such a relatively narrow range of and .
Preferably, one would like to find more such systems in wider orbital periods in the Kepler data where (i) the white dwarf masses are expected to deviate more substantially from 0.2 , and (ii) the relation should work more robustly because of the higher likelihood of purely stable mass transfer. In this regard, we note that Maxted et al. (2011) have recently discovered a possibly related system with an orbital period of only 0.688 days, and a highly bloated white dwarf of , as well as several other similar very short period systems (P. Maxted, private communication). The relation is difficult to quantify theoretically at these very short periods (of 1 day) because it is not yet clear whether these systems result from stable mass transfer, a common envelope scenario, or some less catastrophic form of unstable mass transfer. The relation works only when mass is removed from the progenitor of the white dwarf in a reasonably controlled way such that the core-mass–radius relation, on which it is based, remains valid.
- Allard et al. (2007) Allard, F., Allard, N. F., Homeier, D., Kielkopf, J., McCaughrean, M. J., & Spiegelman, F. 2007, A&A, 474, L21
- Allard et al. (2003) Allard, F., Guillot, T., Ludwig, H.-G., Hauschildt, P. H., Schweitzer, A., Alexander, D. R., & Ferguson, J. W. 2003, in IAU Symposium, Vol. 211, Brown Dwarfs, ed. E. Martín, 325–+
- Allard et al. (2010) Allard, F., Homeier, D., & Freytag, B. 2010, ArXiv e-prints
- Avni & Bahcall (1975) Avni, Y., & Bahcall, J. N. 1975, ApJ, 197, 675
- Beech (1985) Beech, M. 1985, Ap&SS, 117, 69
- Bianchi et al. (2007) Bianchi, L., et al. 2007, ApJS, 173, 659
- Borucki et al. (2011) Borucki, W. J., et al. 2011, ApJ, 736, 19
- Brown et al. (2011) Brown, T. M., Latham, D. W., Everett, M. E., & Esquerdo, G. A. 2011, arXiv:1102.0342
- Carter et al. (2011) Carter, J. A., Rappaport, S., & Fabrycky, D. 2011, ApJ, 728, 139
- Che et al. (2011) Che, X., et al. 2011, ApJ, 732, 68
- Covey et al. (2007) Covey, K. R., et al. 2007, AJ, 134, 2398
- Cox (2000) Cox, A. N. 2000, Allen’s astrophysical quantities, ed. Cox, A. N.
- de Kool (1990) de Kool, M. 1990, ApJ, 358, 189
- Dewi & Tauris (2000) Dewi, J. D. M., & Tauris, T. M. 2000, A&A, 360, 1043
- Di Stefano (2011) Di Stefano, R. 2011, AJ, 141, 142
- Driebe et al. (1999) Driebe, T., Blöcker, T., Schönberner, D., & Herwig, F. 1999, A&A, 350, 89
- Farmer & Agol (2003) Farmer, A. J., & Agol, E. 2003, ApJ, 592, 1151
- Fukugita et al. (1996) Fukugita, M., Ichikawa, T., Gunn, J. E., Doi, M., Shimasaku, K., & Schneider, D. P. 1996, AJ, 111, 1748
- Gelman & Rubin (1992) Gelman, A., & Rubin, D. B. 1992, Statistical Science, 7, 457
- Gies et al. (2008) Gies, D. R., et al. 2008, ApJ, 682, L117
- Gray et al. (2001) Gray, R. O., Graham, P. W., & Hoyt, S. R. 2001, AJ, 121, 2159
- Gregory (2005) Gregory, P. C. 2005, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with ‘Mathematica’ Support, ed. Gregory, P. C. (Cambridge University Press)
- Han et al. (1994) Han, Z., Podsiadlowski, P., & Eggleton, P. P. 1994, MNRAS, 270, 121
- Hendry & Mochnacki (1992) Hendry, P. D., & Mochnacki, S. W. 1992, ApJ, 388, 603
- Kippenhahn & Meyer-Hofmeister (1977) Kippenhahn, R., & Meyer-Hofmeister, E. 1977, A&A, 54, 539
- Kopal (1959) Kopal, Z. 1959, Close binary systems, ed. Kopal, Z.
- Lin et al. (2011) Lin, J., Rappaport, S., Podsiadlowski, P., Nelson, L., Paxton, B., & Todorov, P. 2011, ApJ, 732, 70
- Loeb & Gaudi (2003) Loeb, A., & Gaudi, B. S. 2003, ApJ, 588, L117
- Lucy (1967) Lucy, L. B. 1967, ZAp, 65, 89
- Maxted et al. (2011) Maxted, P. F. L., et al. 2011, ArXiv e-prints 1107.4986v2
- McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22
- Morris (1985) Morris, S. L. 1985, ApJ, 295, 143
- Neckel (2005) Neckel, H. 2005, Sol. Phys., 229, 13
- Nelson et al. (2004) Nelson, L. A., Dubeau, E., & MacCannell, K. A. 2004, ApJ, 616, 1124
- Orosz & Hauschildt (2000) Orosz, J. A., & Hauschildt, P. H. 2000, A&A, 364, 265
- Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F. 2011, ApJS, 192, 3
- Podsiadlowski et al. (2003) Podsiadlowski, P., Rappaport, S., & Han, Z. 2003, MNRAS, 341, 385
- Podsiadlowski et al. (2002) Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107
- Prša & Zwitter (2005) Prša, A., & Zwitter, T. 2005, ApJ, 628, 426
- Prša et al. (2011) Prša, A., et al. 2011, AJ, 141, 83
- Rappaport et al. (2009) Rappaport, S., Podsiadlowski, P., & Horev, I. 2009, ApJ, 698, 666
- Rappaport et al. (1995) Rappaport, S., Podsiadlowski, P., Joss, P. C., Di Stefano, R., & Han, Z. 1995, MNRAS, 273, 731
- Rey et al. (2007) Rey, S.-C., et al. 2007, ApJS, 173, 643
- Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15
- Rowe et al. (2010) Rowe, J. F., et al. 2010, ApJ, 713, L150
- Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
- Shporer et al. (2011) Shporer, A., Brown, T., Mazeh, T., & Zucker, S. 2011, ArXiv e-prints
- Sing (2010) Sing, D. K. 2010, A&A, 510, A21+
- Skrutskie et al. (2006) Skrutskie, M. F., et al. 2006, AJ, 131, 1163
- Tauris & Dewi (2001) Tauris, T. M., & Dewi, J. D. M. 2001, A&A, 369, 170
- van Kerkwijk et al. (2010) van Kerkwijk, M. H., Rappaport, S. A., Breton, R. P., Justham, S., Podsiadlowski, P., & Han, Z. 2010, ApJ, 715, 51
- van Rensbergen et al. (2010) van Rensbergen, W., De Greve, J. P., Mennekens, N., Jansen, K., & De Loore, C. 2010, A&A, 510, A13+
- van Rensbergen et al. (2011) van Rensbergen, W., de Greve, J. P., Mennekens, N., Jansen, K., & de Loore, C. 2011, A&A, 528, A16+
- Vennes et al. (2011) Vennes, S., Kawka, A., & Németh, P. 2011, MNRAS, 410, 2095
- von Zeipel (1924) von Zeipel, H. 1924, MNRAS, 84, 665
- Walsh (2004) Walsh, B. 2004, Markov Chain Monte Carlo and Gibbs Sampling
- Wilson & Devinney (1971) Wilson, R. E., & Devinney, E. J. 1971, ApJ, 166, 605
Appendix A Icarus: Binary Light Curve Synthesis Code
a.1. Stellar Surface Grid
The underlying stellar surface grid is constructed using the triangular tessellation of a unit sphere that was inspired by the GDDSYN synthesis code (Hendry & Mochnacki, 1992). The base vertices are initially generated from the primitives of an icosahedron. Each triangular surface is then sub-divided into four triangles by creating a new vertex at the midpoint of each edge. The new vertices are finally projected back onto the unit sphere by re-normalizing them to unity. This sub-division process can be repeated iteratively in order to obtain the desired number of surface elements. At a sub-division level (where is the base icosahedron, and is one sub-division further), the number of surface elements is , while the number of independent vertices is . The coordinates of each face is calculated as the centroid of the triangle.
There are several advantages to using the algorithm above. First, a triangular tessellation does not suffer from the uneven surface sampling problem common to simple meshes derived from an equally spaced grid in spherical coordinates. To circumvent this problem, other synthesis codes, like PHOEBE, use a variable sampling strategy with fewer points near the poles. The subdivision algorithm of a triangular tessellation is relatively easy to implement and, using an icosahedron as the primitive, yields triangular faces that are organized on a regular tile of hexagons and pentagons, and that have fairly equal areas. Triangular faces are by definition convex and hence are simpler to work with when it comes to dealing with transits and eclipses. For reasons that will become obvious later we also keep track of the associativities: each surface is associated to three vertices, and each vertex is associated to either 5 or 6 faces.
a.2. Equipotential Surface
As in the ELC code, our stellar surface is defined by the gravitational equipotential equation from Avni & Bahcall (1975), which takes into account the effects of stellar rotation:
where and is the distance of a point measured from the barycenter of the star that is being modelled (labeled 1) and from the barycenter of its companion (labeled 2) in units of orbital separation, is the mass ratio, and is the co-rotation factor expressed as the ratio of the stellar to orbital frequency. We work in a coordinate system having its origin at the barycenter of star 1, in which the -axis points towards star 2’s barycenter, the -axis is along the orbital angular momentum, and the -axis is along the orbital plane orthogonal to the and axes. Our code is currently limited to circular orbits but could easily be extended to elliptical orbits as in PHOEBE, for example.
Initially, the position of the L1 point, , along the x-axis is found by identifying the saddle-point of the above equation between stars 1 and 2. Then, as explained in Orosz & Hauschildt (2000), the potential corresponding to the “nose” of the star (the inner point towards star 2) is defined by specifying a filling factor , which is expressed as a fraction of the L1 distance. From there, the radius of star 1 can be evaluated in the direction of every vertex and face of the tessellated surface.
The effective surface gravity is also calculated for each face using the equation as well as the components of the normal to the local surface.
a.3. Surface Temperature and Gravity Darkening
We assign a fiducial temperature to each face that takes into account the gravity darkening according to the equation:
where typically varies between 0.08 and 0.25 depending on whether the star has a convective or a radiative envelope, respectively (Lucy 1967; von Zeipel 1924; for empirical constraints, see Che et al. 2011).
a.4. Irradiation from a Companion
We treat the effect of irradiation from a companion as a source of energy that increases the temperature of the star while preserving thermal equilibrium (i.e., no thermal inversion that gives rises to emission lines). In this case, irradiation changes the temperature distribution of the star as following:
(Beech, 1985) where is the angle between the surface normal and the direction to the companion (i.e., star 2), and is the irradiation temperature which describes the temperature increase. Note that working with is more convenient than using an irradiation luminosity and a bolometric albedo since empirically one effectively measures the back and front temperatures of a star using multi-band light curves. One can convert the irradiation temperature into an irradiation luminosity after the fact using , where is the Stefan-Boltzmann constant and the orbital separation, and work out the stellar albedo if the companion’s luminosity is known such that . If one prefers, it is also possible to fix an albedo and perform the calculation the other way around in order to provide an irradiation temperature to the code.
Our treatment of irradiation is different than that of PHOEBE and GDDSYN, which calculate the reflection effect of each surface element due to the other star in an iterative way. While their technique is certainly more accurate, it comes at the cost of a substantially larger computational burden and only becomes more relevant in the case of contact binaries, where one needs to account for shadowing effects due to the bridge of matter in order to be self-consistent. Moreover, the fact that heat might be partly redistributed over the surface in a way that is non-trivial to calculate from first principles helps justify our use of this approximation.
a.5. Surface Area
We pre-calculate the area of each surface triangle on the unit sphere. Once the equipotential equation has been solved we rescale each pre-calculated value by its proper radius in order to obtain the effective area of the surface elements. By this means, the integrated flux simply corresponds to a discrete summation over the visible surface.
a.6. Atmosphere Models
Once the stellar grid has been constructed, and the effective temperature and surface gravity determined for all the faces, we can evaluate the flux perceived by an observer located at a given orbital inclination and orbital phase. The backend of the code that returns the flux for each face can be swapped between a blackbody or an atmosphere grid.
For the purpose of this paper, we have worked with the BTSettl atmosphere models of Allard et al. (2003, 2007, 2010) which are available online on the Phoenix web simulator121212Available at http://phoenix.ens-lyon.fr/simulator/. We used these atmosphere models to cover a grid over the range 1000-15000 K in temperature and 3.0-4.5 in . Since only integrated spectra were available, we used the empirical limb darkening relationship of Neckel (2005).
Prior to performing the modelling of the KOI 1224 data, we integrated the spectral grid over the Kepler passband. Such integration can be performed for any photometric filter that one wishes to use to model the light curves. However, our code can also handle working with a full set of spectral data. In which case it is possible to model full orbital-resolved spectroscopic data. When this is the case, the spectrum of each face is Doppler shifted at the appropriate velocity and hence the rotational broadening, orbital Doppler shift, and displacement of the light-center away from the barycenter due to irradiation naturally arise from our model spectra. In principle, the Rossiter-McLaughlin effect (Rossiter, 1924; McLaughlin, 1924) (as well as its Doppler boosting analogue (van Kerkwijk et al., 2010; Shporer et al., 2011)) would also appear from spectra modelled in combination with the transit calculation described below.
a.7. Doppler Boosting
Doppler boosting is another feature implemented in our code. For most cases, the amplitude of the effect is negligible, and hence it is turned off to speed up the light curve calculation. For the KOI 1224 data, however, Doppler boosting is clearly detected and hence it was included in the following way.
As described in van Kerkwijk et al. (2010), Doppler boosting can be written as a modulation of the observed flux:
where is the orbital phase measured at the inferior conjunction of the star, is the orbital velocity, is the speed of light, and is a factor that depends on the spectrum of the source and the observing passband. Based on our analytical estimate of the primary’s temperature and surface gravity (see §3.1), we evaluated the coefficient at temperatures of 6200 and 6500 K for of 3.5 and 4.0. The coefficient is a rather smooth function of the temperature of surface gravity, and hence we simply performed a bilinear interpolation to fetch the coefficient to apply to each face of the stellar grid before the total flux was integrated. In the future, we are planning to include the exact calculation at each temperature and value.
a.8. Eclipse and Transit Calculation
Eclipse and transit calculation are also included in our synthesis code. Each time that the stellar surface grid is recomputed with new input parameters, we perform several tests to evaluate whether partial or total occultations are possible as well as the orbital phase range within which they can occur. This allows for our code to revert to algorithms optimized for each situation.
At an orbital phase and an inclination , the projected distance between the barycenter of the stars is given by , where at the inferior conjunction of the primary.
First, we test whether any kind of occultation can occur. We consider the maximum dimension of each star in order to obtain a conservative estimate. Hence, for any orbital range where , one star might be occultated by the other. If this is the case, a second test aims to check if a full occultation occurs, and, if so, what orbital range does it cover. If , star 1 will definitely get fully eclipsed in the orbital range that meets the requirement . The same test is also performed with the indices reversed.
With these criteria identified, the flux calculation can be classified into three categories: (1) no occultation, (2) partial occultation, and (3) full occultation. Case 3 is the simplest as the star does not contribute to the total observed flux and hence no calculation is required. Case 1 is also relatively simple as every face that is not beyond the stellar limb will contribute; a standard flux integration over the visible surface of each star is performed. Case 2 is more complex and we explain our algorithm below.
a.8.1 Partial Visibility
Evaluating the partial visibility of a stellar surface occultated by an other is a computationally intensive task that can be tackled in several ways. Our algorithm is meant to provide a balance between accuracy, execution speed, and ease of implementation.
In short summary, our algorithm aims to calculate the fraction of a surface element that is visible to an observer. The first step consists in the determination of the outline of the star located in front (here we choose star 2 for the explanation). We make the approximation that the orbital range within which the occultation of star 1 occurs is small enough that the outline at the onset of the occultation is the same as at conjunction. This simplifies the calculation considerably since, when the orbital inclination is sufficiently close to edge-on131313A more accurate outline could be calculated but our approximation is good enough in the case of KOI 1224, the limb of star 2 corresponds to its outline in the plane at , which is easy to compute. In this fashion, we obtain a lookup table of values, where and are the polar coordinates of the limb of star 2 in the sky plane calculated using its barycenter as the origin.
In the second phase, we compute the projected position of star 1’s vertices, and , with respect to star 2’s barycenter. If a vertex of star 1 lies within the limb of star 2 – that is – then it is hidden. To calculate the partial visibility of a surface element, we make the approximation the visible fraction of the surface is equal to the fraction of its vertices that are not hidden. Since each individual vertex is associated with either five or six faces, we can use the vertex-to-face associativity to assign the partial visibility weighting factor to all the faces without having to perform redundant calculations.
Our partial visibility algorithm relies on a simple projected distance determination for all the vertices of the visible face. When such a vertex is occulted, an extra weighting factor is added to five or six array elements corresponding to the associated surface elements. The computational cost of this technique remains small compared to a more accurate calculation such as that performed by GDDSYN. Moreover it can easily be parallelized.
The main caveat of our algorithm lies in the limited spatial resolution of the grid faces; the partial visibility goes in increments. Hence, the light curve can display important jitter if the chosen grid resolution is too coarse. A notable symptomatic case is that of stars having very unequal relative sizes and for which the projected area of the smaller star becomes comparable to the area of a surface elements of the larger star. The obvious solution relies on using a higher resolution for the larger star, though this comes at the expense of an increased computational burden.
To work around these problems, we have designed a method that takes advantage of the sub-division tessellation algorithm. Since the sub-division to a finer level adds vertices at the mid-points of the surface sides, it is easy to keep track of the associations of these vertices with the parent surface elements at the lower resolution. Using such a nested sub-division yields partial visibilities that can be computed to a much better precision than the original step-size, increasing to for one extra sub-division and for two. This algorithm comes with a relatively small computing footprint: since the surface properties (temperature and gravity) vary smoothly, it is sufficient to obtain an accurate partial visibility weighting for a low-resolution grid using the sub-division. In other words, there is little if no gain in calculating the emergent intensity for all the surface elements at the higher resolution.