The formulation of Einstein's equations in the NQS gauge is
presented in Ref. [15] and the associated gauge freedom arising from (*u*,
*r*) dependent rotation and boosts of the unit sphere is discussed
in Ref. [16]. As in the PITT code, the main equations involve integrating a
hierarchy of hypersurface equations along the radial null
geodesics extending from the inner boundary to null infinity. In
the NQS gauge the source terms for these radial ODE's are rather
simple when the unknowns are chosen to be the connection
coefficients. However, as a price to pay for this simplicity,
after the radial integrations are performed on each null
hypersurface a first order elliptic equation must be solved on
each
cross-section to reconstruct the underlying metric.

Both codes require the ability to handle tensor fields and their derivatives on the sphere. Spherical coordinates and spherical harmonics are natural analytic tools for the description of radiation but their implementation in computational work requires dealing with the impossibility of smoothly covering the sphere with a single coordinate grid. Polar coordinate singularities in axisymmetric systems can be regularized by standard tricks. In the absence of symmetry, these techniques do not generalize and would be especially prohibitive to develop for tensor fields.

A crucial ingredient of the PITT code is the
*eth*
-module [72] which incorporates a computational version of the
Newman-Penrose eth-formalism [107]. The eth-module covers the sphere with two overlapping
stereographic coordinate grids (North and South). It provides
everywhere regular, second order accurate, finite difference
expressions for tensor fields on the sphere and their covariant
derivatives. The eth-calculus simplifies the underlying
equations, avoids spurious coordinate singularities and allows
accurate differentiation of tensor fields on the sphere in a
computationally efficient and clean way. Its main weakness is the
numerical noise introduced by interpolations (4th order accurate)
between the North and South patches. For parabolic or elliptic
equations on the sphere, the finite difference approach of the
eth-calculus would be less efficient than a spectral approach,
but no parabolic or elliptic equations appear in the Bondi-Sachs
evolution scheme.

The Canberra code handles fields on the sphere by means of a 3-fold representation [20]: (i) discretized functions on a spherical grid uniformly spaced in standard coordinates, (ii) fast-Fourier transforms with respect to (based upon a smooth map of the torus onto the sphere), and (iii) spectral decomposition of scalar, vector and tensor fields in terms of spin-weighted spherical harmonics. The grid values are used in carrying out nonlinear algebraic operations; the Fourier representation is used to calculate -derivatives; and the spherical harmonic representation is used to solve global problems, such as the solution of the first order elliptic equation for the reconstruction of the metric, whose unique solution requires pinning down the gauge freedom. The sizes of the grid and of the Fourier and spherical harmonic representations are coordinated. In practice, the spherical harmonic expansion is carried out to 15th order in but the resulting coefficients must then be projected into the subspace in order to avoid inconsistencies between the spherical harmonic and the grid and Fourier representations.

The Canberra code solves the null hypersurface equations by
combining an 8th order Runge-Kutta integration with a convolution
spline to interpolate field values. The radial grid points are
dynamically positioned to approximate ingoing null geodesics, a
technique originally due to Goldwirth and Piran [67] to avoid inaccuracies near a horizon resulting from a uniform
*r*
-grid due to the coordinate singularity in the case of a
stationary horizon. The time evolution uses the method of lines
with a 4th order Runge-Kutta integrator, which introduces further
high frequency filtering.

**Canberra code**
: Direct stability analysis of the underlying finite difference
equations is impractical because of the extensive mix of spectral
techniques, higher order methods and splines. Similarly, there is
no clear-cut CFL limit on the code, although stability tests show
that there is a limit on the time step. The damping of high
frequency modes due to the implicit filtering would be expected
to suppress numerical instability, but the stability of the
Canberra code is nevertheless subject to two
qualifications [18,
19,
21]: (i) At late times (less than 100
*M*), the evolution terminates as it approaches an event horizon,
apparently because of a breakdown of the NQS gauge condition,
although an analysis of how and why this should occur has not yet
been given. (ii) Numerical instabilities arise from dynamic inner
boundary conditions and restrict the inner boundary to a fixed
Schwarzschild horizon. Tests in the extreme nonlinear regime were
not reported.

