Overview

I will create a model of white dwarf stellar structure which produces a family of equilibrium configurations that is in turn used to explore the fascinating mass-radius relationship in white dwarfs. I will also use the model to find the Chandrasekhar limit, which gives a maximum limit to the possible mass of a stable white dwarf star.


Background

The fate of 99% of all stars, including our own Sun, is to end their life cycles as white dwarf stars. White dwarfs are essentially stellar core remnants, usually arising as the feature that is left over after a red giant has shed its outer layers to form a planetary nebula. At this point in its life cycle, the star has used up all of its fuel and is no longer able to carry out the nuclear fusion reactions which had dominated its energy generation during its time on the main sequence. The material of the star is now a highly compressed collection of nuclei (usually carbon and oxygen ions) and electrons, forming an electron degenerate gas.
Due to their huge masses, all stars are subject to a large self-gravitational force which, if unopposed, would cause the star to collapse in upon itself. Main sequence stars, such as our Sun, use thermal pressure from the nuclear fusion reactions in their core to push against this inward gravitational force.

Because a white dwarf star is simply a core remnant leftover after a main sequence star has exhausted its fuel, nuclear fusion is no longer occurring. So in order to remain structurally stable, the white dwarf must have another source of pressure to push back against gravitational collapse.

This pressure comes from the fact that the the gas is made of fermions which are subject the Pauli exclusion principle. This principle states that no two fermions can occupy the same energy state. So as the electrons are packed more tightly together in the compressed stellar remnant which becomes the white dwarf, the lower energy levels quickly become filled and some electrons are forced into higher energy states. This higher energy imparts a momentum to the electrons which in turn leads to a pressure, termed electron degeneracy pressure. The gravitational collapse of the star is halted by this degeneracy pressure once the pressure and gravity are in equilibrium, and a hydrostatically stable white dwarf is formed.


The Model

Here I model white dwarf structure using a dimensionless set of equations which, after being numerically solved, can be converted back into meaningful parameters. I will choose radius as the independent variable and evaluate mass and density for a given core density value.

My model makes a few assumptions:

Assumption no. 1: The composition of the white dwarf is ionized 12C. Assumption no. 2: The white dwarf is spherical and non-rotating. Assumption no. 3: Effects due to the presence of magnetic fields are to be neglected. Assumption no. 4: The star is in hydrostatic equilibrium.

Due to assumption no. 4, the inward gravitational force acting upon each mass element must be equal to the outward force of the pressure acting on that mass element. Therefore, at any mass element in the star,

(1) \(F_g = F_p\)

From this starting point, one can then go through a lengthy process of deriving differential equations of state for the mass and density with respect to radius. (I'll refrain from posting the details in order to not get too physics-y, but I'm happy to add them in if anyone would like to see the full derivation.) In order to simplify the calculations (and because these kinds of units as fractional relations to solar constants are quite conventional in astrophysics), I will change to dimensionless variables, \(\overline{r}\), \(\overline{m}\), and \(\overline{\rho}\):

(2) \(\overline{r} = \frac{r}{R_\bigodot}\)

(3) \(\overline{m} = \frac{m}{M_\bigodot}\)

(4) \(\overline{\rho} = \frac{\rho}{\rho_\bigodot}\)

where \(R_\bigodot\), \(M_\bigodot\), and \(\rho_\bigodot\) represent the radius, mass, and core density of the Sun, thus making our dimensionless variables fractional percentages of the solar analog for each parameter.

The equations of state can then be expressed as:

(5) \(\frac{d\overline{m}}{d\overline{r}} = C_m\overline{r}^2\overline{\rho}\)

(6) \(\frac{d\overline{\rho}}{d\overline{r}} = \frac{C_\rho \overline{m} \overline{\rho}}{\overline{r}^2\Upsilon(\overline{\rho})}\)

where

(7) \(C_m = \frac{4\pi R_\bigodot ^3\rho_c}{M_\bigodot}\)

(8) \(C_\rho = \frac{-Gm_HM_\bigodot}{0.5m_ec^2R_\bigodot}\)

(9) \(\Upsilon(\overline{\rho}) = \frac{K_F^2\rho_c^{\frac{2}{3}}\overline{\rho}^{\frac{2}{3}}}{3(1+K_F^2\rho_c^{\frac{2}{3}}\overline{\rho}^{\frac{2}{3}})^{\frac{1}{2}}}\)

(10) \(K_F = (\frac{3\times0.5}{8\pi m_H})^{\frac{1}{3}}\frac{h}{m_ec}\)

Initial conditions for the radius and mass at the core of the star are:

(11) \(r_c = 0\)

(12) \(m_c = 0\)

which, using (2) and (3), become

(13) \(\overline{r}_c = 0\)

(14) \(\overline{m}_c = 0\)

An initial condition for the density at the core of the star is also needed and will be collected as an input value by the Python program. So,

