notes-lecture-HSE.org

400A - Hydrostatic equilibrium

Materials: chapter 2 of Onno Pols’ lecture notes, chapters 1 & 2 of Kippenhahn’s book.

The first two equations of stellar structure

So far we have mostly been discussing astronomical facts on stars (e.g., use of populations to understand phenomena too slow to observe, the CMD/HRD diagrams, various types of binaries and how to use observables to derive masses), albeit anticipating some physical modeling (e.g., some basics of black-body radiation to introduce T_eff). In this lecture, we will move on to proper astro physics of stars: we want to derive a set of equations to model the parts of the stars that we cannot see, the interior structure of stars and how they evolve in time. This will take a few lectures to construct, and it will require several approximations - I will try to highlight these when we encounter them.

These equations express the conservation of mass, momentum, and energy, plus equations that describe the changes in chemical composition, and an equation that relates thermodynamic quantities to each other. These are typically expressed in a specific form for stellar evolution applications, to make the problem of determining the internal structure of a star and how it evolves tractable.

Note that deriving these equations will take us away from directly observable phenomena, and push us into theoretical modeling based on physical understanding that needs to be indirectly validated on observations.

Before we get into arguments specific to stars, let’s review briefly hydrodynamical coordinate systems that we will be using to describe the star.

Lagrangian and Eulerian coordinates

In hydrodynamics one is typically interested in describing fields of variables that take one value for each point in the fluid (e.g., density, velocity, pressure, temperature, etc..). Let’s call a generic field Ψ (=ρ, v, P, T, etc…) and discuss the two ways we can use to described it.

Eulerian description

In a Eulerian description the fields are described as a function of position in space (and time): Ψ ≡ Ψ(x, y, z, t) in Cartesian coordinates or Ψ(r, θ, \varphi, t) in spherical-polar coordinates (of course, one could use any other spatial coordinate system as well, the point is the field is expressed as a function of a fixed location in space-time).

This is like looking at the flow from on top of a bridge at a fixed position, describing the variables at that specific position.

For a star assuming spherical symmetry (as we already did to relate the stellar flux with the luminosity and obtain L=4π R²σ T_eff⁴), we would have no dependence on the angular variables ∂_θ Ψ ≡ ∂_\varphiΨ ≡ 0 and we could express Ψ\equivΨ(r, t).

Lagrangian description

In a Lagrangian description the fields are described as a function of the fluid element, following the fluid element itself: Ψ ≡\Psi(fluid element, t).

This is like hopping on a boat that moves with the fluid, and describing the properties of the fluid as the boat moves with it.

In a Lagrangian system of coordinates, the total time derivative can be written using the “chain rule” as:

where the second term on the r.h.s. is the “advective” term, due to the motion of the fluid element.

Because of the extremely large range of spatial scales in stellar evolution (e.g., in the Sun from the center at r≅0 to the outer radius r≅ 10¹¹ cm, and in its future as a red giant 100-1000× larger - even using R_o as unit there are ∼3-5 orders of magnitude of dynamic range without even considering neutron stars and black holes!), compared to the limited range in mass (in the appropriate units 0-1 M_o for the Sun!), typically a Lagrangian description is preferred. This means, for example we can chose as independent coordinate in the star not the radius r, but instead the amount of mass enclosed by a given radius m≡ m(r), thus express Ψ(r,t) =Ψ(r(m), t)→ Ψ(m(r), t).

Derivatives can then be expressed with the chain rule:

This is what is typically done in stellar evolution, and it works because the m(r) function is invertible, and it is in fact monotonically increasing: as r increases, the amount of stellar mass enclosed can also only increase!

N.B.: For some applications, Eulerian and Lagrangian descriptions can also be “combined”, for example with “moving mesh” numerical approaches, see for example AREPO hydrodynamics code, DISCO code.

Spherical symmetry and one dimensional approximation

We have defined a star to be a self-gravitating mass of gas that can produce energy through internal sources. Thus, from the self-gravitating and the gas parts of this definition, at zeroth-order, the dynamical elements (i.e., the forces!) that determine the structure of a star are:

• Gravity: the weight of the gas forming the star itself
• Pressure: the thermal pressure of the gas, sometimes with contribution from the radiation and degeneracy

Gravity is a central force that depends only on the radius (∝ r^-2), and pressure is isotropic. Therefore at zeroth order, we expect stars to be well approximated as spheres. This mathematically means that we can express all the variables that characterize the structure of a star as a function of a single variable, for example the radius from the center of the star. This allows for the calculation of

