In-Class Assignment 5

In-Class Assignment 5#

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

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
import pandas as pd

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

a. Plotting an \(n=1\) Polytrope#

  1. Call the Polytrope Class and obtain the results.

  2. Plot results for \(\theta_{n}(\xi)\) (\(\rho(r)/\rho_{\rm{c}}\)) as a function of the dimensionless radial coordinate \(\xi\).

Hint: You can access the varible of the class by calling the class object followed by .variable. For obtaining \(\xi\) that would be p.xi.

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

## 1 result here

## 2 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. Using Polytrope to determine the degeneracy regine of a WD MESA model#

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

Then,

  1. produce a polytrope of index n between 1 to 4.

  2. Load the MESA profile data for analysis.

  3. compute \(P_{\rm{c}}\) and \(\rho_{\rm{c}}\) using the central/max values from the MESA model.

  4. compute r_n using the central values from the MESA model

  5. use r_n to compute the physical radius (\(r=r_n \xi\)) of the polytrope

  6. Use \(\rho_{\rm{c}}\) from the MESA data to compute the density (\(\rho=\theta_{n}\rho_{\rm{c}}\))

  7. Load the radius \(r\) and density \(\rho\) from the MESA data

  8. plot (and label) them on the same plot as the MESA model from 0 to \(r\approx 1\times 10^9\) cm.

  9. try different values of \(n\) to determine if the WD is in the non-relativstic or relativistic region.

  10. Describe in a few sentences your results.

## 1 results here
# produce a polytrope here
G = 6.67e-8
msun = 2e33    # g
rsun = 6.96e10 # cm

## 2-3 results here
#co_wd_profile = pd.read_csv('data/co_wd_profile.data',sep=r'\s+',header=4)
#co_wd_profile.columns

# 4-6 results here - compute r_n here HKT 7.25, given above for reference

# 7 results here - it will be in units of Rsun so you will need to convert to cm !

# 8 results here - plot the converted polytrope data and the raw MESA data

Explain in words, which index you found to match best, what this states about the material, and why does the profile differ near the surface?

c results here.

Contents