- Linearized waves on a Minkowski background in null cone coordinates.
- Boost and rotation symmetric solutions [23].
- Schwarzschild in rotating coordinates.
- Polarization symmetry of nonlinear twist-free axisymmetric waveforms.
- Robinson-Trautman waveforms from perturbed Schwarzschild black holes.
- Nonlinear Robinson-Trautman waveforms utilizing an independently computed solution of the Robinson-Trautman equation.

**Canberra code**
: The complexity of the algorithm and NQS gauge makes it
problematic to establish accuracy by direct means. Exact
solutions do not provide an effective convergence check because
the Schwarzschild solution is trivial in the NQS gauge and other
known solutions in this gauge require dynamic inner boundary
conditions which destabilize the present version of the code.
Convergence to linearized solutions is a possible check but has
not yet been performed. Instead indirect tests by means of
geometric consistency and partial convergence tests are used to
calibrate accuracy. The consistency tests are based on the
constraint equations, which are not enforced during null
evolution except at the inner boundary. The balance between mass
loss and radiation flux through
is a global version of these constraints. No appreciable growth
of the constraints is noticeable until within 5
*M*
of the final breakdown of the code. In weak field tests where
angular resolution does not dominate the error, partial
convergence tests based upon varying the radial grid size show
the 8th order convergence in the shear expected from the
Runge-Kutta integration and splines. When the radial source of
error is small, reduced error with smaller time step can also be
discerned.

In practical runs, the major source of inaccuracy is the spherical harmonic resolution, currently restricted to by hardware limitations. Truncation of the spherical harmonic expansion has the effect of modifying the equations to a system for which the constraints are no longer satisfied. The relative error in the constraints is for strong field simulations [17].

The Pittsburgh group studied the backscattered waveform (described by the Bondi news function) as a function of incoming pulse amplitude. The computational eth module smoothly handled the complicated time dependent transformation between the non-inertial computational frame at and the inertial (Bondi) frame necessary to obtain the standard ``plus'' and ``cross'' polarization modes. In the perturbative regime, the news corresponds to the backscattering of the incoming pulse off the effective Schwarzschild potential. When the energy of the pulse is no larger than the central Schwarzschild mass, the backscattered waveform still depends roughly linearly on the amplitude of the incoming pulse. However, for very high amplitudes the waveform behaves quite differently. Its amplitude is greater than that predicted by linear scaling and its shape drastically changes and exhibits extra oscillations. In this very high amplitude case, the mass of the system is completely dominated by the incoming pulse, which essentially backscatters off itself in a nonlinear way.

The Canberra code was used to study the change in Bondi mass
due to radiation [17]. The Hawking mass
is first calculated as a function of radius and retarded time,
with the Bondi mass
then obtained by taking the limit
. The limit has good numerical behavior. For a strong initial
pulse with
angular dependence, in a run from
*u*
=0 to
*u*
=70 (in units where the interior Schwarzschild mass is 1), the
Bondi mass drops from 1.8 to 1.00002, showing that almost half of
the initial energy of the system was backscattered and that a
surprisingly negligible amount falls into the black hole. The
Bondi mass decreases monotonically in time, as theoretically
necessary, but its rate of change exhibits an interesting pulsing
behavior whose time scale cannot be obviously explained in terms
of quasinormal oscillations. The Bondi mass loss formula is
confirmed with relative error of less than
. This is impressive accuracy considering the potential sources
of numerical error introduced by taking the limit of the Hawking
mass. The Canberra group also studied the effect of initial pulse
amplitude on the waveform of the backscattered radiation but did
not extend their study to the very high amplitude regime in which
qualitatively interesting nonlinear effects occur.

The code essentially runs ``forever'' even when the worldtube
wobbles with respect to the black hole to produce artificial
periodic time dependence. An initially distorted, wobbling black
hole was evolved for a time of 60,000
*M*, longer by orders of magnitude than stability permits for any
other existing 3D black hole code. This exceptional performance
opens a promising new approach to handle the inner boundary
condition for Cauchy evolution of black holes by the matching
methods reviewed in Sec.
4
.

Characteristic Evolution and Matching
Jeffrey Winicour
http://www.livingreviews.org/lrr-2001-3
© Max-Planck-Gesellschaft. ISSN 1433-8351 Problems/Comments to livrev@aei-potsdam.mpg.de |