Milankovitch Cycles of Terrestrial Planets in Binary Star Systems
Abstract
The habitability of planets in binary star systems depends not only on the radiation environment created by the two stars, but also on the perturbations to planetary orbits and rotation produced by the gravitational field of the binary and neighbouring planets. Habitable planets in binaries may therefore experience significant perturbations in orbit and spin. The direct effects of orbital resonances and secular evolution on the climate of binary planets remain largely unconsidered.
We present latitudinal energy balance modelling of exoplanet climates with direct coupling to an N Body integrator and an obliquity evolution model. This allows us to simultaneously investigate the thermal and dynamical evolution of planets orbiting binary stars, and discover gravitoclimatic oscillations on dynamical and secular timescales.
We investigate the Kepler47 and Alpha Centauri systems as archetypes of P and S type binary systems respectively. In the first case, Earthlike planets would experience rapid Milankovitch cycles (of order 1000 years) in eccentricity, obliquity and precession, inducing temperature oscillations of similar periods (modulated by other planets in the system). These secular temperature variations have amplitudes similar to those induced on the much shorter timescale of the binary period.
In the Alpha Centauri system, the influence of the secondary produces eccentricity variations on 15,000 year timescales. This produces climate oscillations of similar strength to the variation on the orbital timescale of the binary. Phase drifts between eccentricity and obliquity oscillations creates further cycles that are of order 100,000 years in duration, which are further modulated by neighbouring planets.
keywords:
astrobiology, methods:numerical, planets and satellites: general1 Introduction
Approximately half of all solar type stars reside in binary systems (Duquennoy & Mayor, 1991; Raghavan et al., 2010). Recent exoplanet detections have shown that planet formation in these systems is possible. Planets can orbit one of the stars in the socalled S type configuration, such as Cephei (Hatzes et al., 2003) HD41004b (Zucker et al., 2004) and GJ86b (Queloz et al., 2000). If the binary semimajor axis is sufficiently small, then the planet can orbit the system centre of mass in the circumbinary or P type configuration. Planets in this configuration were first detected around postmain sequence stars, in particular the binary pulsar B16026 (Thorsett et al., 1993; Sigurdsson et al., 2003). The Kepler space telescope has been pivotal in detecting circumbinary planets orbiting main sequence stars, such as Kepler16 (Doyle et al., 2011), Kepler34 and Kepler35 (Welsh et al., 2012), and Kepler47 (Orosz et al., 2012).
Planets in binary systems are sufficiently common that we should consider their habitability seriously. As of July 2016, 112 exoplanets have been detected in binary star systems
However, theoretical modelling indicates that the dynamical landscape of the binary significantly affects the planet formation process, both for Stype (Wiegert & Holman, 1997; Quintana et al., 2002, 2007; Thébault et al., 2008, 2009; Xie et al., 2010; Rafikov & Silsbee, 2014b, a) and Ptype systems (Doolin & Blundell, 2011; Rafikov, 2013; G. Martin et al., 2013; Marzari et al., 2013; Dunhill & Alexander, 2013; Meschiari, 2014; Silsbee & Rafikov, 2015). Therefore, when considering the prospects for habitable worlds in the Milky Way, one must take care to consider the effects that companion stars will have on the thermal and gravitational evolution of planets and moons.
The habitable zone (HZ) concept (Huang, 1959; Hart, 1979) is often employed to determine whether a detected exoplanet might be expected to be conducive to surface liquid water (that is, if its mass and atmospheric composition allow it). Initially calculated for the single star case using 1D radiative transfer modelling of the layers of an Earthlike atmosphere (Kasting et al., 1993), this quickly establishes a range of orbital distances that produce clement planetary conditions. Over time, line radiative transfer models have been refined, leading to improved estimates of the inner and outer habitable zone edges (Kopparapu et al., 2013; Kopparapu et al., 2014).
In the case of multiple star systems, the presence and motion of extra sources of gravity and radiation have two important effects:

The morphology and location of the system’s HZ changes with time, and

