Numerical relativity
1797708
223700062
2008-07-05T09:43:35Z
Atonita
1925670
Expanded the computational difficulties section.
'''Numerical relativity''' is one of the branches of [[general relativity]] that uses numerical methods and algorithms to solve and analyze problems. To this end, [[supercomputers]] are often employed to study [[black holes]], [[gravitational waves]], [[neutron stars]] and many other phenomena governed by [[Albert Einstein|Einstein's]] [[General theory of relativity|Theory of General Relativity]]. A current active research in numerical relativity is the simulation of relativistic binaries and their associated gravitational waves, other branches are also quite active.
==Overview==
A primary goal of numerical relativity is to study [[spacetime]]s whose exact form is not known. The spacetimes so found computational can either be fully [[dynamical]], [[stationary spacetime|stationary]] or [[static spacetime|static]] and may contain matter fields or vacuum. In the case of stationary and static solutions, numerical methods may also be used to study the stability of the equilibrium spacetimes. In the case of dynamical spacetimes, the problem may be divided into the initial value problem and the evolution, each requiring different methods.
Numerical relativity is applied to many areas, such as [[Physical cosmology|cosmological]] [[Model (abstract)|models]], [[critical phenomena]], [[Perturbation (astronomy)|perturbed]] [[black hole]]s and [[neutron star]]s, and the [[coalescence (meteorology)|coalescence of black hole]]s and neutron stars, for example. In any of these cases, Einstein's equations can be formulated in several ways that allow us to evolve the dynamics. While [[Cauchy]] methods have received a majority of the attention, characteristic and [[Regge calculus]] based methods have also been used. All of these methods begin with a snapshot of the [[gravitational field]]s on some [[hypersurface]], the initial data, and evolve these data to neighboring hypersurfaces. [http://relativity.livingreviews.org/Articles/lrr-2000-5/]
Like all problems in numerical analysis, careful attention is paid to the [[numerical stability|stability]] and [[Numerical_ordinary_differential_equations#Convergence|convergence]] of the numerical solutions. In this line, much attention is paid to the [[gauge fixing|gauge conditions]], coordinates, and various formulations of the Einstein equations and the effect they have on the ability to produce accurate numerical solutions.
Numerical relativity research is distinct from work on [[Classical field theory|classical field theories]] as many techniques implemented in these areas are inapplicable in relativity. Many facets are however shared with large scale problems in other computational sciences like [[computational fluid dynamics]], electromagnetics, and solid mechanics. Numerical relativists often work with applied mathematicians and draw insight from [[numerical analysis]], [[scientific computing|scientific computation]], [[partial differential equation]]s, and [[geometry]] among other mathematical areas of specialization.
==The "3+1" Approach to Numerical Relativity==
Most active research in NR has involved a 3+1, or "space<math>\times</math>time" decomposition of Einstein's equations. Here the spacetime of interest is sliced up into a stack (or "[[foliation]]") of spacelike three-dimensional hypersurfaces. That is, cross-sections are taken of the spacetime, so that each point in the spacetime lies on exactly one slice. Successive slices then represent--loosely speaking--different instants in time. More specifically, the hypersurfaces are parametrized by a coordinate usually denoted by ''t'', which can be interpreted as a universal coordinate time. (Note that this time does not necessarily coincide with the [[proper time]] of any observer.)
If we now restrict attention to a single hypersurface, we see that it can be treated just like any other geometric space. For example, the space has some [[Riemannian metric]]. This three-dimensional metric is usually written as <math>\gamma_{ij}</math>, so that the [[line element]] is given by <br />
<math> dl^2 = \gamma_{ij} dx^i dx^j\ . </math> <br />
Here, the coordinates <math>x^i</math> are just some coordinates placed on the hypersurface, in some smooth way, and the indices ''i'' and ''j'' run over the three dimensions of the hypersurface.
Next, we need some way to relate the successive hypersurfaces. We imagine a set of observers moving through time along the direction normal to the hypersurface. If we choose one such observer, and denote the proper time by <math>\tau</math>, we can relate it to the time coordinate ''t'' by way of the "lapse function" <math>\alpha</math>: <br />
<math>d\tau = \alpha(t, x^k) dt\ .</math> <br />
Similarly, we need to relate the position of the observer on one slice to its position on the next slice, using coordinates. If the change in time coordinate is ''dt'', and the change in spatial coordinates is <math>dx^i</math>, we can relate the two as: <br />
<math>dx^i = -\beta^i(t,x^k)\, dt\ .</math> <br />
The functions <math>\beta^i</math> are known as the "shift vector".
The way in which we foliate a spacetime is certainly not unique, as there is clearly no need (other than smoothness) for the spatial coordinates on neighboring slices to be related to each other uniquely. This means that the functions <math>\alpha(t,x^k)</math> and <math>\beta^i(t,x^k)</math> can be freely chosen. These two functions determine our choice of coordinate system, and are known as "[[gauge freedom|gauge]] functions" for this reason.
The metric of the full spacetime can now be written in the form <br />
<math>
ds^2
=
\left( -\alpha^2 + \beta^i \beta^j \gamma_{ij} \right) dt^2
+
2\beta^i \gamma_{ij} dt dx^j + \gamma_{ij} dx^i dx^j\ .
</math> <br />
In computations for General Relativity, the desired "solution" is just the metric, given as a function of the coordinates. From the above equation, then, we see that objective is to find <math>\gamma_{ij}(t,x^k)</math>.
The mathematics of the 3+1 split of Einstein's equations was laid out in a series of papers by Arnowitt, Deser and Misner ([[ADM formalism|ADM]]) in the early 1960s; the form actually used for most simulation efforts until the mid-1990s was described by York (1979).<ref>York, J W "Kinematics and Dynamics of General Relativity" from "Sources of Gravitational Radiation: Proceedings of the Battelle Seattle Workshop" ed. Larry Smarr, pp. 83-126 (Cambridge University Press 1979)</ref> ADM showed that data could be specified on an initial hypersurface, and Einstein's equations then show how that initial data evolves onto later hypersurfaces. More specifically, ADM showed that Einstein's equations separate into four "constraint" equations, and a set of "evolution" equations.
The four constraint equations are equations which only refer to quantities defined on each individual hypersurface, and must be satisfied on each hypersurface--including the initial slice. The evolution equations then describe how those quantities change in time. (In technical terms, ADM showed that Einstein's equations form a constrained [[hyperbolic system]].)
==Computational Difficulties==
In addition to the normal difficulties associated with the development of stable and accurate numerical methods, there are some particular problems inherent to numerical relativity.
It can be shown that if the constraint equations are satisfied initially, and the variables are evolved exactly according to the evolution equations, then the constraints will remain satisfied at all later times. In theory, this means that it is only necessary to specify initial data which satisfy the constraints and evolve according to the evolution equations to have a correct solution to Einstein's equations. Unfortunately, discretising the problem produces complications. In an analogous fashion to numerical methods attempting to solve the [[Maxwell equations]] where care must be taken to preserve the divergence free condition of the magnetic field, numerical methods which attempt to solve Einstein's equations must take care to keep the magnitude of the constraint violation under control or the solution loses validity. To this end, [[hyperbolic divergence cleaning]] or [[constrained transport]] schemes are often used to evolve the Einstein equations. In addition, much work is done in fomulating constraint conserving boundary conditions at the analytic level.
A primary goal of numerical relativity is the description of black hole spacetimes, this is problematic from a numerical standpoint because of the difficulties associated with representing a physical singularity on a computational domain. Two techniques are commonly used to treat black holes computationally: excision methods, whereby the singular region is removed and treated by an apparent horizon inner boundary and puncture methods, whereby the singularity is essentially factored out.
==Methodology==
===Elliptic Equations===
In the case of producing stationary spacetimes, static spacetimes or producing initial data, three dimensional [[boundary value problem|elliptic equations]] must be solved. In some cases, symmetry may effectively reduce the dimensionality of the problem. The [[well-posed problem|well posedness]] of these problems depend crucially on the specification of the boundary values. It is common to suppose that isolated systems in general relativity are [[asymptotically flat spacetime|asymptotically flat]], which gives [[dirichlet boundary conditions]] if compact coordinates are used or a [[Robin boundary condition]] which suitably approximates asymptotic flatness if the domain is truncated for computational simplicity. In the case of spacetimes containing a black hole, an inner surface is specified as an apparent horizon which yields a [[Robin boundary condition]].
The resulting partial differential equations are then solved using [[finite difference methods]], [[finite element methods]], [[spectral methods]] or any other numerical technique for solving elliptic equations. Other specific problems in general relativity also require the solution of elliptic equations.
===Evolution===
In general, an evolution equation expresses the time derivative of a variable in terms of spatial derivatives of that variable and other variables. Implementing these equations on a computer, then, requires methods to take spatial derivatives of the functions (spatial differentiation), and use the resulting time derivative to evolve the variables (time stepping).
====Spatial Differentiation====
For ''D'' spatial dimensions in a simulation where physical phenomena are reduced to a set of difference equations (restricting attention to [[Finite difference|finite-difference method]]s) with a corresponding ''D''-dimensional array of ''N'' "mesh points", the [[CPU]] time is of [[Big O notation|order]] <math>\mathcal{O}(N^{D+1})</math> and the [[Memory (computers)|memory]] required is of order <math>\mathcal{O}(N^{D})</math>. Typically, the accuracy of finite-difference methods improves with the number of points as <math>N^2</math> (second-order), <math>N^4</math> (fourth-order), or some times <math>N^6</math> (sixth-order).
An alternative technique is to approximate each variable in a spectral series. Examples of spectral series include [[Fourier series]], expansion in terms of [[Chebyshev polynomials]], and [[spherical harmonics]]. These approximations can then be differentiated exactly. The advantage of this technique is the speed with which its accuracy improves (its [[convergence]] rate); the accuracy of spectral methods typically improves as <math>e^N</math>. Compared to <math>N^2</math> or <math>N^4</math> for finite differencing, this is a vast improvement.
The disadvantage of spectral methods is their sensitivity to smoothness of the functions they are used with. Especially in situations involving matter, shocks and discontinuities can arise easily. In those cases, spectral methods are not worth the trouble, and finite differencing methods are preferred. Nonetheless, many simulations of solutions of the vacuum Einstein equations can be significantly improved by using spectral methods. The great majority of Numerical Relativity research groups currently use finite differencing methods, for reasons of simplicity.
====Time Stepping====
The most common techniques for time stepping are Iterated [[Crank-Nicolson method|Crank-Nicolson]] and [[Runge-Kutta]]. Iterated Crank-Nicolson is a second-order convergent method, which means that its errors decrease with time step <math>\Delta t</math> as <math>1/(\Delta t)^2</math>. Runge-Kutta methods are somewhat more general. The "classic" Runge-Kutta algorithm is fourth-order convergent, but related algorithms can be made arbitrarily convergent.
There is a penalty, however, for moving to higher-order methods. These methods require more computations, and therefore run more slowly. The method used in a given simulation depends on balancing the accuracy with the time required. Many simulations which use second-order finite differencing, for example, use Iterated Crank-Nicolson time stepping because higher accuracy is not always needed.
==See also==
*[[Mathematics of general relativity]]
*[[Post-Newtonian expansion]]
==Notes==
{{reflist}}
==External links==
* [http://relativity.livingreviews.org/Articles/lrr-2000-5/ Initial Data for Numerical Relativity] — A review article which includes a nice (technical) discussion of numerical relativity.
* [http://www.emis.ams.org/journals/LRG/Articles/lrr-2003-3/node19.html Rotating Stars in Relativity] — A (technical) review article about rotating stars, with a section on numerical relativity applications.
* [http://www.black-holes.org/numrel1.html A Relativity Tutorial at Caltech] — A basic introduction to concepts of Numerical Relativity.
*[http://xstructure.inr.ac.ru/x-bin/theme2.py?arxiv=gr-qc&level=2&index1=12 Numerical relativity on arxiv.org]
[[Category:Computational physics]]
[[Category:Mathematical methods in general relativity]]