How to Calculate Planetary Equilibrium Temperature in Binary systems

Question

I'm working on a program that generates the basics stats of a terrestrial planets in binary star systems. I'm not the best at this kind of math, so I'm having trouble calculating the planetary equilibrium temperature.

The formula I'm currently referencing was posted in a thread from almost 4 years ago: How to calculate the expected surface temperature of a planet.

It seems reasonably accurate since (if memory serves me) using it yields a good answer for Earth and Mars. However this is for a system with one radiating body not a system with two radiating bodies (binary stars).

My instinct is to treat the two stars as one big star and add the wattage $(w/m^2)$ the planet receives from the various stars. So total watts = watts star 1 + watts star 2. That seems reasonable to me, at least for the purposes of calculating planetary temperatures.

I've already done work to calculate the watts received by the planets from the various sources. I've used this equation, subbing in values for both stars.

$$ F(L/S^2) $$

• L = Luminosity (in fraction of Solar), stars 1 & 2 (respectively): 0.0159 & 0.7758
• S = Semimajor Axis (in AU), stars 1 & 2 (respectively): 0.1228 & 3.6000
• F = Earth Mean Solar Flux $(w/m^2)$ 1362

So the above calculation yields a Mean Solar Flux for the example planet of $1344w/m^2$ from Star 1 and $92w/m^2$ from Star 2. Adding those together I get a mean global wattage of $1436w/m^2$.

However, I'm unsure how to modify the original equation to utilize the above calculate solar flux (ideally without having to reference the original luminosity and semimajor axis variables). What is the correct way to alter the formula?

Any help would be very much appreciated!

Answer 1

The power radiated from the planet, assuming it is at the same temperature across its entire surface:

$$P_{\mathrm{rad}}=\epsilon\sigma T_{\mathrm{eq}}^4 \cdot 4\pi R_{\mathrm{pl}}^2$$

Where $\epsilon$ is the emissivity (to match the formula in the question, set this to 1), $\sigma$ is the Stefan-Boltzmann constant, $T_{\mathrm{eq}}$ is the equilibrium temperature, and $R_{\mathrm{pl}}$ is the radius of the planet.

The power absorbed by the planet from the stars, assuming it is sufficiently far away from each star that there is negligible illumination from that star beyond 90° from the substellar point:

$$P_{\mathrm{abs}}=\sum_i{f_i\left(1 - A_i\right) \cdot \pi R_{\mathrm{pl}}^2}$$

Where $f_i$ is the flux from the $i$th star, and $A_i$ is the planetary albedo with respect to the $i$th star (this may be different because the reflectivity will be wavelength-dependent, and the stellar temperature affects which wavelengths are being emitted).

So equate radiated and absorbed power, and rearrange:

$$T_{\mathrm{eq}} = \left[ \frac{\sum_i{f_i\left(1-A_i\right)}}{4\epsilon\sigma} \right]^{1/4}$$

The limitation here is that the fluxes are going to be time-variable due to the planet's distance changes from the stars (same applies to a planet in an eccentric orbit around a single star), and it is non-trivial to figure out how to average this out. The real planet will take time to heat up and cool down, and in that time the distances and fluxes will change.

So take the numbers you get from this "equilibrium temperature" calculation with a pinch of salt.

---

As requested, here are some example calculations. First off, the sanity check: Earth. Taking 1361 W/m² as the value of the solar constant and using an albedo of 0.3 and an emissivity of 1, the temperature works out as

$$T_{\mathrm{eq}} = \left[\frac{1361\ \mathrm{W\!\cdot\!m^{-2}\,\times \left(1-0.3\right)}}{4 \times 1 \times \left(5.6704 \times 10^{-8}\ \mathrm{W\!\cdot m^{-2}\!\cdot\!K^{-4}} \right)}\right]^{1/4} \approx 255\ \mathrm{K} $$

A bit chilly, but not too bad: the greenhouse effect keeps the actual Earth warmer, which could be taken into account by using an emissivity $\epsilon < 1$.

Now for the example in your question, with fluxes of 1344 W/m² and 92 W/m². Let's also suppose that while the first star is Sunlike (so I will keep the 0.3 value for the albedo), the second star is cooler than the Sun. This star emits more of its light in the red where the planet is less reflective. I will use an albedo of 0.25 for this star.

The calculation works out as follows (dropping the units in the calculation step for the sake of space):

$$T_\mathrm{eq}=\left[\frac{1344\times\left(1-0.3\right) + 92\times\left(1-0.25\right)}{4\times1\times\left(5.6704\times10^{-8}\right)}\right]^{1/4} \approx 258\ \mathrm{K}$$

Hope that helps!