Introduction
Extrasolar planets (or exoplanets) show a rich variety of orbital architectures, challenging the planet formation theory. Their properties force astrophysicists and astrobiologists to re-examine their assumptions regarding the growth and evolution of potentially habitable planets in the Milky Way.
Typically, the astrophysical ‘habitability’ of a world is determined by the radiative flux received from an external radiation source. This flux is used as input for climate models, which simulate the response of a planet of Earth mass – with similar composition and atmosphere – to this radiation. This allows the construction of a habitable zone (HZ), a region surrounding the radiation source that should be capable of hosting planets of the Earth mass and composition with liquid water on their surface.
This HZ concept pre-dates the detection of exoplanets by several decades (Huang Reference Huang1959; Hart Reference Hart1979), culminating in the seminal one-dimensional (1D) radiative transfer calculations of Kasting, Whitmire & Reynolds (Reference Kasting, Whitmire and Reynolds1993), which defined inner and outer radial boundaries of a spherically symmetric HZ, based on the star's luminosity and effective temperature. These calculations would form the basis for a number of parametrizations of the HZ (Underwood et al. Reference Underwood, Jones and Sleep2003; Selsis et al. Reference Selsis, Kasting, Levrard, Paillet, Ribas and Delfosse2007; Kaltenegger & Sasselov Reference Kaltenegger and Sasselov2011). Recently, Kopparapu et al. (Reference Kopparapu2013) returned to Kasting et al's original calculations, updating the atmospheric absorption models to refine the inner and outer boundaries of the HZ, extending them (in the conservative case) to 0.99 and 1.7 AU, respectively; although these boundaries apply only to Earth mass planets – increasing the mass tends to move the inner boundary closer to the star (Kopparapu et al. Reference Kopparapu, Ramirez, SchottelKotte, Kasting, Domagal-Goldman and Eymet2014).
The HZ model has been strongly challenged since the first detection of an exoplanet in a main sequence star system (Mayor & Queloz Reference Mayor and Queloz1995). For example, the orbital eccentricity of exoplanets is often much higher than exhibited in the Solar System, allowing planets to move in and out of the HZ over the course of an orbit. The habitability of a world then becomes a function of the length of time it spends inside the zone (Williams & Pollard Reference Williams and Pollard2002; Kane & Gelino Reference Kane and Gelino2012; Kane et al. Reference Kane, Ciardi, Gelino and von Braun2012).
The advent of the Kepler Space Telescope has allowed, for the first time, meaningful statistical discussion of how frequently stars in the Milky Way possess planets in the HZ (Dressing & Charbonneau Reference Dressing and Charbonneau2013; Kasting et al. Reference Kasting, Kopparapu, Ramirez and Harman2013; Petigura et al. Reference Petigura, Howard and Marcy2013). These analyses are generally sensitive to properties such as stellar mass, and the precise definition of the inner and outer HZ boundaries, which are in turn sensitive to properties such as the planet's surface liquid water fraction and atmospheric composition. Nevertheless, the above authors find the frequency of ‘Earth-like’ planets around other stars to range from about 15% to as high as 50%.
One of Kepler's achievements which deserves to note is the first detection of planetary systems orbiting binary main sequence stars. These possess planets in the orbit of the centre of mass of a binary system. They are commonly referred to as P-type or circumbinary systems, in contrast to S-type systems where planets orbit a single component of the binary. The first system to be detected was Kepler-16b (Doyle et al. Reference Doyle2011), followed quickly by Kepler-34b and Kepler-35b (Welsh et al. Reference Welsh2012) and two planets orbiting the star Kepler-47 (Orosz et al. Reference Orosz2012). Orosz et al. calculated that Kepler-47c was within the HZ of its multi-star system, although this initial calculation assumed that the radiative flux arrived from the system's centre of mass and not from the individual stars.
As the HZ is sensitive to the flux received at the planet's surface, the HZ of circumbinary systems can differ significantly from the equivalent single-star HZ, due to: (a) the extra source of radiation from the companion star, and (b) the extra effective sink of radiation that may come from eclipses of one star by the other. Several groups have investigated (a) in detail, revising the HZ parametrization of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993) for a multiple star system. Kane & Hinkel (Reference Kane and Hinkel2013) combined the radiation from both stars to approximate a single blackbody flux received at the planet, with the resulting effective temperature used to determine the standard HZ limits at all points in the system. Haghighipour & Kaltenegger (Reference Haghighipour and Kaltenegger2013) take a similar approach, using the previous HZ limits modified to account for a second star, where the flux from each star is spectrally weighted. Most recently, Cuntz (Reference Cuntz2014) presents a general analytic solution to solving for HZs in both S- and P-type binary systems. All approaches result in a HZ that deviates from the standard annular form for single-star systems, with the deviation being time-dependent, and a function of the stars’ orbit around the centre of mass. Forgan (Reference Forgan2014) used 1D energy balance modelling to determine the HZ as a function of planetary orbital parameters (semi-major axis, eccentricity). By carrying out simulations rather than relying on fitting functions, the effects of eclipses could be more naturally incorporated into the model. These varying approaches to HZ calculation agree that, for example, Kepler-47c would be considered habitable if it were not relatively massive – consequently, it may host habitable moons (cf. Quarles et al. Reference Quarles, Musielak and Cuntz2012; Forgan & Kipping Reference Forgan and Kipping2013; Heller & Zuluaga Reference Heller and Zuluaga2013; Hinkel & Kane Reference Hinkel and Kane2013; Kipping et al. Reference Kipping, Forgan, Hartman, Nesvorný, Bakos, Schmitt and Buchhave2013; Forgan & Yotov Reference Forgan and Yotov2014).
Generally speaking, these calculations have been 1D, tracking the radiative response of terrestrial atmospheres as a function of atmospheric depth or latitude. If one wishes to consider the potential nature of biomes on exoplanets, then it is important to measure both the longitudinal and latitudinal distributions of flux. Even without the presence of a second star, the relationship between a planet's orbital period and rotation period can prove crucial in the distribution of flux, depending on the planet's obliquity and eccentricity (Dobrovolskis Reference Dobrovolskis2007, Reference Dobrovolskis2009, Reference Dobrovolskis2013), and must be considered if we are to investigate for example where photosynthetic organisms might reside (Brown et al. Reference Brown, Mead, Forgan, Raven and Cockell2014).
Photosynthetic organisms present important biosignatures in exoplanet atmospheres (Wolstencroft & Raven Reference Wolstencroft and Raven2002; Raven & Cockell Reference Raven and Cockell2006; Seager et al. Reference Seager, Bains and Hu2013a, Reference Seager, Bains and Hub) that are potentially detectable by future spectroscopy missions (although see Livengood et al. Reference Livengood2011; Brandt & Spiegel Reference Brandt and Spiegel2014; Rein et al. Reference Rein, Fujii and Spiegel2014 for the challenges involved in such an endeavour). The distribution of photosynthetic potential on an exoplanet's surface will be an important model parameter for analysing future spectroscopic data of Earth-like planets for evidence of biological activity.
In this paper, we investigate both the latitudinal and longitudinal distributions of flux on the surface of an Earth-like planet in a circumbinary system. We use orbital parameters from Kepler-16 and Kepler-47 as an illustration of the flux patterns already present in the HZ of known exoplanet systems. We also investigate the time spent in darkness as a function of latitude and longitude, which will prove to be important when we consider attributes of photosynthetic organisms, such as photoperiodism.
In the ‘Method’ section, we outline the method by which we obtain two-dimensional (2D) surface flux and darkness patterns; in the ‘Results’ section, we display results obtained for putative Earth-like planets orbiting Kepler-16 and Kepler-47. In the ‘Discussion’ section, we discuss the potential for photosynthesis given the calculated flux and darkness patterns, and the dependence of these patterns on factors such as the binary's orbital parameters, and the relative phase of the planetary and binary orbits.
Method
System setup
The system is composed of two stars, with masses M 1 and M 2 and luminosities L 1 and L 2, orbiting each other with semi-major axis a bin and eccentricity e bin, with orbital period P orb,bin. The planet (with mass M p ) performs a Keplerian orbit about the centre of mass of the system, with semi-major axis a p and eccentricity e p , and orbital period P orb,p . The orbit may be inclined by an angle i p to the binary plane, and the obliquity of the planet's rotation relative to the binary orbital plane is δ p . At any time during the simulation, the positions of the stars and planets are given by r i (t), where i = 1,2, p, and the relative vector
$${{\bf r}_{ij}} = {{\bf r}_i} - {{\bf r}_j}.$$
The surface of the planet is divided into cells of longitude Λ and latitude Φ. We define the longitude of noon Λnoon,i as the longitude on the planet's surface at which star i is at its maximum height, where i = 1,2 (see Figure 1):
$${{\rm \Lambda} _{{\rm noon},i}} = {\tan ^{ - 1}}\left( {\displaystyle{{{{\bf r}_{{\bf pi}}}.{\hat {\bf y}}} \over {{{\bf r}_{{\bf pi}}}.{ \hat {\bf x}}}}} \right).$$
The planet rotates with a period P spin, and hence the apparent longitude of a fixed location on the planet surface (for a given star) varies as
$${{\rm \Lambda} _s}(t) = {\rm \Lambda} (t = 0) - {{\rm \Lambda} _{{\rm noon},i}} + \displaystyle{{2\pi t} \over {{P_{{\rm spin}}}}}.$$
We fix P spin = 24 h for all runs. Resonances between the spin and orbital periods are clearly extremely influential on the flux distribution, and both the dynamics and habitability of planets in spin–orbit resonances around single stars have been investigated in detail (Dobrovolskis Reference Dobrovolskis2007, Reference Dobrovolskis2009; Brown et al. Reference Brown, Mead, Forgan, Raven and Cockell2014). However, it is unclear what spin–orbit resonances are available to a circumbinary planetary system, if at all (see the Discussion section).
Calculating the received flux
The Cartesian surface normal vector for a given longitude/latitude cell, n s (Λ i (t),Φ), is constructed in the customary fashion using spherical polar co-ordinates:
$$\eqalign{& {n_{s,x}} = \sin {\rm \Phi} \cos {\rm \;} {{\rm \Lambda} _i}, \cr & {n_{s,y}} = \sin {\rm \Phi} \sin {\rm \;} {{\rm \Lambda} _i}, \cr & {n_{s,z}} = \cos {\rm \Phi,}} $$
and then rotated along the x-axis to generate the correct obliquity:
$${\bf n} \to {{\bf R}_x}{\bf n}{\rm,} $$
where R x is the standard x-axis rotation matrix. The flux from star s hitting this latitude/longitude cell is:
$${F_s} = \displaystyle{{{L_s}} \over {4\pi {{\vert {{{\bf r}_{{\bf ps}}}} \vert}^2}}}{{\hat {\bf r}}_{{\bf ps}}}{\rm.} {{\bf n}_{\bf s}},\,\,{\rm for}\, \,{{\hat {\bf r}}_{{\bf ps}}}{\rm.} {{\bf n}_{\bf s}} \gt 0.$$
If
${\hat {\bf r}_{{\bf ps}}}$
.n
s
⩽ 0, then the latitude/longitude cell is facing away from the star, and should therefore be in its night side.
Calculating eclipses and periods of darkness
The total flux received by a latitude/longitude cell is simply
$${F_{{\rm tot}}}({\rm \Lambda}, {\rm \Phi}, t) = {F_1}({\rm \Lambda}, {\rm \Phi}, t) + {F_2}({\rm \Lambda}, {\rm \Phi}, t).$$
If one star eclipses the other, then the received flux must be modified. If star i eclipses star j, then star j has its flux reduced by a factor A ij , which is the area of intersection of two circles corresponding to the stars i and j.
We measure the flux received at every longitude/latitude as a function of time, as well as averaging the total flux received over the course of the simulation:
$$\bar F({\rm \Lambda}, {\rm \Phi} ) = \displaystyle{{\mathop \int \nolimits^{{}} {F_{{\rm tot}}}(\Lambda, \Phi, t^{\prime})\,dt^{\prime}} \over {\mathop \int \nolimits^{{}}{dt^{\prime}}}}.$$
This integral is evaluated as a sum with dt→Δt = 30 min. We also integrate the total time spent in darkness D int over the course of the simulation. We define darkness at any latitude/longitude as any time interval over which neither star is visible (i.e. F tot(t) = 0). We shall see that this quantity is not as straightforward as it is in single-star systems.
Results
Kepler-16b
The Kepler-16 system consists of two relatively low-mass stars, of masses M 1 = 0.6897 M ⊙ and M 2 = 0.2026 M ⊙, orbiting each other with a semi-major axis of a bin = 0.224 AU and eccentricity e bin = 0.15944. These parameters were taken from Doyle et al. (Reference Doyle2011), and we use their values for stellar radii and luminosity also.
Kepler-16b is a 0.3 M Jup exoplanet, with a semi-major axis of 0.7 AU and eccentricity 0.069, with an orbital inclination within 0.4° of the binary plane. Kepler-16 provides an example of a system with two low-mass stars that will eclipse each other on a relatively long timescale (as the orbital period of the binary is 41 days), and Kepler-16b's orbit is just outside the HZ (Haghighipour & Kaltenegger Reference Haghighipour and Kaltenegger2013; Kane & Hinkel Reference Hinkel and Kane2013; Forgan Reference Forgan2014). We will therefore replace Kepler-16b with an Earth-like planet to assess how flux is distributed in this particular scenario.
We also perform runs without the second star, where the planet orbits the primary with the same parameters. We will use these data to identify features that are unique to circumbinary systems. We note here that orbits at the semi-major axis of Kepler-16b are stable only for low eccentricities (see the Discussion section). Some of the patterns displayed here are therefore of limited relevance to Kepler-16b in particular, but we can consider them as indicative of patterns displayed in as yet undiscovered circumbinary planetary systems where the binary has a relatively small mass ratio, while permitting stable eccentric orbits.
Flux patterns
Figure 2 shows
$\bar F$
on an Earth-like planet with Kepler-16b's parameters, where we run the simulation for one orbital period. We fix the planet's semi-major axis at 0.7 AU, and vary the eccentricity between e
p
= 0, 0.2, 0.5, and the obliquity between δ
p
= 0°, 30°, 60°. The longitudes of periapsis for both the planet and binary orbits are aligned, which we note will have an effect on the resulting surface patterns.
The simulation setup. Stars 1 and 2 orbit the centre of mass with position vectors r 1 and r 2, and the planet orbits with position vector r p . For illustration, the longitude of noon for star 2 is also shown.
Flux patterns (averaged over the planetary orbital period) for a circumbinary planet in orbit around the Kepler-16 binary, with a p = 0.7 AU, and e p = 0 (top row), e p = 0.2 (middle row) and e p = 0.5 (bottom row). The three columns have obliquities δ p = 0°, 30°, 60°, respectively.
Generally speaking, the features produced by increasing eccentricity are similar to those in the single-star case: the flux peaks at around the equator, in this case at longitudes around 90°. However, in the single-star case this peak appears at the substellar point of the planet when it is at periastron. Circumbinary planets will present two substellar points at periastron – if the obliquity is small, the substellar points both appear on the equator, at slightly separated longitudes, resulting in the flux being smeared over larger angular area than in the single-star case.
If the obliquity is high, the flux is more uniformly distributed, as it is in the single-star case. However, the effect of eclipses becomes increasingly important. The orbital period of the planet is around 200 days, around 5 full orbits of the binary. Therefore, one half of the orbit will experience an extra eclipse period compared with the other, over the course of an orbital period. At high obliquity, this means one polar region will receive less flux than the other – in this case, the North pole receives less flux than the South. This is a function of the relative phases of the binary orbit and the planetary orbit – in future years this pattern will be reversed, and the North will receive more flux than the South.
This may seem surprising initially, as a back-of-the-envelope calculation would suggest that the eclipse timescale is rather short:
$${t_{{\rm eclipse}}} = \displaystyle{{{R_1} + {R_2}} \over {{v_{{\rm orbit}}}}} \approx 3\,{\rm h},$$
where R 1 and R 2 are the stars’ physical radii, and v orbit is the binary's orbital velocity. However, this equation is valid only for a static observer, and the planet moves at its own orbital velocity relative to the centre of mass of the system. Therefore the eclipse timescale should be written
$${t_{{\rm eclipse}}} = \displaystyle{{{R_1} + {R_2}} \over {{v_{{\rm orbit,relative}}}}},$$
where v orbit,relative now describes the orbital velocity of the binary in the frame where the planet is at rest. Figure 3 shows how the flux in a single longitude/latitude cell varies with time for a planet orbiting Kepler-16 at a p = 0.7 AU, with zero eccentricity and obliquity (the flux is measured at the planet's equator along the prime meridian).
The variation of received flux with time, for a planet orbiting the Kepler-16 binary at a p = 0.7 AU, with eccentricity e p = 0 and obliquity δ p = 0°. The flux is measured at the planet's equator, along the prime meridian. Curves are plotted for flux from the primary and secondary, as well as the total received flux.
We can see that the total flux varies according to several distinct timescales. Firstly, there is the standard day/night cycle, which operates at the highest frequency. As this example has e p = 0, δ p = 0°, this frequency is fixed. The motions of the binary around the centre of mass add an amplitude modulation to this signal, with a period equal to the binary's period. Finally, the effect of eclipses reduces or completely erases the flux at regular intervals – we can see that the eclipse timescale is indeed much longer than a few hours. As the secondary has a smaller radius than the primary, the primary is only partially eclipsed by the secondary, but the secondary is totally eclipsed by the primary.
Darkness patterns
For a planet in a circular orbit around a single star, the total amount of time spent in darkness, for any given longitude and latitude, will be half of the orbital periodFootnote 1 . If the planet's orbit is eccentric, then the changing orbital velocity will ensure that for non-zero obliquity, the planet's northern hemisphere can experience a longer winter than the southern hemisphere (for example).
If the planet is in a binary system, then for a given location on its surface to be in darkness, both stars must not be in the sky. Figure 4 shows the total time spent in darkness for planets with obliquity δ p = 30° orbiting the Kepler-16 binary. We also run (but do not show here) simulations with the secondary removed, to assess our ability to measure uniform darkness times.
Darkness patterns for (left column) a circumbinary planet in orbit around the Kepler-16 binary, each with a p = 0.7 AU, obliquity δ p = 30°, and (top row) e p = 0, (middle row) e p = 0.2 and (bottom row) e p = 0.5.
In these single-star runs, for circular orbits, the darkness time in the single-star case is quite uniform, to within the numerical limits of the simulation (i.e. for Δt ≈ 30 min). The range in measured darkness time is a few hours for the circular orbit case (this value is so large due to difficulties measuring the darkness period at the polar circles of the planet). As the eccentricity increases, it is clear that some longitudes receive less daylight than others, and this is a function of orbital phase. The limits introduced by a finite time step introduces problems resolving along the latitude of the polar circles – however, it is clear that for a single-star system there is no evidence of any latitudinal dependence.
In the two-star case (Fig. 4), the difference is stark. The overall time spent in darkness is shorter – this is a natural consequence of having a second radiation source. However, we can see a clear latitudinal differentiation, demarcated by the polar circles. These are regions where darkness can occur in protracted, uninterrupted intervals as well as part of the night/day cycle. As the planet orbits the centre of mass of the system, the motion of the stars can be sufficient to hide them from view when the higher latitudes would normally begin their exit from winter, if their radiation source resided at the local centre of mass. This introduces a slight delay on when daylight returns to these regions, adding extra darkness time. Even in the case of a circular orbit, the difference in darkness time between the polar and equator regions can be as large as 8 days!
Kepler-47c
Kepler-47 is composed of a G star and an M star, with masses M 1 = 1.043 M ⊙ and M 2 = 0.362 M ⊙, respectively. Their orbit has semi-major axis a bin = 0.0836 AU and eccentricity e bin = 0.0234. Other properties are taken from Orosz et al. (Reference Orosz2012).
Kepler-47 is the first circumbinary planetary system with multiple planet detections. At an orbital semi-major axis of 0.989 AU, with an eccentricity constrained to be <0.4, Kepler-47c is in the HZ of this system (Haghighipour & Kaltenegger Reference Haghighipour and Kaltenegger2013; Kane & Hinkel Reference Hinkel and Kane2013; Forgan Reference Forgan2014), whereas Kepler-47b is too close to the binary (and subsequently too hot).
Again, Kepler-47c is not an Earth-like planet, but for our purposes we will assume an Earth-like planet exists at the same semi-major axis as Kepler-47c. This binary system has a relatively large binary mass ratio, with extremely frequent eclipses (the binary orbital period is approximately 7 days).
Flux patterns
Figure 5 shows the time integrated flux for planets orbiting the Kepler-47c binary at a p = 0.989 AU. We again vary the eccentricity between e p = 0, 0.2, 0.5 (top, middle, bottom rows), and the obliquity between δ p = 0°, 30°, 60° (left, middle, right columns).
Flux patterns (averaged over the planetary orbital period) for a circumbinary planet in orbit around the Kepler-47 binary, with a p = 0.989 AU, and e p = 0 (top row), e p = 0.2 (middle row) and e p = 0.5 (bottom row). The three columns have obliquities δ p = 0°, 30°, 60°, respectively.
The low obliquity cases are very similar to those for Kepler-16b. Again, the peak flux is spread over a wider range of longitudes than would be seen for a single star – however, the relatively small a bin places the two substellar points closer together on the planet's surface, and so the smearing of flux is less evident. Also, the secondary luminosity is a good deal smaller than the primary luminosity, so it contributes much less to the total surface flux, relatively speaking.
For high obliquities, we do not see the poles receiving significantly different levels of flux. The eclipse timescale is significantly smaller than that of Kepler-16b – as a result, for circular orbits the planet will experience around 40 primary eclipse events per orbit (see Fig. 6). Consequently, both the northern and southern hemispheres will experience very similar numbers of eclipses during their summer periods. Also, the large mass ratio of the binary ensures that eclipses of the secondary by the primary have an almost negligible effect on the incoming flux.
The variation of received flux with time, for a planet orbiting the Kepler-47 binary at a p = 0.989 AU, with eccentricity e p = 0 and obliquity δ p = 0°. The flux is measured at the planet's equator, along the prime meridian. Curves are plotted for flux from the primary and secondary, as well as the total received flux.
Figure 6 demonstrates the minimal effect of the second star on the flux received at a given latitude/longitude. Eclipses of the primary cause a drop in flux of around 15%, but the event is brief, lasting approximately 3 h.
Darkness patterns
Figure 7 shows the time spent in darkness for a planet with obliquity δ p = 30°, orbiting the Kepler-47 binary. Again, we see that the single-star case (not shown here) has relatively uniform darkness time. The simulation timestep is kept the same as for Kepler-16, so with Kepler-47c's increased orbital period, the simulation can carry out significantly more timesteps, improving its ability to measure uniform darkness maps.
Darkness patterns for a circumbinary planet in orbit around the Kepler-47 binary, each with a p = 0.989 AU, obliquity δ p = 30°, and (top row) e p = 0, (middle row) e p = 0.2 and (bottom row) e p = 0.5.
A circular, zero obliquity orbit around the Kepler-47 binary has a range in darkness time of about 4 days, increasing to about 5 days as the eccentricity is increased. Again, different darkness regimes are separated by the polar circles. This is for the same reasons as seen for Kepler-16 – during prolonged darkness periods in the winter, the orbital motions of the stars introduces a slight delay into the first sunrise of the spring.
We see a slight longitudinal dependence appearing as e p is increased: comparing Λ = 120° to Λ = 300° shows a difference in darkness time of about 1 day. This is two to three times larger than the range seen in single-star cases where we expect uniform darkness maps, so we are confident that this effect is not merely numerical.
Discussion
The effect of orbital inclination
Figure 8 shows dependence of both flux and darkness patterns on orbital inclination for the Kepler-47 system. As a bin and e bin are both relatively low, the binary sufficiently resembles an extended single star that increasing i p is equivalent to increasing δ p , making the flux and darkness patterns difficult to distinguish from a simulation with high obliquity orbiting in the binary plane (i p = 0).
The effect of increasing orbital inclination. Flux patterns (left column) and darkness patterns (right column) for a circumbinary planet in orbit around the Kepler-47 binary with obliquity δ p = 30°, eccentricity e p = 0, and (top row) orbital inclination i p = 0, (middle row) i p = 15° and (bottom row) i p = 30°.
Kepler-16 presents a larger a bin and more equal mass ratio. Figure 9 shows that there are some small differences in the presented flux and darkness features compared to Kepler-47. The flux is distributed over a larger latitudinal range for i p = 15°, and the eclipse timing effects which gave the South pole more flux for δ p = 60°, i p = 0° become apparent again for δ p = 30°, i p = 30°.
The effect of increasing orbital inclination. Flux patterns (left column) and darkness patterns (right column) for a circumbinary planet in orbit around the Kepler-16 binary with obliquity δ p = 30°, eccentricity e p = 0, and (top row) orbital inclination i p = 0, (middle row) i p = 15° and (bottom row) i p = 30°.
Sensitivity to the binary orbital phase
Many of the flux patterns we have studied show strong longitudinal dependences, in particular where the substellar points are during the planet's periastron passage. As P orb,bin and P orb,p are generally incommensurate, we should investigate how the relative orbital phases affect the resulting flux pattern. The simplest probe of this effect is to change the longitude of periastron of the planet's orbit, set at 0° throughout. As planet periastron occurs at t = 0, the relative orbital phase is changed as a result.
Figure 10 displays the flux patterns observed in the Kepler-16 system (for e p = 0.5, δ p = 60°), as the longitude of periapsis is increased. As the relative orbital phase changes, the strong peak in flux at southern latitudes weakens, and finally disappears. This sensitivity to phase will have important consequences for fluxes averaged over more than one orbit, as we will see in the following section.
The effect of relative orbital phases. Flux patterns for a circumbinary planet in orbit around the Kepler-16 binary with obliquity δ p = 60°, eccentricity e p = 0.5, and (top) longitude of periapsis of 0°, (middle row) longitude of periapsis of 5° and (bottom) longitude of periapsis of 45°.
The persistence of patterns over many orbits
The patterns seen thus far relate to one full orbit of the planet around the centre of mass. As we have seen that the underlying cause of these patterns are in the relationship between the planet's orbital phase, the binary's orbital phase and their other respective orbital parameters, will such unusual patterns persist over many orbits, or will they ‘average out’ into a morphology more consistent with a single-star system?
Generally speaking, the flux patterns average out to a more uniform distribution, as the difference between the binary and planet's orbital longitudes precesses, allowing the stars’ substellar points to move across all longitudes. The darkness patterns, however, remain distinct. Figure 11 compares the darkness patterns over 1 planetary orbital period to 100 planetary orbital periods for two cases. The first is the Kepler-16b system, with e p = 0.2 and δ p = 30°; the second is the Kepler-47c system with e p = 0.5, δ p = 30°.
The persistence of patterns with time. Top row: Darkness patterns for the Kepler-16b system, with e p = 0.2 and δ p = 30°. Bottom row: Darkness patterns for the Kepler-47 system, with e p = 0.5 and δ p = 30°. The left column displays results from one orbital period of the planet; the right column displays results from 100 orbital periods of the planet.
Orbital stability and the potential for spin–orbit resonances
Until now, we have considered circumbinary planets with the same rotation period as the Earth, P spin = 24 h, but this is clearly not generally true, even for ‘Earth-like’ planets around single stars. As a consequence of total angular momentum conservation, tidal interactions between a host star and a planet can transfer spin angular momentum into orbital angular momentum. Typically, the long-term tidal evolution of exoplanets at close distances to their star results in circular orbits with synchronous rotation, where P spin÷ P orb,p , with obliquities generally locked into a Cassini state corresponding either to high or low values (Peale Reference Peale1969; Henrard & Murigande Reference Henrard and Murigande1987). If the initial orbit is eccentric, the planet can be captured into an m:n spin–orbit resonance, where mP spin÷ nP orb,p , with m, n integers.
Mercury is currently locked in a 3 : 2 spin–orbit resonance with the Sun, thanks to its modest eccentricity of 0.206. It can be shown that this capture was probable (Correia & Laskar Reference Correia and Laskar2004; Dobrovolskis Reference Dobrovolskis2007) and stable over long timescales. As such, spin–orbit resonances such as 3 : 2 may be a common feature of habitable planets around M stars. As has been previously mentioned, the effects of resonances between P spin and P orb,p on flux patterns have been considered in detail for single-star systems (Dobrovolskis Reference Dobrovolskis2007, Reference Dobrovolskis2009, Reference Dobrovolskis2013; Brown et al. Reference Brown, Mead, Forgan, Raven and Cockell2014).
To date, there has not been any in-depth dynamical study of spin–orbit resonances for circumbinary planets. Some authors have considered the tidal interactions between the binary stars, and the resulting impact on hazardous flare activity (Mason et al. Reference Mason, Zuluaga, Clark and Cuartas-Restrepo2013), but none consider the spin–orbit evolution of the planet itself. This may be a reasonable lapse, as the calculations for single stars suggest that spin–orbit evolution such as that seen for Mercury (and proposed for other exoplanets) only occurs if a p is sufficiently small.
In circumbinary planetary systems, there is a lower limit for a p , underneath which planetary orbits are no longer dynamically stable (Holman & Wiegert Reference Holman and Wiegert1999)
$$\eqalign{{a_p} \gt {a_{{\rm dyn}}} &= {a_{\rm bin}}\left( {1.6 + 5.1{e_{{\rm bin}}} + 4.12{\mu _{\rm bin}} - 2.22e_{\rm bin}^2 - 4.27{\mu _{{\rm bin}}}{e_{{\rm bin}}} }\right. \cr &\quad \left.{ - 5.09\mu _{{\rm bin}}^2 + 4.61\mu _{{\rm bin}}^2 e_{{\rm bin}}^2} \right),}$$
where
${\mu _{{\rm bin}}} = \displaystyle{{{M_2}} \over {{M_1} + {M_2}}}$
. The Kepler-47 system can be thought of as a solar-type system perturbed by the presence of the M-type companion, so it may be reasonable to assume that tidal interaction calculations made for the Solar System, while inaccurate for this case, will give a sense of broad trends. In the absence of a detailed calculation of the combined tidal interaction between both stars and any planet in the system, we can compare a
dyn to the orbital semi-major axis of Mercury, which is 0.381 AU, and consider whether in the first instance Mercury's orbit would be stable, and leave the details of the interaction to later work. In the case of Kepler-47, a
dyn = 0.202 AU, so without the benefit of detailed calculation, we should consider the possibility that some planets could enter a spin–orbit resonance. However, as in the case of the Solar System, it is likely that planets susceptible to spin–orbit resonances in the Kepler-47 system would also lie outside the habitable zone.
For Kepler-16, a dyn = 0.645 AU, and most authors agree this means most of, if not all, of its HZ is orbitally unstable (Liu et al. Reference Liu, Zhang and Zhou2013). As we have already stated, eccentric orbits in the Kepler-16 system with Kepler-16b's semi-major axis are not stable, which we confirm via N-body simulation (Fig. 12). Given the low mass of both stars, it seems that planets cannot orbit sufficiently close to the binary to undergo significant tidal interactions, but once again we stress that this is at best an educated guess, and requires further work.
Orbital stability of planets orbiting the Kepler-16 binary under N-body integration. We show the semi-major axis (left), and eccentricity (right) for planets with initial eccentricities of 0 (red), 0.2 (green) and 0.5 (blue). The e = 0.5 planet is ejected from the system after ~2 years, whereas the e = 0.2 planet is ejected after several hundred years.
We should also note that circumbinary planets generally display rapid apsidal precession, behaviour that is confirmed both analytically and by N-body simulation by other authors (Doyle et al. Reference Doyle2011; Orosz et al. Reference Orosz2012; Leung & Lee Reference Leung and Lee2013). This precession will have important consequences for planetary seasons, and may also produce secular resonances, with corresponding effects on much longer timescales.
Prospects for photosynthesis
It seems clear that the patterns of both flux and darkness measured in this work could lead to patterns of rhythmicity that are unusual on the Earth (although we can appeal to some terrestrial analogues). The changing spectral quality reaching the planetary surface caused by the changing dominance of the radiation from the primary and secondary could potentially allow for a selection pressure that causes organisms to specialize in using particular regions of the spectrum, particularly as it shifts from short to longer wavelengths.
This spectral niche differentiation is seen on the Earth, across scales similar to that observed at different depths in the Earth's oceans (Raven Reference Raven2007; Stomp et al. Reference Stomp, Huisman, Stal and Matthijs2007). O'Malley-James et al. (Reference O'Malley-James, Raven, Cockell and Greaves2012) considered photosynthesis on an Earth-like planet orbiting in the HZ of a close binary system comprising a G and an M star. That paper emphasizes the possibility of a greater significance of photosynthetic organisms using radiation at wavelengths out to about 1000 nm than is the case on Earth, although the radiation environment of the Earth-like planet orbiting a G–M binary is dominated by the G star. Figure 13 shows the spectral flux density reaching the surface of a planet in the Kepler-16 system, orbiting at zero eccentricity and zero obliquity. In the left plot, the primary is close to its maximum altitude in the sky (primary midday), and in the right plot it is close to the horizon (primary sunset). Note that the secondary flux is increasing as the primary flux is decreasing, and exceeds the primary flux at infrared wavelengths during and after primary sunset. The effects of such changes in infrared flux on photosynthesis depends on the extent to which global photosynthesis depends on photon absorption in the 700–1000 nm range; this is in turn a function of the global occurrence of infrared-absorbers and the depth of infrared-absorbing water overlying them (Wolstencroft & Raven Reference Wolstencroft and Raven2002; Stomp et al. Reference Stomp, Huisman, Stal and Matthijs2007; O'Malley-James et al. Reference O'Malley-James, Raven, Cockell and Greaves2012).
Spectral variation over the course of a (primary) day in the Kepler-16 system. The left-hand plot shows the spectral flux density when the primary is at maximum altitude (i.e. midday), and the right-hand plot shows the same for when the primary is setting.
The changing intensity of light and darkness on the planets discussed is analogous to the onset of lunar day, in which an additional light cycle with a lower frequency is superposed on the diurnal day/night cycle. On the Earth, it has been found that animals can synchronize their activities to a great diversity of natural, geophysical rhythms, for example diurnal (24 h), tidal (12.4 h), semilunar (14.8 days), lunar (29.5 days), annual, seasonal or photoperiodic (365 days), as well as longer cycles involving prime numbers of years (13 or 17) of reproduction related to predator satiation of the hatching organisms (Sota et al. Reference Sota, Yamamoto, Cooley, Hill, Simon and Yoshimura2013). There is an evolutionary advantage of responding to stimuli which correlate with recurring environmental conditions. The lunar cycle is known to cause an enormous number of biological responses (Endres & Schad Reference Endres and Schad2002), for example the timing of breeding in marine organisms. This is despite the moon being relatively faint (some 400 000 times fainter than the Sun in the optical). It is true that some of these biological cycles are driven by tides, but we argue that biological rhythms driven by light cycles are likely to be present in organism behaviour of inhabited circumbinary planets, especially as the secondary light source is significantly brighter. Given that the periodicity of light cycles in circumbinary systems can be both greater than and less than the planet's rotational period, we should expect that in biospheres of circumbinary planets, selection pressures will exist for biological cycles tuned to a variety of different rhythms, with the principal cycles linked to the period of the binary.
Conclusions
We have calculated time-dependent surface maps of both flux and darkness on planets in circumbinary systems. Using Kepler-16 and Kepler-47 as archetypes, we model the flux arriving from both stars as a function of planet latitude and longitude, given an initial Keplerian orbit with fixed elements and a fixed spin obliquity.
We identify patterns in both the flux and darkness that are unique to binary star systems. These patterns exist as both changing spatial distributions and temporal fluctuations. The spectral quality of radiation that might be considered photosynthetically useful also varies with time and surface position. All the above variations have periods, amplitudes and phases that depend on the relative orbital phase between the stellar orbit and the planetary orbit, with timescales both much larger than the planetary orbital period and much smaller.
It is clear that if climate modelling of Earth-like circumbinary planets is to go beyond simple 1D approximations (Forgan Reference Forgan2014) towards more sophisticated 3D global circulation models (e.g. Shields et al. Reference Shields, Bitz, Meadows, Joshi and Robinson2014; Yang et al. Reference Yang, Boué, Fabrycky and Abbot2014) then these surface flux and darkness patterns will play a crucial role in determining both atmospheric physics and oceanic circulation (and the subsequent interactions thereof).
With such a wide variety of forcing timescales for photosynthesis, we conclude that any inhabited circumbinary planet will produce a biosphere rich in rhythms and cycles, determined by non-trivial relationships between the planet's orbit, the binary orbit and the location of biomes on the planetary surface.
Acknowledgements
DF and CC gratefully acknowledge support from STFC grant ST/J001422/1. AM acknowledges the support of a STFC studentship. The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. The University of Dundee is a charitable body registered in Scotland, with registration number SC 15096.