• Q: can you think of cases where a star may not be spherical?

Mass conservation

Let’s consider the amount of mass in a parcel of stellar gas. This will depend on the local gas density $ρ(r, t)$ (or equivalently in the Lagrangian formalism $ρ(m, t)$!) and the amount of volume in the shell

where dA is the element base area, and dr its radial thickness. We can integrate over the base to get the parcel to be a spherical shell

where r is the radius of the shell, therefore

In principle gas could also flow in/out of the shell at a rate determined by the inflow/outflow velocity such that in a time interval dt an amount $-ρ v dA dt$ flows out (for $v>0$, the quantity is negative) or in ($v<0$). Again integrating over $dA$:

This is the complete mass continuity equation in spherical symmetry. From this complete form we can take the partial derivatives w.r.t. $r$ (at fixed $t$) and $t$ (at fixed $r$):

We can also derive the first equation above w.r.t. $t$ and the second one w.r.t $r$, and demand the two forms are the same. Since $r$ and $t$ are the independent variables here (i.e., $∂ r/∂ t = 0$) we obtain:

with $∂_θ ≡ ∂_φ ≡ 0$ for the last one, that is the typical form of the mass continuity equation in spherical symmetry.

To turn these equations in the more typical form for stellar structure, just take the first equation in \ref{eq:mass_continuity_rt} and express it with $m$ as independent variable:

where the partial derivatives become total derivatives in a static situation (where by definition $∂_t = 0$, which is also why we don’t typically focus on the second equation in \ref{eq:mass_continuity_rt} - by the end of this lecture we will be able to discuss whether this is an acceptable approximation). This is the first stellar structure equation that expresses mass conservation, and it depends on a yet unknown variable, the gas density $ρ$.

Momentum conservation and hydrostatic equilibrium

Consider the equation of motion of a parcel of stellar gas, $F = dp/dt = ma$ (for constant $m$), or often more conveniently in fluid dynamics, work per unit volume with $f = dF/dV$ and thus $f=ρ a$ with $ρ = dm/dV$ and $dV=dAdr ⇒ V=∫ dAdr$ the volume. Let’s start by writing down explicitly the forces that we think are important for an isolated, non-rotating, non-magnetic star.

Gravity

Since by definition a star is a self-gravitating body (N.B.: so is a planet, that’s not the whole definition of a star!), we want to include the gravitational force on the l.h.s. of our f=ρ a equation. This can be obtained as the gradient of the gravitational potential Φ which is a solution of the Poisson equation:

where the second form assumes already spherical symmetry. Note how this equation does not make the problem worse: we have a new variable $Φ$ but the r.h.s. only depends on the density $ρ$ which is already appearing in Eq. \ref{eq:mass_conservation}.

We can introduce the gravitational acceleration $g = - ∇\Phi$, which in spherical symmetry only has a non-zero radial component ⇒ $g = - dΦ/dr$ which from Newton’s theory of gravity we know to be

where $m≡ m(r)$ is the mass enclosed within a certain radius $r$, which we already encountered. Thanks to the spherical symmetry assumption, we don’t even need to really solve Poisson’s equation to make a stellar model! The gravitational force acting on a spherical shell of mass $dm = 4π r²ρ dr$ is thus just $-gdm = -Gmdm/r²$, or per unit volume $f_\mathrm{grav} = -gρ = -Gmρ/r²$, where the minus sign is to explicitly indicate that this force points towards the center of the star.

Pressure gradient

The other contribution we need to include in our $f = ρ a$ equation is from the pressure. We could already use dimensional analysis to guess in what form pressure can enter the l.h.s. of the equation:

This suggests that the pressure divided an appropriate length scale has the right dimension to enter the force per unit volume f. This in turn suggests that maybe what we need is the pressure gradient!

Let’s now have a slightly more formal look. Because of spherical symmetry, the pressure in the horizontal direction (which in stellar context always means in the plane orthogonal to the radial direction) is perfectly balanced, and the pressure only depends on the radius $P≡ P(r) (≡ P(r(m)) ≡ P(m))$.

The net force per unit area on each side of a spherical shell of gas of thickness dr is $P(r)$ at the inner boundary and $P(r+dr)$ at the outer boundary. Therefore, $dF_\mathrm{press} = P(r)dA - P(r+dr)dA ≅ dP/dr dA$ where we used $P(r+dr)≅ P(r)+(dP/dr)dr$. Now using $dm = ρ drdA$ and dividing by $dV = drdA$ we finally obtain $f_\mathrm{press} = - dP/dr$.

Combining the two