(15) \(\rho_c = \rho_{input}\)

which, using (4), becomes

(16) \(\overline{\rho}_c = \frac{\rho_{input}}{\rho_\bigodot}\)

Euler's method can then be used in Python to numerically integrate the set of differential equations of state. I chose a step size of 0.000000001, which corresponds to 0.0000001% of \(1 R_\bigodot\). When I converted this back into meters, it looks like this step size ended up requiring a span of ~26 thousand to ~37 million steps for the models I ran, depending upon the \(\rho_c\) input. Wow! That certainly explains why the runtime was so long for lower input values of \(\rho_c\) (which ended up corresponding to larger radii). Aside from using a small step size to make up for the fact that Euler's formula is only of order one, part of the reason I used such a small value for \(h\) is due to the issue which follows below.

There is a problem with the initial conditions. Notice that in (6), \(\overline{r}\) appears in the denominator. So we can't start Euler's method with \(\overline{r} = 0\). To get around this hiccup, we instead start \(r\) one step size away at \(h\). We then assume that the density along this short distance \(h\) is constant (hence needing an especially small \(h\) value) and, keeping the initial value density as \(\rho_{input}\), we can derive an approximated initial value for the mass at one step size away from the core of the white dwarf. Our new initial conditions, \(\overline{r}_0\), \(\overline{m}_0\), and \(\overline{\rho}_0\), for the Euler's method algorithm become:

(17) \(\overline{r}_0 = h\)

(18) \(\overline{m}_0 = 4\pi \overline{r}_0^2 \overline{\rho}_0\)

(19) \(\overline{\rho}_0 = \frac{\rho_{input}}{\rho_\bigodot}\)

The numerical integration can now be performed. I ran the program with 45 different core density inputs, with each run generating the mass-radius pair which would keep the star in hydrostatic equilibrium given the core density that I entered for that particular run. The core densities I used ranged from \(10^6\)kg/m3 to \(10^{12}\)kg/m3. I entered the resulting data into Fathom, producing a graph which shows the possible mass-radius configurations for a white dwarf star:
The units of mass in the graph are in fractions of 1 solar mass (\(M_\bigodot\)), and the units of radius are in fractions of 1 Earth radius (\(R_\bigoplus \)).

I extended our Search program to handle 4 parameters in order to fit a curve to this data. The family of curves \(y = a((\frac{1}{x})^b - x^c)^d\), where \(a, b, c > 0\) and \(0 \lt d \lt 1\) provide a template for the curve fit, as can be seen here in Desmos:
Running the 4D Search program to solve \(y = a((\frac{1}{x})^b - x^c)^d\) for \(a\), \(b\), \(c\), and \(d\) produced:

(20) \(radius = 1.14((\frac{1.45615}{mass})^{0.74} - (\frac{mass}{1.45615})^{0.96})^{0.48}\)

where the mass is in units of fractional Earth mass and the radius is in units of fractional solar radius.

Plotting the function on Figure (1) gives a solid fit for the curve:
Two interesting characteristics stand out on this graph. First, the slope throughout the curve is negative. This means that the more massive the white dwarf, the smaller its radius! While this is counterintuitive to our everyday experiences with matter, this phenomenon is a natural result of the degenerate nature of the gas which makes up a white dwarf star. The effect arises because higher mass leads the electrons within the star to fill higher and higher energy states (the electrons being subject to the Pauli exclusion principle). This in turn forces the high energy electrons to have high momentum and thus raises the pressure of the gas. So unlike in a normal, non-degenerate gas, the pressure within a white dwarf is actually correlated to the particle density and is independent of the temperature. The degeneracy pressure is needed to balance against the inward self-gravitating force, so if higher mass gives higher pressure, then an equilibrium configuration would require a white dwarf with high mass to have a high gravitational force and therefore a smaller radius since the radius is inversely proportional to the self-gravitation force. Similarly, less mass leads to a lower pressure and so equilibrium is found at a higher radius since a lower gravitational force is required for the balancing act.

The second point of interest in Figure (3) is that the curve eventually hits the mass axis: there is a limit to how massive a white dwarf can be. This point is called the Chandrasekhar mass limit and marks the mass upper bound at which the degeneracy pressure of the Fermi gas is inadequate to support the star against its self-gravitating force. So if the mass reaches above this limit, the star will collapse and become either a neutron star or a black hole. The limit that I found with my model is \(~1.56M_\bigodot\). This is consistent with the accepted values for the Chandrasekhar limit, which I found to be between \(1.4M_\bigodot\) and \(1.6M_\bigodot\), depending upon the source I was reading.

As a final note, a quick glance at the graph shows that a white dwarf with the same amount of mass as the Sun has about the same radius as the Earth. So picture the entire Sun being packed into a sphere the size of the Earth. And more massive white dwarfs are packed into even smaller radii. These stars are extremely dense!