After seeing in science fiction too many planets with two suns that look too much like a geocentric system, I'm trying for my own amusement to understand if it's really possible to have a planet with two suns that can sustain life and how different will it be from our planet.
Since finding a 3-body solution that has stable orbits and stable temperatures for the planet is really hard, I decided to narrow the choices and cheat a little bit. I looked at another 3-body system that we all know and is stable enough to last a few billion years: Sun-Earth-Moon. To scale it up, the idea is to have a planet that goes around a red dwarf (Star 1), that in turn orbits around a blue giant (Star 2). To be realistic, the planet should be tidally locked to the red dwarf (and that should simplify the calculations too). To simplify the situation even more, I imagined the planet to be Earth-like: same mass, density, albedo, composition, tilt, etc.
Let $\alpha$ be the latitude (0 at the equator, 90° at north pole), $\beta$ the longitude (0 at the hot pole, where the red dwarf shines perpendicular to the ground) and $\delta = 23.5^{\circ} \sin \left( \frac{2 \pi}{\tau_2} t \right)$ the tilt of the planet with respect to Star 2, where $\tau$ is the orbital period.
Where Star 1 shines it has a fixed angle with the azimuth $$\sin(\gamma_1) = \cos(\alpha) \cos(\beta)$$ while Star 2 has a variable angle $$\sin(\gamma_2(t)) = \sin(\alpha) \sin(\delta) - \cos(\alpha) \cos(\delta) \cos \left( \frac{2 \pi}{\tau_1} t - \beta \right)$$ where $t = 0$ means it's midnight at $\beta = 0$. Since Star 2 is a blue giant, the period is big enough that the time dependency of $\delta$ is negligible and can be treated as a constant for daily variations.
Sunrise and sunset times can be calculated by requesting that $\sin(\gamma_2) = 0$:
- If $\delta = 0$, then $t_{sr} = \tau_1 \left( \frac{1}{4} + \frac{\beta}{2 \pi} \right)$ and $t_{ss} = \tau_1 \left( \frac{3}{4} + \frac{\beta}{2 \pi} \right)$.
- If $\delta > 0$ then:
- If $- \frac{1}{\tan(\delta)} \le \tan(\alpha) \le \frac{1}{\tan(\delta)}$ then $t_{sr} = \frac{\tau_1}{2 \pi} \arccos(\tan(\alpha) \tan(\delta)) + \frac{\tau_1 \beta}{2 \pi}$ and $t_{ss} = \tau_1 - t_{sr}$
- If $\tan(\alpha) > \frac{1}{\tan(\delta)}$ then Star 2 always shines
- If $\tan(\alpha) < -\frac{1}{\tan(\delta)}$ then Star 2 never shines
- If$\delta < 0$ then: as above, except the second and third case are exchanged
I can evaluate a mean temperature by looking at the stellar constants $I_1$ and $I_2$ and comparing their sum with our solar constant to obtain a mean temperature for a certain latitude and longitude. This helps a little, but the problem is that it's a good estimate only if $\tau_1$ is similar enough to Earth's period. A mean temperature on Earth means that minimum and maximum temperatures are usually less than 10 K from the mean temperature, but if $\tau_1$ is bigger the planet has more time to absorb heat during the "day" and release heat during the "night", widening the difference.
My second approach was to set up a differential equation:
- The total energy is $\mbox{d}E_{tot} = c m\,\mbox{d}T$, where $c$ is the specific heat and $m$ the mass where the heat is stored.
- The incoming energy is $\mbox{d}E_{in} = a (1-A) (I_1 \sin(\gamma_1) + I_2 \sin(\gamma_2 (t))) \mbox{d}t$, where $a$ is an area and $A$ the planet's albedo
- The outgoing energy is $\mbox{d}E_{out} = \sigma a T^4 \mbox{d}t$ (Stefan-Boltzmann law)
The resulting equation is: $$\frac{\mbox{d}T}{\mbox{d}t} = \frac{a (1-A)}{c m} (I_1 \sin(\gamma_1) + I_2 \sin(\gamma_2 (t))) - \frac{\sigma a}{c m} T^4$$
where $I_1 = 0$ on the darker side and $I_2= 0$ when Star 2 is not visible. This means that this equation is really 4 different equations.
This approach has two great problems: First the parameters $c$, $m$, $a$ seems to be easy to define when talking about a whole planet, but quite hard when we analyze only a small portion. Second, and more important, the equation seems to be unsolvable at least when Star 2 is visible. In general it's a Chini differential equation (see here). When both stars are not visible it becomes a Bernoulli equation and the solution is easy to find; when only Star 1 is visible the Chini invarant is constant ($C=0$ to be precise), so there is a precise solution to that case too. When Star 2 is visible, however, there seems to be no way to find a solution. So, after this wall of text, my question is:
- Is there a way to calculate the explicit solution to the equation?
- Alternatively, is there a better approach to the problem that can lead to a solution? To summarize, I'm interested in calculating minimum and maximum temperature for a given latitude, longitude, tilt and orbital period $\tau_1$. The temperature every instant is just a bonus, not a necessary feature.
Edit: I've done some computer simulations and it seems to work, but it needs some calibrations due to the not-so meaningful parameters and greenhouse effect, that acts only on the emitted thermal radiation. The model still doesn't take into account the heat redistribution due to air and water flows, but this is probably way too complicated to put in the model. A few plots for Earth ( 1 star, $\tau = 1$ day) can be seen in this gallery. Blue line indicates sunrise, red line sunset, green line is the instantaneous temperature, yellow line is the average temperature and the arrow field should give an idea of the slopes the differential equation produces, but it always seems less inclined than the green curve.
The plot is produced with a periodicity condition: the algorithm starts with a temperature of 300 K and cycles through a few days until the temperature at the beginning of the day matches the temperature at the end of the day. The seasonal variation is slow enough that this should be a good approximation.
