Solution of the Tolman-Oppenheimer-Volkov equations

The class o2scl::tov_solve provides a solution to the Tolman-Oppenheimer-Volkov (TOV) equations given an equation of state, provided as an object of type o2scl::eos_tov. These classes are particularly useful for static neutron star structure: given any equation of state one can calculate the mass vs. radius curve and the properties of any star of a given mass.

Mathematical background:

In units where $ c=1 $, the most general static and spherically symmetric metric is of the form

\[ ds^2 = - e^{2 \Phi(r)} d t^2 + e^{2 \Lambda(r)} d r^2 + r^2 d \theta^2 + r^2 \sin^2 \theta~d \phi^2 \]

where $ \theta $ is the polar angle and $ \phi $ is the azimuthal angle. Often we will not write explicitly the radial dependence for many of the quantities defined below, i.e. $ \Phi \equiv \Phi(r) $.

This leads to the TOV equation (i.e. Einstein's equations for a static spherically symmetric object)

\[ \frac{d P}{d r} = - \frac{G \varepsilon m}{r^2} \left( 1+\frac{P}{\varepsilon}\right) \left( 1+\frac{4 \pi P r^3}{m} \right) \left( 1-\frac{2 G m}{r}\right)^{-1} \]

where $r$ is the radial coordinate, $m(r)$ is the gravitational mass enclosed within a radius $r$, and $\varepsilon(r)$ and $P(r)$ are the energy density and pressure at $r$, and $G$ is the gravitational constant. The mass enclosed is related to the energy density through

\[ \frac{d m}{d r} = 4 \pi r^2 \varepsilon \]

and these two differential equations can be solved simultaneously given an equation of state, $ P(\varepsilon) $. The total gravitational mass is given by

\[ M = \int_0^R 4 \pi r^2 \varepsilon d r \]

The boundary conditions are $m(r=0)=0$ and the condition $P(r=R)=0$ for some fixed radius $R$. These boundary conditions give a one-dimensional family solutions to the TOV equations as a function of the central pressure. Each central pressure implies a gravitational mass, $ M $ and radius $ R $ and thus defines a mass-radius curve.

The metric function $ \Lambda $ is

\[ e^{2 \Lambda} = \left( 1-\frac{2 G m}{r}\right)^{-1} \]

The other metric function, $\Phi(r)$ is sometimes referred to as the gravitational potential. In vacuum above the star, it is

\[ e^{2 \Phi} = \left( 1-\frac{2 G M}{r}\right) \]

and inside the star it is determined by

\[ \frac{d \Phi}{d r} = - \frac{1}{\varepsilon} \frac{ d P}{d r} \left(1+\frac{P}{\varepsilon}\right)^{-1} = \frac{G m}{r^2} \left( 1+\frac{4 \pi P r^3}{m} \right) \left( 1-\frac{2 G m}{r}\right)^{-1} \]

Note that this can be rewritten as

\[ -d \Phi = \frac{d P}{P+\varepsilon} \, . \]

If the neutron star is at zero temperature and there is only one conserved charge, (i.e. baryon number), then

\[ -d \Phi = \frac{d P}{\mu n} = \frac{d \mu}{\mu} \]

and this implies that $ \mu e^{\Phi} $ is everywhere constant in the star. Alternatively,

\[ -d \Phi = \frac{d P}{h} \]

where $ h $ is the enthalpy density. Since $ d h = T d s + d P $, then if the entropy is everywhere a constant we also have

\[ -d \Phi = \frac{d h}{h} \]

and $ h e^{\Phi} $ is everywhere constant (even if there is more than one conserved charge).

The proper boundary condition for the gravitational potential is

\[ \Phi(r=R) = \frac{1}{2} \ln \left(1-\frac{2 G M}{R} \right) \]

which ensures that $ \Phi(r) $ matches the metric above in vacuum. Since the expression for $ d\Phi/dr $ is independent of $ \Phi $, the differential equation can be solved for an arbitrary value of $ \Phi(r=0) $ and then shifted afterwards to obtain the correct boundary condition at $ r=R $ .

The surface gravity is defined to be

\[ g = \frac{G m}{r^2}\left(1-\frac{2 G m}{r}\right)^{-1/2} \]

which is computed in units of inverse kilometers, and the redshift is defined to be

\[ z = \left(1-\frac{2 G m}{r}\right)^{-1/2} - 1 \]

which is unitless.

The baryonic mass is typically defined by

\[ M_B = \int_0^R 4 \pi r^2 n_B m_B \left(1-\frac{2 G m}{r}\right)^{-1/2} d r \]

where $ n_B(r) $ is the baryon number density at radius $ r $ and $ m_B $ is the mass one baryon (taken to be the mass of the proton by default and stored in o2scl::tov_solve::baryon_mass). If the EOS specifies the baryon density (i.e. if o2scl::eos_tov::baryon_column is true), then this class will compute the associated baryonic mass for you.

In the case of slow rigid rotation with angular velocity $ \Omega $, the moment of inertia is

\[ I = \frac{8 \pi}{3} \int_0^R dr~r^4\left(\varepsilon+P\right) \left(\frac{\bar{\omega}}{\Omega}\right) e^{\Lambda-\Phi} = \frac{8 \pi}{3} \int_0^R dr~r^4\left(\varepsilon+P\right) \left(\frac{\bar{\omega}}{\Omega}\right) \left(1-\frac{2 G m}{r}\right)^{-1/2} e^{-\Phi} \]

where $ \omega(r) $ is the rotation rate of the inertial frame, $ \Omega $ is the angular velocity in the fluid frame, and $ \bar{\omega}(r) \equiv \Omega - \omega(r) $ is the angular velocity of a fluid element at infinity. The function $ \bar{\omega}(r) $ is the solution of

\[ \frac{d}{dr} \left( r^4 j \frac{d \bar{\omega}}{dr}\right) + 4 r^3 \frac{d j}{dr} \bar{\omega} = 0 \]

where the function $ j(r) $ is defined by

\[ j = e^{-\Lambda-\Phi} = \left( 1-\frac{2 G m}{r} \right)^{1/2} e^{-\Phi} \, . \]

Note that $ j(r=R) = 1 $. The boundary conditions for $ \bar{\omega} $ are $ d \bar{\omega}/dr = 0 $ at $ r=0 $ and

\[ \bar{\omega}(R) = \Omega - \left(\frac{R}{3}\right) \left(\frac{d \bar{\omega}}{dr}\right)_{r=R} \, . \]

One can use the TOV equation to rewrite the moment of inertia as

\[ I= \left(\frac{d \bar{\omega}}{dr}\right)_{r=R} \frac{R^4}{6 G \Omega} \, . \]

The star's angular momentum is just $ J = I \Omega $.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).