We have now an explicit form for the two most important forces in a (isolated, non-rotating, non-magnetic) star $f = f_grav + f_pres = -gρ - dP/dr ≡ ρ a$.

Since stars don’t change that much on short timescales (we will see exceptions later, and define relevant timescales too), we can assume that overall the acceleration a of each parcel of gas is zero in most cases, that is $a=0$. Stars are generally in hydrostatic equilibrium. In this case the conservation of momentum becomes

or changing to have $m$ has the independent variable, to have a Lagrangian treatment:

and thus

which is the second stellar structure equation that expresses the fact that the gravitational pull of the stellar gas is compensated by the pressure gradient inside the star. This also means that it is the gravity of the star that imposes the pressure stratification of the star and ultimately its structure: the pressure in each layer is just what is needed to support the weight of the layer(s) above. And finally, the fact that $dP/dm<0$, that is the pressure decreases as the enclose mass increases, or equivalently, the pressure increases towards the center (smaller radii, smaller amount of enclosed mass) makes sense, if the gradient has to compensate the gravitational pull.

All of stellar evolution can be though as gas re-arranging itself to fight against gravity, delaying gravitational collapse.

N.B.: The hydrostatic equilibrium equation can also be obtained starting from the Navier-Stokes equation assuming no viscosity (the microscopic viscosity is generally negligible in stars).

N.B.: We have implicitly assumed that the star we model is sufficiently far away from anything else that there are no external forces. This may not hold in a binary system, for which in the Euler equation there will be other terms, such as the gravity of the binary companion, and tidal forces arising from its presence. While these are important, they often affect most directly only the outer layers of a star (that can be significantly tidally distorted), and can maybe be neglected further in the interior.

N.B.: Similarly, we have neglected rotation, which also breaks the spherical symmetry by adding in the reference frame co-rotating with the star non-inertial forces (centrifugal, Euler, and Coriolis). However, the centrifugal foce depends on the distance from the rotation axis (r sin(θ)) and thus mostly impacts the outer layers of the stars and is less critical in the inner regions (though not at all always negligible!). The Euler force depends on dω/dt with ω rotation rate, so it is typically negligible on evolutionary timescales for the star. The Coriolis force depends on ω × v, so it does not affect static gas (but it does have important effects if there are velocities, think for examples hurricanes in the Earth atmosphere). Moreover, rotation can interplay with many hydrodynamical and secular instabilities affecting the stellar gas in ways that are only roughly approximated in models. All these complications introduced by rotation and how to model them in stellar evolution are still active research topics (see for example Maeder & Meynet 2000 and Heger et al. 2000).

N.B.: Finally, we have neglected also the impact of magnetic fields on the stellar gas. We know that stellar magnetism exists from observational phenomena such as stellar flares, seeing Zeeman splitting in stellar spectra, etc. We should also expect magnetic fields theoretically, because stars are giant balls of ionized gas. However, the global dynamical impact of magnetic field should be small in most cases, given the success of stellar evolution in explaining many observations neglecting them. Stellar magnetism (and it’s important interplay with rotation) is also an active field of research both observationally and theoretically.

Equations Eq. \ref{eq:mass_conservation} and \ref{eq:HSE} are two differential equations, that under the assumption of spherical symmetry are ordinary differential equations ($∂_r → d/dr$), for the function $m≡ m(r)$ that depend on $P$, $ρ$. We thus have three variables $(m, P, ρ)$ and two equations: we cannot yet solve for the structure of a star. We will close the system of equations (meaning, obtain as many equations as variables, so we can solve for the stellar structure) later in the course.

Estimate for the central pressure

A first estimate for the central pressure can be obtained substituting the local gradient with the difference from surface to the core across the entire mass of the star $dP/dm → (P_surface - P_center)/M ≅ -P_center/M$, where we also use $P$ increases inwards and thus it is legitimate to expect $P_center\gg P_surface$. Then, on the l.h.s. of Eq. \ref{eq:HSE}, we should take as estimates some fraction of the total mass M and radius R. For the sake of simplicity, let’s take the fraction to be 1 and drop the 4π:

Plugging in the numbers for the Sun this gives $P_center≅ 10¹⁶ dyne cm^-2≅ 10¹⁰$ atmospheres. Although this a is very imprecise estimate, it already gives the idea that the pressure in the center of the Sun must be extremely high. See Onno Pols chapter 2 for more precise estimates and lower bounds.

Dynamical timescale estimates

Let’s say that the star was not in hydrostatic equilibrium, but still spherically symmetric. Returning to the general form for the momentum conservation $f = ρ a ≡ ρ ∂²r/∂ t²$ we have

