An Apple Pie From Scratch, Part II Supplement: Simulating Star Evolution with Starpasta

 

Simulated evolution of the sun using Starpasta

As useful as the various estimation formulas for stars in the main post are, there are some obvious gaps; we’ve only accounted for main-sequence stars, and we only really managed to get good estimates for the evolution on the main-sequence for stars of 1 solar mass or less. We also didn’t account for metallicity at all, which can pretty significantly impact the lifetime and evolution.

So, for those who want to get really specific about stellar evolution, I’ve been working on an implementation of the formulas of Hurley et al. 2000, “Comprehensive analytic formulae for stellar evolution as a function of mass and metallicity”, which gives an algorithm for estimating the key features—luminosity, radius (and thus effective temperature), lifetime, mass accounting for loss to solar wind—over the full evolution of stars over a broad range of mass and metallicity. I've called my version Starpasta.

This is not, strictly speaking, a model of star evolution, as it does not directly represent the internal mechanisms that drive the evolution of stars. Instead, models of star evolution were run and then fits were created—formulae that approximate the important features of those models. These have been designed to achieve no more than 5% error from the models across most of the parameter space represented. However, the paper is a bit old and the models it was based on have been updated since then. I’ve included a couple updates from more recent papers (cited on the github repository) regarding the later stages of evolution for high-mass stars, but I cannot guarantee the results are fully accurate by today's modelling standards.

I also haven’t included any of the proscriptions in those paper for modelling the evolution of very close-orbiting binaries that might exchange mass with each other; I might be interested in exploring that in the future, but it would be a fair bit of work, so don’t hold your breath.

There’s also a decent possibility that there might still be bugs in the code that I haven’t identified. The actual math is fairly straightforward, mostly just arithmetic and algebra, but there are a huge number of these formulae which refer to each other in complex ways. I’ve tried to test out the code across a broad sample of stars, but there’s still a possibility of bugs arising for some specific cases—please do contact me if any arise.

Using Starpasta

Starpasta is available on this github in 2 forms: First, a python script, starpasta.py, which should be usable if Python 3.6 or later is installed with the numpy package. You can look through the script if you like and make any alterations you might desire. The second is a standalone package containing starpasta.exe and other files required for it to run, which doesn’t need python installed but isn’t editable.

Running either version should produce a command prompt window. First, it will prompt you for the stars initial mass in solar masses (multiples of our sun’s mass). The results should be reasonably accurate for the whole evolution of stars between 0.8 and 150 solar masses, but there’s no specific restriction on what you can put in here. At least some researchers have extended it up to 300 solar masses, which is probably around the maximum mass for stars in this metallicity range. For smaller stars, the results should be reasonably accurate for the main sequence down to around 0.5 solar masses and at least the early main sequence down to the minimum star mass of around 0.08 solar masses, but the later stages of evolution sometimes deviate from the results of better models, especially the evolution into a blue dwarf for stars below 0.25 solar masses; though given the long lifetimes of these stars, these aren’t issues if you’re only interested in stars that could exist in the current universe. For stars below 0.8 solar masses, you’ll be given the option to include the full evolution or stop at the end of the main sequence.

Then, you’ll be prompted for the star’s metallicity, the portion of the star’s mass composed of elements other than hydrogen and helium. This is around 0.02 for the sun and the code is designed for an intended range of 0.0001 to 0.03. Again, there’s nothing stopping you from inputting values outside that range but accuracy is not guaranteed and there’s a decent chance it’ll crash for values far outside the range.

The script will then run, usually not taking more than a few seconds, and show you a quick summary of the star’s evolution as it does. Once it says it’s complete and has written out the data, you can press “enter” to close the prompt.

The script will record data about the star as it runs and then, when complete, write it all out to a .csv spreadsheet named based on the star metallicity and mass. Each row of the spreadsheet corresponds to one timestep in the simulation, and from left to right the columns are:

• Timestep number.
• The evolutionary stage of the star, represented as a number (see later in this post for the meaning of these acronyms and an explanation of the stages): 0, end of the MS for stars <0.7 solar masses; 1, MS; 2, HG; 3, GB; 4, CHeB; 5, EAGB; 6, TPAGB; 7, HeMS; 8, HeHG; 9, HeGB; 10, HeWD; 11, COWD; 12, ONeWD; 13, NS; 14, BH; 15, no remnant.
• The star’s age in millions of years.
• The star’s mass in solar masses.
• The mass of the star’s helium core in solar masses (including all cores of heavier elements inside it in later stages, even when all helium is exhausted).
• The mass of the star’s carbon/oxygen core in solar masses (again, including all heavier elements).
• The rate of mass loss in solar masses per million years.
• The star’s luminosity in solar luminosities.
• The star’s radius in solar radii.
• The radius of the helium core in solar radii.
• The star’s effective temperature in Kelvin.
• The optimistic inner edge of the habitable zone in AU (based off Kopparapu et al. 2014; these may be inaccurate for effective temperatures over 7200 K).
• The conservative inner edge of the habitable zone in AU.
• The conservative outer edge of the habitable zone in AU.
• The optimistic outer edge of the habitable zone in AU.

