PROBLEM SET #5

For the problems below, you need to write programs using C language. Turn in your own source programs written independently (by email) together with reports.

1.

Lane-Emden equation: A polytrope refers to a self-gravitating spherical object in hydrostatic equilibrium. Its density distribution is described by $\Theta (\xi)$ that satisfies the Lane-Emden equation

(1)
\begin{align} {{1} \over {\xi^2}}{{d} \over {d \xi}} \left( \xi^2 {{d \Theta} \over {d \xi}} \right) = - \Theta^n ( \xi ) \end{align}

where $\xi$ denotes the dimensionless radius and $n$ is the polytropic index. The proper boundary conditions are $\Theta = 1$ and $d \Theta / d \xi = 0$ at $\xi = 0$.

(a)

One can seek for a power series solution

(2)
\begin{align} \Theta(\xi) = \sum^\infty_{m=0} a_m \xi^m \end{align}

of equation (1) for arbitrary index $n$, where $a_m$ ’s are coefficients. Find the values of $a_0, a_1, a_2, a_3,$ and $a_4$ in terms of $n$.

(b)

Write a program to solve equation (1) for arbitrary $n$. Check your codes for $n = 0, 1, 5$ by comparing your numerical results with the analytic solutions:

(3)
\begin{align} \Theta_{n=0} = 1 - {1 \over 6} \xi^2, \quad \Theta_{n=1} = {{\sin \xi} \over \xi}, \quad \Theta_{n=5} = \left( 1 + {1 \over 3} \xi^2 \right)^{-{1 \over 2}} \end{align}

(c)

Make a table that shows the values of the first zeros, $\xi_1$ , and $\left. d\Theta / d \xi \right\vert_{\xi_1}$ $n = 0, 0.5, 1.0, 1.5, \cdots , 4.0, 4.5, 4.9.$ Your answers should be accurate to 6 digits at least.

2.

Friedmann Universe: In cosmology, the cosmic expansion is usually described by the scale factor $a(t)$. By definition, $a = 0$ at the time of the Big Bang (i.e., $t = 0$), and $a = 1$ at the present time $t_0$. The scale factor is related to the redshift $z$ through

(4)
\begin{align} a = {{1} \over {1 + z}}, \end{align}

and also to the Hubble parameter

(5)
\begin{align} H = {{\dot{a}} \over {a}}, \end{align}

where the dot represents a time derivative. The Friedmann equation that governs the expansion of our Universe is given by

(6)
\begin{align} H^2 (z) = {H_0}^2 \left[ \Omega_R (1 + z)^4 + \Omega_M (1 + z)^3 + \Omega_k (1 + z)^2 + \Omega_\Lambda \right], \end{align}

where $H_0 = 71\ \mathrm{km\ s^{−1}\ Mpc^{−1}}$ is the Hubble constant at the present time, $\Omega_R, \Omega_M, \Omega_k,$ , and $\Omega_\Lambda$ denote the $\Omega$ parameters for radiation, matter, curvature, and dark energy, respectively, with the condition that $\Omega_R + \Omega_M + \Omega_k + \Omega_\Lambda = 1.$

(a)

Suppose a fictitious, mass-dominated universe with $\Omega_R = \Omega_k = \Omega_\Lambda = 0$ and $\Omega_M = 1.$ Solve equation (6) numerically, and compare your results with the analytic solution

(7)
\begin{align} a(t) = \left( {{3 H_0} \over {2}} t \right)^{2 \over 3}. \end{align}

Note that the age of the universe is $t_0 = 2/(3H_0).$

(b)

The current estimates of the $\Omega$ parameters in our Universe are ΩR = 3 × 10−5 , Ωk = 0 (flat universe), ΩM = 0.27 and ΩΛ = 0.73. Solve equation (6) numerically, and plot a as a function of time $t$. What is the current age of the Universe in Gyr?

(c)

When did the Universe switch from decelerating to accelerating expansion? What were the values of the the redshift and Hubble parameter at this time?

3.

Leap-frog Integrator: Consider a binary system consisting of two stars with masses $m_1 = 1$ and $m_2 = 0.5$ placed in the x-y plane. Initially$(t = 0)$, the two stars are located at $r_1 = (x_1 , y_1 ) = (−0.5, 0)$ and $r_2 = (1, 0),$ and have velocities of $v_1 = (0.01, 0.05)$ and $v_2 = (0.02, 0.2).$ The two stars orbit with each other due to the mutual gravity. The relevant equation of motion is

(8)
\begin{align} \ddot{\vec{r}_i} = - m_j {{\vec{r}_i - \vec{r}_j} \over {\left\vert \vec{r}_i - \vec{r}_j \right\vert^3 }}, \quad (i \ne j). \end{align}

for $i, j = 1$ or $2.$

(a)

Integrate equation (8) from $t = 0$ to $t = 50$ using the Leap-frog integrator, and plot the orbits of the two stars in the x-y plane. You need to choose a small enough dt for accurate orbit calculations.

(b)

Indicate the motion of the center of mass in the figure you draw in Part (a).

(c)

Plot the total energy defined by

(9)
\begin{align} E = {1 \over 2} m_1 {\vec{v}_1}^2 + {1 \over 2} m_2 {\vec{v}_2}^2 - {1 \over 2} {{m_1 m_2} \over { \left\vert \vec{r}_1 - \vec{r}_2 \right\vert }}. \end{align}

over $t = 0 − 1000$. Comment on the accuracy of the calculated orbits in terms of the energy conservation.