Next semester is approaching fast, and I haven’t had time to dabble with code much so far this summer. So, in order not to take on more than I can chew, I decided I wanted to revisit an assignment from the Computational Physics course. Using numerical integration, we had to find the mass-radius relationship of white dwarf stars. In the assignment, we had to implement the “standard” fourth-order Runge-Kutta scheme, known as RK4. So I thought it would be fun to try to extend it by using adaptive Runge-Kutta schemes instead. In this first part, I’ll provide some of the background and theory needed. The next part will contain the implementation and simulation results.
White dwarf stars
I won’t perform the derivations of the equations used, as that is not the main focus of this post. However a short introduction to white dwarf stars might be in order. A star such as our sun survive by burning the hydrogen in their core. Immense gravitational forces, caused by the matter of the star itself, squeezes the core so much that hydrogen fuses into helium. This nuclear fusion releases heat which provides an outwards pressure. This pressure counteracts the gravitational force, leading to a fine balance.
When the fuel runs out, there is nothing to stop gravity. The star gets crushed under its own weight, so to speak. The star will continue to shrink until the atoms making up the star are squashed together and quantum mechanical effects become dominant. What happens next is weird. The electrons try to crowd together as best as they can, however like people in a bus, only one electron can occupy a specific quantum state (bus seat) at time. This is called the Pauli exclusion principle. When all the states are occupied, the other electrons have nowhere to go (they don’t like to stand). This creates a kind of force, preventing the star from collapsing further.
Unless the star is much more massive than our sun, then there’s not enough mass to overcome this force, and the star will remain a small but massive ember.
From forces to equations
For this project, we want to find the radius of the white dwarf by considering the forces inside the white dwarf. As mentioned, I’m basing this on an earlier assignment which can be found here. It contains a brief derivation of the equations. Assuming the star is in an equilibrium, the pressure must counterbalance the gravitational force. The gravitational force depends on the mass, which again depends on the density of the star. We end up with the following set of coupled first-order ordinary differential equations \[\frac{d\rho}{dr} = -\left(\frac{dP}{d\rho}\right)^{-1}\frac{Gm}{r^2}\rho\] \[\frac{dm}{dr} = 4\pi r^2 \rho.\] Here \(\rho = \rho(r)\) is the mass density of a small volume at a distance \(r\) from the center of the star, \(m = m(r)\) is the integrated (total) mass within a radius \(r\), \(P\) is the pressure and \(G\) is the graviational constant. The radius \(R\) of the star is found when \(\rho(R) = 0\), giving the star a mass \(M\) of \(M = m(R)\). These equations relate the change in density and the change in mass.
Since several quantities involved are either very large or very small, the equations should be written on dimensionless form before they can be implemented. Otherwise one coulde experience precision problems due to the limits of the floating-point representiation used by computers. On dimensionless form the equations read \[\label{eq:drhodr}\frac{d\bar{\rho}}{d\bar{r}} = -\frac{\bar{m}}{\gamma}\frac{\bar{\rho}}{\bar{r}^2}\] \[\label{eq:dmdr} \frac{d\bar{m}}{d\bar{r}} = \bar{r}^2 \bar{\rho},\] where \(r = R_0\bar{r}\), \(m = M_0\bar{m}\), \(\rho = \rho_0\bar{\rho}\). \(\gamma\) is given as \[\gamma(x) = \frac{x^2}{3\sqrt{1+x^2}},\] where \(x = (\rho / \rho_0)^{1/3} = \bar{\rho}^{1/3}\). The constants \(R_0\), \(M_0\) and \(\rho_0\) depend on \(Y_e\), the number of electrons per nucleon. This is an input parameter, indicating which element the white dwarf consists of. For iron (\(^{56}\text{Fe}\)) this is \(26/56\). The constants are given as follows \[R_0 = 7.72 \times 10^6 Y_e \text{ m}\] \[M_0 = 5.67 \times 10^{30} Y_e \text{ kg}\] \[\rho_0 = 9.79 \times 10^8 Y_e^{-1} \text{ kg m}^{-3}.\]
General implementation issues
When integrating the above equations, we start at the center of the star and integrate outwards until we reach the stopping condition \(\rho(r) = 0\). Due to the nature of floating point numbers we should instead to stop when \(\rho(r) = \epsilon\), where \(\epsilon\) is a small positive number. I used \(\epsilon = 10^{-9}\) in my simulations. When \(r = 0\), ie the center of the star, \(\bar{m}(r) = 0\) and the density \(\bar{\rho}(0)\) is given by \(\bar{\rho}_c\), the core density. This is an input parameter which we’ll vary.
Unfortunately we can’t start integrating at the core just like that. The reason is that \(\ref{eq:drhodr}\) is not defined when \(r = 0\). This means that we need to start the integration at some small initial radius \(h\). In order to do so we need to approximate \(\bar{\rho}(h)\) and \(\bar{m}(h)\). Using a backwards Euler scheme, we can find an approximation of the initial values. From \(\ref{eq:drhodr}\) and \(\ref{eq:dmdr}\) we get the following scheme \[\label{eq:be_rho}\frac{\bar{\rho}(h)-\bar{\rho}(0)}{h} = -\frac{\bar{m}(h)}{\gamma(x(h))}\frac{\bar{\rho}(h)}{h}\]\[\label{eq:be_m}\frac{\bar{m}(h) – \bar{m}(0)}{h} = h^2\bar{\rho}(h).\] If we rewrite these equations we get \[\label{eq:be_rho}\frac{\bar{\rho}(h)-\bar{\rho}(0)}{h} = -\frac{\bar{m}(h)}{\gamma(x(h))}\frac{\bar{\rho}(h)}{h}\] \[\label{eq:be_m}\frac{\bar{m}(h) – \bar{m}(0)}{h} = h^2\bar{\rho}(h).\] If we rewrite these equations we get \[\label{eq:be_rho2}\left(1 + \frac{\bar{m}(h)}{\gamma(x(h))}\right)\bar{\rho}(h) – \bar{\rho}(0) = 0\] \[\label{eq:be_m2}\bar{m}(h) – h^3\bar{\rho}(h) – \bar{m}(0) = 0.\]We see that \(\ref{eq:be_m2}\) implies that \[\label{eq:be_mh}\bar{m}(h) = h^3\bar{\rho}(h) + \bar{m}(0),\]so we can insert this into \(\ref{eq:be_rho2}\) and using \(\bar{\rho}(0) = \bar{\rho}_c\) and \(\bar{m}(0) = 0\) we get \[\label{eq:be_rho3}\left(1 + h^3\frac{\bar{\rho}(h)}{\gamma(x(h))}\right)\bar{\rho}(h) – \bar{\rho}_c = 0.\] From this we can see that if we select a very small initial radius \(h\) then \(\bar{\rho}_c\) is actually a very good approximation to \(\bar{\rho}(h)\). Thus we set \(\bar{\rho}(h) \approx \bar{\rho}_c\). We then insert this approximation into \(\ref{eq:be_mh}\) to find the missing \(\bar{m}(h)\), and we can start the integration!
There remains only one small issue: \(\gamma(x)\) is not defined for negative values of \(\bar{\rho}\) (we’re not dealing with imaginary stars here 🙂 ). Even though we stop the integration when \(\bar{\rho}\) is close to zero, we run the risk of overstepping during the integration step. In my original implementation I used \(\text{max}(\bar{\rho}, 10^{-12})\) instead of \(\bar{\rho}\) when calculating \(\gamma\).
Integrating the equations
At this point, all that remains is to implement an integrator and off you go. In the assignment we had to implement and use the common fourth-order Runge-Kutta scheme, RK4. As mentioned I wanted to try to implement an adaptive RK scheme. In the next post I’ll hopefully provide a working implementation. As an inbetween snack, I leave you with this beautiful image of the white dwarf in NGC 2440 (small dot in center), surrounded by its ghostly remains.
NGC 2440 - Credits: NASA, ESA, and The Hubble Heritage Team (STScI/AURA)
One Response to Adaptive Runge-Kutta methods – Part 1