To help you interpret all that data, I’ve created the starpasta_out spreadsheet. To use it, copy the data from the Starpasta output and paste it into the “Raw Data” tab of starpasta_out (make sure to paste it exactly so that it starts in the upper left corner, and delete any data in the tab before replacing it with data from another output).

You can then go to the “Results” tab and use the tools there to browse through the data:

The Output Summary box gives you a quick summary of the star’s evolution between stages and the duration of each one in millions of years.

Single Output lets you specify a specific age of the star and see the parameters at that particular point in time (in addition to the outputs listed above, the spreadsheet also calculates the peak wavelength of the star’s light and the distance at which a planet would receive the same light as Earth). These values are all linearly interpolated between the 2 nearest timesteps, which should be a decent estimate because the Starpasta code tries to ensure that timesteps are always short enough that none of these values change too much (excepting special events like supernovae).

The other boxes allow you to control graphs of the star’s whole evolution—both a Hertzprung–Russel diagram on the right and a user-defined graph below it, both of which display different stages in a star’s lifetime as different colors. Output Controls allows you to control which stages of the star’s lifetime are represented, or alternatively pick a start and end time to the data shown (this is necessarily rounded to the next timestep for the start time or the previous timestep for the end time). Starpasta will track the evolution of remnants (white dwarfs, neutron stars, black holes) for 10 trillion years, so you may need to exclude them or pick an appropriate end times to get reasonable graphs over time of the star’s evolution during its fusing lifetime.

Graph Control below that will allow you to pick two parameters from the list, which will be charted against each other in the graph to its right. The graph will tend to automatically start its axes at 0, so if you want a more detailed graph you may have to adjust the settings. For those that haven’t used Excel, right-click the graph, select “Format Chart Area…”, select the horizontal or vertical axis in the dropdown menu in the sidebar that appears, select the icon that looks like a bar graph, open “Axis Options”, and adjust the bounds there; other programs probably have a similar process.

Stellar Evolution

Now that we’ve got that settled, it’s probably worth going over the different stages of stellar evolution represented by these estimates. I’ve already mentioned several of these in the main post but I’ll get a bit little more detailed here. Note that all the mass boundaries I mention are approximate and will vary by metallicity, as you can discover using Starpasta, and I’m mostly referring to the initial mass of the star at the start of the main sequence, before any mass loss (and usually at sunlike metallicity).

First, the star begins on the main sequence (MS), where it will spend most of its fusing lifetime. Low-mass stars may have a significant pre-MS stage, but that’s not currently modeled. Stars will get a tad brighter and larger over the MS and either hotter for low-mass stars or cooler for high-mass stars, but not by much. As the MS comes to a close, stars over ~1.1 solar masses will become suddenly brighter, smaller, and hotter, as they run short of fuel and the outer layers contract, a behavior sometimes called the hook.

A distinct helium core forms at the end of the MS and hydrogen fusion now continues in a shell around the core. During the brief Hertzsprung gap (HG) stage (a.k.a. the subgiant stage), this pushes the outer layers outwards and the star rapidly becomes larger and cooler while maintaining about the same luminosity.

Stars of less than ~12 solar masses then enter the first giant branch (GB), during which the star becomes gradually brighter, larger, and cooler as the helium core becomes denser and more massive and the hydrogen-fusing shell grows around it, pushing out the outer layers more and more. Stars below ~0.5 solar masses will remain in the GB until their hydrogen fuel is exhausted, leaving a helium white dwarf (HeWD) (or they’ll become blue dwarfs if small enough instead, which isn’t modeled in this code, with the same end result).

For more massive stars, the helium core eventually becomes dense and hot enough to start fusing helium to carbon and oxygen. This core helium burning (CHeB) stage varies between stars: Below ~2 solar masses, helium fusion begins suddenly in an event called the helium flash; this "flash" will be absorbed by the outer layers and not cause any immediate change in its outward appearance, but the star abruptly transitions to a much dimmer, smaller, and hotter horizontal branch star over the course of a few thousand years (the code treats this as instantaneous, so you’ll sometimes see a sudden jump between the GB and CHeB) and then remains there as helium fusion continues in the core and hydrogen fusion in a shell around it. Between ~2 and ~12 solar masses, helium fusion begins earlier and more gradually; the star more gradually dims and shrinks, then enters a period called the blue loop, becoming first much smaller and hotter (moreso for more massive stars) and then expanding and cooling again. Stars above ~12 solar masses begin fusing helium during the HG, skipping the GB entirely and simply continuing to expand and cool as they rapidly burn through their helium cores, becoming a red supergiant.

Evolution across the HR diagram of a typical sun-mass star (left), with a short HG, long GB, and rapid (helium flash) transition the the horizontal branch; and a 5-solar-masses star, with a long HG, short GB, and gradual transition onto the blue loop. Lithopsian, Wikimedia

