April 20, 2017

Abstract

This thorn solves the Tolman-Oppenheimer-Volkov equations of hydrostatic equilibrium for a spherically symmetric static star.

The Tolman-Oppenheimer-Volkoﬀ solution is a static perfect ﬂuid “star”. It is frequently used as a test of relativistic hydro codes. Here it is intended for use without evolving the matter terms. This provides a compact strong ﬁeld solution which is static but does not contain singularities.

The equations for a TOV star [1, 2, 3] are usually derived in Schwarzschild coordinates. In these coordinates, the metric can be brought into the form

$$d{s}^{2}=-{e}^{2\varphi}d{t}^{2}+{\left(1-\frac{2m}{r}\right)}^{-1}d{r}^{2}+{r}^{2}d{\Omega}^{2}\phantom{\rule{1em}{0ex}}.$$ | (1) |

This thorn is based on the notes of Thomas Baumgarte [4] that have been partially included in this documentation. However the notation for the ﬂuid quantities follows [5]. Here we are assuming that the stress energy tensor is given by

$${T}^{\mu \nu}=\left(\mu +P\right){u}^{\mu}{u}^{\nu}+P{g}^{\mu \nu},$$ | (2) |

where $\mu $ is the total energy, $P$ the pressure, ${u}^{\mu}$ the ﬂuid four velocity, $\rho $ the rest-mass density, $\mathit{\epsilon}$ the speciﬁc internal energy, and

$$\begin{array}{rcll}\mu & =& \rho \left(1+\mathit{\epsilon}\right),& \text{(3)}\text{}\text{}\\ P& =& \left(\Gamma -1\right)\rho \mathit{\epsilon},& \text{(4)}\text{}\text{}\\ P& =& K{\rho}^{\Gamma}.& \text{(5)}\text{}\text{}\end{array}$$

This enforces a polytropic equation of state. We note that in Cactus the units are $c=G={M}_{\odot}=1$.

The equations to give the initial data are solved (as usual) in the Schwarzschild-like coordinates with the areal radius labelled $r$. The equations of the relativistic hydrostatic equilibrium are

$$\begin{array}{rcll}\frac{dP}{dr}& =& -\left(\mu +P\right)\frac{m+4\pi {r}^{3}P}{r\left(r-2m\right)},& \text{(6)}\text{}\text{}\\ \frac{dm}{dr}& =& 4\pi {r}^{2}\mu ,& \text{(7)}\text{}\text{}\\ \frac{d\varphi}{dr}& =& \frac{m+4\pi {r}^{3}P}{r\left(r-2m\right)}\phantom{\rule{1em}{0ex}}.& \text{(8)}\text{}\text{}\\ & & & \text{(9)}\text{}\text{}\end{array}$$

Here $m$ is the gravitational mass inside the sphere radius $r$, and $\varphi $ the logarithm of the lapse. Once the integration is done for the interior of the star we match to the exterior (see below). In the exterior we have

$$\begin{array}{lll}\hfill P& =TOV\text{\_}atmosphere,\phantom{\rule{2em}{0ex}}& \hfill \text{(10)}\\ \hfill m& =M,\phantom{\rule{2em}{0ex}}& \hfill \text{(11)}\\ \hfill \varphi & =\frac{1}{2}log\left(1-2M\u2215r\right).\phantom{\rule{2em}{0ex}}& \hfill \text{(12)}\end{array}$$In order to impose initial data in cartesian coordinates, we want to transform this solution to isotropic coordinates, in which the metric takes the form

$$d{s}^{2}=-{e}^{2\varphi}d{t}^{2}+{e}^{2\Lambda}\left(d{\stackrel{\u0304}{r}}^{2}+{\stackrel{\u0304}{r}}^{2}d{\Omega}^{2}\right)\phantom{\rule{1em}{0ex}}.$$ | (13) |

Here $\stackrel{\u0304}{r}$ denotes the isotropic radius. Matching the two metrics, one obviously ﬁnds

