Equations of state for white dwarf stars

Question

I am dealing with some model describing the given composition of a white dwarf made with photons, degenerate electrons, positrons and ions. From that, I can get tabulated values of densities ranging from small to large numbers and corresponding to it the total pressure and total energy. Now I'm interested in getting the parameters like density, mass and pressure with respect to the radius at some given initial conditions. I can do it using two coupled ordinary differential equations, i.e.: $$ \frac{\partial m}{\partial r}=4\pi r^2\rho(r)\\ \frac{\partial P}{\partial r}=-\frac{Gm(r)\rho(r)}{r^2} $$ This can be solved, for example, using the forward Euler method which will give us as follows: $$ m(i)=m(i-1)+4\pi r^2 \rho(i-1)\Delta r\\ P(i)=P(i-1)-\frac{Gm(i-1)\rho(i-1)}{r^2}\Delta r $$ Firstly, I want to implement the following initial conditions i.e. $$ m(0)=0 \ \ ; \ \ \rho(0)=\rho_{central} $$ I'm going outward with my calculations. Using my tabulated values of pressure and density I can easily calculate for small $\Delta r$ how the mass and pressure change with respect to the radius breaking the computations when the density reaches a very small value (the surface of the WD). Having all of these I can determine what is the mass of the WD for a given central density.

Secondly, I want to use different initial conditions to get, in the end, the central density for a given mass, i.e. $$ m(R)=M \ \ ; \ \ \rho(R)=0 \ \ or/and \ \ P(R)=0 $$ It means that I have to do my calculations inward. In that case, I have to modify the above result of the Euler method to get: $$ m(i-1)=m(i)-4\pi r^2 \rho(i)\Delta r\\ P(i-1)=P(i)+\frac{Gm(i)\rho(i)}{r^2}\Delta r $$ And this is the moment I got stuck for several weeks. If I take my initial conditions which are the total mass, density and pressure are equal to zero at the surface then the only parameter that stops me from proceeding with these computations is the total radius. And here I need some help or hint. How can I perform these calculations inward when I don't know this radius? Is there maybe some way of transformation of getting rid of this unknown? Or perhaps I can do computations outward even though I've got surface boundary conditions?

Answer 1

You've described a boundary value problem: you're trying to solve a differential equation where different conditions apply at different boundaries. These are distinguished from initial value problems, where you're given all the conditions you need at some point (usually in time), from which you can solve the equations. The method you describe of starting from, say, the centre and integrating outwards would solve an initial value problem but, as you've realised, you need something else to solve a boundary value problem.

The simplest approach I'm aware of is the shooting method. Conceptually, what you do is

• guess the missing condition at one end (e.g. the central density),
• solve the equations outwards using whatever method¹ until you reach the surface condition (e.g. P(R)=0) and
• then compare your result with the condition you haven't considered yet (the total mass M(R)).

You can then iterate on central density by using a root-finding algorithm on the total mass less the target mass to find the value of the central density that works.

Put differently, the outward integration gives you a function that looks something like $M(\rho_c)$ and you do a root-find on $M(\rho_c)-M_\mathrm{target}$.

¹ You mentioned Euler, which is a fine place to start.