[List and description as of September 2001; written by the Scientific Organizing Committee and extracted from the proposal to NSF for partial funding of the Visitors Program.]
For the case of data representing binaries in nearly circular orbits, thin-sandwich methods [5] have been used for binary neutron star data [6][7] and in initial trials with binary black hole configurations [8,9]. It seems clear that these methods, as applied to black hole data, can be extended to allow for general, astrophysically realistic configurations [10]. An important question is: how relevant are waveforms from binaries in elliptic orbits? If these configurations are important, then new methods must be developed.
Astrophysically realistic initial data must be based on well motivated choices for the freely-specifiable portions of the data. We expect that post-Newtonian methods will be critical in this area. While first steps have been taken in providing the free data, e.g. [11][12], there are still significant problems to overcome in specifying the boundary and gauge conditions. It seems likely that non-maximal slicings will be important, particularly for the black hole case where horizon penetrating gauges are needed. Little is known about the well-posedness of the constraint equations for non-vanishing mean extrinsic curvature.
Finally, efficient numerical methods are important for solving the coupled sets of elliptic equations. Parameter space surveys are often needed to locate appropriate solutions and to fully understand the content of the initial data [3][13]. These are particularly important for locating the innermost stable circular orbit. Finite-difference [14,15], pseudo-spectral [16][17], and finite-element [18] methods have all been used in the construction of initial data and there is no consensus as to which methods are most appropriate and efficient.
Radiation reaction effects are likely to circularize the binary orbit, and cause it to slowly shrink while approaching the ISCO. This suggests a quasi-equilibrium approach, in which the orbit is approximated as quasi-circular, and the orbital decay is treated as a small correction. The inspiral can then be modeled as a sequence of quasi-circular binary models at varying separation. Evolving the gravitational fields in the presence of these predetermined background models yields the gravitational wave signal and luminosity, and hence the rate at which the system loses energy. The luminosity determines the inspiral rate, and from that the separation as a function of time, and accordingly the entire gravitational wave train.
Various approaches can be taken to solve various aspects of the problem. Duez, Baumgarte and Shapiro [20] and Duez et. al. [21] inserted models of corotating and irrotational binary neutron stars into a fully relativistic evolution code to construct the quasi-adiabatic gravitational wavetrain (see also Yo, Baumgarte and Shapiro [22]). Shibata and Uryu [23] have adopted a perturbational approach to simulate the gravitational wave emission from binary neutron stars. Other related approaches have been suggested by Brady, Creighton and Thorne [19], Laguna [24], Whelan and Romano [25] and Whelan, Krivan and Price [26].
So far, all applications of such quasi-adiabatic evolution schemes have been for binary neutron stars. The primary goal of this workshop is therefore to develop methods that are suitable for binary black holes and to solve the "intermediate binary black hole problem". We also hope that simulations of the late binary neutron star inspiral can be improved in various ways, including improvements of the background models and the wave extraction. Related issues are also being considered in Topics #1 and #6 of this workshop.
The goal of this topic is to develop new formulations of Einstein's equations and to evaluate their suitability for numerical implementation. There are two main issues to consider: analytic properties of the equations and numerical properties of particular implementations. These issues are closely related, and ideally, a suitable analytic formulation should simplify the construction of an acceptable numerical implementation. But this is far from being guaranteed, and additional mathematical insight is required to obtain the desired numerical properties. To provide this insight is the most important aspect of this topic.
Key questions are: What variables are appropriate? Should one perform a trace or a conformal decomposition? Should one use a first or second order form of the equations? Is it essential for the set of evolution equations to be strongly hyperbolic? Can one identify unstable modes analytically? Can one identify instabilities caused by particular numerical implementations? Do finite element methods offer significant advantages? Can one obtain better satisfaction of the constraints, by adding driver terms to the equations or by the choice of a particular formulation? Can geometric (symplectic or time-reversible) integrators used in constrained Hamiltonian systems be designed to enforce the constraints? What are the advantages of formulations different from the 3+1 approach (characteristic formulations, conformal field equations)? What gauge choices (see also #4) and boundary conditions (see also #5) are appropriate for different formulations?
Gauge choices involve relating the coordinates to geometric conditions on the dynamic variables, and fall into several categories: 1) Algebraic gauges, where the gauge variables are obtained from the dynamical variables by algebraic prescriptions; examples are the harmonic [49,29] and "1+log" [50] slicings. 2) Conditions that involve solving evolution equations for the gauge variables, such as the Bona-Masso [30] and the "driver" conditions [51]. 3) Choices that involve the solution of elliptic equations, such as maximal slicing [52,53] and the minimal distortion shift [53]. Algebraic and evolution-type conditions are easy to implement numerically and fast to solve, but are prone to developing coordinate singularities. Elliptic conditions are more robust but computationally costly, and they also require appropriate boundary conditions [54]. In choosing a good gauge several considerations are involved. How does the gauge influence the stability and accuracy of numerical simulations? How are coordinates well adapted to a specific problem chosen? How are gauges found that make the system time-independent at late times? How are coordinate singularities avoided?
A.- Outer Boundaries: If the domain of an evolution code does not extend to infinity, one may impose an approximate boundary condition (e.g., an outgoing-wave condition) at an artificial boundary, or match to an exterior evolution computed using a perturbative [55] or a characteristic formulation [56,57]. The last option allows evolution out to null infinity, so it does not introduce errors caused by artificial boundaries. However, Cauchy-characteristic matching has yet to be successful in 3D simulations; explicit use of the constraints at the matching boundary [58] may prove helpful.
B.- Inner boundary: One may introduce boundaries inside black hole horizons to excise their singular interiors and evolve only the non-singular region [59,60,61,62]. For hyperbolic equations, causality ensures that the excised interior cannot affect the exterior, so (in principle) no boundary condition need be imposed at the excision boundary. Numerical schemes have been constructed in which no excision boundary condition is required [48,44]; however, for numerical reasons one may wish to impose an artificial boundary condition inside the horizon [63,43]. Special attention is required if the excised regions are chosen to move through the numerical grid [64]. Furthermore, for elliptic equations (e.g., constrained evolution), assigning boundary values is particularly delicate [43] as they can influence physics outside the horizon.
C.- Physical and grid boundaries: Special treatment might also be needed at physical boundaries like the surface of stars, and boundaries between distinct numerical grids that represent multiple coordinate patches or are used for adaptive mesh refinement techniques.
A.- Standard methods: The basic idea is to assume that, at finite distance from the source on a ``far-away'' time-like surface, the spacetime is well approximated by perturbations of a Kerr spacetime since away from the strong field region one expects only low amplitude gravitational waves propagating on a black hole background [55]. (The inherent error in this method is diminished as the outer boundary is moved further away, assuming spurious reflections on the boundaries do not significantly contaminate the evolution.) These methods are based on the gauge invariant formalisms of Moncrief [65] and Teukolsky [66].
B.- The Lazarus approach: Another natural `extraction' boundary occurs on a `space-like' interface defined by the time beyond which non-linearities no longer contribute significantly to the evolution. The `Lazarus project' [67,68] takes the following approach to providing a full numerical-close limit interface: once the system has evolved towards a perturbative regime (e.g., two coalescing black holes form a distorted one), one extracts the relevant gravitational wave data, and evolves them on the appropriate black hole background. This approach has produced the first merger waveforms for the coalescence of two black holes starting from an estimate of the innermost stable circular orbit. Further work is needed to improve the accuracy of these waveforms[68].
C.- Cauchy-characteristic matching: This approach combines a Cauchy formulation with a characteristic implementation of Einstein's equations in the exterior region of the Cauchy domain reaching to future null infinity [56]. As a result, the equations are solved through the whole spacetime. For this strategy to work, appropriate communication between the regions treated by the Cauchy and characteristic approaches is required[57,69]. Although this technique has been successfully applied in a variety of cases, it has not yet been shown to work in 3D. At least part of the problem appears to lie in communicating all dynamical variables without actively enforcing the constraints [58]. Preliminary investigations of a revised matching strategy in the linearized regime has provided correct results but further work is needed to fully realize the promises of this approach; namely the possibility of accurately providing the waveforms but more generally producing a complete description of the spacetime.
In the coming years, as interferometric gravitational wave detectors take observational data, it will become very important to predict reliable and accurate gravitational waveforms from a variety of astrophysical sources. For these templates to be useful in gravitational wave data analysis, increased dialogue between numerical relativists and gravitational wave data analysts is urgently needed. This dialogue will help define the parameters of these template catalogues; namely: accuracy required, parameter space sampling, data format, etc. The program will dedicate time to spur this dialogue and expedite the transition from numerical simulation results to useful information for gravitational wave data analysis.
In the mergers of compact binaries (neutron star/neutron star & black hole/neutron star), hydrodynamic forces and microphysics become important once the binary components are close enough to induce tidal disruption. Because these mergers are one of the most promising sources of gravitational waves [70] and are possible sources of r-process nucleosynthesis [71] and short-duration gamma-ray bursts [72,73], they have been studied in some detail, with a variety of approximations (see [72,74,75,71,76,77] and references therein). For obtaining accurate templates of gravitational waveforms, it is likely that simulations using fully general relativistic gravity and realistic equations of state and microphysical processes will be needed.
The collapses which form neutron stars and black holes are of importance in many facets of astronomy (supernovae, hypernovae, gamma-ray bursts, sites of nucleosynthesis, and the formation of neutron stars and black holes), and are studied in great detail. Although general relativistic effects can be treated in the weak field limit for collapses which form neutron stars, they can significantly alter the dynamics [78,79] of the supernova explosion as can deviations from spherical symmetry[80,81]. Stellar collapses to black holes (collapsars, hypernovae, collapse of super-massive stars) require strong field general relativity. Most of the current work including general relativity has assumed spherical symmetry[82,78,83].
Neutron stars may be born rapidly rotating or they may be spun up by accretion. In both cases, the rotation rate may be sufficiently high to produce instabilities (e.g. bar-mode, r-mode) which may remove angular momentum and also produce detectable amounts of gravitational waves [84,81,85][86,87].
The goal of this focus topic is to identify these gravitational wave sources and the numerical techniques needed to study them. We must consider the entire evolution when determining the appropriate models for these simulations. For example, at first, a point-mass approximation is sufficient for modelling neutron star binary mergers, but once a merged object is formed, strong-field general relativity is required. If the hydrodynamics-dominated phase is not embedded into the appropriate general relativistic surroundings, it may be impossible to connect these two regimes. All of these problems require high resolution and large scales, demanding lagrangian and adaptive mesh techniques coupled to efficient parallelization algorithms. We will conduct a systematic survey of these sources and incorporate terascale computing resources to construct detailed models including gravitational waveforms. These simulations will be used first to predict measurements, but ultimately we hope to solve the reverse problem and use measurements to gain a better understanding of the astrophysical objects.
While the particular form of the discrete equations depends on both the chosen formulation of the equations as well as the space and time discretizations, a number of challenging computational issues arise which are common to most of the possible discrete formulations [43]. The goal of this topic is to study a number of the most pressing of these common issues, and to develop new techniques for addressing them. Three primary issues will be addressed through a number of subtopics: (1) adaptive numerical methods and supporting techniques for solving hyperbolic and elliptic systems, (2) datastructures and algorithms that support adaptivity on parallel platforms, and (3) supporting software infrastructure, development methodologies and tools (e.g. visualization, archiving, monitoring/steering, etc.).
Several questions involving the use of adaptive methods in numerical relativity will be addressed. For example, how do spectral and adaptive finite element methods compare to (locally uniform) finite difference (AMR) methods? Do spectral methods offer advantages over other methods, and can they be used effectively with adaptive techniques? Do new finite element approaches [88,89] show promise for strong field hyperbolic systems, and are these methods superior to other techniques currently in use? Can error estimates be established for the evolution equations, similar to those recently established for the constraint equations [90]? What types of error indicators should be used to drive adaptive mesh refinement algorithms?
A number of questions on algorithm and datastructure support will also be addressed. What computational infrastructure is necessary to support both structured (finite difference/AMR) and unstructured (finite element) mesh refinement frameworks, on sequential, parallel and distributed computers, and on computational ``Grids''? What types of (distributed) dynamic/adaptive datastructures are necessary for supporting these adaptive algorithms [91,92]? What dynamic load-balancing and scheduling algorithms work well with such datastructures [93]? Which distributed elliptic solvers are most effective in the adaptive mesh setting?
Regarding software and development infrastructure: can existing parallel software frameworks, such as GrACE/DAGH, Cactus and PETsc, be leveraged for adaptive refinement? Can a collaborative, robust, stable and reliable code base be developed which the numerical relativity community can collectively develop and leverage? Can advanced component-based development models for scientific software[94] and software engineering techniques for collaborative application development be leveraged and effectively applied? How can the large data sets produced by 3D simulations of Einstein's equations be managed to vizualise the outcome both in real time and in a postprocessing phase? Can interactive monitoring/steering systems [95,96] be used to provide runtime diagonostics and to increase the overall utility of the simulations?