II WORKSHOP :NOVA FÍSICA NO ESPAÇO
Samuel Rocha de Oliveira
SOME NUMERICAL
RELATIVITY METHODS FOR GRAZING BLACK HOLES USING HIERARCHICAL 3+1 DECOMPOSITION
OF EINSTEIN'S EQUATIONS
Samuel
Rocha de Oliveira*
Depto.
de Matemática Aplicada, UNICAMP
Abstract: We use a specific 3+1 decomposition of Einstein
equations for numerical implementation. The decomposition keeps the volume
element under control, allows analytical extraction of some singularities and
the imposition of the boundary conditions can be consistent with the
constraints. The coordinates does not avoid singularities though. So we have to
extract them either analytically or numerically with the assumption that an
event horizon hide both the singularities and the procedure to deal with them.
The resulting equations are numerically solved in the following way: I) The
initial data is constructed solving the nonlinear constraint equations using
analytical exact solutions superpositions as initial guess and relaxation of
the elliptical equations using finite difference algebraic equations. The
constraints are not solved further more. They are used as measure of accuracy
during the evolution. The outer boundary conditions for the initial data are
exact solutions of supposed configurations. II) For the lapse function we chose
a few variations of "1+log" slice. III) For the shift vectors we use
minimal distortion conditions on the factored metric. This expensive choice
results in 3 more elliptical equations which have to be solved in every time
step but it keeps the codes running longer. Steps II and III fix a set of
observers along the evolution. IV) With respect to those observers the
evolution equations are solved using second order (quasi) implicit finite
difference/volume methods or adaptive method of lines. We use Richardson's
extrapolations methods to show the second order of our code and to improve the
accuracy. At the outer boundaries we have a combination of standard Newman or
radiation conditions and up to four of the Einstein equations as constraints on
the conditions -- for simplicity we either use an analytical exact solution or
"tangential" linearized Einstein equations. For inner boundary
conditions (at the holes) we damp and freeze the functions close to the
"singularities". We apply these methods for the grazing black holes
problem. Let m the the total ADM like
mass of the isolated system. Our code, nowadays, is running up to about 14m with the constraints L2
violations below 1% for a grid size of the order of m/8 in each spatial direction and a time step equivalent to m/32. After that time spurious mode
prevails. An estimate of the amount of gravitational wave energy emitted
between the time 4m and 10m is below 10-4m. The horizons are not localized in our
simulations yet. The dimensions of our grid for a typical run is 20 ´ 20 ´ 10 ´ m3.
1. Introduction
Let the space-time metric be
where the factored spatial metric is
We further decompose the bi-dimensional metric sIJ as follows:
Thus, we have the following 13 metric functions: N, A,
B, C, W, Y, S, bi, cJ, d.
This decomposition allows, the control of
where the usual coordinates singularities appears.
Thus, one of the roles of the extra functions is to control the volume element
of the space and time. Above and below, we use the indices as follows: Greek
indices, like m ranges the spacetime dimensions (0..3), lower case latin
indices, like i, ranges the space
dimensions (1..3) and the upper case latin indices like I ranges the surface indices (2, 3).
Note also that
Thus, none of the "covectors" bi, cJ, d affect
the volume element. Let us remark that these covectors are defined in the
natural cotangent space of the manifolds (M3,
g), (M2, s) and (M1, 1) respectively.
Therefore there are some topological restrictions on our hierarchical
decomposition.
The "scale" functions N, A, B, C of the "principal" directions have to be
positive definite and bounded. The other functions may vanish but have to be
bounded as well. These choices imply the 4-vectors along each direction 1,2,3
be space-like and it cannot change the behavior to become null or time-like.
Similarly the 4-vector along the direction 0 is always time-like and infinite red-shifts
are avoided. Of course, some of the observers will hit the curvature
singularities in a finite time and have to be discarded. The coordinates and
naked singularities are eliminated from either the computational grid or from
the equations themselves, as much as possible, with the aid of the extra
functions.
Let us call dynamical functions those which the
Einstein equations give second order derivative in time. There are only six of
them.
So the initial boundary value problem for General
Relativity requires the solution for the constraints equations Gtm = kTtm for the
dynamic functions and their time derivatives. Then the evolution equations Gij = k
. At the outer boundary we have to solve the constraints Gnormal = k
for some of the functions and their normal derivatives.
Our general boundary value problem can be cast as the
following set of second order quasi-linear coupled equations (gBV2oqlpdes):
for a list of functions Ue[x] and the
boundary equations:
where (¶nUe)
is the normal derivative of Ue[x]. We use the standard summation on
repeated indexes for the appropriate ranges. The non singular arrays A and V are given functions of the position x; F and B are, in general, non-linear functions
of their arguments.
The algebraic equations come from the differential
equations through a second order accuracy estimate of the integral in a
adjustable finite volume around every grid point of X and ¶X. The
non-linear algebraic equations are then transformed to
linear algebraic equations to be iterated starting from known approximate exact
solutions. For each iteration the linear algebraic equations are solved by
Gauss-Seidel method with 2 term Conjugate Gradient acceleration for the
convergence
(NSPCG and ITPACK packages). For very simple
configurations (up to two functions in 2-d) we use public available elliptic
solvers because they are faster than our present gBV2oqlpdes implementation.
The initial value boundary problem can be written as
the set of second order quasi-linear coupled equations (gIBV2oqlpdes):
, for a list of functions Ue[x, t] having the initial conditions:
and the boundary equations:
Similarly to the previous setting, the non singular
arrays A and V are given functions of x
and t; F and B are, in general,
non-linear functions of their arguments and the index b allows both time and
space differentiation.
At the moment our numerical scheme for the initial
value boundary problem is a fully implicit, second order in time and space,
finite difference method (Crank Nicholson like schemes). For non stiff problems
we allow semi-implicit schemes with two iterations for each time step,
otherwise the tri-diagonal linear problem is combined with a Newton
linearization of the non-linear terms with a Gauss-Seidel solver. The parameter
space is considerably large and no attempts of automatic choices were made.
2. Grazing Black
Holes
The iteration process of the boundary value problem
starts with the exact solution for two black holes (with conical singularities)
in Weyl coordinates. See figures. So we prescribe, as free initial data, all
the metric functions but cJ
and ¶tcJ. These four functions are the unknowns of the
constraint equations.
Figure 1: Front, top and side view of the computational grid
3. Computational
Details
Most of our code contains Fortran 77 standard subroutines,
specially for linear algebra and ordinary differential equations. Several
pieces of the code now are in Fortran 95 (Lahey compiler) and in C++ (gcc
compiler). For several tasks we use Maple 7: (symbolic manipulation, code
generation, grtensor, parameters manipulation etc.). No parallel implementation
yet. Our codes are directly linked to simple graphics software (gnu-plot and
octave) but we intend to use the recent Open DX (based on IBM's Visualization
Data Explorer) which can be linked to Fortran/Linux/X11 software.
Figure 2: Color map of the Hamiltonian constraint violation
References
[1] T. W. Baumgarte and S. L. Shapiro,
"Numerical Relativity and Compact Binaries", gr-qc/02110828 (2002).
[2] L. Lidblom and M. A. Scheel, "Energy
Norms and the Stability of the Einstein Evolution Equations", gr-qc/0206035 (2002).
[3] P. Anninos et. al., "Collision of two
Black Holes", Phys. Rev. Lett. 71, 2851 (1993).
[4] P. Letelier and S. R. Oliveira,
"Superposition of Weyl solutions: the equilibrium forces", Class. Quantum Grav., 15, 412 (1998).
[5] G. Cook e. al., "Three-dimensional
initial data for the collision of two black holes". Phys. Rev. D 47
1471-1490 (1993).