Regions of the system are orbitally unstable
These joint thermaldynamical constraints on habitability have been addressed in a largely decoupled fashion using a variety of analytical and numerical techniques.
The thermal time dependence of the HZ can be evaluated by combining the flux from both stars, taking care to weight each contribution appropriately, and applying the single star constraints to determine whether a particular spatial location would receive flux conducive to surface water. Kane & Hinkel (2013) use the aggregate flux to find a peak wavelength of emission. Assuming the combined spectrum resembles a blackbody, Wien’s Law provides an effective temperature for the total insolation, and hence the limits of Kopparapu et al. (2013) can be applied. This approximation is acceptable for P type systems, where the distance from each star to the planet is similar.
Haghighipour & Kaltenegger (2013) and Kaltenegger & Haghighipour (2013) weight each star’s flux by its effective temperature, and then determine the regions at which this weighted flux matches that of a 1 star at the habitable zone boundaries. This approach is suitable for both S type and P type systems. A detailed analytic solution for calculations of this nature has been undertaken by Cuntz (2014).
Mason et al. (2013) take a similar approach, but they also note that for P type systems, the tidal interaction between primary and secondary can induce rotational synchronisation, which can reduce extreme UV flux and stellar wind pressure, improving conditions in the habitable zone compared to the single star case (see also Zuluaga et al. 2016).
The dynamical constraints on habitability rely heavily on N Body simulation, most prominently the work of Dvorak (Dvorak, 1984, 1986) and Holman & Wiegert (1999). By integrating an ensemble of test particles in a variety of orbits around a binary, regions of dynamical instability can be determined. Holman & Wiegert (1999) used these simulations to develop empirical expressions for a critical orbital semimajor axis, . In the case of a P type system, this represents a minimum value  anything inside is orbitally unstable, as given by the following expression:
(1) 
In the case of an Stype system, represents a maximum value:
(2) 
where is the binary semimajor axis, is the binary orbital eccentricity, and represents the binary mass ratio:
(3) 
The majority of binary habitability calculations rely on the above dynamical constraints. Notable exceptions include Eggl et al. (2012)’s use of Fast Lyapunov Indicators for chaos detection, which yield slightly smaller values of for S type systems (PilatLohinger & Dvorak, 2002), and Jaime et al. (2014)’s use of invariant loops to discover nonintersecting orbits (Pichardo et al., 2005). There is a good deal of research into spinorbit alignments of extrasolar planets under the influence of inclined stellar companions (e.g. Anderson et al. 2016), but this work rarely pertains to terrestrial planet habitability. On the other hand, the evolution of planetary rotation period has been studied intently with regards to habitability of planets in single star systems (e.g. Bolmont et al., 2014; Brown et al., 2014; CuartasRestrepo et al., 2016).
All the above approaches to determining habitability in binary systems rely on an initial 1D calculation of the atmosphere’s response to radiative flux, where the key dimension is atmospheric depth. Equally, 1D approaches can consider the latitudinal variation of flux on a planet’s surface, giving rise to the socalled latitudinal energy balance models or LEBMs, which have been used both in the single star case (Spiegel et al., 2008; Dressing et al., 2010; Vladilo et al., 2013) and for multiple stars (Forgan, 2012, 2014). These are better suited to capture processes that depend on atmospheric circulation, such as the snowball effect arising from icealbedo feedback (Pierrehumbert, 2005; Tajika, 2008), which is likely to occur in systems where the orbits undergo Milankovitch cycles and other secular evolution (Spiegel et al., 2010).
However, all these approaches typically decouple the thermal from the dynamical. The orbital constraints on the HZ are considered separately from the radiative transfer calculations. While they are eventually combined, the binary habitable zones that are constructed do not incorporate the effects of coupled gravitothermal perturbations. Indeed, Holman & Wiegert (1999) admit that their empirical limits on semimajor axis ignore the potential for stable resonances inside the instability region, as well as unstable resonances in stable regions (cf Chavez et al. 2014). It is likely that planets on stable orbits in binary systems will experience relatively strong orbital element evolution. For example, circumbinary planets can undergo rapid precession of periapsis, which affects their ability to be detected via transit (Kostov et al., 2014; Welsh et al., 2015). Presumably the spin evolution of planets in this situation can proceed with similar rapidity. Crucially, climate systems are nonlinear, and can alter their state on very short timescales compared to the planet’s orbital period.
In this work, we consider coupled gravitothermal perturbations on the climate of exoplanets in binary systems. To do so, we present a LEBM directly coupled to an NBody integrator and an obliquity evolution model. We use this combined code to investigate the spinorbitalclimate dynamics of putative planets in two archetypal binary systems: the Ptype system Kepler47, a multiplanet circumbinary system which possesses one exoplanet inside the habitable zone (Orosz et al., 2012); and Alpha Centauri, the nearest star system to the Sun, an S type binary system which was thought to possess a short period, Earthmass exoplanet (Dumusque
et al., 2012)
2 Method
2.1 Latitudinal Energy Balance Modelling
Typically, LEBMs solve the following diffusion equation:
(4) 
Where is the surface temperature, is the effective heat capacity of the atmosphere, is the insolation flux, is the IR cooling and is the albedo. In the above equation, , , and are functions of (either explicitly, as is, or implicitly through ). The latitude appears through . This equation is evolved with the boundary condition at the poles (where ), and requires the assumption that the planet rotates rapidly relative to its orbital period. Our implementation of the LEBM follows that of Spiegel et al. (2008), and has been used previously in studying the climate evolution of planets in binary systems on timescales of order a few hundred years (Forgan, 2012, 2014). In our approach, we consider a given latitude to be habitable if its temperature resides within K, i.e. that surface water is liquid.
The diffusion coefficient determines the efficiency of heat redistribution across latitudes. Its value is defined such that a fiducial Earthlike planet, rotating with period 1 day, orbiting at 1 au around a star of , produces the correct average temperature profile (see e.g. Spiegel et al. 2008; Vladilo et al. 2013). If the planet’s rotation is more rapid, the Coriolis effect will inhibit latitudinal heat transport (see Farrell 1990):
(5) 
where is the rotational angular velocity of the planet, and is the rotational angular velocity of the Earth. This is a necessarily simple expression, but can be made more rigorous through including terms for atmospheric pressure and mean molecular weight (e.g. Williams & Kasting 1997, but see also Vladilo et al. 2013’s attempts to introduce a latitudinal dependence to to mimic the Hadley convective cells on Earth). Beyond this, full global circulation modelling is needed to explore the effects of rotation (Del Genio, 1993, 1996).
As in previous work, we solve equation (4) using an explicit forward time, centre space finite difference algorithm. A global timestep is used, with standard constraint
(6) 
The atmospheric heat capacity , is a function of the planet’s surface ocean fraction and how much of that is frozen, :
(7) 
where . The heat capacities of land, ocean and ice covered areas are
The infrared cooling function is
(8) 
with the optical depth of the atmosphere given as
(9) 
The albedo function is
(10) 
This correctly reproduces the icealbedo feedback phenomenon, which allows a rapid nonlinear increase in albedo as the ice coverage increases.
At any instant, for a single star, the insolation received at a given latitude at an orbital distance is
(11) 
where is the bolometric flux received from the star at a distance of 1 AU, and is the zenith angle:
(12) 
(13) 
The solar hour angle is , and is the solar declination, which is calculated by computing the scalar product of the spinaxis vector and the planetstar separation vector . We obtain the spinaxis vector by rotation of the angular momentum vector in the xaxis by , followed by a rotation around the axis defined by the angular momentum vector by , the axial precession angle (or longitude of winter solstice).
Our rapid rotation assumption requires that we use diurnally averaged quantities, so we also diurnally average :
(14) 
We do this by integrating over the sunlit part of the day, i.e. , where is the radian halfday length at a given latitude. Multiplying by (as if a latitude is illuminated for a full rotation) gives the total diurnal insolation as
(15) 
The radian half day length is calculated as
(16) 
The total insolation is a simple linear combination of the contributions from both stars. If one star is eclipsed by the other, then we set its contribution to to zero. We ensure that the simulation can accurately model an eclipse by adding an extra timestep criterion, ensuring that the transit’s duration will not be less than ten timesteps.
We fix the parameters of the model to those of the Earth: the initial obliquity is set to 23.5 degrees, and the ocean fraction . The rotation period of the body is 1 day. It is important to note that altering these parameters will alter the strength of climate fluctuations, especially if orbits are eccentric. Indeed, Forgan (2012) showed that reducing the planet’s ocean fraction can significantly boost temperature fluctuations in Stype binary systems with fixed orbits, and that increasing obliquity while holding other parameters fixed typically increases the average temperature of the planet. The following results should be considered with these facts in mind.
2.2 The NBody Model
The dynamical evolution of the system utilises a standard 4thorder Hermite integrator with an adaptive shared timestep. We calculate this N Body timestep for all bodies, , by finding the minimum value of :
(17) 
Here, represents the magnitude of the body’s acceleration, and are the magnitudes of the first, second and third derivatives of the acceleration of particle respectively, and is a tunable parameter which we set to 0.002. This is a fairly strict timestep condition, and as such the error in angular momentum is typically one in or better throughout.
2.3 Obliquity Evolution
We adopt the obliquity evolution model of Laskar (1986a, b), developed for the Solar System and subsequently used for putative exoplanet systems (Armstrong et al., 2004; Armstrong et al., 2014a). In this paradigm, the evolution of the obliquity and precession are functions of the inclination variables
(18)  
(19) 
Where is the inclination, and is the longitude of the ascending node. The obliquity and precession evolve according to the following:
(20)  
(21) 
and are all functions of and :
(22)  
(23)  
(24) 
Note that these terms ensure increases in inclination mediate changes in obliquity. Equivalently, if the inclination of a planet’s orbit is increased, the obliquity decreases, as the angle between the orbital plane and the fundamental plane defined by the planet’s spin axis decreases (see Figure 1 of Armstrong et al. 2014a).
That being said, the spin axis of the planet can change regardless of the inclination, due to either direct torques from the star () or from the relativistic precession term . Laskar (1986a) give the direct torque from a single host star as
(25) 
Where is the dynamical ellipticity (i.e. the nonsphericity) of the planet (which we set equal to 0.00328005 for the remainder of this work),
(26) 
and (where the units of must be selected to be appropriate for comparison with ). For a single star, the relativistic precession is
(27) 
where
(28) 
The mean motion can be determined by considering in the context of Kepler’s third law:
(29) 
In this work, we make the following assumptions about these equations in their use for binary stars. In the S type case, we assume that direct torques and precession is generated by the host star only. The secondary can influence the obliquity only through modification of the planet’s orbital elements .
In the case of a P type system, we assume that the torques from both stars coadd. The planet’s orbital elements relative to the system centre of mass are employed in both cases for simplicity. Given the distance of both stars from the centre of mass is small relative to the planet’s semimajor axis, this seems a reasonable assumption (although we do note the need for further investigation of this problem, see Discussion).
2.4 Coupling the Models
To couple the LEBM to the N Body integrator and obliquity evolution model, we elect the simplest route, by forcing all systems to evolve according to a shared timestep. In practice, this means comparing the LEBM and N Body timesteps, i.e.
(30) 
Typically the obliquity evolution timestep is much larger than the other two. This does limit the code’s efficacy when evolving systems with either short dynamical timescales, or short thermal timescales. In the case of a fiducial EarthSun model, we are able to evolve the coupled LEBMN Body system with similar runtime to a LEBM using fixed Keplerian orbits. We will see that in the S type configuration, the addition of N Body physics makes little appreciable difference to computational speed. However, in the P type configuration, the short dynamical timescale of the binary increases the runtime significantly. This could be alleviated by other timestepping approaches, which we address in the Discussion.
We emphasise that correctly resolving the LEBM is crucial  it is a nonlinear system, with positive feedback mechanisms that can operate rapidly compared to the system’s spinorbit dynamical time. It is this property that requires the models to be fully coupled in order to truly understand the climate of planets in dynamically rich systems over secular timescales.
We have tested the N Body integrator and obliquity evolution model against the results of Armstrong et al. (2014a) (their System 1), and find a good match for their orbital elements and spin parameters. In a companion paper (Forgan and Mead, in prep) we test the spinorbitclimate evolution of the Earth under the influence of the Solar System planets, and find that appropriate Milankovitch cycles in the planet’s spinorbit parameters do indeed arise.
3 Results
We now apply our combined model to the two archetypal P and S type binary systems. We will be comparing runs with obliquity evolution switched on and off to investigate what climate features are due to either orbital or spin evolution.
3.1 Kepler47
Setup
The Kepler47 system contains a 1.043 star and an 0.362 star orbiting each other with a period of around 7.5 days. We adopt the orbital parameters of Orosz et al. (2012), with a semimajor axis of 0.0836 AU and eccentricity 0.0234, and assume that the stars’ luminosities are determined by standard main sequence relations.
Kepler47c orbits inside the circumbinary habitable zone at 0.989 AU, with an eccentricity upper limit of 0.41. As we are using the Kepler47 system as an archetype for terrestrial habitability in P type systems, we replace Kepler47c with an Earth mass planet orbiting at the same semimajor axis, and investigate both low and high eccentricity orbits. Kepler47b orbits interior to Kepler47c with a semimajor axis of 0.2956 AU with eccentricity 0.034, and period 49.5 days. We investigate the climate of our terrestrial planet both with and without Kepler47b’s presence.
Zero Eccentricity, Without Kepler47b
Figure 1 shows the orbital evolution of a terrestrial planet orbiting the Kepler47 binary at AU with zero eccentricity and an initial inclination of 0.5 relative to the binary plane. We run the simulation for 10,000 years, with sufficiently high snapshot frequency that the orbital period of the binary (0.0205 years) is well resolved. The planet’s orbit is relatively stable, undergoing small eccentricity and inclination variations of around 800 and 400 year periods respectively (note also that the argument of periapsis precesses on a similar timescale).
In the case where the obliquity is fixed, the planet’s climate settles to a stable state, with mean temperatures fluctuating by around 0.1 K (top row of Figure 2). We can see in the periodogram for fixed obliquity that the major contribution to temperature fluctuation is seasonal variation over the orbital period of 0.829 years (and its harmonics at of the period), closely followed by a contribution at the binary period of 0.0205 years as the relative insolation from each object varies. Finally, we see a significantly weaker contribution from eccentricity variation at 800 years. There are no low order mean motion resonances between the binary and planet period  the system is closest to a 80:2 resonance. There is no evidence of such a resonance in the temperature data, which would result in a peak at approximately 1.66 years in the periodogram.
In the case where obliquity is allowed to vary (bottom row of Figure 2), we can immediately detect climatic variations from inspecting the maximum, mean and minimum temperature curves. The presence of an extra peak at around 400 years in the temperature periodogram (bottom right of Figure 2) shows that the inclination is forcing similar variations in obliquity and precession angle (Figure 1). Generally speaking, the planet’s climate now shows a richer set of resonant features in the periodogram with periods greater than that of the orbital periods in play.
Zero Eccentricity, with Kepler47b
The previous section has shown that single planets in P type systems will undergo secular evolution quite similar to that of Milankovitch cycles (albeit at a much reduced timescale). We now add Kepler47b to the system (with zero eccentricity and inclination) to gauge what effect neighbouring planets might have on the secular evolution of circumbinary habitable climates.
Figure 3 shows the orbital evolution of the Kepler47c substitute. Comparing to the previous section (Figure 1), we see that the eccentricity variation has not changed much, but the inclination variation has decreased its period by a factor of roughly two. Interestingly, no such changes are seen in the obliquity and precession evolution, indicating that stellar torques are presumably dominant.
The periodograms for both cases (Figure 4) show little change in the climate by adding a neighbour planet. The periodograms show no signs of Kepler47b’s influence at its orbital period of 0.1355 years. The features seen at 0.1355 years with obliquity evolution exist in the previous run without Kepler47b. The planets are not in mean motion resonance  they are closest to a 49:8 mean motion resonance, which would indicate a peak at approximately 6.63 years, which is not seen in either case.
High Eccentricity, no b
We now remove Kepler47b from the system, and increase the eccentricity of our habitable planet to 0.4. The dynamical evolution (Figure 5) is more rapid, with small eccentricity and inclination oscillations about the original value with a period of around 550 years, and similar obliquity and precession evolution. Note the amplitude modulation of the inclination, which coincides with peak eccentricity.
Naturally, the climate of the body experiences stronger temperature oscillations even with obliquity switched off (top row of Figure 6). The periodogram shows greater importance for the seasonal variation, as well as the eccentricity variation peak at 550 years. As the planet and binary are not in mean motion resonance, the contribution of the binary to the planet’s eccentricity periodogram is smeared between 0.02 and 0.03 years due to the planet’s increased eccentricity. Note that this increased eccentricity raises the maximum temperature beyond the runaway greenhouse limit of 340K. The runaway greenhouse effect is not modelled by the LEBM, and we should be careful when making statements about this configuration’s habitability. Some weak modes appear around the planet’s orbital period of 0.829 years, but their origin is unclear  presumably they are linked to the precession of the planet’s periapsis relative to that of the binary.
Allowing obliquity to vary allows other oscillations to assume greater importance. Indeed, the variations caused by binary motion are close to negligible in this case, especially compared to variations in the yeardecade range.
3.2 Alpha Centauri B
Setup
The Alpha Centauri system is in fact a hierarchical triple system, with Alpha Centauri A and B orbiting each other at 23.4 AU with eccentricity 0.5179. We neglect the third component, Proxima Centauri, as it orbits at great distance and is of sufficiently low mass (Wertheimer & Laughlin, 2006). We consider Cen B as the host star for a planetary system.
The stellar masses are , , and their luminosities are and respectively (Thevenin et al., 2002). This modifies the location of the habitable zone as was previously measured by Forgan (2012), as they used main sequence relations for the luminosity.
We do not model the presence of Cen Bb, as its 3 day orbit would place it extremely close to Cen B, and hence is unlikely to produce a significant perturbation on any planets within the habitable zone. Instead, we place a single Earthlike planet in the system near the outer edge of the habitable zone, on a circular orbit at AU , where the effects of Cen A are maximal. To ensure obliquity evolution occurs, we give our planet a small inclination of relative to the binary plane.
However, we do wish to consider the relative strength of Milankovitch cycles resulting from the binary compared to those induced by neighbouring planets (cf Figure 8 of AndradeInes & Michtchenko 2014). We attempt to maximise this effect by running another set of models with a second Earthmass body orbiting in 3:2 resonance with our habitable world (with a zero inclination orbit).
Single Planet Runs
Figure 7 shows the dynamical evolution of the planet around Cen B. The initially zero eccentricity is forced to a maximum of 0.05 on a cycle of approximately 14,500 years. The obliquity and precession evolve with a slightly longer period, resulting in the eccentricity and obliquity cycles drifting in and out of phase.
This phase drift results in markedly different climate evolution of the body, compared to the case where obliquity is held fixed (Figure 8). In the fixed obliquity case, the eccentricity cycle induces a temperature oscillation of approximately 2K (to add to the radiative oscillation of 5K due to the changing proximity of Cen A). The periodogram shows the two dominant oscillation modes at 79.9 and 14,500 years. Their strength is indicated by the strength of their subsequent harmonics, which can be seen down to the tenth level!
A quite different picture emerges if obliquity evolution is activated (bottom row of Figure 8). The temperature oscillations are now modulated by the phase drift between eccentricity and obliquity, which is periodic over 200,000 year timescales. When the two cycles are in phase, we see the largest temperature oscillations (e.g. at 200,000 years).
Adding a planet in 3:2 mean motion resonance
We now consider joint planetarybinary Milankovitch cycles by adding an Earth mass planet on a circular orbit at 0.9293 AU, placing it in 3:2 mean motion resonance with the habitable planet. Test runs with Cen A absent show the additional planet induces regular eccentricity oscillations in the habitable planet with amplitude of approximately 0.01, and a period of approximately 500 years. Incidentally, the absence of Cen A would also place both planets outside the habitable zone.
With Cen A present, the combination of stellar and planetary forcings produces eccentricity oscillations of maximum amplitude 0.08 (left panel of Figure 9) and with a mix of dominant periods, as opposed to the distinct 14,500 year period observed in the single planet case. The inclination varies with a period of approximately 30,000 years, with a distinctive shift in mean inclination of around 0.001 radians (i.e. 0.05). The obliquity and precession continue to evolve at close to the eccentricity oscillation period, but the amplitude of their oscillations varies on approximately twice this timescale.
The uniform temperature evolution cycles seen in Figure 8 are now more confused with the addition of a neighbour planet (Figure 10). With obliquity evolution switched off (top row), the extra structure introduced into the eccentricity and inclination oscillations leaves an imprint on the temperature curves. This can be seen in its periodogram (top right panel of Figure 10), which shows a relatively weak feature at the perturbing planet’s orbital period, and at the resonant period of twice the perturber’s period (or equivalently, three times the habitable planet’s period). The perturbations induced by the additional planet produce temperature variations of up to 2K compared to the single planet case.
With obliquity evolution turned on (bottom panel), the eccentricity/obliquity relationship seen in the previous case is preserved, resulting in phase drift between the two oscillations. However, the extra structure in the eccentricity oscillation prevents the smooth amplitude modulation of temperature that we saw in the bottom right panel of Figure 8. It is broadly present, but heavily modified by the presence of the neighbouring planet. The periodogram still reveals weak signals at the perturbing planet’s period, and the strong peak feature at approximately 14,500 years is now split in two. There is also a significant increase in signal for periods of order 1001000 years.
Additional giant planets in a system like this might be expected to produce even larger excursions from circular orbits and stronger Milankovitch cycling. Given that planet formation models disfavour the creation of Jupiter mass bodies in this system (Xie et al., 2010) and are ruled out by observations of the Cen system, at least at periods less than 1 year (Endl et al., 2001; Dumusque et al., 2012; Demory et al., 2015) this is not a particular concern. But, one might imagine that undetected Neptune mass bodies could be present in this system on relatively long period orbits, and such bodies would be responsible for longer period Milankovitch cycles similar to that of Earth’s.
4 Discussion
4.1 Limitations of the Model
LEBM modelling is by its definition a compromise between the granularity of a climate simulation and computational expediency. This compromise is stretched further by the coupling of the NBody integrator and obliquity evolution. We have adopted a very simple coupling where both the NBody and LEBM components are constrained to follow the same global timestep.
This timestep system works extremely well for systems where the dynamical timescale is relatively long, such as the S type binary systems. In this scenario, the system timestep is limited only by the LEBM, and as such we can run simulations with similar wallclock times as that of a LEBM using fixed Keplerian orbits. However, in the P type scenario, the dynamical timescale is relatively short, and the system is limited by the N Body timestep required to resolve the binary.
There are several possible strategies for mitigating this timestep issue. The most straightforward solution is to adopt a nonshared timestep for the NBody component, allowing some of the bodies to possess shorter N Body timesteps. This would reduce the computational load of evolving all the bodies (and the LEBM) at what can be very short timesteps. Another solution would require the interpolation of body motions (in the case where the LEBM timestep is small compared to the N Body timestep), but this would likely produce only marginal gains in speed. Perhaps the best solution for P type systems would be chain regularisation of the tight binary orbit (Mikkola & Aarseth, 1990, 1993).
Aside from the new challenges arising from the adoption of the NBody integrator, there are the usual limitations that many LEBMs are subject to. Our implementation of the LEBM is among the most simple available which can still broadly reproduce the seasonal temperature profiles of a fiducial Earth model. The principal advantage of this simplicity is its ease of interpretation, but we must acknowledge that more advanced models may produce features we cannot.
For example, we do not model the carbonatesilicate (CS) cycle, which moderates fluctuations in atmospheric temperature by increasing and reducing the partial pressure of carbon dioxide. The timescale on which we expect levels to vary depends on the planet’s geochemical properties, especially its ocean circulation. For Earthlike planets, the equilbriation timescale of is approximately half a million years (Williams & Kasting, 1997) which is far shorter than the Milankovitch cycles experienced by the planetary bodies in this analysis. However, our understanding of the CS cycle is rooted firmly in our understanding of the Earth, which orbits a single star. It remains unclear whether a planet in a binary star system would possess a similar equilibriation timescale, even if the planet was effectively identical to the Earth.
While we have taken the first steps towards coupling celestial dynamics and LEBM climate modelling here, there are still several steps ahead of us. For example, the tidal interactions between bodies will also modify orbits of habitable worlds, in particular reducing their eccentricity and modifying their rotational period (Bolmont et al., 2014; Cunha et al., 2014). While this is unlikely to be an issue for the orbital configurations adopted in this analysis, it remains the case that while the tidal interactions between the binary stars is well characterised (e.g. Mason et al. 2013; Zuluaga et al. 2016), the tidal evolution of planets in P type systems remains relatively unexplored.
Also to be explored in full are the obliquity variations felt by planets in binary systems. We have adopted a set of equations designed for a single star planetary system, and assumed they are valid when there are two stars present. In effect we have assumed that in S type systems, the secondary’s direct tidal torque on planetary spin is negligible, and that in P type systems the direct torques always coadd. Is this always the case? More investigation is needed.
We should also note that the strength of Milankovitch cycles measured by the LEBM will be an underestimate. Tests conducted using Solar system parameters (Forgan & Mead, in prep.) give Milankovitch cycles for the Earth that are an order of magnitude smaller in temperature variation than observed in paleoclimate data (Zachos et al., 2001; Lisiecki & Raymo, 2005). Paradoxically, stochastic EBMs, with additional random noise, can enhance periodic variations through the phenomenon of stochastic resonance (Imkeller, 2001; Benzi, 2010). Obliquity variation does produce a much richer set of temperature variations on decadal timescales, which may be forced into stochastic resonant behaviour under appropriate circumstances. Future investigations should consider adding a random noise term to the LEBM equation to permit this behaviour.
4.2 Implications for Habitability
So what have we gained by this coupling of N Body and LEBM integrators? Initially, we are able to confirm that in general, the decoupled approach of considering the radiative and gravitational perturbations separately is broadly acceptable.
Previous work in this field is not invalidated by our results, but it makes explicit some general principles that are already known implicitly. Firstly, the habitable zone of a planetary system is defined by more than where the radiation sources are in the system. The gravitational sources are equally important. We know this on Earth thanks to our understanding of Milankovitch cycles, and the Earth’s orbital and spin cycles are relatively weak when compared to measured cycles for Earthlike planets in typical exoplanet system configurations around a single star (Spiegel et al. 2010, Forgan & Mead, in prep.).
Secondly, the habitable zone of binary systems is even more sensitive to the gravitational field than single star systems. This is already demonstrated implicitly by the NBody simulations of orbital stability discussed in the Introduction. Our results clearly identify the effect of orbital and spin stability on climate. We show that relatively strong Milankovitch cycles exist in binary systems, even if there is only one planet present. The periods of these cycles are in general shorter than that of single star systems, but of similar amplitudes. Even on short timescales, the radiative perturbations induced over the orbital period of the binary are detectable in the mean temperature of the planet.
Thirdly, the circadian rhythms of life on planets in binary systems will be forced to adapt to the rhythms present in the binary system, as is evidenced by analogous studies of lunar photoperiodism in terrestrial organisms (O’MalleyJames et al. 2012; Forgan et al. 2015 and references within). Temperature fluctuations of several K on timescales ranging from less than a year to almost a century (depending on whether the system is P or S type) is likely to produce significant fluctuations in surface coverage of biomes. The rapid Milankovitch cycles are likely to play a stronger role also. More sophisticated climate models coupled to NBody physics (for example, 3D global circulation models) may show potential for more, shorter Ice Ages, and briefer interglacial periods. The presence of such rapid changes to environmental selection pressure will have an indelible effect on the evolution of organisms in binary planetary systems. Future work should build on recent attempts to produce 3D General Circulation Models of circumbinary planets (cf May & Rauscher 2016), incorporating the systems’ gravitational evolution to determine these effects in detail.
5 Conclusions
We have investigated Milankovitch cycles both circumbinary (P type) and distant binary (S type) systems, using Kepler47 and Centauri as archetypes. To do this, we coupled a 1D latitudinal energy balance climate model (LEBM) with an NBody integrator to follow the orbital evolution, and an obliquity evolution algorithm to study the spinaxis evolution.
We find that the combined spinorbitradiative perturbations induced by a companion star on a habitable planet produce Milankovitch cycles for both types of binary system, even when other planets are not present. Periodogram analysis identifies both dynamical and secular oscillations in the mean temperature of planets in these systems, over a variety of short and long periods, as well as the presence of radiative perturbations directly linked to the period of the binary. The strength of these oscillations is sensitive to the orbital configuration of the system. The relative phase between eccentricity, precession and obliquity cycles is important, just as it is for the Earth.
In general, we find these Milankovitch cycles are significantly shorter than comparable cycles on the Earth (in some cases shorter than 1000 years), although the amplitude of the changes they produce in the planets’ orbital elements are comparable to those experienced by Earth. This work demonstrates the need to consider joint dynamicsclimate simulations of habitable worlds in binary systems, if we are to truly assess the potential for the birth and growth of biospheres on worlds with two suns.
Acknowledgments
DHF gratefully acknowledges support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC2011ADG. This work relied on the compute resources of the St Andrews MHD Cluster. The author thanks both Nader Haghighipour and James Gilmore for insightful comments on an early version of this manuscript. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. The code used in this work is now available open source as OBERON, which can be downloaded at github.com/dh4gan/oberon.
Footnotes
 pagerange: Milankovitch Cycles of Terrestrial Planets in Binary Star Systems–References
 pubyear:
 http://www.univie.ac.at/adg/schwarz/multiple.html
 This detection is no longer considered to be credible by some groups, due to concerns with how stellar activity is filtered out of radial velocity data (Hatzes, 2013). Recent attempts to detect Cen Bb via transit show a null result (Demory et al., 2015), and reanalysis of the radial velocity data suggests that Cen Bb does not exist (Rajpaul et al., 2016).
