Home

Analytic Insolation Calculator for Circumbinary Planets

Project

This project provides a time resolved insolation calculator for circumbinary planets. The relativistic three-body problem (star-star-planet) is solved analytically for planets on initially circular orbits. It is not necessary to propagate orbits numerically. For any moment in time t, the state vectors of the system as well as the incident stellar flux in predefined spectral bands (insolation) can be retrieved instantaneously.

The module twostars_m.f90 contains the analytic orbit propagator (Georgakarakos & Eggl 2015) that describes the orbit evolution of circumbinary planets as well as routines to evaluate the amount and spectral distribution of the combined stellar flux on an extended, rotating planet. Stellar eclipses are simulated as well. In Popp & Eggl (2017) this module was coupled with a (proprietary) general circulation model (GCM) to investigate the impact of variable insolation conditions on the climate of Earth-like circumbinary planets.

The planet’s orbit tends to evolve rapidly and significantly in such systems. Yet, its time evolution of can be modeled analytically given certain assumptions (Georgakarakos & Eggl, 2015). We assume the system to be co-planar, i.e. all bodies move in the same plane perpendicular to the system’s angular momentum vector. The variations in the semimajor axes of the binary and planetary orbits are small for dynamically stable orbits far from resonance. Hence, changes in the semimajor axes are neglected. Furthermore, we assume that the eccentricity vector of the binary star evolves due to perturbations of the outer body and General Relativity. The initial eccentricity of the planet shall be small. Then, the orbital evolution of the system is determined by

the system’s masses: m0; m1; m2,
the semimajor axes of the binary star’s and the planet’s orbit: a1; a2,
the eccentricity of the binary star’s orbit: e1,
and the initial positions (mean longitudes) of all bodies on their respective orbits: M1; M2.

For details on the analytic theory see Georgakarakos & Eggl (2015).

Stellar Radiative Properties

In the code the following parameters are used to evaluate the radiative properties of a star such as its luminosity and the spectral distribution of its light:

the stellar radius,
the star's effective temperature,
the stellar metallicity,
the stellar surface gravity / mass.

In order to allow for a self-consistent treatment of these quantities we offer the subroutine insolation1. The latter requires the input of stellar masses and metallicities and completes the other values based on a least squares interpolation of observational data of stars on the Zero Age Main Sequence (ZAMS) (Tout et al., 1996). As the fit underestimates the luminosity values for the Sun, we decided to include the subroutine insolation2 where explicit values for the stellar radii and effective temperatures can be entered independently from the stars’ masses. Additionally, we assume that the stars’ energy output is constant in time and that the stellar components of the binary star are distant enough to evolve independently on the ZAMS. In other words, we exclude systems that fill their mutual Roche-radii as well as systems that exhibit non-negligible mass transfer or common envelope evolution. Light travel time corrections for the planetary insolation have been neglected as well. The stellar spectra are modeled via effective temperature based Planck-curves without emission or absorption lines. The spectral band flux is calculated using an iterative method developed by Widger Jr & Woodall (1976). This can be changed by modifying the subroutine intplanckwn.

Stellar Eclipses

The stars may also block each other’s light in a co-planar configuration. Stellar eclipses effect are accounted for. Eclipses scale with the stars’ luminosity and covered surface ratios.

Planetary Spin

Planets are approximated as rigid spheres. As such, effects due to tides or the planet’s rotational oblateness are neglected. As a consequence, the rotation period and the direction of the planet’s spin vector remain constant.

How-to?

Initialize orbital and physical parameters of the system. Choose number and width of the desired spectral bands.
The routine checkconf should be called before the start of the simulation, to ensure that all variables have been correctly initialized and that the desired initial conditions constitute a dynamically stable system compatible with the limitations of the analytical propagator.
call insolation1 or insolation2 to calculate the top-of-the-atmosphere insolation at the chosen epoch in the desired spectral bands.

The user can choose between two subroutines, insolation1 or insolation2. Both deliver the same output. Subroutine insolation2 requires a consistent manual choice of stellar radiative properties, such as mass, radius and effective temperature, whereas insolation1 handles this automatically by choosing proper luminosities and effective temperatures on the basis of given stellar mass and metallicity values. Since the results are based on fits, they do not reproduce exact values for the Sun. Should the user require exact values for stars with known properties such as the Sun, it is recommended to use the routine insolation2.

In the subroutine insolation3, stellar spectra are taken from input files input_spectra1.csv and input_spectra2.csv, instead of being calculated as blackbody spectra. Input spectra can be modified by updating these files. This subroutine outputs the weighted stellar spectra and weighted combined spectra at each time step, as well as the set of time-dependent weights. The weights do not depend on the spectra, so if only the weight factors are needed any input spectra (or spectra of one band with flux of 1 unit) can be used to provide a more general output.

Common input variables for both subroutines are:

Input Variables	Code	Units	Description
Ω	spin(1)	[deg]	Initial planetary obliquity, i.e. the angle between orbit angular momentum (z-axis) and the system's x-axis. The latter is pointing along the initial eccentricity vector of the binary star's orbit.
α	spin(2)	[deg]	Initial angle of precession of the planet's spin axis, i.e. angle between the system's x-axis and the planetary spin vector projected onto the orbital xy plane.
P	spin (3)	[days]	Rotation period of the planet.
η	spin(4)	[deg]	Initial angle between x-axis and zero-meridian measured along the planet’s equator.
M1, M2	man(1:2)	[deg]	Initial positions (mean longitudes) of the binary star and the planet on their respective orbits.
t	t	[days]	Epoch of the system. At t=0 the binary star has a mean longitude of M1, the planet M2, and the planet's orbital eccentricity is e2=0.
m0, m1, m2	mass(1:3)	[Msun]	Masses of the primary star, secondary star, and the planet.
a1, a2	semia(1:2)	[au]	Semimajor axes of the binary (inner) and planetary (outer) orbit. 0<a1<a2
e1	e1	[]	eccentricity of binary star's orbit 0 <= e1 < 1; the planet starts out on a circular orbit

Input variables specific for the subroutine insolation1 are:

Input Variables insolation1	Code	Units	Description
Z1, Z2	met(1:2)	[]	Metallicity of each star (0.0001-0.03); The solar value is 0.02.

Input variables specific for the subroutine insolation2 are:

Input Variables insolation2	Code	Units	Description
Teff1, Teff2	teff(1:2)	[K]	Effective temperatures of the stars.
R1, R2	rs(1:2)	[m]	Stellar (mean) radii.

Input variables specific for the subroutine insolation3 are:

Input Variables insolation3	Code	Units	Description
R1, R2	rs(1:2)	[m]	Stellar (mean) radii.
Spectrum1, Spectrum2	spectra(1:2)	[photons m^-2 s^-1 A^-1]	Stellar spectra

Output variables common to all routines are:

Output Variables	Code	Units	Description
φ	phi(1:2)	[deg]	Hour angles between zero meridian and stars A and B measured along the equator.
δ	decl(1:2)	[deg]	Declinations of star A and star B with respect to the planet’s equator.
Planetary Insolation	bndflux(1:2,1:nbands)	[W m^-2]	Top of the atmosphere flux on the planet for each of the two stars in each spectral band.

Output variables specific for the subroutine insolation3 are:

|Output Variables insolation3 Code Units Description |------------------|---------------|-----------|-------------------------------------------------------------| |Distances | dr1, dr2 | [au] | Distance from planet to each star | |Spectral weights | Weights(1:2) | [] | Weight factors to scale each stellar spectrum due to distance to stars and eclipsed area of stars at given time step. |

Figures 1 and 2 provide a visualization of the setup.

circumbinary planet

Figure 1: The circumbinary planet setup projected onto the orbit plane. The quantities m0; m1; m2 and a1, a2 denote the system masses and semimajor axes of the binary and the planetary orbit, respectively. The x-axis is aligned with the initial direction of the pericenter of the binary. Please note that the angles φ and δ are projections. The actual angles are measured along the planet’s equator, not the orbital plane (ecliptic).

circumbinary planet

Figure 2: Projection onto the xz-plane. The system is coplanar, and the z-axis is aligned with the orbital angular momentum axis.

Validation Test Cases

In order to demonstrate the use of the module twostars_m, a test-program is provided: twostars2 works with the subroutine insolation2. The program simulates a simple time propagation of the system. The evolution of planetary insolation, declinations δ and the hour angles η is written in a .csv file. Figures 4-6 contain the results of several test cases. They feature stars with effective temperatures of 5777 [K] and a radius of 1 [Rsun] for G-stars and 5000 [K] and 0.9 [Rsun] for K-stars. The input parameters used to generate Figure 4 as well as the corresponding program output are provided in the "validation" folder.

circumbinary planet

Figure 4: An Earth-like planet orbits a binary with one sunlike G-star and a smaller K-star. The stellar flux arriving at the planet is calulated in a single spectral band from 1-10000 [1/cm]. The binary star features an elliptic orbit (e1=0.5). Significant changes in insolation occur within days. The sudden and drops in insolation are due to stellar eclipses. Two dominant frequencies appear in the insolation. The higher frequency (around 1=20 [1/D]) is due to the stellar motion. The lower frequency (around 1=280 [1/D]) stems from the evolution of the planet’s orbit.

circumbinary planet Figure 5: The planet orbits two sun-like stars. Stellar eclipses are complete for both stars (spikes reach down to the abscissa), since they are of identical size and luminosity. In contrast, the smaller K-class star in Figure 4 is not able to block all the light of the larger G-class star.

circumbinary planet Figure 5: The time evolution of the declination for each star for a planet with Earth-like obliquity orbiting a binary star at an initial distance of 1 au.