# The Resilience of Habitable Climates Around Circumbinary Stars

Eric Wolf, Jacob Haqq-misra, Ravi Kopparapu, Thomas Fauchez, William

Welsh, Stephen Kane, Veselin Kostov, Siegfried Eggl

► To cite this version:

Eric Wolf, Jacob Haqq-misra, Ravi Kopparapu, Thomas Fauchez, William Welsh, et al.. The Resilience of Habitable Climates Around Circumbinary Stars. Journal of Geophysical Research. Planets, 2020, 125 (9), 10.1029/2020JE006576. [hal-03250313](https://hal.science/hal-03250313)

HAL Id: hal-03250313

https://hal.science/hal-03250313

Submitted on 5 Jun 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Planets

## RESEARCH ARTICLE

10.1029/2020JE006576

### Special Section:

Exoplanets: The Nexus of Astronomy and Geoscience

### Key Points:

* Earth-like planets in circumbinary systems are resilient against climate catastrophe, despite significant time-dependent variations in the incident stellar flux

* While ocean temperatures are not appreciably affected, land surface temperatures can exhibit an additional mode of circumbinary-seasonal variability

* Circumbinary systems should be considered as viable hosts for habitable worlds

### Correspondence to:

E. T. Wolf,
eric.wolf@colorado.edu

### Citation:

Wolf, E. T., Haqq-Misra, J., Kopparapu, R., Fauchez, T. J., Welsh, W. F., Kane, S. R., et al. (2020). The resilience of habitable climates around circumbinary stars. *Journal of Geophysical Research: Planets, 125*, e2020JE006576. https://doi.org/10.1029/2020JE006576

Received 15 JUN 2020
Accepted 18 AUG 2020
Accepted article online 25 AUG 2020

### Author Contributions:

**Conceptualization:** Eric T. Wolf, Jacob Haqq-Misra, Ravi Kopparapu, William F. Welsh, Stephen R. Kane, Veselin Kostov

**Data curation:** Thomas J. Fauchez

**Formal analysis:** Eric T. Wolf, Thomas J. Fauchez

**Funding acquisition:** Jacob Haqq-Misra

(continued)

©2020. The Authors.
This is an open access article under the
terms of the Creative Commons
Attribution-NonCommercial-NoDerivs
License, which permits use and distri-
bution in any medium, provided the
original work is properly cited, the use
is non-commercial and no modifica-
tions or adaptations are made.

# The Resilience of Habitable Climates Around Circumbinary Stars

<figure><img src="image_2.png" /></figure>

Eric T. Wolf¹,²,³ ID, Jacob Haqq-Misra²,⁴ ID, Ravi Kopparapu²,⁴,⁵ ID, Thomas J. Fauchez³,⁵,⁶ ID, William F. Welsh⁷ ID, Stephen R. Kane⁸ ID, Veselin Kostov⁵,⁹ ID, and Siegfried Eggl¹⁰,¹¹,¹² ID

¹Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO, USA, ²NASA NExSS Virtual Planetary Laboratory, Seattle, WA, USA, ³NASA Goddard Sellers Exoplanet Environments Collaboration, Greenbelt, MD, USA, ⁴Blue Marble Space Institute of Science, Seattle, WA, USA, ⁵NASA Goddard Space Flight Center, Greenbelt, MD, USA, ⁶Goddard Earth Sciences Technology and Research (GESTAR), Universities Space Research Association, Columbia, MD, USA, ⁷Department of Astronomy, San Diego State University, San Diego, CA, USA, ⁸Department of Earth and Planetary Sciences, University of California, Riverside, CA, USA, ⁹SETI Institute, Mountain View, CA, USA, ¹⁰Vera C. Rubin Observatory / DiRAC Institute, Department of Astronomy, University of Washington, Seattle, WA, USA, ¹¹IMCCE Observatoire de Paris, Université Lille 1, UPMC, Paris, France, ¹²Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, Urbana, IL, USA

**Abstract** Here we use a 3-D climate system model to study the habitability of Earth-like planets orbiting in circumbinary systems. In the most extreme cases, Earth-like planets in circumbinary systems could experience variations in the incident stellar flux of up to ~50% on ~100-day timescales. However, we find that Earth-like planets, having abundant surface liquid water, are generally effective at buffering against these time-dependent changes in the stellar irradiation due to the high thermal inertia of oceans compared with the relatively short periods of circumbinary-driven variations in the received stellar flux. Ocean surface temperatures exhibit little to no variation in time; however, land surfaces can experience modest changes in temperature, thus exhibiting an additional mode of climate variability driven by the circumbinary variations. Still, meaningful oscillations in surface temperatures are only found for circumbinary system architectures featuring the largest physically possible amplitudes in the stellar flux variation. In the most extreme cases, an Earth-like planet could experience circumbinary-driven variations in the global mean land surface temperature of up to ~5 K, and variations of local daytime maximum temperatures of up to ~12 K on seasonal timescales, while the global mean ocean temperatures vary by less than ~2 K. Such seasonal temperature swings over land areas could potentially pose adaptability challenges for extant life. Still, habitable planets in circumbinary systems appear to be remarkably resilient against circumbinary-driven climate variations and can avoid any true climate catastrophes.

**Plain Language Summary** A circumbinary system is one where a planet orbits around two stars simultaneously, resulting in large and rapid changes to both the stellar energy distribution and the total stellar energy received by the planet. Planets in circumbinary systems can experience significant variations in their received stellar energy over ~10- to ~100-day timescales. However, habitable planets in circumbinary systems are highly effective at buffering these changes in the received stellar energy due to the large thermal inertia of the oceans. Ocean surface temperatures exhibit little to no variation in time; however, land surfaces can experience modest changes in temperature. Even for the most extreme variations in stellar flux possible for circumbinary systems, Earth-like planets could avoid a permanent climate catastrophe, such as a runaway greenhouse or glaciation.

## 1. Introduction

Planets that orbit two stars at once are a common trope of science fiction. These are known as circumbinary systems, or Planet (P)-type systems (Dvorak, 1984). Indeed, stellar population statistics indicate that perhaps half of Sun-like stars exist in binary and higher multiple systems (Raghavan et al., 2010), and planets in circumbinary systems have already been discovered in transit light curves (Doyle et al., 2011; Kostov et al., 2014, 2020; Orosz et al., 2019; Welsh et al., 2012; Welsh & Orosz, 2018) despite the inherent difficulty in their observation due to significant transit timing variations (Martin, 2018). A number of previous works have considered the climate and habitability of planets in circumbinary systems using the results from
**Investigation:** Eric T. Wolf, Jacob Haqq-Misra, Ravi Kopparapu

**Methodology:** Eric T. Wolf, Jacob Haqq-Misra, William F. Welsh, Stephen R. Kane

**Project administration:** Jacob Haqq-Misra

**Software:** Eric T. Wolf, Siegfried Eggl

**Visualization:** Eric T. Wolf

**Writing - original draft:** Eric T. Wolf, Jacob Haqq-Misra

**Writing - review & editing:** Eric T. Wolf, Jacob Haqq-Misra, Ravi Kopparapu, Thomas J. Fauchez, William F. Welsh, Stephen R. Kane, Veselin Kostov, Siegfried Eggl

one-dimensional climate models. These have included both analytic calculations that weight the relative spectral contribution of two stars into the total irradiance received by planets (Haghighipour & Kaltenegger, 2013; Kane & Hinkel, 2012; Moorman et al., 2020; Wang & Cuntz, 2019a; Wang & Cuntz, 2019b), 1-D radiative-convective climate calculations that prognostically blend stellar spectra within their radiative transfer calculations (Cukier et al., 2019), and time-marching 1-D energy balance models that explicitly account for the time-dependent changes to the total incident stellar radiation (Forgan, 2014; Haqq-Misra et al., 2019; May & Rauscher, 2016). Additionally, two previous studies have explored the climate and atmospheres of planets in circumbinary systems using 3-D general circulation models (GCMs). May and Rauscher (2016) used an idealized GCM to study Kepler-47 b, assuming it to be a Neptune-like world. Popp and Eggl (2017) used a 3-D terrestrial planet climate system model to study theoretical Earth-sized and completely ocean covered planets in the Kepler-35 system. In both of these studies, it was shown that global mean temperatures are not significantly affected by periodic variations in the stellar instellation in these specific circumbinary systems. The high thermal inertia of both a thick Neptunian atmosphere, and of a completely ocean covered planet surface, are each highly effective at buffering time-dependent oscillations in the total stellar irradiance, resulting in variations in global mean temperatures of ~1.0% or less in these cases. However, both the Kepler-47 and Kepler-35 systems feature binary periods and amplitudes of stellar flux variation that are relatively small compared to what is physically possible.

To drive more meaningful variations in circumbinary planetary climate, the magnitude and period of stellar irradiance changes experienced by the planet must be maximized, and the thermal inertia of the surface/atmosphere system must be minimized (Haqq-Misra et al., 2019). The circumbinary periods used in May and Rauscher (2016) and Popp and Eggl (2017) are 7.448 and 22 Earth days respectively, and the normalized amplitude of variation in the incident stellar flux due to the circumbinary architecture were ~20% and ~15%, respectively. Here, we consider the normalized amplitude of the stellar flux variation as the maximum received flux minus the minimum received flux, divided by the time-averaged flux received by the planet over many orbits. Note that eclipsing binaries can cause drops in stellar irradiance by up to nearly ~50% during the eclipse if the stars are of equal radius and on a short period; however, eclipses last only a few hours, and thus, their time-integrated effect on the total stellar irradiance is actually quite small. Thus, eclipses are not generally important for driving changes to planetary climate nor weather. The scenarios studied by May and Rauscher (2016) and Popp and Eggl (2017), while based on known systems, represent circumbinary system architectures where the variation in stellar forcing is relatively modest compared to what is physically possible and dynamically stable. Using a simplified orbital calculation, which assumed circular orbits for all bodies, Haqq-Misra et al. (2019) found that for a planet receiving a time-averaged stellar flux equal to that of the modern Earth (1,360 W m⁻²) in a maximal scenario would experience a change in stellar flux of ~55% over a synodic period of ~150 Earth days. Note that this result is based on the orbital dynamical stability criteria of Holman and Wiegert (1999), who argued that such planets can remain dynamically stable in circumbinary systems so long as their orbital period is at least ~3 times larger than the circumbinary orbital period (see also Quarles et al., 2018). The latitudinal energy balance model results of Haqq-Misra et al. (2019) also showed that land surface characteristics and planetary obliquity can play key roles in the climate response to circumbinary forcings. Land surfaces heat up and cool down much more rapidly than does a deep optically thick atmosphere or an ocean surface, a consideration that was overlooked in previous GCM works. The coupling of obliquity-driven seasons like on Earth, with circumbinary-driven seasons due to the variation in received stellar flux, has the potential to generate more extreme local weather conditions when the warm seasons coincide.

While the global mean temperature is often the first consideration in whether or not a planet may be deemed habitable, seasonal and local extremes could pose challenges for the evolution and maintenance of surface life as we know it. Here, we revisit the problem of climate and habitability for terrestrial planets orbiting in circumbinary systems by using a GCM developed for Earth-like planets, while assuming modern Earth continental configurations, atmospheric composition, rotation rate, and time-averaged received stellar instellation (1,360 W m⁻²). First, we examine the full-phase space of potential circumbinary orbital configurations using the analytic orbital calculations of Georgakarakos and Eggl (2015) to identify system architectures that cause the most extreme variations in incident stellar flux received by the planet. Then we assess the impact on global, regional, and seasonal climates of these Earth-like planets in circumbinary systems using a 3-D climate system model. Ultimately, our results show that for physically plausible systems, land
<table><caption>**Table 1**
Stellar Characteristics Assumed in This Study</caption><thead><tr><th>#</th><th>Stellar type</th><th>Mass (M<sub>☉</sub>)</th><th>T<sub>eff</sub> (K)</th><th>Radius (R<sub>☉</sub>)</th><th>Luminosity (L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>1</td><td>G2V</td><td>1.0</td><td>5714</td><td>0.999</td><td>0.958</td></tr><tr><td>2</td><td>K0V</td><td>0.85</td><td>5229</td><td>0.815</td><td>0.448</td></tr><tr><td>3</td><td>M0V</td><td>0.55</td><td>3824</td><td>0.556</td><td>0.0578</td></tr><tr><td>4</td><td>M5V</td><td>0.15</td><td>3070</td><td>0.194</td><td>0.00298</td></tr></tbody></table>

Note. All circumbinary systems studied are a combination of these stars.

surface temperatures can exhibit circumbinary-driven seasonality, but global climate remains resilient against climatic catastrophe and habitability remains possible.

## 2. Methods

### 2.1. Description of the Circumbinary Treatment

To represent the variation in stellar irradiance received by planets in circumbinary systems, we use the analytical solutions for coplanar planets in circumbinary systems of Georgakarakos and Eggl (2015).

This method was also used to describe the Kepler-35 system in the 3-D climate simulations of Popp and Eggl (2017). We use the Georgakarakos and Eggl (2015) analytical circumbinary calculations first independently to explore what circumbinary system architectures are possible (section 3), and then fully integrated into our 3-D climate system model to explore the climatic effects for specific cases (section 4). Note we have made some minor modifications to the interfacing routine in order to suit our particular operational needs; however, the circumbinary system solutions are unaltered. As inputs we define the masses of the circumbinary stars and planet, the separation of the circumbinary stars, and the desired time-mean incident stellar flux received by the planet (1,360 W m⁻²). The modified code then solves for the planetary orbit that meets these conditions, if such an orbit is dynamically stable. Cases with unstable orbital dynamics are rejected. The Georgakarakos and Eggl (2015) model solves for the planetary orbit, and provides the incident stellar flux received by the planet from each star individually as a function of time, which includes the effect of eclipsing binaries, dual evolving zenith angles, and planetary eccentricity. The planet is assumed to start out on a circular orbit which evolves over time to an equilibrium state. The resulting short periodic and secular changes in the planetary eccentricity are tracked. Whether planets form on circular orbits or orbits with forced eccentricity is still a matter of debate (Mardling, 2007; Silsbee & Rafikov, 2014). Circumbinary planets on forced orbits tend to experience less variation in instellation (Eggl, 2018). In this work we have chosen initially circular orbits, as we aim to explore how planetary climates react to extreme variability. To describe the stars, we use mass, luminosity, radius, and effective temperature relationships derived from Pecaut et al. (2012) and Pecaut & Mamajek (2013). The choice of mass, luminosity, and radius relationship underpin the circumbinary solutions of Georgakarakos and Eggl (2015), which affect orbital properties, stellar fluxes at the planet, and magnitude of eclipse effects respectively. Note that using our derived mass, luminosity, radius, and effective temperature relationships, a 1 M☉ star has an effective temperature (T<sub>eff</sub>) of 5714 K, a radius of 0.99 R☉, and a luminosity of 0.958 L☉. Thus, our nominal solar twin is slightly cooler, smaller, and dimmer than our Sun. For this study, we have considered four different stars in shown in Table 1; G2V (1 M☉, 5714 K), K0V (0.85 M☉, 5229 K), M0V (0.55 M☉, 3824 K), and M5V (0.15 M☉, 3070 K).

### 2.2. Description of the 3-D GCM

Here we use a 3-D climate system model to simulate Earth-like planets in a variety of dynamically stable circumbinary systems. We use the Community Earth System Model Version 1.2.1 from the National Center for Atmospheric Research (Neale et al., 2010), along with the ExoCAM modeling package (https://github.com/storyofthewolf/ExoCAM), which facilitates simulations of extrasolar planets. ExoCAM is a modification to the Community Atmosphere Model version 4. As noted in section 2.1, we have further modified ExoCAM to incorporate the analytic circumbinary orbital system solutions of Georgakarakos and Eggl (2015), which drive the variations in stellar flux received by planets in circumbinary systems. While ExoCAM continues to evolve, the exact circumbinary version of ExoCAM used in this study has been frozen and made publicly available via Github and Zenodo (Wolf, 2020). For this study the 3-D model was further developed so that the light from each star is individually resolved in the radiative transfer computations, which accounts for stellar spectra and zenith angle of each star respectively. Thus, we are explicitly resolving the changes in received stellar flux driven by both stars simultaneously and self-consistently. In its present form, the circumbinary treatment increases the computational cost of our 3-D model by about a factor of 1.5 compared to the nominal single-star model due to increased burden in the radiation scheme. To model the stellar energy distribution, we use BT-Settl spectra for each star in the circumbinary pair (Allard et al., 2007) assuming [Fe/H] = 0 and log(g) = 5, and interpolated to desired T<sub>eff</sub> (Table 1). In all cases we assume that the planets receive a stellar flux of 1,360 W m⁻² averaged over many orbits (including the effect of eclipses). In this
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The stellar flux amplitudes possible for dynamically stable solutions for Earth-like planets in circumbinary systems derived from the Georgakarakos and Eggl (2015) model. The stellar flux amplitude is the percent change in the received stellar flux by the planet relative to the time-mean value (1,360 W m⁻²) and is indicated by contours and the color bar. White space indicates regions in our phase space where no stable solutions exist. Shown are solutions for circumbinary pairs having eccentricities of 0, 0.2, 0.4, and 0.6. The primary star is always assumed to be 1 M₀, while the secondary star mass is given on the x axis in units of stellar masses. The separation of the binary stars is given on the y axis in units of AU.</figcaption></figure>

study, one of the circumbinary stars is always assumed to be a Sun-like star with a mass of 1 M☉. The addition of a secondary star means that the combined intrinsic luminosity of the circumbinary pair is greater than our Sun alone, and thus, our simulated planets' semimajor axes and orbital periods both must be larger than Earth's in order to ensure the time-averaged received stellar flux remains 1,360 W m⁻².

We assume Earth-like planets that have the same mass, radius, surface gravity, continental configuration, and diurnal period as the Earth. The planet's absolute rotation period (i.e., the sidereal period) is modified slightly in each case to ensure that the diurnal period remains fixed at 24 hr, despite changes to the planet's orbital period in each case, but these variations are generally negligible and unlikely to yield any discernible effects on climate (Charnay et al., 2013; Wolf & Toon, 2014). We focus primarily on cases where the planet has 0° obliquity in order to isolate the effects driven by the circumbinary orbital configuration, but sensitivity tests are also conducted with a 23.5° obliquity. The planets are assumed to have a total atmospheric pressure of 1 bar, with 400 ppm of CO₂ and the remainder as N₂. We do not consider O₂, O₃, nor the effects of photo-chemistry in this study. The primary climatological contribution of O₂ is through its contribution to the total atmospheric pressure, which causes scattering and pressure broadening. While O₂ could have a significant effect on stratospheric temperatures, replacing O₂ with N₂ has a negligible effect on the surface climate (Wolf & Toon, 2013), which is our primary interest in this paper. Water vapor, liquid water clouds, and ice water clouds are variable and advected constituents in the model.

We assume the modern Earth continental configuration in all simulations. For treatment of the ground surface we use the Community Land Model (CLM) Version 4 (Oleson et al., 2010). CLM4 is a fully active land model, with multiple land surface, soil, and vegetation types simulated, along with active ground hydrology. The distribution of land surface, soil, and plant coverages are taken to be identical to the modern Earth. The ocean is treated as a 50-m-deep slab ocean. Implied ocean heat transport is included via heat flux
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The instellation synodic period in Earth days derived from the Georgakarakos and Eggl (2015) model. The instellation synodic period is the period of variation in circumbinary-driven instellation as viewed from the reference frame of the planet, and is indicated by contours and the color bar. The plot axis and conventions are identical to Figure 1.</figcaption></figure>

convergence terms to approximate that of the modern Earth (Bitz et al., 2012), where cold water is transported equatorward via eastern boundary currents (e.g., the Humboldt current and the California current), and warm water is transported poleward via western boundary currents (e.g., the Kuroshio current and the Gulf Stream). We have removed seasonal dependencies from the heat flux convergence terms and instead use the annual mean implied transport at all times. Thus, our ocean heat transport assumptions differ slightly from the modern Earth on a seasonal basis. Simulations are run with 4° × 5° horizontal resolution and 40 vertical layers extending from the surface up to ~1-mbar pressures, using a finite-volume dynamical core (Lin & Rood, 1996). Land, ocean, snow, and ice albedos are treated with two bands, visible and near infrared split at 0.77 µm. Open-ocean albedos are held fixed at 0.06 at all wavelengths. Sea ice albedos in visible (near-IR) are given as 0.68 (0.3), and snow albedos are given as 0.91 (0.63). Land albedos depend on the assumed surface type, surface color, and soil saturation. Permanent glacier ice has been removed from Antarctica and Greenland, however snow can accumulate and cover these regions. The specific heat capacity of sea water is 3,996 J kg⁻¹ K⁻¹, and the volumetric heat capacity of land surfaces varies from between 2 × 10⁶ and 2.5 × 10⁶ J m⁻³ K⁻¹ depending on the surface type (i.e., soil types and vegetation types). Surface types are assumed to be identical to the modern Earth. All simulations are initiated from a modern Earth climate state, but with the stratospheric inversion due to ozone heating removed. Simulations are considered equilibrated when the 5-year running mean top-of-atmosphere energy balance is less than ~1 W m⁻². Generally, radiative balance is first achieved in about 30 Earth years of simulation, and simulations are run for 60 Earth years in total.

### 3. Analytic Solutions for Earth-Like Worlds in Circumbinary Systems

First, we examine which specific circumbinary system architectures can harbor dynamically stable planets at modern Earth instellation levels, and which architectures can drive the largest magnitude and longest period changes to the incident stellar flux. To do so, we conduct an array of off-line calculations using our modified version of the Georgakarakos and Eggl (2015) model, solving for planets receiving 1,360 W m⁻² of stellar flux on time-average. Figures 1–3 describe the steady state properties of Earth-like planets in circumbinary
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The planet eccentricity derived from the Georgakarakos and Eggl (2015) model, indicated by contours and the color bar. The plot axis and conventions are identical to Figure 1.</figcaption></figure>

systems considering the stellar flux amplitude, the instellation synodic period, and the planet eccentricity respectively. In all figures, we assume that the primary star has a mass of 1 M☉, with the secondary star mass given on the x axis, and the separation of the binary stars in AU given on the y axis. Note we always refer to the lesser mass star as the secondary star. The time mean stellar flux is 1,360 W m⁻² in all cases. Binary eccentricities of 0, 0.2, 0.4, and 0.6 are shown and labeled in each panel. The white spaces at the top of each panel are regions in phase space where no dynamically stable solution exists for our Earth-like planets as described.

Figure 1 shows the stellar flux amplitude (%), defined as the maximum received flux by the planet at any time, minus the minimum received stellar flux, divided by the time-average flux (1,360 W m⁻²) and multiplied 100 to yield values in percentage. In this calculation we have excluded the effect of eclipses in the stellar flux amplitude, because the precipitous drop in stellar radiation during a binary eclipse only lasts for a few hours, and results in only a small time-integrated reduction in the total flux received by the planet with no apparent effect on planetary climate in our full climate calculations described in section 4. In Figure 1, note that the greatest stellar flux amplitude variation occurs when the secondary star has a mass of about half (~0.55 M☉) of the primary star. In this case, assuming the binary pair has zero eccentricity, a maximum stellar flux amplitude of ~50% can occur when the binary separation is 0.4 AU. The change in the stellar flux received by the planet is dominated by the motion of the primary star relative to the planet due to the circumbinary orbit. If the binary separation is held constant, a smaller and thus dimmer secondary star would result in reduced motion of the primary while contributing very little flux to the planet itself. The result is a reduction in the stellar flux amplitude. In the case of a larger and thus brighter secondary star, the secondary star can contribute meaningfully to the stellar flux received by the planet. When the secondary star moves toward its closest position to the planet, its stellar flux can begin to offset the reduction in stellar flux from the primary star at its furthest position from the planet. Intermediate mass secondary stars thus are a sweet spot, where they can induce significant enough motion of the primary star, without being too bright at as to offset stellar flux changes during the circumbinary orbit.

Naturally, if the binary separation is varied, larger separations yield greater stellar motions and correspondingly larger stellar flux amplitudes. However, if the binary separation becomes too wide, no dynamically
stable planet can exist in an orbit where it would receive 1,360 W m⁻² of time-averaged flux. Although stable planets could still exist on larger orbits, they would necessarily receive less incident stellar flux than the present-day Earth. Likewise, as the eccentricity of the circumbinary pair is increased, the phase space where an Earth-like planet can remain dynamically stable is reduced. The elongated orbits of eccentric circumbinary pairs can effectively scatter away surrounding planets, fewer permutations of stable planet orbits exist in our phase space, and the maximum possible amplitudes of stellar flux variation are at most ~35%.

As noted, the time period of incident stellar flux variation is also critical to the problem of circumbinary-driven climate variability. Figure 2 shows the instellation synodic period, which we define as the period of variation in the circumbinary-driven instellation when viewed in the reference frame of the planet. The instellation synodic period is greater than the binary orbital period, due to the motion of the planet relative to the circumbinary pair, much like the diurnal period of a planet is greater than its side-real period. For circumbinary pairs with zero eccentricity and separations greater than 0.3 AU, the instellation synodic period can exceed 75 days and reach as high as ~200 days in certain corner cases of our phase space. However, for higher eccentricity binaries, the instellation synodic period generally remains less than about 50 days, due to the constriction in the binary semimajor axes possible which permit dynamically stable Earth-like planets. In Figure 3, we show the planet eccentricities, which are direct output of the Georgakarakos and Eggl (2015) model. Planet eccentricities remain relatively small (<0.03) in all cases studied; however, the highest planet eccentricities predicted cluster around the more extreme circumbinary architectures with the widest binary separations. The planet's eccentricity introduces another mode of variation in the received stellar flux by the planet, and as we will see next, in some cases the stellar flux variations due to planetary eccentricity can sometimes outweigh that driven by the circumbinary pair.

In Figure 4 we show the time-dependent behavior of the global mean combined stellar flux from both stars impinging on our Earth-like planets, which includes the effect of binary eclipses, manifested as vertical straight lines, and the effect of planet eccentricities. The orange lines indicate the instantaneous stellar flux received by the planet tabulated once per hour. The black lines indicate the 30-day running average of the incident stellar flux. We present 12 different time-varying stellar flux patterns in Figure 4, corresponding to the system descriptions shown in Table 2. Figure 4 is labeled according to the cases listed in Table 2. Case 1 is our control case, where the planet orbits a single Sun-like star and the planet's orbital eccentricity is 0. Here, naturally, the global mean stellar flux received by the planet is constant in time, with a value of 340 W m⁻² (i.e., 1,360 W m⁻² divided by 4, accounting for diurnal and geometric effects). All other cases shown are for a variety of permutations of physically plausible circumbinary systems. While we are most interested in the cases where the variation in the stellar flux received by the planet is of the largest amplitude and longest period, we also show some cases where the variations in the instellation synodic period and flux amplitude are small. Note that periodic variations on order ~100 days or less are driven by the motion of the circumbinaries, while periodic variations on order of ~400 days are driven by the eccentricity of the planets. Recall, that the planet eccentricity in these circumbinary systems is not a free parameter but is a prediction outputted from Georgakarakos and Eggl (2015) model.

Case 7 illustrates a scenario of maximum amplitude stellar flux variation. This case features a G2V-M0V star pairing with a binary separation of 0.4 AU, a stellar flux amplitude of ~51%, and an instellation synodic period 97.3 days. In this case the amplitude of the stellar flux variation is dominated by the motion of the larger, brighter G2V star, which is ~18 times more luminous than the secondary M0V star (Table 1). Note also the subtle variation in the stellar flux maxima in Case 7, driven by the nonzero eccentricity of the planet. Case 8 shows the identical star pairing (G2V-M0V), but with the binary separation cut in half to 0.2 AU. In this case the stellar flux amplitude is reduced to ~25%, while the instellation synodic period is reduced to 28.7 days. In this case the 30-day running average variation in incident stellar flux is now dominated by the planet's eccentricity instead of the circumbinary effect. If we consider a smaller, dimmer, secondary star (M5V), the stellar flux amplitudes decrease further, because the motion of the larger, brighter relative to the center of mass decreases for a constant binary separation (e.g., Cases 9 and 10 vs. Cases 7 and 8, respectively).

Cases 2, 3, and 4 feature two equal mass stars (G2V-G2V pairings). In the case of two equal mass, equally luminous stars, the effective period of stellar flux variation at the planet is half that of the instellation synodic period. With each star being identical, the change in incident flux from each star is precisely balanced as each star approaches/recedes from the planet in tandem and at an equal velocity. The maximum incident
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Variations in the stellar flux received by our test planet for the 12 scenarios described in Table 2. Case 1 is our single-star control case, and the time-mean stellar flux in all cases is 1,360 W m⁻². The orange lines are the instantaneous flux received by the planet tabulated once per hour, while the black lines are the 30-day running average. Note that the vertical orange lines that fall to low stellar fluxes are caused by binary eclipses.</figcaption></figure>

flux occurs whenever one of the two stars is closest to the planet. The minimum incident flux occurs when the binary stars are at quadrature relative to the planet. For one complete orbit of the binary pair relative to the planet, thus, there are two cycles of instellation that occur. Similar as before, as the binary separation between the pair is reduced, the stellar flux amplitude and instellation synodic period decrease. Note in Case 2, and Case 3, that again the planet eccentricity imparts a significant signal on the total received instellation, rivaling the changes due to the circumbinary. Case 4 highlights when the effect of binary eclipsing is maximized. Note that the frequency of eclipses washes out the Case 4 panel in Figure 4. Here, the circumbinary stars orbit each other every 8.2 days, and eclipses occur every 8.3 days. Because the stars are of equal size, each eclipse reduces the total flux received by the planet by ~46%. However, due to the brevity of eclipses, still, the net reduction in the stellar flux is only 1.2% compared to an idealized architecture where no eclipses occur. Still, the instantaneous effects of eclipses are lost in the noise of weather and diurnal cycling, with eclipses causing a reduction in the surface downwelling stellar flux on par with that of a moderately cloudy day. However, for all other cases studied, eclipses are infrequent and the net flux reductions or no more than a few tenths of a percent.

When the stellar masses are different by less than about 20% (Case 5 and Case 6), an irregular pattern of stellar instellation results due to the 1/R² dependence of radiative flux combined with the ~M⁴ relationship of stellar luminosities, where R is distance between each star and planet respectively. Naturally, the stellar instellation is greatest when the larger star is closest to the planet; however, the lesser star can still contribute a meaningful fraction of the total irradiance (~10%). As the larger star moves away from the planet, the total instellation decreases. The stellar flux reaches a minimum at quadrature, but then begins to increase as the
<table><caption>**Table 2**

Circumbinary System Characteristics Used in This Study</caption><thead><tr><th>#</th><th>Stellar pair</th><th>Binary separation (AU)</th><th>Binary eccentricity</th><th>Binary period (days)</th><th>Planet period (days)</th><th>Planet eccentricity</th><th>Insolation synodic period (days)</th><th>Planet semimajor axis (AU)</th><th>Stellar flux amplitude (%)</th><th>Net eclipse flux reduction (%)</th></tr></thead><tbody><tr><td>1</td><td>G2V</td><td>—</td><td>0</td><td>—</td><td>353.9</td><td>0</td><td>—</td><td>0.979</td><td>—</td><td>—</td></tr><tr><td>2</td><td>G2V-G2V</td><td>0.55</td><td>0</td><td>105.3</td><td>449.8</td><td>2.17e-2</td><td>137.5</td><td>1.448</td><td>16.06</td><td>0.174</td></tr><tr><td>3</td><td>G2V-G2V</td><td>0.3</td><td>0</td><td>42.4</td><td>428.9</td><td>1.48e-2</td><td>47.1</td><td>1.402</td><td>6.77</td><td>0.372</td></tr><tr><td>4</td><td>G2V-G2V</td><td>0.1</td><td>0</td><td>8.2</td><td>418.1</td><td>8.54e-5</td><td>8.3</td><td>1.378</td><td>0.89</td><td>1.224</td></tr><tr><td>5</td><td>G2V-K0V</td><td>0.45</td><td>0</td><td>81.1</td><td>368.3</td><td>5.00e-3</td><td>103.9</td><td>1.235</td><td>28.35</td><td>0.179</td></tr><tr><td>6</td><td>G2V-K0V</td><td>0.25</td><td>0</td><td>33.6</td><td>353.0</td><td>1.08e-2</td><td>37.1</td><td>1.200</td><td>14.63</td><td>0.371</td></tr><tr><td>7</td><td>G2V-M0V</td><td>0.4</td><td>0</td><td>74.2</td><td>313.2</td><td>5.85e-3</td><td>97.3</td><td>1.045</td><td>51.29</td><td>0.099</td></tr><tr><td>8</td><td>G2V-M0V</td><td>0.2</td><td>0</td><td>26.2</td><td>300.8</td><td>3.52e-3</td><td>28.7</td><td>1.017</td><td>25.80</td><td>0.227</td></tr><tr><td>9</td><td>G2V-M5V</td><td>0.45</td><td>0</td><td>102.8</td><td>340.5</td><td>2.18e-2</td><td>147.3</td><td>0.999</td><td>25.24</td><td>0.012</td></tr><tr><td>10</td><td>G2V-M5V</td><td>0.25</td><td>0</td><td>42.6</td><td>333.8</td><td>4.85e-3</td><td>48.8</td><td>0.987</td><td>14.55</td><td>0.022</td></tr><tr><td>11</td><td>G2V-G2V</td><td>0.35</td><td>0.6</td><td>53.5</td><td>439.4</td><td>7.74e-2</td><td>60.9</td><td>1.425</td><td>32.81</td><td>0.305</td></tr><tr><td>12</td><td>G2V-K0V</td><td>0.3</td><td>0.4</td><td>44.1</td><td>358.7</td><td>4.34e-2</td><td>50.3</td><td>1.213</td><td>28.80</td><td>0.298</td></tr></tbody></table>

lesser star moves closer to the planet. There is a secondary maximum that occurs when the secondary star is closest to the planet, but because the secondary is dimmer, this stellar flux maximum is lesser. At this phase, the brighter star still contributes to the largest fraction of stellar radiation, but now the secondary contributes up to 40%.

Finally, we also illustrate two scenarios (Cases 11 and 12) where the circumbinary pair has significant eccentricity. Note that again we have cherry-picked system architectures where the variation in stellar flux at the planet is largest. Here, complex interactions occur between star and planet positions that yield complex time-dependent patterns of stellar instellation. Note that in both Cases 11 and 12, the variation in incident flux is dominated by the planet's eccentricity, with the circumbinary signal laid over top. Note that Cases 11 and 12 are corner cases, located on the stability line in Figures 1–3. Cases 11 and 12 have planet eccentricities of 0.0774 and 0.0434, respectively, which are both considerably larger than that of the modern Earth (0.0167), enhancing eccentricity driven stellar flux variations over the course of its year (Figure 4). Note that all other cases have planet eccentricities which are on par with or considerably less than that of Earth (Table 2). While the stellar flux variations in Cases 11 and 12 (and others) can be dominated by the planet eccentricity rather than the circumbinary motions, recall that the nonzero eccentricity of these planets is a direct consequence of the complex orbital dynamical interactions present in circumbinary systems and thus must be accounted for in any self-consistent circumbinary planets residing in the habitable zone (Georgakarakos & Eggl, 2015; Popp & Eggl, 2017).

## 4. Three-Dimensional Climate Modeling Results

Next, we conducted a variety of 3-D climate simulations to illustrate the possible changes to climate and weather for Earth-similar planets in circumbinary systems. Table 2 shows a list of the circumbinary system architectures that were explored, along with relevant statistics that describe the systems. We explore 12 unique system configurations, with stellar flux evolution shown in Figure 4. For our primary set of simulations, we consider Earth-like planets with zero obliquity. Recall that in all cases the time-averaged stellar flux received by the planet is 1,360 W m⁻². Figure 5 illustrates surface temperature variations as a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The temperatures variations corresponding to our twelve system configurations, described in Table 2 and with stellar flux variations shown in Figure 4. We show the instantaneous maximum local temperature (red), the land surface mean (green), and the ocean mean temperature. Black lines indicated the 30-day running mean, while the colors show the instantaneous values at a 1-hr time cadence.</figcaption></figure>

function of time for Earth-like planets with zero obliquity. Shown are the maximum local temperature taken from anywhere on the planet (red), the global mean ocean temperatures, and the global mean land surface mean temperature (green). The solid black lines indicate the 30-day running mean temperature for each, while the colors show the instantaneous values sampled once every hour. Note that the maximum local temperature is defined as the maximum temperature measured at any location on the planet at a given time, and all temperatures shown are surface temperatures. Recall that Case 1 is the single-star control case. Due to the significant thermal inertia of the oceans, the instantaneous variations in ocean temperatures are virtually indistinguishable from the 30-day running mean value, and even under the most extreme circumbinary architectures ocean mean temperatures do not vary more than a few degrees kelvin, whether driven by circumbinary forcings (e.g., Case 7 and Case 9) or via planet eccentricity effects (Case 11 and Case 12). Thus, for ocean dominated planets in the habitable zones of circumbinary stars, even in the most extreme physically plausible cases, dramatic global changes in climate do not occur due to the thermal inertia of the oceans, combined with the fact that instellation synodic periods (Figure 2) are generally no more than ~150 Earth days. Thus, the ocean does not have sufficient time to catastrophically cool or warm. This result agrees with the results of Popp and Eggl (2017) and Way and Georgakarakos (2017). For evaluating the global mean climate of habitable circumbinary planets, the mean-flux approximation appears to hold true (Bolmont et al., 2016).

Greater variations in the land surface temperatures compared to ocean temperatures are evident in Figure 5, which agrees with the energy balance modeling of Haqq-Misra et al. (2019). The relatively low thermal inertia of land surfaces allows continents to heat up and cool off much more quickly and thus signals of the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The time evolution of latitudinal temperature distributions from four different longitudes around the planet for Case 7. Plotted are the maximum daily temperatures recorded at each location over each 24-hr period. Ocean surface temperatures remain generally stable. Tropical and subtropical land surfaces experience the most extreme circumbinary-driven warming episodes, while higher latitude continents fair marginally better.</figcaption></figure>

circumbinary forcing and the planet eccentricity are more apparent. Variations to land surface mean and local maximum temperatures are evident for several cases. Cases 5, 7, and 9 indicate a relatively strong circumbinary signal in land surface temperatures, implying a circumbinary pseudoseasonality is induced in these systems. In our most extreme case, Case 7, the 30-day running mean land surface temperatures may vary by ~5 K between maximum and minimum instellation phases of the circumbinary motion, while the 30-day running mean maximum local temperatures may vary by as much as ~12 K. If we consider the instantaneous maximum temperatures, in Case 7 maximum instellation phases may drive local land surface temperatures in the tropics to reach ~345 K, whereas in our single-star control simulation, maximum local temperatures generally hover near 330 K. Note that these extreme local instantaneous temperatures are also seen in Case 11 and Case 12, although these temperature maximums are driven primarily by the planet eccentricity and not from the circumbinary forcing. Still, note that we are highlighting cases of maximum plausible circumbinary-driven variations in stellar flux. Even in such cases of maximal stellar forcing, the heat stress thresholds for human habitability based on wet bulb temperatures (Sherwood & Huber, 2010) are not exceeded during warm periods because the maximum temperatures indicated occur over regions of low relative humidity (i.e., subtropical desert areas). For more moderate system architectures (e.g., Cases 3, 4, 6, and 8) the effects of local climate are muted or nearly indistinguishable from the single-star control case.

In Figure 6, we examine the time evolution of the latitudinal distribution of the maximum daily surface temperature for Case 7. The maximum daily surface temperature is defined as the highest instantaneous temperature attained at a given location during a 24-hr period. Previous circumbinary climate studies have constructed similar time evolution and latitude plots based on the results of 1-D energy balance models
<figure><img src="image_9.png" /><figcaption>**Figure 7.** This illustrates that there is a subtle phase offset between cloud water (middle row) and the variations in the land surface temperatures (top row). While the TOA albedo (bottom row), naturally remains in phase with the cloud water, clouds and albedo subtly lag the temperature maximum. Still, elevated clouds and albedo during warm phases act as a negative feedback, stabilizing the climate against circumbinary-driven effects.</figcaption></figure>

(Haqq-Misra et al., 2019; May & Rauscher, 2016). However, advancing from these previous works, here we also have important longitudinal contrasts in the surface properties (i.e., land/ocean contrasts), which must be illustrated to yield a more complete view of the circumbinary-driven temperature variations that are occurring. Thus, in Figure 6, we plot the maximum daily surface temperatures as a function of time at four different longitudes; 180°, which transects the center of the Pacific Ocean; 105°W, which intersects Canada and the western United States; 20°E, which crosses central Europe and Africa; and finally, 60°W, which crosses the heart of South America.

Reinforcing the results from Figure 4, surface temperatures across the Pacific Ocean show very little variation over time, even in this case of maximum circumbinary-driven instellation change. As expected, maximum daily surface temperatures swing more greatly over the land surface areas. The most extreme surface temperature variations occur in the tropics and subtropics. In this scenario, Africa and South America, which predominantly occupy the tropics and subtropics, experience the most extreme circumbinary-driven periodic heating events. North America and Europe fare generally better, with the amplitude of the daily maximum surface temperature changes being relatively muted. Considering Case 7 with zero-obliquity, life forms would be best suited to migrate poleward to seek more stable seasonal climates.

In Figure 7 we illustrate the land surface mean temperature, along with the global mean cloud water column, and the TOA albedo, for Case 7 and Case 9 over several instellation synodic periods. We see that there is a phase offset between the land surface temperatures, and the cloud water and albedo respectively. Cloud water and albedo remain in phase due to the strong contribution of clouds to modulating the planetary albedo. However, the maximum in cloud water and the planetary albedo lag behind the land surface temperature maximum. Indeed, the warmer the air is, the more water vapor it can hold before condensing into cloud droplets. As the planet temperature warms, more water vapor is evaporated into the atmosphere. When temperatures begin to cool, after the instellation maximum has passed, the cooling temperatures decrease the ability of the atmosphere to hold water vapor, and so water condenses to form clouds. As this process is not instantaneous a slight offset is visible between the minimum cloud water and maximum temperature. Still, the general trend of rising cloud water and albedo at and just past the instellation peak represent a negative feedback at work that acts to mute circumbinary-driven climate changes.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The effect of nonzero obliquity. Recall Case 1 is the single-star control case, while Case 7 exhibits the most extreme circumbinary variation. For all but the most extreme circumbinary scenarios, seasonal effects from nonzero obliquity coupled with asymmetric north-south land distributions dominate over the circumbinary effects.</figcaption></figure>

We have also conducted sensitivity tests considering planets with a 23.5° obliquity. In the majority of our cases, climate changes driven by obliquity seasons completely swamp any circumbinary-driven climate signals. In Figure 8 we show the maximum local, ocean mean, and land surface mean temperatures for Case 1 (single-star control) and Case 7 (our most extreme circumbinary architecture), while assuming the planet has a 23.5° obliquity. Note that for the control case, seasonal variations in land surface temperatures are a result of the asymmetric distribution of continents on Earth. Most of the land mass is located in the Northern Hemisphere, land surfaces have a lower thermal inertia, and thus, the Northern Hemisphere summer exhibits warmer temperatures. When nonzero obliquity is coupled with an extreme circumbinary forcing, temperature variations can be either strengthened or weakened depending on the alignment of the eccentricity seasons with the circumbinary variations. When Northern Hemisphere summer aligns with a circumbinary maximum instellation phase, then land surface temperatures can be more extreme. However, if conversely the circumbinary minimum occurs during the Northern Hemisphere summer, climatological effects are muted. Because the planet orbital periods do not divide evenly into the circumbinary periods, there may be longer term trends not explored here as the mutual planet and circumbinary orbits may only precisely repeat on many decades to century scale timescales.

**5. Discussion and Conclusions**

Our results indicate that there is no climate catastrophe that precludes the habitability of Earth-like planets in the HZs of circumbinary systems. Earth-like planets in dynamically stable orbits in circumbinary systems are as likely to remain as habitable as planets in single-star systems. Circumbinary-driven variations in surface climate are insignificant in all but the most extreme system architectures. In the most extreme physically plausible circumbinary system architectures, land surface temperatures can exhibit an additional mode of variability with predictable periods of more intense heat waves; however, the thermal inertia of habitable planets (e.g., ocean dominated) is such that radical climate changes are prevented. Aspects such as the stochastic nature of planet formation, the acquisition and availability of volatiles, the orbital dynamical stability of planets, and the appropriate balance of time-averaged stellar flux and atmospheric greenhouse effect, are all more important arbiters of habitability than are periodic variations in stellar instellation found in circumbinary systems.
There are several caveats that should come with this conclusion. Here, we have limited our study to Earth-twin planets, and we have not explored the effects of other planetary parameter variations such as planetary rotation period, orbital eccentricity, surface pressure, and atmospheric composition. For instance, desiccated planets with exceedingly thin atmospheres and no oceans would be expected to show greater degree of temperature variations compared to the Earth-like planets studied here. Conversely, planets with thicker atmospheres, or larger fractions of ocean compared to land, should experience even less climate variability than is shown here (May & Rauscher, 2016). Results from slow rotating planets indicate that their climate sensitivity to changing stellar fluxes is much less than for rapidly rotating planets (Yang et al., 2013; Way et al., 2018), and thus, we do not expect changes in the planetary rotation rate to significantly change the resiliency of habitable climates to circumbinary forcings. Furthermore, larger assumed planetary eccentricity or obliquity would swamp any circumbinary signals, as is already evident from Figure 5 (Cases 11 and 12) and Figure 8. Alternatively, precise orbital resonances between eccentricity, obliquity, and circumbinary seasonality could exist that lead to more extreme climate variations in specific cases. Other possibilities exist, such as the circumbinary stellar flux variations lining up precisely with modes of internal climate oscillations, such as the Madden-Julian Oscillation, causing either an amplification or damping of internal variability depending on their relative phases. If we find that HZ planets are as common in circumbinary systems as in single-star systems, we can conclude that they have an equal chance of hosting life. Considering that perhaps half of Sun-like stars exist in binary and higher multiple systems, the resilience of habitable worlds in these systems may increase the chances of life-hosting planets existing in our galaxy.

The prospect of finding habitable Earth-like planets in circumbinary systems raises further speculation as to how a circumbinary biosphere might differ from a single-star one. Several studies have suggested that the trajectory of biological evolution on Earth has been driven in part by orbital Milankovitch cycles and the resulting changes in regional instellation and glaciation (Bennett, 1990; Jansson & Dynesius, 2002; Lister, 2004). The “seasonality hypothesis” even suggests that the evolution of traits such as animal body size could be connected with environmental pressures associated with climate cycles (Troost et al., 2009). Such ideas suggest that the effect of circumbinary seasons, as well as circumbinary Milankovitch forcing (Forgan, 2016), could similarly act as an evolutionary selection mechanism. Spectral biosignatures on circumbinary planets could likewise exhibit unique features that depend upon such variations; for example, the observability of a planet’s photosynthetic “red edge” (Kiang et al., 2007; Seager et al., 2005; Turnbull et al., 2006) could depend upon the phase of its circumbinary seasons. Seasonal and orbital variations provide a selection pressure for life on Earth, so we can reasonably expect that such pressures should apply to any life that develops on circumbinary planets.

We can even speculate further as to the implications of habitable circumbinary planets for the development of technology. DeVito (2013) discusses the connection between human mathematical systems and our physical spatial environment; for example, the model of natural numbers could have been extrapolated from observing the consistency of the day-night cycle on Earth. The period of instellation experienced by circumbinary planets poses a more complex pattern than a single-star system: The inference of natural numbers from the setting stars might be more difficult in such a system, although perhaps such added complexity would enhance the evolution of technology. The “proportional evolutionary time” hypothesis (Haqq-Misra, 2019) suggests that the evolutionary timescale for complex life depends upon the availability of free energy between 200 and 1,000 nm. The incident radiation for an Earth-like circumbinary planet with different spectral host stars (e.g., a G-dwarf primary and K-dwarf secondary) will show a similar time-mean irradiance within this wavelength range compared to a planet orbiting a single-star at an equivalent orbital distance. Thus, the proportional evolutionary time hypothesis suggests that circumbinary systems with a primary star of a given mass *M* are at least as promising candidates as single-star systems of stellar mass *M* to search for evidence of complex life or technology. Finally, the possibility that life, and even technology, could exist around circumbinary systems raises an important philosophical question: “Why do we find ourselves around a single star instead of a binary pair?” If binary systems are as likely to support life, then this might give us reason to suspect that the circumstances of life on Earth are less typical than we might otherwise assume (e.g., Haqq-Misra et al., 2018). Conversely, binary systems may be atypical in ways not explored in this manuscript, precluding the formation of dynamical stable and volatile rich worlds. The answer to these questions, and many of our other speculations, will require detection and spectral characterization of a
statistically significant sample of circumbinary planets in order to begin estimating the prevalence of bio-
signatures and technosignatures in such systems.

Conflict of Interest

No competing financial interests exist.

Data Availability Statement

Computing resources for this project were provided by the NASA Center for Climate Simulation and used
the Discover supercomputer at Goddard Space Flight Center (http://www.nccs.nasa.gov). The core model
used in this study, CESM1.2.1, is freely available to all at http://www.cesm.ucar.edu/models/cesm1.2/ via
the National Center for Atmospheric Research, Boulder, CO. The specific code changes we have made to
the model to enable modeling circumbinary planets, along with necessary initial files for implementation
into CESM1.2.1, have been frozen and released on the lead authors Github page. See https://github.com/
storyofthewolf/ExoCAM/releases/tag/circumbinary_jgrplanets_frozen for the circumbinary release, or on
Zenodo at http://doi.org/10.5281/zenodo.3968971 (Wolf, 2020). For more information on exoplanet model-
ing with CESM1.2.1, see https://github.com/storyofthewolf/ExoCAM/ and https://github.com/storyofthe-
wolf/ExoRT websites. The primary model outputs used in this study for analysis are multiyear high time
cadence outputs and have been made publicly available on Zenodo at http://doi.org/10.5281/
zenodo.3964961 (Wolf et al., 2020a), http://doi.org/10.5281/zenodo.3972021 (Wolf et al., 2020b), and also
available at the lead authors data repository through The Internet Archive (see https://archive.org/details/
Circumbinary_Climates).

Acknowledgments

The authors gratefully acknowledge funding from NASA Habitable Worlds program under Award 80NSSC17K0741. J. H., E. T. W., and R. K. K. also acknowledge funding from the Virtual Planetary Laboratory under NASA Award 80NSSC18K0829, and R. K. K. and V. K. acknowledge funding from the Sellers Exoplanet Environments Collaboration. W. F. W. gratefully thanks John Hood Jr. for his support of exoplanet research at SDSU. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NASA. We thank Dorian Abbot and Tad Komacek for constructive reviews which improved the quality and clarity of this manuscript.

References

Allard, F., Allard, N. F., Homeier, D., Kielkopf, J., McCaughrean, M. J., & Spiegelman, F. (2007). K-H₂ quasi-molecular absorption detected in the T-dwarf *ε* Indi Ba. *Astronomy & Astrophysics*, *474*(2), L21–L24. https://doi.org/10.1051/0004-6361:20078362

Bennett, K. D. (1990). Milankovitch cycles and their effects on species in ecological and evolutionary time. *Paleobiology*, *16*(1), 11–21.
https://doi.org/10.1017/S0094837300009684

Bitz, C. M., Shell, K. M., Gent, P. R., Bailey, D. A., Danabasoglu, G., Armour, K. C., et al. (2012). Climate sensitivity of the community climate system model, version 4. *Journal of Climate*, *25*(9), 3053–3070. https://doi.org/10.1175/JCLI-D-11-00290.1

Bolmont, E., Libert, A.-S., Leconte, J., & Selsis, F. (2016). Habitability of planets on eccentric orbits: The limits of the mean flux approximation. *Astronomy & Astrophysics*, *591*, A106. https://doi.org/10.1051/0004-6361/201628073

Charnay, B., Forget, F., Wordsworth, R., Leconte, J., & Millour, E. (2013). Exploring the faint young Sun problem and the possible climates of the Archean Earth with a 3-D GCM. *Journal of Geophysical Research: Atmospheres*, *118*, 10,414–10,431. https://doi.org/10.1002/jgrd.50808

Cukier, W., Kopparpau, R. K., Kane, S. R., Welsh, W., Wolf, E., Kostov, V., & Haqq-Misra, J. (2019). Habitable zone boundaries for circumbinary planets. *PASP*, *131*(1006), 124402 (5pp). https://doi.org/10.1088/1538-3873/ab50cb

DeVito, C. L. (2013). *Science*, *SETI*, and *mathematics*. New York, NY: Berghahn Books.

Doyle, L. R., Carter, J. A., Fabrycky, D. C., Slawson, R. W., Howell, S. B., Winn, J. N., et al. (2011). Kepler-16: A transiting circumbinary planet. *Science*, *333*(6049), 1602–1606. https://doi.org/10.1126/science.1210923

Dvorak, R. (1984). Numerical experiments on planetary orbits in double stars. In *The stability of planetary systems* (pp. 369–378). Dordrecht: Springer. https://doi.org/10.1007/978-94-009-5331-4_31

Eggl, S. (2018). Habitability of planets in binary star systems. In *Handbook of exoplanets* (Chap. 141, pp. 3041–3067). Switzerland: Springer international publishing AG. ISBN 978-3-319-55332-0

Forgan, D. (2014). Assessing circumbinary habitable zones using latitudinal energy balance modelling. *Monthly Notices of the Royal Astronomical Society*, *437*(2), 1352–1361. https://doi.org/10.1093/mnras/stt1964

Forgan, D. (2016). Milankovitch cycles of terrestrial planets in binary star systems. *Monthly Notices of the Royal Astronomical Society*, *463*(3), 2768–2780. https://doi.org/10.1093/mnras/stw2098

Georgakarakos, N., & Eggl, S. (2015). Analytic orbit propagation for transiting circumbinary planets. *The Astrophysical Journal*, *802*(2), 94 (16pp). https://doi.org/10.1088/0004-637X/802/2/94

Haghighipour, N., & Kaltenegger, L. (2013). Calculating the habitable zone of binary star systems. II. *P*-type binaries. *The Astrophysical Journal*, *777*(2), 166.

Haqq-Misra, J. (2019). Does the evolution of complex life depend on the stellar spectral energy distribution? *Astrobiology*, *19*(10), 1292–1299. https://doi.org/10.1089/ast.2018.1946

Haqq-Misra, J., Kopparapu, R. K., & Wolf, E. T. (2018). Why do we find ourselves around a yellow star instead of a red star? *International Journal of Astrobiology*, *17*(1), 77–86. https://doi.org/10.1017/S1473550417000118

Haqq-Misra, J., Wolf, E. T., Welsh, W. F., Kopparapu, R. K., Kostov, V., & Kane, S. R. (2019). Constraining the magnitude of climate extremes from time-varying instellation on a circumbinary terrestrial planet. *Journal of Geophysical Research: Planets*, *124*, 3231–3243. https://doi.org/10.1029/2019JE006222

Holman, M. J., & Wiegert, P. A. (1999). Long-term stability of planets in binary systems. *The Astronomical Journal*, *117*(1), 621–628. https://doi.org/10.1086/300695

Jansson, R., & Dynesius, M. (2002). The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution. *Annual Review of Ecology and Systematics*, *33*(1), 741–777. https://doi.org/10.1146/annurev.ecolsys.33.010802.150520
Kane, S. R., & Hinkel, N. R. (2012). On the habitable zones of circumbinary planetary systems. *The Astrophysical Journal*, **762**(1), 7.

Kiang, N. Y., Siefert, J., & Blankenship, R. E. (2007). Spectral signatures of photosynthesis. *I. Review of Earth organisms*. *Astrobiology*, *7*, 222–251.

Kostov, V. B., McCullough, P. R., Carter, J. A., Deleuil, M., Díaz, R. F., Fabrycky, D. C., et al. (2014). Kepler-413b: A slightly misaligned, Neptune-size transiting circumbinary planet. *The Astrophysical Journal*, **784**(1), 14. https://doi.org/10.1088/0004-637X/784/1/14

Kostov, V. B., Orosz, J. A., Feinstein, A. D., Welsh, W. F., Cukier, W., Haghighipour, N., et al. (2020). TOI-1338: TESS' first transiting circumbinary planet. *The Astronomical Journal*, **159**(6), 253 (26pp). https://doi.org/10.3847/1538-3881/ab8a48

Lin, S. J., & Rood, R. B. (1996). Multidimensional flux-form semi-Lagrangian transport schemes. *Monthly Weather Review*, **124**(9), 2046–2070. https://doi.org/10.1175/1520-0493(1996)124<2046:MFFSLT>2.0.CO;2

Lister, A. M. (2004). The impact of quaternary ice ages on mammalian evolution. *Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences*, **359**(1442), 221–241.

Mardling, R. A. (2007). Long-term tidal evolution of short-period planets with companions. *Monthly Notices of the Royal Astronomical Society*, **382**(4), 1768–1790. https://doi.org/10.1111/j.1365-2966.2007.12500.x

Martin, D. V. (2018). Populations of planets in multiple star systems. In H. Deeg, & J. Belmonte (Eds.), *Handbook of exoplanets* (Chap. 100, pp. 2035–2060). Cham: Springer. https://doi.org/10.1007/978-3-319-55333-7_156

May, E., & Rauscher, E. (2016). Examining tatoome: Atmospheric models of Neptune-like circumbinary planets. *The Astrophysical Journal*, **826**(2), 225. https://doi.org/10.3847/0004-637X/826/2/225

Moorman, S. Y., Wang, Z., & Cuntz, M. (2020). Dead zones of classical habitability in stellar binary systems. *Astrophysics and Space Science*, **365**(1), 1–13.

Neale, R. B., et al. (2010). Description of the NCAR Community atmosphere model (CAM 4.0). *NCAR/TN-486+STR NCAR TECHNICAL NOTE*.

Oleson, K. W., et al. (2010). Technical description of version 4.0 of the Community Land Model (CLM). *NCAR/TN-478+STR NCAR TECHNICAL NOTE*.

Orosz, J. A., Welsh, W. F., Haghighipour, N., Quarles, B., Short, D. R., Mills, S. M., et al. (2019). Discovery of a third transiting planet in the kepler-47 circumbinary system. *The Astronomical Journal*, **157**(5), 174. https://doi.org/10.3847/1538-3881/ab0ca0

Pecaut, M. J., & Mamajek, E. E. (2013). Intrinsic colors, temperatures, and bolometric corrections of pre-main-sequence stars. *The Astrophysical Journal Supplement Series*, **208**(1), 9. https://doi.org/10.1088/0067-0049/208/1/9

Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. (2012). A revised age for upper Scorpius and the star formation history among the *f*-type members of the Scorpius-Centaurus ob Association. *The Astrophysical Journal*, **746**(2), 154. https://doi.org/10.1088/0004-637X/746/2/154

Popp, M., & Eggl, S. (2017). Climate variations on Earth-like circumbinary planets. *Nature Communications*, **8**(1), 14957. https://doi.org/10.1038/ncomms14957

Quarles, B., Satyal, S., Kostov, V., Kaib, N., & Haghighipour, N. (2018). Stability limits of Circumbinary planets: Is there a pile-up in the *Kepler* CBPs? *The Astrophysical Journal*, **856**(2), 150. https://doi.org/10.3847/1538-4357/aab264

Raghavan, D., McAlister, H. A., Henry, T. J., Latham, D. W., Marcy, G. W., Mason, B. D., et al. (2010). A survey of stellar families: Multiplicity of solar-type stars. *The Astrophysical Journal Supplement Series*, **190**(1), 1–42. https://doi.org/10.1088/0067-0049/190/1/1

Seager, S., Turner, E. L., Schafer, J., & Ford, E. B. (2005). Vegetation's red edge: A possible spectroscopic biosignature of extraterrestrial plants. *Astrobiology*, **5**(3), 372–390. https://doi.org/10.1089/ast.2005.5.372

Sherwood, S. C., & Huber, M. (2010). An adaptability limit to climate change due to heat stress. *PNAS*, **107**(21), 9552–9555. https://doi.org/10.1073/pnas.0913352107

Silsbee, K., & Rafikov, R. R. (2014). Planet formation in binaries: Dynamics of planetesimals perturbed by the eccentric protoplanetary disk and the secondary. *The Astrophysical Journal*, **798**(2), 71. https://doi.org/10.1088/0004-637X/798/2/71

Troost, T. A., van Dam, J. A., Kooi, B. W., & Tuenter, E. (2009). Seasonality, climate cycles and body size evolution. *Mathematical Modelling of Natural Phenomena*, **4**(6), 135–155. https://doi.org/10.1051/mmnp/20094605

Turnbull, M. C., Traub, W. A., Jucks, K. W., Woolf, N. J., Meyer, M. R., Gorlova, N., et al. (2006). Spectrum of a habitable world: Earthshine in the near-infrared. *The Astrophysical Journal*, **644**(1), 551–559. https://doi.org/10.1086/503322

Wang, Z., & Cuntz, M. (2019a). S-type and p-type habitability in stellar binary systems: A comprehensive approach. III. Results for Mars, Earth, and super-Earth planets. *The Astrophysical Journal*, **873**(2), 113.

Wang, Z., & Cuntz, M. (2019b). S-type and p-type habitable zones of stellar binary systems: Effect of the planetary mass. *Research Notes of the AAS*, **3**(5), 70. https://doi.org/10.3847/2515-5172/ab1ea3

Way, M. J., & Georgakarakos, N. (2017). Effects of variable eccentricity on the climate of an Earth-like world. *The Astrophysical Journal Letters*, **835**(1), L1. https://doi.org/10.3847/2041-8213/835/1/L1

Way, M. J., Del Genio, A. D., Aleinov, I., Clune, T. L., Kelley, M., & Kiang, N. Y. (2018). Climates of Warm Earth-like Planets. I. 3D Model Simulations. *The Astrophysical Journal Supplement Series*, **239**(2), 1–22. https://doi.org/10.3847/1538-4365/aae9e1

Welsh, W. F., & Orosz, J. A. (2018). Two suns in the sky: The Kepler circumbinary planets. In J. D. Hans, & J. A. Belmonte (Eds.), *Handbook of exoplanets, springer living reference work* (pp. 2749–2768). Cham, Switzerland: Springer International Publishing. https://doi.org/10.1007/978-3-319-55333-7_34

Welsh, W. F., Orosz, J. A., Carter, J. A., Fabrycky, D. C., Ford, E. B., Lissauer, J. J., et al. (2012). Transiting circumbinary planets Kepler-34 b and Kepler-35 b. *Nature*, **481**(7382), 475–479. https://doi.org/10.1038/nature10768

Wolf, E. T. (2020). Storyofthewolf/ExoCAM: circumbinary_jgrplanets_frozen (version circumbinary_jgrplanets_frozen). Zenodo. http://doi.org/10.5281/zenodo.3968971

Wolf, E. T., Haqq-Misra, J. H., Kopparapu, R., Fauchez, T. J., Welsh, W. F., Kane, S. R., et al. (2020a). The resilience of habitable climates around circumbinary stars: 3D climate model data part 1 [data set]. Zenodo. http://doi.org/10.5281/zenodo.3964961

Wolf, E. T., Haqq-Misra, J. H., Kopparapu, R., Fauchez, T. J., Welsh, W. F., Kane, S. R., et al. (2020b). The resilience of habitable climates around circumbinary stars: 3D climate model data part 2 [data set]. Zenodo. http://doi.org/10.5281/zenodo.3972021

Wolf, E. T., & Toon, O. B. (2013). Hospitable Archean climates simulated by a general circulation model. *Astrobiology*, **13**(7), 1–18. https://doi.org/10.1089/ast.2012.0936

Wolf, E. T., & Toon, O. B. (2014). Controls on the Archean limate system investigated with a global climate model. *Astrobiology*, **14**(3), 241–253. https://doi.org/10.1089/ast.2013.1112

Yang, J., Cowan, N. B., & Abbot, D. S. (2013). Stabilizing Cloud Feedback Dramatically Expands the Habitable Zone of Tidally Locked Planets. *The Astrophysical Journal Letters*, **771**(L45). https://doi.org/10.1088/2041-8205/771/2/L45