In any case, the helium fuel in the core is eventually exhausted, and the stars brighten, expand, and cool again. When a distinct carbon/oxygen core forms and helium fusion in the core stops, the star enters the asymptotic giant branch. During the early asymptotic giant branch (EAGB), helium fusion continues in a shell around the CO core, with no hydrogen fusion. Much as with the GB, the star brightens, expands, and cools as the CO core becomes denser and more massive. Often, the star becomes even larger than during the GB, with a radius hundreds of times that during the MS.

Eventually, the helium fuel is completely exhausted. The star briefly dims, contracts, and warms again, until hydrogen fusion begins again in a shell around the CO core, starting the thermally pulsating asymptotic giant branch (TPAGB). The star passes through a regular cycle of pulses: hydrogen shell fusion builds up a shell of helium around the CO core over periods of ~10,000–100,000 years, until enough helium has gathered to rapidly fuse; the star suddenly brightens, expands, and cools for several hundred years; but hydrogen fusion stops, the helium is eventually exhausted, and the star dims, contracts, and warms again until hydrogen fusion resumes. The formulae used here don’t show these individual pulses, but do show an average brightening, expansion, and cooling.

Detailed evolution of a sun-mass star during the TPAGB. Lithopsian, Wikimedia

Alternatively, metal-rich stars over ~25 solar masses may completely lose their outer layers of hydrogen due to solar wind during the CHeB or EAGB, leaving the exposed core as a small but bright and extremely hot naked helium star (a.k.a. a Wolf-Rayet star). These will pass through a main sequence (HeMS) of core helium fusion and then Hertzsprung gap (HeHG) and giant branch (HeGB) of shell helium fusion around a CO core, similar to regular stars but much briefer—skipping the HeMS and sometimes HeHG if the hydrogen layers were lost in the EAGB.

Even if this doesn’t occur, stronger solar winds will lead to ever faster mass loss in later stages. Most stars will eventually lose their outer layers during the TPAGB, briefly becoming bright and extremely hot as the hydrogen shell thins (these are also considered Wolf-Rayet stars, though the code counts this period within the TPAGB for convenience). The ejected hydrogen (with some heavier elements mixed in) forms a planetary nebula around the star, lit various colors by the intense UV light it puts out (it doesn’t actually have anything to do with planets, it’s an old name). But the star then rapidly dims as the remaining hydrogen and helium is exhausted and fusion ceases. For stars with an initial mass (pre-mass-loss) below ~7 solar masses, the remaining CO core remains as a small carbon/oxygen white dwarf (COWD), which gradually dims and cools thereafter.

Stars with initial masses of ~7­–8.3 solar masses will grow CO cores massive enough (1.08 solar masses) to fuse carbon to neon and magnesium. For stars below ~7.6 solar masses, this results in an oxygen/neon/magnesium white dwarf (ONeWD) when the outer layers are lost. But for more massive stars, the ONe core will reach a critical mass (1.38 solar masses) when the core becomes dense enough for electron capture to convert magnesium to neutrons, allowing the core to rapidly contract; the sudden collapse and then rebound of the outer layers (aided by rapid fusion around the core) produces an electron-capture supernova, blowing away the outer layers and leaving the core as an extremely dense neutron star (NS).

Stars with initial masses over ~8.3 solar masses will instead form cores hot and dense enough to fuse neon and oxygen, which supports the core against collapse. Fusion continues to ever heavier elements until eventually producing iron and nickel. Past this point, further fusion will consume energy, not produce it, so the outward push of energy holding up the star’s outer layers is lost and the core collapses. In most cases, this will result in a supernova. For stars below ~21 masses, the outer layers will be lost completely, usually leaving the core as a NS; for higher initial masses, some of the material initially ejected by the supernova will fall back down on the core, which usually results in it reaching a high enough mass (~2.5 solar masses) to collapse further into a black hole (BH). Especially massive (~40 solar masses) and somewhat metal-poor stars may instead undergo direct collapse, with the entire remaining mass of the star promptly collapsing into a BH, save for around 1/10 which is ejected as a burst of neutrinos.

Even more massive and metal-poor stars that grow helium cores above ~45 solar masses (which corresponds to a varied range of initial masses due to nuances in core growth and mass loss rates) may develop the right conditions in their extremely hot cores to spontaneously produce pairs of electrons and positrons. This effectively drains some of the energy holding up the star’s outer layers, causing some of the core to collapse inwards. But the sudden increase in pressure and temperature then triggers rapid fusion: For stars with helium cores below ~65 solar masses, the result is a pair-instability pulsation supernova that strips away the star’s outer layers but leaves the core mostly intact; the star may experience multiple pulses until the core is too small to produce another, and then ultimately collapses into a BH (these pulses aren’t simulated here, only their impact on the final BH mass). For stars with helium cores of ~65-135 solar masses, the resulting pair-instability supernova is so intense that it blows apart the star completely, leaving no remnant. Even more massive stars, with helium cores over ~135 solar masses, will collapse too forcefully to be stopped and again undergo direct collapse to a BH.

For reference, here's a quick overview of all these possible evolutionary stages:

And that’s about it for today. Have fun playing around with Starpasta, and do let me know if any issues arise.

Buy me a cup of tea (on Patreon)