$$\begin{array}{lll}\hfill {r}^{2}& ={e}^{2\Lambda}{\stackrel{\u0304}{r}}^{2}\phantom{\rule{1em}{0ex}},\phantom{\rule{2em}{0ex}}& \hfill \text{(14)}\\ \hfill {\left(1-\frac{2m}{r}\right)}^{-1}d{r}^{2}& ={e}^{2\Lambda}d{\stackrel{\u0304}{r}}^{2}\phantom{\rule{1em}{0ex}}.\phantom{\rule{2em}{0ex}}& \hfill \text{(15)}\end{array}$$As a result, we have an additional diﬀerential equation to solve in order to have $\stackrel{\u0304}{r}\left(r\right)$, that is

$$\frac{d\left(log\left(\stackrel{\u0304}{r}\u2215r\right)\right)}{\partial r}=\frac{{r}^{1\u22152}-{\left(r-2m\right)}^{1\u22152}}{r{\left(r-2m\right)}^{1\u22152}}\phantom{\rule{1em}{0ex}}.$$ | (16) |

Given such a solution, the missing metric potential is simply given by

$${e}^{\Lambda}=\frac{r}{\stackrel{\u0304}{r}}\phantom{\rule{1em}{0ex}}.$$ | (17) |

In the following section we concentrate on solving Eq. (16) in the exterior and in the interior of the star.

Then, given these one-dimensional data we interpolate to get data on the three-dimensional Cactus grid; that is, we interpolate on the three dimensional r given by the x, y, z variables the physical hydro and spacetime quantities that are function of the isotropic radius $\stackrel{\u0304}{r}$ computed above. Only linear interpolation is used. This avoids problems at the surface of the star, and does not cause problems if the number of points in the one dimensional array is suﬃcient ($1\times 1{0}^{5}$ is the default, which should be suﬃcient for medium-sized grids).

In the exterior of the star, $r>R$, the mass $M\equiv m\left(R\right)$ is constant, and Eq. (16) can be solved analytically up to a constant of integration. Fixing this constant such that $r$ and $\stackrel{\u0304}{r}$ agree at inﬁnity, we ﬁnd

$$\stackrel{\u0304}{r}=\frac{1}{2}\left(\sqrt{{r}^{2}-2Mr}+r-M\right)\phantom{\rule{1em}{0ex}},$$ | (18) |

or, solving for $r$ [cfr. Exercise 31.7 of MTW [3]]

$$r=\stackrel{\u0304}{r}{\left(1+\frac{M}{2\stackrel{\u0304}{r}}\right)}^{2}\phantom{\rule{1em}{0ex}}.$$ | (19) |

The metric potential as a function of $\stackrel{\u0304}{r}$ is obviously

$${e}^{2\Lambda}={\left(1+\frac{M}{2\stackrel{\u0304}{r}}\right)}^{2}\phantom{\rule{1em}{0ex}}.$$ | (20) |

In the interior, Eq, (16) can not be integrated analytically, because $m$ is now a function of $r$. Instead, we have to integrate

$${\int}_{0}^{\stackrel{\u0304}{r}}\frac{d\stackrel{\u0304}{r}}{\stackrel{\u0304}{r}}={\int}_{0}^{r}{\left(1-\frac{2m}{r}\right)}^{-1\u22152}\frac{dr}{r}\phantom{\rule{1em}{0ex}}.$$ | (21) |

The left hand side can be integrated analytically, and has a singular point at $\stackrel{\u0304}{r}=0$. The right hand side cannot be integrated analytically, but will also be singular at $r=0$, which poses problems when trying to integrate the equation numerically. We therefore rewrite the right hand side by adding and subracting a term $1\u2215r$, which yields

Since $m\sim {r}^{3}$ close to the origin, the ﬁrst term on the right hand side is now regular and the second one can be integrated analytically. As a result, we ﬁnd