where since P decreases outwards, $dP/dr<0$, so the first term on the l.h.s. pushes outwards (positive radial acceleration), while gravity pulls inward, as one would expect.

Normally, for a star, we expect these two terms to balance each other, but what happens if we turn one off?

Explosion timescale

Let’s turn off gravity, setting $g = - Gm/r² → 0$! To estimate how long it takes for the pressure gradient to push the gas out to a radius comparable to the radius of the star we can do the following rough substitution in the dynamical equation above:

• $∂² r → R$ (outer radius of the star)
• $∂ t² → τ_mathrm{expl}²$ (what we want to estimate)
• $dP/dr→ P_\mathrm{avg}/R$ with $P_\mathrm{avg}$ some averaged pressure in the star
• $ρ → ρ_\mathrm{avg}$ some averaged density of the star

and we obtain:

where, if we interpret $P$ and $ρ$ as some average values throughout the star the sound speed $c_s² = P/ρ$ appears! The timescale for the internal pressure of the stellar gas to “blow up” the star is of the order of the sound-crossing time.

Free fall timescale

Almost by definition, this is how the star would collapse if there were no forces other than gravity, so let’s turn off the pressure gradient $dP/dr→0$. Then, as above:

• $∂² r → R$ (outer radius of the star)
• $∂ t² → τ_ff²$ (what we want to estimate)
• $m → M$ (total mass)

we get:

with $ρ_\mathrm{avg} = 3M/(4π R³)$ average density of the star. Note that here we have been very loose with the π factors and averages.

• Q: you all have estimated the Sun’s mean density, calculate the Sun free fall time now. Does the Sun vary on this timescale? Do you think this justifies our assumption of hydrostatic equilibrium?
• Q: are stars in hydrostatic equilibrium? How do we know observationally?

Introduction to MESA_web

We will discuss in detail stellar evolution codes, numerical strategies for solving the stellar structure equations, and what goes on in MESA/MESA-web. For now I just want to introduce this tool and show you how you can obtain numerical stellar models.

• Description of Input
• Submission website
• Example output:
  1. Download the zip file from the email you receive when the calculation is done
  2. Unzip the file, the content has a *.mp4 video with the evolution of some quantities (depending on the star you asked, it may be very short), an input.txt file that reminds you of what you put into MESA-web, the trimmed_history.data and a few profile*.data, and a profiles.index that contains a map of which profile*.data maps to which “model number” (i.e., timestep of the code).
  3. You can inspect the txt, list, and *.data files using your text editor.

     The trimmed_history.data contains in each column global variables of the star (e.g., surface luminosity, outer radius, etc.) and each row correspond to a specific timestep. This is what you can use to plot, for example, an Herzsprung-Russell diagram using the columns log_L and log_Teff.

     The profile*.data files contain each a snapshot of the internal structure of the star you simulated at fixed time, so each column corresponds to a quantity that takes different values at different locations in the star (e.g., Lagrangian mass coordinate, density, pressure, opacity). Each row corresponds to a “mesh point”, that is a discretized spatial coordinate (we will see later what the full set of equations is and how codes like MESA solve them).

     Refer to the MESA-web output page for a full description of the output.

  4. If you want you can use the python module =mesa_web.py= provided by MESA-web to read the output in the *.data, but remember these are just plain text, so you can also write your own.

Homework

• Calculate the Keplerian period of a point mass orbiting at the surface of a star of mass M and radius R and compare it to the free fall timescale of the star.
• Calculate the free fall timescale for the Sun, for a Red Supergiant with M=10M_o and R=1000R_o and a White Dwarf with M=1M_o and R=1000 km, and a Neutron star with M=1M_o and radius R=10km. Compare also their average densities.
• Skim MESA-web paper by Fields et al. 2022.
• Using MESA-web make a 1 M_o star until age 4.5×10⁹ years (a very rough model of the Sun as it is today!). Hint: look a the Custom Stopping condition menu for possible definitions of “end of the simulation”. Plot m(r), make sure to label your axes properly (including units!). Are there other variables with a qualitatively similar behavior that one could use as independent coordinate for the stellar structure? Try to make other plots to find some, and explain what is the mathematical property that allows to use m(r) and or any other variable you found as a coordinate.
• With the model above, check the central pressure of the star (you can also plot P(m) and P(r), or look at the final frame in the movie made by MESA-web for you) and compare it with the estimate above and the one provided in Onno Pols’ lecture notes.
• Check also the outer luminosity: is it the value you expected?