stars are complex, dynamic objects. they are the “atoms” of galaxies, the building blocks of our cosmos, the parents of uncountable solar systems like and unlike our own. we’ll never perfectly describe how they work, but it turns out that we can accurately model certain kinds of stars with a fairly simple set of equations. you can check out the python code i wrote here, on my github.
What do I mean by certain “kinds” of star? Stars, for the most part, are big balls of gas. They are so massive, their cores become so dense and hot, that they fuse hydrogen into helium, which produces energy that they radiate away as light. How massive a star is determines how dense and hot its core is, which determines how much energy it produces, (which controls how puffy the star is, which controls the … you get the picture). So, when I’m speaking about different types of stars in this post, I’m talking about a sequence in mass for stars of the same age (things can change a lot as a star grows older, but that’s for another blog). For instance, stars with masses between 0.3 and 1.3 times the mass of the sun behave differently than stars with masses between 1.3 and 3.
Based on this view of stars, its natural to imagine that a star’s mass might be the “independent variable” in a set of equations describing its behavior. For a population of stars, the mass would describe how different stars look compared to one another. Within a single star of a given total mass, the “enclosed mass” $m$ controls the gravitational force pulling inwards towards the core, and therefore the density at that location within the star, and therefore the structure of the star generally. Let’s take a look at a set of equations we might construct to describe the “structure” of a star.
$$\frac{\partial r}{\partial m} = \frac{1}{4\pi r^2 \rho},$$
$$\frac{\partial P}{\partial m} = -\frac{G m}{4\pi r^4},$$
$$\frac{\partial l}{\partial m}= \varepsilon,$$
$$\frac{\partial T}{\partial m} = -\frac{GmT}{4\pi r^4 P} \nabla,$$
$$P = \frac{\rho kT}{\mu}+\frac{1}{3}aT^4,$$
See, for example, Hansen, Kawaler & Trimble (2004) or Chandrasekhar (1939) for more.
Let’s talk through what these differential equations mean.
$$\frac{\partial r}{\partial m} = \frac{1}{4\pi r^2 \rho},$$
This equation is called the “mass transport” or “mass continuity” equation. Its a deceptively simple equation, related to density (which is $\rho$). When you take the integral of this equation from the center of the star ($r=0$) to the surface ($r=R_\star$), you sum up all the mass, and the integral gives you the total mass of the star $M_\star$. This equation may seem frivolous, but it makes sure that mass is not created or destroyed (something we know to be true) by the other differential equations.
$$\frac{\partial P}{\partial m} = -\frac{G m}{4\pi r^4},$$
This is a fun one! We call this an equation of “hydrostatic equilibrium.” This equation is a claim that force due to gravity (that $Gm$ term on the right) is balanced by the outwards pressure. That is, the star isn’t collapsing, and there is something pressing outwards to keep the star in the same basic shape. But what is supplying that outwards balancing pressure? For that, we need an “equation of state” that is dependent on what states of matter apply to the specific composition of gas in the star. This equation,
$$P = \frac{\rho kT}{\mu}+\frac{1}{3}aT^4,$$
tells us that the pressure, temperature, and density within our star are related by the kind of matter that comprise the star. In this case, we assume that the pressure support is due to the intrinsic pressure of an ideal gas (so we are making an ideal gas assumption), the $\frac{\rho kT}{\mu}$ term, and also due to the radiation pressure, the $\frac{1}{3}aT^4$ term (all the photons of light being generated by the star’s nuclear fusion impart some outwards flowing energy and excite the gas, producing pressure support).
$$\frac{\partial l}{\partial m}= \varepsilon,$$
This equation is deceptively simple, mostly because we’ve kick the can down the road on what actually generates energy within a star, and simply call whatever total energy is generated in a star over time $\varepsilon$. The equation itself states that the change in the amount of light produced within a portion of the star is proportional to the energy generation, which makes sense because light is just energy in the form of a little particle we call a photon. The actual value of $\varepsilon$ depends on correctly calculating all the nuclear and quantum physics that go into predicting the outcomes of nuclear fusion. We can take the results of those calculations and plug them into our equations here to determine how much light the star is producing.
$$\frac{\partial T}{\partial m} = -\frac{GmT}{4\pi r^4 P} \nabla,$$
Last, but certainly not least, this equation describes how the temperature of a star changes from the core to the surface. It depends on another can we’ve kicked down the road, $\nabla$ which is called the “temperature gradient.” The temperature gradient depends on how heat is transported in a star and is a combination of radiation (the flow of light through the star) and convection (the bubbling of pockets of gas with different heats moving up and down through the star). The heat transport in a star changes with the total mass of the star, and the specific location within the star, since at certain densities and temperatures, it is easier or harder to push light through gas. When it becomes too difficult to push light through the gas, convection will take over, and move globs of gas up and down to transport the energy. For instance, the core of our Sun is radiative, but the outer region, its “envelope,” is convective.
In my code, I included a few good guesses for $\varepsilon$ and $\nabla$ for stars with masses around 1-5 times the mass of the sun. For stars of this mass, it is safe to assume $\varepsilon$ comes from nuclear fusion from the proton-proton chain and the CNO-cycle, two methods for converting 4 hydrogen atoms into 1 helium atom. Similarly, it is safe to assume $\nabla$ is usually radiative, but when it is convective, the convection is a simple adiabatic process that we can describe with a few physical constants thanks to our other assumptions, like the ideal gas assumption.
To solve our system of differential equations, we need to use an integrator, like the Runge-Kutta method, to compute the change in our equations over our independent variable $m$. We also need to know where to start our integrator, we call these “boundary conditions.” We know that at the center of the star, $m,r,l=0$ and at the surface $m=M_\star, r=R_\star, l=L_\star$, and $P,T\rightarrow0$. Even still, we don’t know what the pressure and temperature at the center, $P_c$ or $T_c$ are, and we don’t know what the final values for $R_\star$ and $L_\star$ will be. Luckily we have two sets of boundary conditions, and two sets of unknown initial conditions. We can integrate inwards from the surface using one set of boundary conditions, and a reasonable guess for the unknowns, and then do the same integrating outwards from the surface with the other set of boundary conditions, and another reasonable guess. If the inwards and outwards integrations agree, we know our guesses were the correct ones. So, we’re going to need to test different values for these guesses, to find the ones that make our inwards and outwards integrations meet.
What’s the final result? We have a code that will repeatedly integrate our equations from either boundary until our guesses for the core temperature and pressure agree with our guesses for the surface luminosity and radius. The final result is that we can choose a star of a certain total mass, say, $M_\star = 3~M_\odot$, and produce a graph like the one below, which shows how the variables of the star’s structure (the density, pressure, temperature, and luminosity) change as you move from the center ($m=0$) to the surface ($m=M_\star$). We’ve built our own star!
Using the code I wrote, we can compare stars of different masses. We can also take a peek under the hood, and see where these stars are convective vs radiative (when the blue line is below the red line, the star is radiative). And, we can see what proportion of the total energy of that star is being produced by the two nuclear reactions I mentioned.
Pretty neat, right? With only a handful of equations and some python code we can simulate what’s happening inside a star. Until next time, clear skies.
Leave a Reply