Replacing the lower limits ($r=\stackrel{\u0304}{r}=0$) temporarily with ${r}_{0}$ and ${\stackrel{\u0304}{r}}_{0}$, we can integrate the right hand side and ﬁnd

or

$$\stackrel{\u0304}{r}=Crexp\left[{\int}_{0}^{r}\frac{1-{\left(1-2m\u2215r\right)}^{1\u22152}}{r{\left(1-2m\u2215r\right)}^{1\u22152}}dr\right]\phantom{\rule{1em}{0ex}}.$$ | (25) |

Here the constant of integration $C$ is related to the ratio ${\stackrel{\u0304}{r}}_{0}\u2215{r}_{0}$ evaluated at the origin (which is perfectly regular). It can be chosen such that the interior solution matches the exterior solution at the surface of the star. This requirement implies

In this respect, let us recall how we do initial data for the system of equations at $r=0$. Given a value of the central density ${\rho}_{c}$ (in Cactus units) we pose ${P}_{0}=K{\rho}_{c}^{\Gamma}$, ${m}_{0}=0$, ${\varphi}_{0}=0$ and ${\stackrel{\u0304}{r}}_{0}=TOV\text{\_}Tiny$, ${r}_{0}=TOV\text{\_}Tiny$. The TOV_Tiny number is hardwired into the code to avoid divide by zero errors; it is $1{0}^{-20}$. Also the default parameters will give the TOV star used for the long term evolutions in [6]. That is, a nonrotating $N=1$ ($\gamma =1+1\u2215N=2$) polytropic star with gravitational mass $M=1.4{M}_{\odot}$, circumferential radius $R=14.15$km, central rest-mass density ${\rho}_{c}=1.28\times 1{0}^{-3}$ and $K=100$.

To use this thorn to provide initial data for the ADMBase variables $\alpha $, $\beta $, $g$ and $K$ just activate the thorn and set ADMBase:initial_data = ‘‘TOV’’.

There are two ways of coupling the matter sources to the thorn that evolves the Einstein equations. One is to use the CalcTmunu interface. This will give the components of the stress energy tensor pointwise across the grid. For an example of this, see thorn ADM in CactusEinstein.

To use the CalcTmunu interface you should

- put the lines
friend: ADMCoupling

USES INCLUDE: CalcTmunu.inc

USES INCLUDE: CalcTmunu_temps.inc

USES INCLUDE: CalcTmunu_rfr.incin your interface.ccl

- In any routine requiring the matter terms, put
- #include ‘‘CalcTmunu_temps.inc’’ in the variable declarations
- declare CCTK_REALs Ttt, Ttx, Tty, Ttz, Txx, Txy, Txz, Tyy, Tyz, Tzz.
- Inside an i,j,k loop put #include ‘‘CalcTmunu.inc’’. This must be a Fortran routine (We could probably ﬁx this if requested).

- You then use the real numbers Ttt etc. as the stress energy tensor at a point.

As an alternative you can use the grid functions StressEnergytt, StressEnergytx, etc. directly to have the stress energy tensor over the entire grid. To do this you just need the line friend: ADMCoupling in your interface.ccl. Although this seems much simpler, you will now only get the contributions from the TOVSolver thorn. If you want to use other matter sources, most of the current thorns (CosmologicalConstant, the hydro code, the scalar ﬁeld code) all use the CalcTmunu interface.

You also have the possibility to use a parameter GRHydrotovsolver::TOV_Separation to obtain a spacetime consisting of one
TOV-system for $x>0$ and a
second (similar) for $x<0$.
This parameter sets the separation of the centers of two neutron stars, has to be positive and should be larger
than twice the radius of one star.

Be aware that the spacetime obtained by this is no physical spacetime and no solution of Einsteins Equations
and therefore an IVP-run has to follow. This parameter was only introduced for testing purposes of the
IVP-Solver and should only be considered as such. There would be better (and also easy) ways to obtain initial
data for two TOVs than that.

[1] R. C. Tolman, Phys. Rev. 55, 364 (1939).