References
 Anderson K. R., Storch N. I., Lai D., 2016, MNRAS, 456, 3671
 AndradeInes E., Michtchenko T. A., 2014, MNRAS, 444, 2167
 Armstrong J. C., Leovy C. B., Quinn T., 2004, Icarus, 171, 255
 Armstrong J. C., Barnes R., DomagalGoldman S., Breiner J., Quinn T. R., Meadows V. S., 2014a, Astrobiology, 14, 277
 Armstrong D. J., Osborn H. P., Brown D. J. A., Faedi F., Gomez Maqueo Chew Y., Martin D. V., Pollacco D., Udry S., 2014b, MNRAS, 444, 1873
 Benzi R., 2010, Nonlinear Processes in Geophysics, 17, 431
 Bolmont E., Raymond S. N., von Paris P., Selsis F., Hersant F., Quintana E. V., Barclay T., 2014, ApJ, 793, 3
 Brown S., Mead A., Forgan D., Raven J., Cockell C., 2014, International Journal of Astrobiology, 13, 279
 Chavez C. E., Georgakarakos N., Prodan S., ReyesRuiz M., Aceves H., Betancourt F., PerezTijerina E., 2014, MNRAS, 446, 1283
 CuartasRestrepo P., Melita M., Zuluaga J., Portilla B., Sucerquia M., Miloni O., 2016, MNRAS, pp submitted, eprint arXiv:1606.07546
 Cunha D., Correia A. C., Laskar J., 2014, International Journal of Astrobiology, 14, 233
 Cuntz M., 2014, ApJ, 780, 14
 Del Genio A., 1993, Icarus, 101, 1
 Del Genio A., 1996, Icarus, 120, 332
 Demory B.O., et al., 2015, MNRAS, p. in press
 Desidera S., Barbieri M., 2007, A&A, 462, 345
 Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656
 Doyle L. R., et al., 2011, Science (New York, N.Y.), 333, 1602
 Dressing C. D., Spiegel D. S., Scharf C. A., Menou K., Raymond S. N., 2010, ApJ, 721, 1295
 Dumusque X., et al., 2012, Nature, 491, 207
 Dunhill A. C., Alexander R. D., 2013, MNRAS, 435, 2328
 Duquennoy A., Mayor M., 1991, A&A, 248, 485
 Dvorak R., 1984, Celestial Mechanics, 34, 369
 Dvorak R., 1986, A&A, 167, 379
 Eggl S., PilatLohinger E., Georgakarakos N., Gyergyovits M., Funk B., 2012, ApJ, 752, 74
 Endl M., Kï¿½rster M., Els S., Hatzes A. P., Cochran W. D., 2001, Astronomy and Astrophysics, 374, 675
 Farrell B. F., 1990, Journal of Atmospheric Sciences, 47, 2986
 Forgan D., 2012, MNRAS, 422, 1241
 Forgan D., 2014, MNRAS, 437, 1352
 Forgan D. H., Mead A., Cockell C. S., Raven J. A., 2015, International Journal of Astrobiology, 14, 465
 G. Martin R., Armitage P. J., Alexander R. D., 2013, ApJ, 773, 74
 Haghighipour N., Kaltenegger L., 2013, ApJ, 777, 166
 Hart M., 1979, Icarus, 37, 351
 Hatzes A. P., 2013, ApJ, 770, 133
 Hatzes A. P., Cochran W. D., Endl M., McArthur B., Paulson D. B., Walker G. A. H., Campbell B., Yang S., 2003, ApJ, 599, 1383
 Holman M. J., Wiegert P. A., 1999, The Astronomical Journal, 117, 621
 Huang S.S., 1959, PASP, 71, 421
 Imkeller P., 2001, in Imkeller P., Von Storch J.S., eds, , Stochastic Climate Models. Birkhäuser Basel, Basel, pp 213–240, doi:10.1007/9783034882873_9, http://link.springer.com/10.1007/9783034882873{_}9
 Jaime L. G., Aguilar L., Pichardo B., 2014, MNRAS, 443, 260
 Kaltenegger L., Haghighipour N., 2013, ApJ, 777, 165
 Kane S. R., Hinkel N. R., 2013, ApJ, 762, 7
 Kasting J., Whitmire D., Reynolds R., 1993, Icarus, 101, 108
 Kopparapu R. K., et al., 2013, ApJ, 765, 131
 Kopparapu R. K., Ramirez R. M., SchottelKotte J., Kasting J. F., DomagalGoldman S., Eymet V., 2014, ApJ, 787, L29
 Kostov V. B., et al., 2014, ApJ, 784, 14
 Laskar J., 1986a, A&A, 157, 59
 Laskar J., 1986b, A&A, 164, 437
 Lisiecki L. E., Raymo M. E., 2005, Paleoceanography, 20, id PA1003
 Marzari F., Thebault P., Scholl H., Picogna G., Baruteau C., 2013, Astronomy & Astrophysics, 553, A71
 Mason P. a., Zuluaga J. I., Clark J. M., CuartasRestrepo P. a., 2013, ApJ, 774, L26
 May E. M., Rauscher E., 2016, The Astrophysical Journal, 826, 225
 Meschiari S., 2014, ApJ, 790, 41
 Mikkola S., Aarseth S. J., 1990, Celestial Mechanics and Dynamical Astronomy (ISSN 09232958), 47, 375
 Mikkola S., Aarseth S. J., 1993, Celestial Mechanics & Dynamical Astronomy, 57, 439
 O’MalleyJames J. T., Raven J. A., Cockell C. S., Greaves J. S., 2012, Astrobiology, 12, 115
 Orosz J. A., et al., 2012, Science (New York, N.Y.), 337, 1511
 Pichardo B., Sparke L. S., Aguilar L. A., 2005, MNRAS, 359, 521
 Pierrehumbert R. T., 2005, Journal of Geophysical Research, 110, D01111
 PilatLohinger E., Dvorak R., 2002, Celestial Mechanics and Dynamical Astronomy, pp 143–153
 Queloz D., et al., 2000, A&A, 354, 99
 Quintana E. V., Lissauer J. J., Chambers J. E., Duncan M. J., 2002, ApJ, 576, 982
 Quintana E. V., Adams F. C., Lissauer J. J., Chambers J. E., 2007, ApJ, 660, 807
 Rafikov R. R., 2013, ApJ, 764, L16
 Rafikov R. R., Silsbee K., 2014a, ApJ, 798, 69
 Rafikov R. R., Silsbee K., 2014b, ApJ, 798, 70
 Raghavan D., et al., 2010, ApJSS, 190, 1
 Rajpaul V., Aigrain S., Roberts S., 2016, MNRAS, 456, L6
 Sigurdsson S., Richer H. B., Hansen B. M., Stairs I. H., Thorsett S. E., 2003, Science (New York, N.Y.), 301, 193
 Silsbee K., Rafikov R. R., 2015, ApJ, 798, 71
 Spiegel D. S., Menou K., Scharf C. A., 2008, ApJ, 681, 1609
 Spiegel D. S., Raymond S. N., Dressing C. D., Scharf C. A., Mitchell J. L., 2010, ApJ, 721, 1308
 Tajika E., 2008, ApJ, 680, L53
 Thébault P., Marzari F., Scholl H., 2008, MNRAS, 388, 1528
 Thébault P., Marzari F., Scholl H., 2009, MNRAS, 393, L21
 Thevenin F., Provost J., Morel P., Berthomieu G., Bouchy F., Carrier F., 2002, Astronomy and Astrophysics, 392, L9
 Thorsett S. E., Arzoumanian Z., Taylor J. H., 1993, ApJ, 412, L33
 Vladilo G., Murante G., Silva L., Provenzale A., Ferri G., Ragazzini G., 2013, ApJ, 767, 65
 Welsh W. F., et al., 2012, Nature, 481, 475
 Welsh W. F., et al., 2015, ApJ, 809, 26
 Wertheimer J. G., Laughlin G., 2006, The Astronomical Journal, 132, 1995
 Wiegert P. A., Holman M. J., 1997, The Astronomical Journal, 113, 1445
 Williams D., Kasting J., 1997, Icarus, 129, 254
 Xie J.W., Zhou J.L., Ge J., 2010, ApJ, 708, 1566
 Zachos J., Pagani M., Sloan L., Thomas E., Billups K., 2001, Science (New York, N.Y.), 292, 686
 Zucker S., Mazeh T., Santos N. C., Udry S., Mayor M., 2004, A&A, 426, 695
 Zuluaga J. I., Mason P. A., CuartasRestrepo P. A., 2016, The Astrophysical Journal, 818, 160