3.3 The Bondi problem

The first characteristic code based upon the original Bondi equations for a twist-free axisymmetric spacetime was constructed by J. Welling in his PhD thesis at Pittsburgh [153Jump To The Next Citation Point] . The spacetime was foliated by a family of null cones, complete with point vertices at which regularity conditions were imposed. The code accurately integrated the hypersurface and evolution equations out to compactified null infinity. This allowed studies of the Bondi mass and radiation flux on the initial null cone, but it could not be used as a practical evolution code because of instabilities.

These instabilities came as a rude shock and led to a retreat to the simpler problem of axisymmetric scalar waves propagating in Minkowski space, with the metric

2 2 2( 2 2 2) ds = − du − 2du dr + r dπœƒ + sin πœƒ dΟ• (9 )
in outgoing null cone coordinates. A null cone code for this problem was constructed using an algorithm based upon Equation (8View Equation), with the angular part of the flat space Laplacian replacing the curvature terms in the integrand on the right hand side. This simple setting allowed one source of instability to be traced to a subtle violation of the CFL condition near the vertices of the cones. In terms of the grid spacing Δx α, the CFL condition in this coordinate system takes the explicit form
Δu-- [ 2 2 ]1βˆ•2 Δr < − 1 + K + (K − 1) − 2K (K − 1)cos Δπœƒ , (10 )
where the coefficient K, which is of order 1, depends on the particular startup procedure adopted for the outward integration. Far from the vertex, the condition (10View Equation) on the time step Δu is quantitatively similar to the CFL condition for a standard Cauchy evolution algorithm in spherical coordinates. But condition (10View Equation) is strongest near the vertex of the cone where (at the equator πœƒ = πβˆ•2) it implies that
Δu < K Δr (Δ πœƒ)2. (11 )
This is in contrast to the analogous requirement
Δu < K Δr Δ πœƒ (12 )
for stable Cauchy evolution near the origin of a spherical coordinate system. The extra power of Δ πœƒ is the price that must be paid near the vertex for the simplicity of a characteristic code. Nevertheless, the enforcement of this condition allowed efficient global simulation of axisymmetric scalar waves. Global studies of backscattering, radiative tail decay, and solitons were carried out for nonlinear axisymmetric waves [153Jump To The Next Citation Point], but three-dimensional simulations extending to the vertices of the cones were impractical at the time on existing machines.

Aware now of the subtleties of the CFL condition near the vertices, the Pittsburgh group returned to the Bondi problem, i.e. to evolve the Bondi metric [54Jump To The Next Citation Point]

( V ) ( ) ds2 = --e2β − U 2r2e2γ du2 + 2e2βdu dr + 2U r2e2γdud πœƒ − r2 e2γdπœƒ2 + e−2γ sin2πœƒ dΟ•2 , (13 ) r
by means of the three hypersurface equations
β = 1r(γ )2, (14 ) ,r 2 [,r ] [ ] ( ) 2 r4e2(γ−β)U,r = 2r2 r2 β- − (sin-πœƒ-γ),rπœƒ + 2γ,rγ,πœƒ , (15 ) ,r r2 ,rπœƒ sin2πœƒ 4 V,r = − 1r4e2(γ− β)(U,r)2 + (r-sinπœƒ-U-),rπœƒ 4 [ 2r2 sin πœƒ ] 2(β−γ) (sin-πœƒβ,πœƒ),πœƒ- 2 +e 1 − sin πœƒ + γ,πœƒπœƒ + 3cot πœƒγ,πœƒ − (β,πœƒ) − 2γ,πœƒ(γ,πœƒ − β,πœƒ) , (16 )
and the evolution equation
{ [ ]} ( U ) (γ U sin πœƒ) 4r(rγ),ur = 2r γ,rV − r2 2γ,πœƒU + sin πœƒ ----- − 2r2--,r-------,πœƒ sinπœƒ ,πœƒ ,r sin πœƒ [ ( ) ] 1-4 2(γ− β) 2 2(β−γ) 2 -β,πœƒ- + 2r e (U,r) + 2e (β,πœƒ) + sin πœƒ sinπœƒ . (17 ) ,πœƒ
The beauty of the Bondi equations is that they form a clean hierarchy. Given γ on an initial null hypersurface, the equations can be integrated radially to determine β, U, V, and γ ,u on the hypersurface (in that order) in terms of integration constants determined by boundary conditions, or smoothness conditions if extended to the vertex of a null cone. The initial data γ is unconstrained except for smoothness conditions. Because γ represents an axisymmetric spin-2 field, it must be π’ͺ (sin2πœƒ) near the poles of the spherical coordinates and must consist of l ≥ 2 spin-2 multipoles.

