The formulation of Einstein's equations in the NQS gauge is presented in Ref.  and the associated gauge freedom arising from (u, r) dependent rotation and boosts of the unit sphere is discussed in Ref. . 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  which incorporates a computational version of the Newman-Penrose eth-formalism . 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 : (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  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.
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 .
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 . 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
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to firstname.lastname@example.org