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

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.