In the computational implementation of this system by the Pittsburgh group [120Jump To The Next Citation Point], the null hypersurfaces were chosen to be complete null cones with nonsingular vertices, which (for simplicity) trace out a geodesic worldline r = 0. The smoothness conditions at the vertices were formulated in local Minkowski coordinates.

The vertices of the cones were not the chief source of difficulty. A null parallelogram marching algorithm, similar to that used in the scalar case, gave rise to another instability that sprang up throughout the grid. In order to reveal the source of this instability, physical considerations suggested looking at the linearized version of the Bondi equations, where they can be related to the wave equation. If this relationship were sufficiently simple, then the scalar wave algorithm could be used as a guide in stabilizing the evolution of γ. A scheme for relating γ to solutions Φ of the wave equation had been formulated in the original paper by Bondi, Metzner, and van der Burgh [54]. However, in that scheme, the relationship of the scalar wave to γ was nonlocal in the angular directions and was not useful for the stability analysis.

A local relationship between γ and solutions of the wave equation was found [120Jump To The Next Citation Point]. This provided a test bed for the null evolution algorithm similar to the Cauchy test bed provided by Teukolsky waves [239]. More critically, it allowed a simple von Neumann linear stability analysis of the finite difference equations, which revealed that the evolution would be unstable if the metric quantity U was evaluated on the grid. For a stable algorithm, the grid points for U must be staggered between the grid points for γ, β, and V. This unexpected feature emphasizes the value of linear stability analysis in formulating stable finite difference approximations.

It led to an axisymmetric code [186Jump To The Next Citation Point120Jump To The Next Citation Point] for the global Bondi problem which ran stably, subject to a CFL condition, throughout the regime in which caustics and horizons did not form. Stability in this regime was verified experimentally by running arbitrary initial data until it radiated away to ℐ+. Also, new exact solutions as well as the linearized null solutions were used to perform extensive convergence tests that established second order accuracy. The code generated a large complement of highly accurate numerical solutions for the class of asymptotically flat, axisymmetric vacuum spacetimes, a class for which no analytic solutions are known. All results of numerical evolutions in this regime were consistent with the theorem of Christodoulou and Klainerman [75] that weak initial data evolve asymptotically to Minkowski space at late time.

An additional global check on accuracy was performed using Bondi’s formula relating mass loss to the time integral of the square of the news function. The Bondi mass loss formula is not one of the equations used in the evolution algorithm but follows from those equations as a consequence of a global integration of the Bianchi identities. Thus it not only furnishes a valuable tool for physical interpretation but it also provides a very important calibration of numerical accuracy and consistency.

An interesting feature of the evolution arises in regard to compactification. By construction, the u-direction is timelike at the origin where it coincides with the worldline traced out by the vertices of the outgoing null cones. But even for weak fields, the u-direction generically becomes spacelike at large distances along an outgoing ray. Geometrically, this reflects the property that ℐ+ is itself a null hypersurface so that all internal directions are spacelike, except for the null generator. For a flat space time, the u-direction picked out at the origin leads to a null evolution direction at + ℐ, but this direction becomes spacelike under a slight deviation from spherical symmetry. Thus the evolution generically becomes “superluminal” near ℐ+. Remarkably, this leads to no adverse numerical effects. This remarkable property apparently arises from the natural way that causality is built into the marching algorithm so that no additional resort to numerical techniques, such as “causal differencing” [78], is necessary.

3.3.1 The conformal-null tetrad approach

Stewart has implemented a characteristic evolution code which handles the Bondi problem by a null tetrad, as opposed to metric, formalism [231Jump To The Next Citation Point]. The geometrical algorithm underlying the evolution scheme, as outlined in [23399], is Friedrich’s [94Jump To The Next Citation Point] conformal-null description of a compactified spacetime in terms of a first order system of partial differential equations. The variables include the metric, the connection, and the curvature, as in a Newman–Penrose formalism, but in addition the conformal factor (necessary for compactification of ℐ) and its gradient. Without assuming any symmetry, there are more than 7 times as many variables as in a metric based null scheme, and the corresponding equations do not decompose into as clean a hierarchy. This disadvantage, compared to the metric approach, is balanced by several advantages:

The code was intended to study gravitational waves from an axisymmetric star. Since only the vacuum equations are evolved, the outgoing radiation from the star is represented by data (Ψ 4 in Newman–Penrose notation) on an ingoing null cone forming the inner boundary of the evolved domain. This inner boundary data is supplemented by Schwarzschild data on the initial outgoing null cone, which models an initially quiescent state of the star. This provides the necessary data for a double-null initial value problem. The evolution would normally break down where the ingoing null hypersurface develops caustics. But by choosing a scenario in which a black hole is formed, it is possible to evolve the entire region exterior to the horizon. An obvious test bed is the Schwarzschild spacetime for which a numerically satisfactory evolution was achieved (although convergence tests were not reported).

Physically interesting results were obtained by choosing data corresponding to an outgoing quadrupole pulse of radiation. By increasing the initial amplitude of the data Ψ4, it was possible to evolve into a regime where the energy loss due to radiation was large enough to drive the total Bondi mass negative. Although such data is too grossly exaggerated to be consistent with an astrophysically realistic source, the formation of a negative mass was an impressive test of the robustness of the code.

3.3.2 Axisymmetric mode coupling

Papadopoulos [183Jump To The Next Citation Point] has carried out an illuminating study of mode mixing by computing the evolution of a pulse emanating outward from an initially Schwarzschild white hole of mass M. The evolution proceeds along a family of ingoing null hypersurfaces with outer boundary at r = 60 M. The evolution is stopped before the pulse hits the outer boundary in order to avoid spurious effects from reflection and the radiation is inferred from data at r = 20 M. Although gauge ambiguities arise in reading off the waveform at a finite radius, the work reveals interesting nonlinear effects: (i) modification of the light cone structure governing the principal part of the equations and hence the propagation of signals; (ii) modulation of the Schwarzschild potential by the introduction of an angular dependent “mass aspect”; and (iii) quadratic and higher order terms in the evolution equations which couple the spherical harmonic modes. A compactified version of this study [256Jump To The Next Citation Point] was later carried out with the 3D PITT code, which confirms these effects as well as new effects which are not present in the axisymmetric case (see Section 4.5 for details).

3.3.3 Spectral approach to the Bondi problem

Oliveira and Rodrigues [81] have presented the first code based upon the Galerkin spectral method to evolve the axisymmetic Bondi problem. The strength of spectral methods is to provide high accuracy relative to computational effort. The spectral decomposition reduces the partial differential evolution system to a system of ordinary differential equations for the spectral coefficients. Several numerical tests were performed to verify stability and convergence, including linearized gravitational waves and the global energy momentum conservation law relating the Bondi mass to the radiated energy flux. The main feature of the Galerkin method is that each basis function is chosen to automatically satisfy the boundary conditions, in this case the regularity conditions on the Bondi variables on the axes of symmetry and at the vertices of the outgoing null cones and the asymptotic flatness condition at infinity. Although + ℐ is not explicitly compactified, the choice of radial basis functions allows verification of the asymptotic relations governing the coefficients of the leading gauge dependent terms of the metric quantities. The results promise a new approach to study spacetimes with gravitational waves. It will be interesting to see whether the approach can be generalized to the full 3D case.

3.3.4 Twisting axisymmetry

The Southampton group, as part of its goal of combining Cauchy and characteristic evolution, has developed a code [85Jump To The Next Citation Point86Jump To The Next Citation Point194Jump To The Next Citation Point] which extends the Bondi problem to full axisymmetry, as described by the general characteristic formalism of Sachs [211]. By dropping the requirement that the rotational Killing vector be twist-free, they were able to include rotational effects, including radiation in the “cross” polarization mode (only the “plus” mode is allowed by twist-free axisymmetry). The null equations and variables were recast into a suitably regularized form to allow compactification of null infinity. Regularization at the vertices or caustics of the null hypersurfaces was not necessary, since they anticipated matching to an interior Cauchy evolution across a finite worldtube.

The code was designed to insure standard Bondi coordinate conditions at infinity, so that the metric has the asymptotically Minkowskian form corresponding to null-spherical coordinates. In order to achieve this, the hypersurface equation for the Bondi metric variable β must be integrated radially inward from infinity, where the integration constant is specified. The evolution of the dynamical variables proceeds radially outward as dictated by causality [194Jump To The Next Citation Point]. This differs from the Pittsburgh code in which all the equations are integrated radially outward, so that the coordinate conditions are determined at the inner boundary and the metric is asymptotically flat but not asymptotically Minkowskian. The Southampton scheme simplifies the formulae for the Bondi news function and mass in terms of the metric. It is anticipated that the inward integration of β causes no numerical problems because this is a gauge choice which does not propagate physical information. However, the code has not yet been subject to convergence and long term stability tests so that these issues cannot be properly assessed at the present time.

The matching of the Southampton axisymmetric code to a Cauchy interior is discussed in Section 5.6.

  Go to previous page Go up Go to next page