[2] J. R. Oppenheimer and G. Volkoﬀ, Physical Review 55, 374 (1939).

[3] C.W. Misner, K.S. Thorn and J.A. Wheeler, Gravitation (Freeman and co. NY, 1973).

[4] T. W. Baumgarte. There is a copy of his notes in this directory:

TOVSolver/doc.

[5] J. A. Font, M. Miller, W. Suen and M. Tobias, Phys. Rev. D61, 044011 (2000).

[6] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla, E. Seidel, N. Stergioulas, W. Suen and M. Tobias, Phys. Rev. D65, 084024 (2002).

pert_amplitude | Scope: private | REAL |

Description: Amplitude of perturbation
| ||

Range | Default: 0.0 | |

*:* | Anything
| |

perturb | Scope: private | BOOLEAN |

Description: Add density perturbation (you should solve the IVP if true)
| ||

Default: no | ||

tov_clear_initial_data | Scope: private | BOOLEAN |

Description: Clear initial data (spacetime)? Default is yes.
| ||

Default: yes | ||

tov_combine_method | Scope: private | KEYWORD |

Description: Which combine method should be used.
| ||

Range | Default: average | |

maximum | Take the maximum of rho and gxx as clue for the rest as clue.
| |

average | Take the average of all available parts.
| |

tov_conformal_ﬂat_three_metric | Scope: private | BOOLEAN |

Description: Use conformal factor to get the 3-metric ﬂat. default is no
| ||

Default: no | ||

tov_dr | Scope: private | REAL |

Description: The spacing in the radial direction on the 1d grid
| ||

Range | Default: 5.e-4 | |

