Equations and Numerical Methods

We wish to study the evolution of a fluid in the background spacetime of a Kerr (rotating) black hole. We adopt Boyer-Lindquist coordinates, $(t,r,\theta,\phi)$, for which the line element has the form,

$${ds}^2=g_{t t}\,{dt}^2+2\,g_{t \phi}\,{dt}\,{d \phi}+g_{r r}......+g_{\theta \theta}\,{d \theta}^2+g_{\phi \phi}\,{d \phi}^2 .$$ (1)

In keeping with Misner, Thorne, & Wheeler (1973), we use the metric signature (-,+,+,+), along with geometrodynamic units where G = c = 1; the black hole mass is unity, M=1. The determinant of the 4-metric is g, and $\sqrt{-g} =\alpha\,\sqrt{\gamma}$ where $\alpha$ is the lapse function, $\alpha=1/\sqrt{-g^{tt}}$, and $\gamma$ is the determinant of the spatial 3-metric.

For a relativistic test fluid described by a density $\rho$, specific internal energy $\epsilon$, and 4-velocity $U^\mu$, we define the transport velocity $V^\mu$ as follows:

$$V^\mu = {U^\mu \over U^t} ,$$ (2)

where $U^t = W/\alpha$, and W is the gravitational redshift factor,

$$W = {1 \over (1-V^\mu\,V_\mu)^{1/2}} .$$ (3)

We also define the momentum,

$$S_\mu = \rho\,h\,W\,U_\mu ,$$ (4)

and auxiliary density and energy functions$D = \rho\,W$ and $E = D\,\epsilon$.

Using these variables the conservation laws can be rewritten into a form suitable for finite differencing. The equation for mass conservation is written

$$\partial_t\,D + {1 \over\sqrt{\gamma}}\,\partial_j\,(D\,\sqrt{\gamma}\,V^j) = 0.$$ (5)

Conservation of the fluid energy-momentum tensor $\nabla_\mu {T}^{\mu \nu} = 0$ yields momentum conservation equations,

$$\partial_t\,S_i + {1 \over \sqrt{\gamma}}\,\partial_j\,(S_i......over 2}\,{S_\mu\,S_\nu \over S^t}\partial_j\,g^{\mu\,\nu}= 0,$$ (6)

and an internal energy conservation equation

$$\partial_t\,E + {1 \over\sqrt{\gamma}}\,\partial_j\,(E\,\sq......\over \sqrt{\gamma}}\,\partial_j\,(W\,\sqrt{\gamma}\,V^j) = 0.$$ (7)

Spatial indices are indicated by roman characters i,j=1,2,3. We assume an ideal gas equation of state $P=\rho\,\epsilon\,(\Gamma-1)$, where$\Gamma$ is the adiabatic exponent. This is the same system of equations as described in HSWa.

The GR hydrodynamics code evolves time-explicit, operator-split, finite difference forms of equations (5)--(7). The algorithm is a three-dimensional generalization of the solver described in HSWb, with some additional modifications. For example, velocity renormalization, $U^\mu\,U_\mu= -1$, which is invoked in the source and transport steps following the solution of the discretized momentum equation (6) is now implemented using the algebraically equivalent condition $S^\mu\,S_\mu = -{(D + \Gamma\,E)}^2$ for improved numerical stability. Several improvements are possible because of the increase in computer power since 1984. Analytically determined metric terms and metric derivatives are now calculated and stored at all needed grid locations rather than averaged. Array operations have also been streamlined wherever possible using Fortran 90 array syntax. The code uses message passing parallelism with a form of domain decomposition, where the global grid is partitioned into subgrids, with each subgrid assigned to a processor. Data on each subgrid is evolved independently during the source and transport phases of a timestep and data on subgrid boundaries is exchanged at the end of each phase through message-passing calls. This results in a highly scalable code that exhibits good speedup over the full range of practically realizable subgrids.

The present simulations have been designed as a follow-on to the results in H91, but unlike the code used there, we do not need to impose equatorial symmetry. The calculations performed on a 64 x 32 x 64 $(r,\theta,\phi)$-grid in H91 correspond to a 64 x 64 x 64 grid here. The ergosphere can be part of the computational domain, but at the event horizon, $r_H = M + \sqrt{M^2-a^2}$, some of the Boyer-Lindquist metric terms are singular, so the inner edge of the radial grid must lie at some point outside the horizon, continuing out to some selected outer grid boundary. Although the horizon itself can never be part of the computational grid, tests indicate that it is possible to come arbitrarily close to the horizon (e.g., r_in = 2.001 M for the Schwarzschild case). But this comes at a steep price: the narrow region where the metric terms diverge rapidly must be covered by a large number of grid zones (in the extreme Kerr limit, the number of zones inside the ergosphere could be as high as the number of zones that lie outside it). In addition to consuming large amounts of memory in a three dimensional simulation, these innermost grid zones can significantly reduce the time step size, greatly increasing the total computational time. The specific location of the inner edge of the grid is determined by a balance between the physical phenomena of interest in a simulation and the memory and performance demands. For the hydrodynamic torus simulations considered here, the region near the horizon is of secondary interest, so the inner grid boundary, r_min, was set at or just inside the static limit on the equator, $r_{min} \ler_{static}{\vert}_{\theta=\pi/2}=2\,M$. In the runs reported here, this ranged from r_in = 1.90 M for the smaller grids to r_in = 2.05 M for the larger grids. The polar grid extends almost to the rotational axes at $\theta=0$ and $\pi$, but the axes themselves are not included. As in non-relativistic codes, the polar axes are problematic due to the shrinking of zone volumes, and the collapse of the inner radial zone face at r = 0. These problems are further compounded by the fact that, here again, some metric components are singular. Torus stability studies do not require the axes, and their omission is not important to the present simulations. Regularized operators will be introduced in future work. The radial and polar grids use logarithmic scaling for all runs in order to concentrate zones near the horizon and the equator. The azimuthal angle $\phi$ spans the full range from 0 to $2\pi$ with equally-spaced grid zones and periodic boundary conditions. Previous studies of the Papaloizou-Pringle instability show that using the full $\phi$ range is important as the most important unstable modes have azimuthal wavenumber m=1.

The time step $\Delta t$ is determined by the extremal light-crossing time for a grid zone, as described in HSWb. This time step size remains fixed for the duration of the simulation since it is a purely geometric result.

The initial state for a simulation is generated from the thick disk equilibrium solution (described in §3) with given input parameters. This initial configuration is then perturbed with small, random enthalpy fluctuations with a maximum amplitude of 1%. These perturbations are the seeds from which the full PPI will develop.

2002-06-05