Milankovitch Cycles in Binary Systems

Milankovitch Cycles of Terrestrial Planets in Binary Star Systems


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 gravito-climatic oscillations on dynamical and secular timescales.

We investigate the Kepler-47 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.

astrobiology, methods:numerical, planets and satellites: general

1 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 so-called 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 post-main sequence stars, in particular the binary pulsar B160-26 (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 Kepler-16 (Doyle et al., 2011), Kepler-34 and Kepler-35 (Welsh et al., 2012), and Kepler-47 (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 systems3, giving an occurrence rate of around 4% (previous estimates on a much smaller exoplanet population by Desidera & Barbieri 2007 placed the fraction of planets in S type systems at 20%). At gas giant masses, the occurrence rate of planets around P type binaries is thought to be similar to that of single stars (Armstrong et al., 2014b).

However, theoretical modelling indicates that the dynamical landscape of the binary significantly affects the planet formation process, both for S-type (Wiegert & Holman, 1997; Quintana et al., 2002, 2007; Thébault et al., 2008, 2009; Xie et al., 2010; Rafikov & Silsbee, 2014b, a) and P-type 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:

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

  2. Regions of the system are orbitally unstable

These joint thermal-dynamical 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:


In the case of an S-type system, represents a maximum value:


where is the binary semimajor axis, is the binary orbital eccentricity, and represents the binary mass ratio:


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 (Pilat-Lohinger & Dvorak, 2002), and Jaime et al. (2014)’s use of invariant loops to discover non-intersecting orbits (Pichardo et al., 2005). There is a good deal of research into spin-orbit 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; Cuartas-Restrepo 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 so-called 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 ice-albedo 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 gravito-thermal 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 gravito-thermal perturbations on the climate of exoplanets in binary systems. To do so, we present a LEBM directly coupled to an N-Body integrator and an obliquity evolution model. We use this combined code to investigate the spin-orbital-climate dynamics of putative planets in two archetypal binary systems: the P-type system Kepler-47, a multi-planet 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, Earth-mass exoplanet (Dumusque et al., 2012)4. By evolving the orbits of the bodies in tandem with the climate, we are able to detect climate variations that are directly linked to the binary’s orbit, and the secular evolution of the planet’s orbit and spin.

In section 2, we describe the LEBM, and how the N Body model is coupled to it. In section 3 we describe the simulation setup and results on dynamical and secular timescales, in section 4 we discuss the implications for habitability, and in section 5 we summarise the work.

2 Method

2.1 Latitudinal Energy Balance Modelling

Typically, LEBMs solve the following diffusion equation:


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):


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


The atmospheric heat capacity , is a function of the planet’s surface ocean fraction and how much of that is frozen, :


where . The heat capacities of land, ocean and ice covered areas are

The infrared cooling function is


with the optical depth of the atmosphere given as


The albedo function is


This correctly reproduces the ice-albedo feedback phenomenon, which allows a rapid non-linear 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


where is the bolometric flux received from the star at a distance of 1 AU, and is the zenith angle:


The solar hour angle is , and is the solar declination, which is calculated by computing the scalar product of the spin-axis vector and the planet-star separation vector . We obtain the spin-axis vector by rotation of the angular momentum vector in the x-axis 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 :


We do this by integrating over the sunlit part of the day, i.e. , where is the radian half-day length at a given latitude. Multiplying by (as if a latitude is illuminated for a full rotation) gives the total diurnal insolation as


The radian half day length is calculated as


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 S-type 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 N-Body Model

The dynamical evolution of the system utilises a standard 4th-order Hermite integrator with an adaptive shared timestep. We calculate this N Body timestep for all bodies, , by finding the minimum value of :


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


Where is the inclination, and is the longitude of the ascending node. The obliquity and precession evolve according to the following:


and are all functions of and :


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


Where is the dynamical ellipticity (i.e. the non-sphericity) of the planet (which we set equal to 0.00328005 for the remainder of this work),


and (where the units of must be selected to be appropriate for comparison with ). For a single star, the relativistic precession is




The mean motion can be determined by considering in the context of Kepler’s third law:


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 co-add. 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 semi-major 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.


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 Earth-Sun model, we are able to evolve the coupled LEBM-N 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 spin-orbit 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 spin-orbit-climate evolution of the Earth under the influence of the Solar System planets, and find that appropriate Milankovitch cycles in the planet’s spin-orbit 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 Kepler-47


The Kepler-47 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 semi-major axis of 0.0836 AU and eccentricity 0.0234, and assume that the stars’ luminosities are determined by standard main sequence relations.

Kepler-47c orbits inside the circumbinary habitable zone at 0.989 AU, with an eccentricity upper limit of 0.41. As we are using the Kepler-47 system as an archetype for terrestrial habitability in P type systems, we replace Kepler-47c with an Earth mass planet orbiting at the same semi-major axis, and investigate both low and high eccentricity orbits. Kepler-47b orbits interior to Kepler-47c 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 Kepler-47b’s presence.

Zero Eccentricity, Without Kepler-47b

Figure 1: The dynamical evolution of the terrestrial planet with Kepler-47c’s semimajor axis, and zero eccentricity. Left, the orbital evolution of the body, as given by its eccentricity and inclination. Right, the spin evolution as given by the obliquity and precession angles.

Figure 2: The climate evolution of the Kepler-47c terrestrial planet, with obliquity evolution switched off (top row) and switched on (bottom row). Left: The global maximum, minimum and mean temperatures on the surface over 10,000 years. Right: Periodograms for the mean temperature. The red dashed lines indicate the planet’s orbital period of 0.829 years, and its harmonics (, … of the period).

Figure 1 shows the orbital evolution of a terrestrial planet orbiting the Kepler-47 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 Kepler-47b

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 Kepler-47b 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: The dynamical evolution of the terrestrial planet with Kepler-47c’s semimajor axis, and zero eccentricity, in the presence of Kepler-47b. Left, the orbital evolution of the body, as given by its eccentricity and inclination. Right, the spin evolution as given by the obliquity and precession angles.

Figure 4: The climate evolution of the Kepler-47c terrestrial planet in the presence of Kepler-47b, with obliquity evolution switched off (top row) and switched on (bottom row). Left: The global maximum, minimum and mean temperatures on the surface over 10,000 years. Right: Periodograms for the mean temperature. The red dashed lines indicate Kepler-47b’s orbital period of 0.1355 years, and its harmonics (, … of the period).

Figure 3 shows the orbital evolution of the Kepler-47c 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 Kepler-47b’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 Kepler-47b. 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 Kepler-47b 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.

Figure 5: The dynamical evolution of the terrestrial planet with Kepler-47c’s semimajor axis, and eccentricity 0.4. Left, the orbital evolution of the body, as given by its eccentricity and inclination. Right, the spin evolution as given by the obliquity and precession angles.

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 year-decade range.

Figure 6: The climate evolution of the Kepler-47c terrestrial planet at high eccentricity, with obliquity evolution switched off (top row) and switched on (bottom row). Left: The global maximum, minimum and mean temperatures on the surface over 10,000 years. Right: Periodograms for the mean temperature. The red dashed lines indicate the planet’s orbital period of 0.829 years, and its harmonics (, … of the period).

3.2 Alpha Centauri B


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 Andrade-Ines & Michtchenko 2014). We attempt to maximise this effect by running another set of models with a second Earth-mass body orbiting in 3:2 resonance with our habitable world (with a zero inclination orbit).