(0.0:* | Greater than 0
| |

tov_enforce_interpolation | Scope: private | BOOLEAN |

Description: Enforce the interpolation of the data onto the Hydro GFs even without tov as
speciﬁed initial data
| ||

Default: no | ||

tov_fake_evolution | Scope: private | INT |

Description: Fake evolution by setting ID at every step
| ||

Range | Default: (none) | |

*:* | anything, 0 as oﬀ (default), everything else as on
| |

tov_fast_interpolation | Scope: private | BOOLEAN |

Description: Use faster interpolation algorithm? Default is yes.
| ||

Default: yes | ||

tov_gamma | Scope: private | REAL |

Description: The polytropic constant in P = K rhoĜamma
| ||

Range | Default: 2.0 | |

1.0: | The physical range at high Lorentz factors is [1,2], but otherwise
higher values of gamma can also be used
| |

tov_k | Scope: private | REAL |

Description: The polytropic constant in P = K rhoĜamma
| ||

Range | Default: 100.0 | |

(0.0:* | Greater than 0
| |

tov_momentum_psi_power | Scope: private | INT |

Description: Power of Psi to be multiplied with Jî for Mom
| ||

Range | Default: (none) | |

*:* | anything, 0 as default
| |

tov_num_radial | Scope: private | INT |

Description: The number of radial points for the ODE integration
| ||

Range | Default: 100000 | |

1:* | Greater than 0
| |

tov_num_tovs | Scope: private | INT |

Description: The number of TOVs
| ||

Range | Default: 1 | |

1:* | Greater than 0
| |

tov_populate_timelevels | Scope: private | INT |

Description: Populate that amount of timelevels
| ||

Range | Default: 1 | |

1:3 | 1 (default) to 3
| |

tov_position_x | Scope: private | REAL |

Description: Position of neutron star, x coordinate
| ||

Range | Default: 0.0 | |

*:* | real
| |

tov_position_y | Scope: private | REAL |

Description: Position of neutron star, y coordinate
| ||

Range | Default: 0.0 | |

*:* | real
| |

tov_position_z | Scope: private | REAL |

Description: Position of neutron star, z coordinate
| ||

Range | Default: 0.0 | |

*:* | real
| |

tov_properposition | Scope: private | BOOLEAN |

Description: For use only with two NSs, atm only handles equal mass
| ||

Default: no | ||

tov_rho_central | Scope: private | REAL |

Description: The central density
| ||

Range | Default: 1e-3 | |

(0.0:* | Must be positive
| |

tov_save_to_dataﬁle | Scope: private | STRING |

Description: Only save data to ﬁle and exit
| ||

Range | Default: (none) | |

.* | Any ﬁlename, not used if empty
| |

tov_solve_for_tovs | Scope: private | INT |

Description: Solve for TOVs even if no TOV initial data was requested?
| ||

Range | Default: 3 | |

0:3 | ”depreciated in favour of TOVSolver::TOV_Enfor
ce_Interpolation”
| |

tov_use_old_initial_data | Scope: private | BOOLEAN |

Description: Take old initial data into account (spacetime)? Default is no.
| ||

Default: no | ||

tov_use_old_matter_initial_data | Scope: private | BOOLEAN |

Description: Use also old matter initial data? Default is no.
| ||

Default: no | ||

tov_velocity_x | Scope: private | REAL |

Description: (ﬁxed) Velocity of neutron star, x coordinate (caution!)
| ||

Range | Default: 0.0 | |

*:* | real
| |

tov_velocity_y | Scope: private | REAL |

Description: (ﬁxed) Velocity of neutron star, y coordinate (caution!)
| ||

Range | Default: 0.0 | |

*:* | real
| |

tov_velocity_z | Scope: private | REAL |

Description: (ﬁxed) Velocity of neutron star, z coordinate (caution!)
| ||

Range | Default: 0.0 | |

*:* | real
| |

conformal_storage | Scope: shared from STATICCONFORMAL | KEYWORD |

Implements:

tovsolver

Inherits:

admbase

hydrobase

constants

staticconformal

Uses header:

constants.h

Provides:

Set_Rho_ADM to

Set_Momentum_Source to

Set_Initial_Guess_for_u to

Rescale_Sources to

This section lists all the variables which are assigned storage by thorn EinsteinInitialData/TOVSolver. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

NONE

CCTK_PARAMCHECK

tov_c_paramcheck

check parameters

Language: | c | |

Options: | global | |

Type: | function | |

CCTK_WRAGH

tov_c_allocatememory

allocate memory for tovsolver_c

Language: | c | |

Options: | global | |

Type: | function | |

CCTK_POSTPOSTINITIAL (conditional)

tov_c_freememory

free memory from tovsolver_c

Language: | c | |

Options: | global | |

Type: | function | |

HydroBase_Initial (conditional)

tov_initial_data

group for the tov initial data

Sync: | admbase::metric | |

admbase::curv | ||

admbase::lapse | ||

admbase::shift | ||

rho | ||

press | ||

eps | ||

vel | ||

w_lorentz | ||

Type: | group | |

TOV_Initial_Data (conditional)

tov_c_integrate_rhs

integrate the 1d equations for the tov star

Language: | c | |

Options: | global | |

Type: | function | |

TOV_Initial_Data (conditional)

tov_write_1d_dataﬁle

save data to ﬁle and exit

After: | tov_c_integrate_rhs | |

Before: | tov_c_exact | |

Language: | c | |

Options: | global | |

Type: | function | |

TOV_Initial_Data (conditional)

tov_set_properpositions

steer ns position parameters according to proper distance

After: | tov_c_integrate_rhs | |

Before: | tov_c_exact | |

Language: | c | |

Options: | global | |

Type: | function | |

TOV_Initial_Data (conditional)

tov_c_exact

set up the 3d quantities for the tov star

After: | tov_c_integrate_rhs | |

Language: | c | |

Type: | function | |

(conditional)

tov_prepare_fake_evolution

prepare for fake evolution

After: | tov_c_exact | |

Language: | c | |

Type: | function | |

MoL_PostStep (conditional)

tov_c_exact

use fake evolution

After: | hydrobase_poststep | |

Language: | c | |

Type: | function | |