In-Class Assignment 13

In-Class Assignment 13#

Credit: Mike Zingale for the Polytrope Routines and Notes.

Comparing Polytropes to MESA models#

Learning Objectives#

explore the relationship between lane-emden solutions and physical variables
identify physics informed polytrope profiles for models of massive stars and white dwarfs

Lane-Emden polytrope equation#

Derivation#

We consider just the equations of hydrostatic equilibrium and mass continuity from our equations of stellar structure.

\[\frac{dP}{dr} = -\rho \frac{GM(r)}{r^2}\]

\[\frac{dM}{dr} = 4\pi r^2 \rho\]

We can close the system via an equation of state of the form: \(P(\rho)\), which we write as:

\[P = K \rho^{1+1/n}\]

where \(n\) is called the polytropic index.

With a bit of algebra, we can combine these equations into a single second-order ODE for density:

\[\left ( \frac{n+1}{n} \right ) \frac{K}{4\pi G} \frac{1}{r^2} \frac{d}{dr} \left ( \frac{r^2}{\rho^{(n-1)/n}} \frac{d\rho}{dr} \right ) = -\rho\]

and then make it dimensionless, but expressing the density in terms of the central density, \(\rho_c\):

\[\rho(r) = \rho_c \theta^n(r)\]

where we note that \(0 \le \theta \le 1\), and a lenghtscale \(\alpha\) such that \(r = \alpha \xi\):

\[\alpha^2 = \frac{(n+1)P_c}{4\pi G\rho_c^2}\]

Giving

\[\frac{1}{\xi^2}\frac{d}{d\xi} \left ( \xi^2 \frac{d\theta}{d\xi} \right ) = -\theta^n\]

The boundary conditions are:

\(\theta(\xi=0) = 1\)
\(d\theta / d\xi |_{\xi=0} = 0\) (\(g = 0\) at the center, so \(dP/dr = d\rho/dr = 0\))

System of first order equations#

We will rewrite this as 2 first order equations, taking \(y = \theta\), \(z = \theta^\prime\):

\[\begin{split} \begin{align} \frac{dy}{d\xi} &=& z \\ \frac{dz}{d\xi} &=& -\frac{2}{\xi}z - y^n \end{align} \end{split}\]

Then we have:

\(y(0) = 1\)
\(z(0) = 0\)

Notice that if we evaluate the system at \(\xi = 0\), then the \(dz/d\xi\) term blows up. We need to use an expansion there.

In the limit \(\xi \rightarrow 0\), the solution takes the form:

\[\theta(\xi) = 1 - \frac{1}{6}\xi^2 + \frac{n}{120} \xi^4 + \ldots\]

which allows us to simply the \(dz/d\xi\) near \(\xi = 0\) as:

\[\frac{dz}{d\xi} \approx -\frac{1}{3}\]

Stopping point#

But we have the problem that we do not know the stopping point, \(\xi_1\).

We estimate the radius, \(\xi_1\) each step and make sure that the next step does not take us past that estimate. This prevents us from having negative \(\theta\) values. Given a point \((\xi_p, \theta(\xi_p))\), and the derivative at that point, \(\theta^\prime_p = d\theta/d\xi |_{\xi_p}\), we can write the equation of a line as:

\[ \theta(\xi) - \theta_p = \theta^\prime_p (\xi - \xi_p) \]

Then we can ask when does \(\theta\) become zero, finding:

\[ \xi = -\frac{\theta_p}{\theta^\prime_p} + \xi_p \]

or in terms of \(y\) and \(z\),

\[ \xi = -\frac{y_p}{z_p} + \xi_p \]

This is our estimate of \(\xi_1\). We then make sure our stepsize \(h\) is small enough that we do not go beyond this estimate.

Implementation#

This is our main class that does the integration. We initialize it with the polytopic index, and then it will integrate the system for us. We can then plot it or get the parameters \(\xi_1\) and \(-\xi_1^2 d\theta/d\xi |_{\xi_1}\)