Single Planet Runs

Figure 7: The dynamical evolution of the terrestrial planet orbiting Cen B. Left, the orbital evolution of the body, as given by its eccentricity. We refrain from plotting the inclination, as its fluctuations are extremely low with no obvious periodic oscillation. Right, the spin evolution as given by the obliquity and precession angles.

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.

Figure 8: The climate evolution of the Cen B terrestrial planet, with obliquity evolution switched off (top row) and switched on (bottom row). Left: The global maximum, minimum and mean temperatures on the surface over 100,000 years (obliquity evolution off) and over approximately 300,000 years (obliquity evolution on). Right: Periodograms for the mean temperature. The red dashed lines indicate the binary’s orbital period of 79 years, and its harmonics (, … of the period).

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 planetary-binary 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.

Figure 9: The dynamical evolution of the terrestrial planet orbiting Cen B. Left, the orbital evolution of the body, as given by its eccentricity and inclination. Right, the spin evolution as given by the obliquity and precession angles.

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 100-1000 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.

Figure 10: The climate evolution of the Cen B terrestrial planet, with obliquity evolution switched off (top row) and switched on (bottom row). Left: The global maximum, minimum and mean temperatures on the surface over 100,000 years (obliquity evolution off) and over approximately 300,000 years (obliquity evolution on). Right: Periodograms for the mean temperature. The red dashed lines indicate the binary’s orbital period of 79 years, and its harmonics (, … of the period).

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 N-Body integrator and obliquity evolution. We have adopted a very simple coupling where both the N-Body 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 non-shared timestep for the N-Body 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 N-Body 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 carbonate-silicate (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 co-add. 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 N-Body 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’Malley-James 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 N-Body 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 Kepler-47 and Centauri as archetypes. To do this, we coupled a 1D latitudinal energy balance climate model (LEBM) with an N-Body integrator to follow the orbital evolution, and an obliquity evolution algorithm to study the spin-axis evolution.

We find that the combined spin-orbit-radiative 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 dynamics-climate 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.


DHF gratefully acknowledges support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. 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


  1. pagerange: Milankovitch Cycles of Terrestrial Planets in Binary Star SystemsReferences
  2. pubyear:
  4. 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 re-analysis of the radial velocity data suggests that Cen Bb does not exist (Rajpaul et al., 2016).


  1. Anderson K. R., Storch N. I., Lai D., 2016, MNRAS, 456, 3671
  2. Andrade-Ines E., Michtchenko T. A., 2014, MNRAS, 444, 2167
  3. Armstrong J. C., Leovy C. B., Quinn T., 2004, Icarus, 171, 255
  4. Armstrong J. C., Barnes R., Domagal-Goldman S., Breiner J., Quinn T. R., Meadows V. S., 2014a, Astrobiology, 14, 277
  5. 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
  6. Benzi R., 2010, Nonlinear Processes in Geophysics, 17, 431
  7. Bolmont E., Raymond S. N., von Paris P., Selsis F., Hersant F., Quintana E. V., Barclay T., 2014, ApJ, 793, 3
  8. Brown S., Mead A., Forgan D., Raven J., Cockell C., 2014, International Journal of Astrobiology, 13, 279
  9. Chavez C. E., Georgakarakos N., Prodan S., Reyes-Ruiz M., Aceves H., Betancourt F., Perez-Tijerina E., 2014, MNRAS, 446, 1283
  10. Cuartas-Restrepo P., Melita M., Zuluaga J., Portilla B., Sucerquia M., Miloni O., 2016, MNRAS, pp submitted, eprint arXiv:1606.07546
  11. Cunha D., Correia A. C., Laskar J., 2014, International Journal of Astrobiology, 14, 233
  12. Cuntz M., 2014, ApJ, 780, 14
  13. Del Genio A., 1993, Icarus, 101, 1
  14. Del Genio A., 1996, Icarus, 120, 332
  15. Demory B.-O., et al., 2015, MNRAS, p. in press
  16. Desidera S., Barbieri M., 2007, A&A, 462, 345
  17. Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656
  18. Doyle L. R., et al., 2011, Science (New York, N.Y.), 333, 1602
  19. Dressing C. D., Spiegel D. S., Scharf C. A., Menou K., Raymond S. N., 2010, ApJ, 721, 1295
  20. Dumusque X., et al., 2012, Nature, 491, 207
  21. Dunhill A. C., Alexander R. D., 2013, MNRAS, 435, 2328
  22. Duquennoy A., Mayor M., 1991, A&A, 248, 485
  23. Dvorak R., 1984, Celestial Mechanics, 34, 369
  24. Dvorak R., 1986, A&A, 167, 379
  25. Eggl S., Pilat-Lohinger E., Georgakarakos N., Gyergyovits M., Funk B., 2012, ApJ, 752, 74
  26. Endl M., K�rster M., Els S., Hatzes A. P., Cochran W. D., 2001, Astronomy and Astrophysics, 374, 675
  27. Farrell B. F., 1990, Journal of Atmospheric Sciences, 47, 2986
  28. Forgan D., 2012, MNRAS, 422, 1241
  29. Forgan D., 2014, MNRAS, 437, 1352
  30. Forgan D. H., Mead A., Cockell C. S., Raven J. A., 2015, International Journal of Astrobiology, 14, 465
  31. G. Martin R., Armitage P. J., Alexander R. D., 2013, ApJ, 773, 74
  32. Haghighipour N., Kaltenegger L., 2013, ApJ, 777, 166
  33. Hart M., 1979, Icarus, 37, 351
  34. Hatzes A. P., 2013, ApJ, 770, 133
  35. 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
  36. Holman M. J., Wiegert P. A., 1999, The Astronomical Journal, 117, 621
  37. Huang S.-S., 1959, PASP, 71, 421
  38. Imkeller P., 2001, in Imkeller P., Von Storch J.-S., eds, , Stochastic Climate Models. Birkhäuser Basel, Basel, pp 213–240, doi:10.1007/978-3-0348-8287-3_9,{_}9
  39. Jaime L. G., Aguilar L., Pichardo B., 2014, MNRAS, 443, 260
  40. Kaltenegger L., Haghighipour N., 2013, ApJ, 777, 165
  41. Kane S. R., Hinkel N. R., 2013, ApJ, 762, 7
  42. Kasting J., Whitmire D., Reynolds R., 1993, Icarus, 101, 108
  43. Kopparapu R. K., et al., 2013, ApJ, 765, 131
  44. Kopparapu R. K., Ramirez R. M., SchottelKotte J., Kasting J. F., Domagal-Goldman S., Eymet V., 2014, ApJ, 787, L29
  45. Kostov V. B., et al., 2014, ApJ, 784, 14
  46. Laskar J., 1986a, A&A, 157, 59
  47. Laskar J., 1986b, A&A, 164, 437
  48. Lisiecki L. E., Raymo M. E., 2005, Paleoceanography, 20, id PA1003
  49. Marzari F., Thebault P., Scholl H., Picogna G., Baruteau C., 2013, Astronomy & Astrophysics, 553, A71
  50. Mason P. a., Zuluaga J. I., Clark J. M., Cuartas-Restrepo P. a., 2013, ApJ, 774, L26
  51. May E. M., Rauscher E., 2016, The Astrophysical Journal, 826, 225
  52. Meschiari S., 2014, ApJ, 790, 41
  53. Mikkola S., Aarseth S. J., 1990, Celestial Mechanics and Dynamical Astronomy (ISSN 0923-2958), 47, 375
  54. Mikkola S., Aarseth S. J., 1993, Celestial Mechanics & Dynamical Astronomy, 57, 439
  55. O’Malley-James J. T., Raven J. A., Cockell C. S., Greaves J. S., 2012, Astrobiology, 12, 115
  56. Orosz J. A., et al., 2012, Science (New York, N.Y.), 337, 1511
  57. Pichardo B., Sparke L. S., Aguilar L. A., 2005, MNRAS, 359, 521
  58. Pierrehumbert R. T., 2005, Journal of Geophysical Research, 110, D01111
  59. Pilat-Lohinger E., Dvorak R., 2002, Celestial Mechanics and Dynamical Astronomy, pp 143–153
  60. Queloz D., et al., 2000, A&A, 354, 99
  61. Quintana E. V., Lissauer J. J., Chambers J. E., Duncan M. J., 2002, ApJ, 576, 982
  62. Quintana E. V., Adams F. C., Lissauer J. J., Chambers J. E., 2007, ApJ, 660, 807
  63. Rafikov R. R., 2013, ApJ, 764, L16
  64. Rafikov R. R., Silsbee K., 2014a, ApJ, 798, 69
  65. Rafikov R. R., Silsbee K., 2014b, ApJ, 798, 70
  66. Raghavan D., et al., 2010, ApJSS, 190, 1
  67. Rajpaul V., Aigrain S., Roberts S., 2016, MNRAS, 456, L6
  68. Sigurdsson S., Richer H. B., Hansen B. M., Stairs I. H., Thorsett S. E., 2003, Science (New York, N.Y.), 301, 193
  69. Silsbee K., Rafikov R. R., 2015, ApJ, 798, 71
  70. Spiegel D. S., Menou K., Scharf C. A., 2008, ApJ, 681, 1609
  71. Spiegel D. S., Raymond S. N., Dressing C. D., Scharf C. A., Mitchell J. L., 2010, ApJ, 721, 1308
  72. Tajika E., 2008, ApJ, 680, L53
  73. Thébault P., Marzari F., Scholl H., 2008, MNRAS, 388, 1528
  74. Thébault P., Marzari F., Scholl H., 2009, MNRAS, 393, L21
  75. Thevenin F., Provost J., Morel P., Berthomieu G., Bouchy F., Carrier F., 2002, Astronomy and Astrophysics, 392, L9
  76. Thorsett S. E., Arzoumanian Z., Taylor J. H., 1993, ApJ, 412, L33
  77. Vladilo G., Murante G., Silva L., Provenzale A., Ferri G., Ragazzini G., 2013, ApJ, 767, 65
  78. Welsh W. F., et al., 2012, Nature, 481, 475
  79. Welsh W. F., et al., 2015, ApJ, 809, 26
  80. Wertheimer J. G., Laughlin G., 2006, The Astronomical Journal, 132, 1995
  81. Wiegert P. A., Holman M. J., 1997, The Astronomical Journal, 113, 1445
  82. Williams D., Kasting J., 1997, Icarus, 129, 254
  83. Xie J.-W., Zhou J.-L., Ge J., 2010, ApJ, 708, 1566
  84. Zachos J., Pagani M., Sloan L., Thomas E., Billups K., 2001, Science (New York, N.Y.), 292, 686
  85. Zucker S., Mazeh T., Santos N. C., Udry S., Mayor M., 2004, A&A, 426, 695
  86. Zuluaga J. I., Mason P. A., Cuartas-Restrepo P. A., 2016, The Astrophysical Journal, 818, 160
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description