import matplotlib.pyplot as plt
import numpy as np

class Polytrope:
    """a polytrope of index n"""
    def __init__(self, n, h0=1.e-2, tol=1.e-12):
        self.n = n
        self.xi = []
        self.theta = []
        self.dtheta_dxi = []
        
        self._integrate(h0, tol)

    def _integrate(self, h0, tol):
        """integrate the Lane-Emden system"""

        # our solution vector q = (y, z)
        q = np.zeros(2, dtype=np.float64)
        xi = 0.0

        h = h0

        # initial conditions
        q[0] = 1.0
        q[1] = 0.0

        while h > tol:
            # 4th order RK integration -- first find the slopes
            k1 = self._rhs(xi, q)
            k2 = self._rhs(xi+0.5*h, q+0.5*h*k1)
            k3 = self._rhs(xi+0.5*h, q+0.5*h*k2)
            k4 = self._rhs(xi+h, q+h*k3)

            # now update the solution to the new xi
            q += (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)
            xi += h

            # set the new stepsize--our systems is always convex
            # (theta'' < 0), so the intersection of theta' with the
            # x-axis will always be a conservative estimate of the
            # radius of the star.  Make sure that the stepsize does
            # not take us past that.
            R_est = xi - q[0]/q[1]

            if xi + h > R_est:
                h = -q[0]/q[1]

            # store the solution:
            self.xi.append(xi)
            self.theta.append(q[0])
            self.dtheta_dxi.append(q[1])

        self.xi = np.array(self.xi)
        self.theta = np.array(self.theta)
        self.dtheta_dxi = np.array(self.dtheta_dxi)

    def _rhs(self, xi, q):
        """ the righthand side of the LE system, q' = f"""

        f = np.zeros_like(q)

        # y' = z
        f[0] = q[1]
        
        # for z', we need to use the expansion if we are at xi = 0,
        # to avoid dividing by 0
        if xi == 0.0:
            f[1] = (2.0/3.0) - q[0]**self.n
        else:
            f[1] = -2.0*q[1]/xi - q[0]**self.n

        return f

    def get_params(self):
        """ return the standard polytrope parameters xi_1,
        and [-xi**2 theta']_{xi_1} """
        xi1 = self.xi[-1]
        p2 = -xi1**2 * self.dtheta_dxi[-1]
        return xi1, p2

    def plot(self):
        """ plot the solution """
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.plot(self.xi, self.theta, label=r"$\theta$")
        ax.plot(self.xi, self.theta**self.n, label=r"$\rho/\rho_c$")
        ax.set_xlabel(r"$\xi$")
        ax.legend(frameon=False)
        return fig

a.#

Plot the results of an \(n=4\) polytrope.

## a result here

From the readings in HKT 7.1-7.2. We found some of the following relations:

Equation 7.20-7.22:

\[ \Large \rho(r) = \rho_{\rm{c}}\theta^{n}(r) \]

\[ \Large P(r) = P_{\rm{c}} \theta^{1+n}(r) \]

\[ \Large P_{\rm{c}} = K \rho^{1+1/n}_{\rm{c}} \]

and

Equation 7.25:

\[ \Large r^{2}_{n} = \frac{(n+1)P_{c}}{4\pi G \rho^{2}_{c}} \]

and

\[ \Large r = r_{n} \xi \]

b.#

Download the following model files locally. These data were produced using the make_co_wd test suite.

\(0.6 M_{\odot}\) WD: co_wd_profile.data;

Then,

produce a polytrope of index n
compute K using the central values from the MESA model (for pandas you can use iloc).
compute r_n using the central values from the MESA model
compute the pressure and density profiles from the above constants and the polytrope solutions
plot them on the same plot as the MESA model - one plot for pressue one for density
try different values of n until a best fit is obtained

## b results here
G = 6.67e-8

c.#

Explain in words, the physics motivating the choice of n in comparing to the white dwarf.

c results here.