Initial data and numerical relativity

20.1 An equation is not yet a prediction#

Chapter 19 found both an expanding and a contracting branch of a cosmological solution. The field equation alone did not choose one: we also needed an initial scale and rate of change. The same need appears without cosmological symmetry. To evolve a general gravitational field, specify its spatial geometry and how that geometry is changing, together with the matter data.

The rate-of-change information is carried by extrinsic curvature, which Chapter 14 introduced as the change of a boundary’s normal direction. We will relate that definition to the evolving spatial metric. Four projections of Einstein’s equation constrain the allowed initial geometry, extrinsic curvature, and matter; they cannot all be chosen independently.

For this chapter set $c=1$, so time and length have the same units. Keep $G_N$ explicit.

20.2 Spatial slices, lapse, and shift#

Choose spacelike hypersurfaces $\Sigma_t$ labeled by a coordinate $t$. Locally, wherever this foliation is valid, write the metric in 3+1 form:

Here are the three ingredients.

Let $n^\mu$ be the future-directed unit normal, $n^\mu n_\mu=-1$. The coordinate time vector decomposes as

The sign of the shift has operational meaning: someone following a normal worldline has $dx^i/dt=-\beta^i$. The coordinates drift relative to that person.

For a flat-spacetime example, start with $ds^2=-dT^2+dX^2+dY^2+dZ^2$ and set $T=2t$, $X=x+vt$, $Y=y$, $Z=z$, with constant dimensionless $v$. The metric becomes

Here $N=2$, $\beta^x=v$, and $\gamma_{ij}=\delta_{ij}$. A normal observer stays at fixed $X$, so $dx/dt=-v$ and $d\tau=2dt$. Lapse changes the clock labeling; shift changes the spatial labeling between slices. This example has no spacetime curvature.

There need not be a convenient global slicing of an arbitrary spacetime. We will meet the additional causal conditions that support one in Chapter 22.

SPACETIME LAB / 09

One step in time. Two different choices.

Lapse carries you along the normal. Shift slides the coordinate grid sideways. Rotate the slices to separate these two pieces.

ExploreDiagram

Interactive geometry is loading.

Drag to orbit · pinch to zoom

Coordinate shift

Negative shiftPositive shift

$$\begin{gathered}\partial_t=\htmlClass{math-observer}{Nn}+\htmlClass{math-transport}{\beta}\\\beta=\beta^i\partial_i\end{gathered}$$

Read the scene. Two spatial dimensions are shown on each slice; the vertical direction represents time schematically. The Euclidean display is not a spacetime distance measurement.

20.3 Extrinsic curvature and the change of spatial geometry#

Extend the spatial projector to spacetime:

It removes a vector’s component normal to the slice. Define extrinsic curvature with the following sign convention:

The derivative asks how the normals change as we move along the slice. The projectors retain the part visible within the slice. A plane has parallel normals; a curved surface generally does not.

For tangent indices this becomes

where $D_i$ is the Levi-Civita derivative of $\gamma_{ij}$. The Lie derivative $\mathcal L$ measures change under the flow of a vector field. In this equation it subtracts change caused merely by sliding the coordinates.

Equivalently,

This is a precise version of “extrinsic curvature contains the velocity of the spatial metric.” The qualification matters: it is velocity relative to the chosen slicing, after correcting for coordinate drift.

Extrinsic-curvature convention. Chapter 14 used the boundary convention $K_{\mathrm{boundary}}=h^{\mu\nu}\nabla_\mu n_\nu$. For the same spacelike hypersurface and the same normal, the present ADM convention gives $K=-K_{\mathrm{boundary}}$. Other textbooks also differ in this choice. Every equation containing an odd number of $K$ factors must be translated consistently. A sign difference here is not a disagreement about expanding universes.

Intrinsic and extrinsic curvature. extrinsic curvature need not indicate spacetime curvature. Curved slices can be drawn inside flat Minkowski spacetime. Intrinsic spatial curvature, extrinsic curvature, and four-dimensional spacetime curvature are related objects, not synonyms.

20.4 Four constraints on the initial data#

Define the matter energy and momentum seen by the normal observers:

$E$ is not automatically the fluid’s rest-frame density $\epsilon$: a moving fluid has extra energy in the normal frame. The minus sign in $j_i$ compensates for the timelike metric sign so that its ordinary local interpretation is momentum density.

The normal-normal projection of Einstein’s equation is the Hamiltonian constraint:

Here $K=\gamma^{ij}K_{ij}$ and ${}^{(3)}R$ is the scalar curvature of the spatial metric. The three mixed normal-spatial projections give the momentum constraints:

These combinations follow by comparing spacetime transport with transport restricted to the slice. To calculate that comparison at a chosen slice, use Gaussian normal coordinates locally: $N=1$, $\beta^i=0$, and $ds^2=-dt^2+\gamma_{ij}(t,x)dx^idx^j$. The normal geodesics construct this chart until they cross. Its connection contains

In the purely spatial Riemann components, the derivative terms and spatial connection products form ${}^{(3)}R_{ijkl}$. The two products with intermediate index 0 remain:

This is the Gauss relation. Contracting with $\gamma^{ik}\gamma^{jl}$ gives ${}^{(3)}R+K^2-K_{ij}K^{ij}$. The same contraction of spacetime curvature equals $2G_{\mu\nu}n^\mu n^\nu$: the normal contributions in $R+2R_{\mu\nu}n^\mu n^\nu$ cancel, leaving just the spatial contraction. Thus

Projecting $G_{\mu\nu}+\Lambda g_{\mu\nu}=8\pi G_NT_{\mu\nu}$ twice along $n$ contributes $-\Lambda$, because $g(n,n)=-1$. Move it across and multiply by two: the $+2\Lambda$ in the constraint follows. The mixed curvature components in this chart give the Codazzi relation,

The ordinary derivatives come from differentiating $\Gamma^0{}_{ij}=-K_{ij}$; the remaining products supply the spatial covariant-derivative corrections. Contracting yields $G_{0i}=R_{0i}=D_iK-D_jK^j{}_i$. Since $T_{0i}=-j_i$ in this normal chart, the mixed field equation gives exactly the momentum constraint above. These projected identities are tensorial, so the result does not depend on having used the convenient chart to derive them.

These equations contain no second time derivative of the geometry. They constrain what can consistently exist on one slice. The 3+1 projection and this extrinsic-curvature convention are developed systematically in Éric Gourgoulhon’s author-written notes on the 3+1 formalism.

Suppose we try to prescribe positive matter density together with an exactly Euclidean spatial metric and $K_{ij}=0$. With $\Lambda=0$, the left side of the Hamiltonian constraint is zero, so it requires $E=0$. To describe that matter, we must change the spatial geometry, its extrinsic curvature, or both.

This is analogous to specifying an electric field with zero divergence everywhere while also inserting a charge. Evolution cannot repair an inconsistent starting point without changing the data.

20.5 Worked check: the Friedmann equation is an initial-data constraint#

For homogeneous, isotropic slices,

Take normal observers comoving with the fluid, so $E=\epsilon$. From the definition,

Therefore $K=-3H$, while $K_{ij}K^{ij}=3H^2$. The difference is $9H^2-3H^2=6H^2$. The Hamiltonian constraint becomes

or

The familiar cosmological equation is the statement that the initial geometry, expansion rate, and matter density fit together. Homogeneity made the constraint algebraic. In a binary-black-hole calculation, solving its spatially varying version is substantial work.

20.6 Evolution and the two physical degrees of freedom#

The spatial Ricci projection contains the time derivative absent from the constraints. In the Gaussian normal chart used above, direct substitution of the same connection gives

For general lapse and shift, the time derivative becomes $N^{-1}(\partial_t-\mathcal L_\beta)K_{ij}$ and the projection has an additional term $-N^{-1}D_iD_jN$. Setting the vacuum Ricci tensor to zero, with $\Lambda=0$, gives

Read it together with the equation for $\partial_t\gamma_{ij}$. Spatial curvature and extrinsic-curvature products govern the next change of $K$; lapse gradients describe the acceleration associated with the chosen normal observers. Matter adds appropriate spatial stress and energy terms.

The Einstein-Hilbert action, after separating its boundary contribution, contains the bulk combination

Here $\gamma=\det(\gamma_{ij})$. One way to track the boundary contribution is the scalar identity

where $a^\mu=n^\nu\nabla_\nu n^\mu$ is the normal observers’ acceleration. Multiplying the last term by $\sqrt{-g}=N\sqrt\gamma$ turns it into an ordinary divergence, as in Chapter 14. Its boundary integral must be treated with the prescribed boundary data. The remaining interior terms give the bulk action displayed above.

There are no independent time derivatives of lapse and shift. In the Hamiltonian formulation they act as multipliers enforcing the constraints, rather than adding propagating gravitational polarizations.

Before counting, give phase space a concrete meaning. For a particle with coordinate $q$ and Lagrangian $L(q,\dot q)$, its conjugate momentum is $p=\partial L/\partial\dot q$. For $L=m\dot q^2/2-V(q)$ this gives $p=m\dot q$. A state requires both position and momentum: $(q,p)$ is one point of a two-dimensional phase space. A field has a coordinate value and its conjugate momentum at each spatial point. Conjugate momentum need not equal mass times velocity in a general Lagrangian; the derivative definition is the rule.

A constraint is an equation restricting the allowed states. Removing a redundant description is a separate operation. As a small model, start with $(q_1,q_2,p_1,p_2)$, impose $p_2=0$, and declare that changing $q_2$ does not change the physical state. The constraint removes one direction and the equivalence removes another, leaving the physical pair $(q_1,p_1)$.

The formal test uses the Poisson bracket, defined for ordinary canonical coordinates by

A constraint is first-class when its bracket with every constraint vanishes on the allowed constraint surface. In the regular canonical formulation of GR, the Hamiltonian and three momentum constraints are first-class and supply the associated gauge redundancy. Establishing their complete bracket algebra is an additional Hamiltonian calculation; we use that result here rather than deriving it from a count of components. Field-theory brackets replace the coordinate derivatives and sum above with functional derivatives and a spatial integral. The reduction to independent canonical variables is developed in Arnowitt, Deser, and Misner’s original account of GR dynamics.

Now count the metric sector, excluding lapse and shift as multipliers. The symmetric $\gamma_{ij}$ has six components. Its six conjugate momenta make twelve phase-space variables per spatial point. The four independent first-class constraints each remove one phase-space direction by imposing an equation and one by identifying gauge-equivalent descriptions:

That means two configuration degrees of freedom, each with its conjugate momentum. In weak gravitational waves, these become the two familiar polarizations. This is a local count in ordinary four-dimensional GR; boundaries, topology, special backgrounds, and matter require additional care.

“Ten minus four equals six” does not perform this count. It subtracts coordinate functions while overlooking the constrained dynamical structure.

20.7 A stable spacetime simulator needs a good coordinate policy#

Begin with a smaller question: can a computed solution look convincing while violating an equation it is supposed to obey? We can test this in the flat dust universe already solved in Chapter 19.

Measure the scale factor relative to its initial value, calling the ratio $A$. Measure elapsed time in units of the initial Hubble time, and call the resulting dimensionless time $t$. Write $V=dA/dt$ for the expansion rate. The dust acceleration equation and initial conditions become

The Friedmann constraint is $\mathcal C=V^2-1/A=0$. Differentiating it gives $\dot{\mathcal C}=2V\dot V+\dot A/A^2=0$: an exact evolution preserves the constraint. A numerical evolution uses finite steps, so this cancellation need not remain exact.

The known solution $A(t)=(1+3t/2)^{2/3}$ gives us two checks. Compare the computed scale factor with the exact curve, then inspect the constraint residual. Halve the step and repeat, keeping the final time fixed. A smaller residual and a more accurate scale factor are related evidence, but they are different measurements.

Numerical relativity lab

Make a universe—and audit the answer

Will a smaller time step repair the apparent motion, the constraint, or both?

• Numerical scale factor
• Exact scale factor

Final scale factor A3.5571

Error against the exact solution-0.10218

Friedmann constraint residual-0.057916

Change in A when step is halved0.051543

Compare Euler with RK4, then reduce the step. A curve can look convincing while its constraint residual reveals an error.

How this is calculated

Dimensionless flat-dust reduction with A(0) = V(0) = 1, zero cosmological constant and exact A(t) = (1 + 3t/2)^(2/3).

$$\dot A=V,\qquad\dot V=-\frac1{2A^2}$$

$$\mathcal C=V^2-\frac1A=0$$

• Homogeneous FLRW dust: Ȧ = V and V̇ = −1/(2A²).
• Choose first-order forward Euler or classical fourth-order Runge–Kutta.
• All refinements reach the same final time; the actual step can be slightly smaller than the requested step.

This evolves a symmetry reduction of Einstein’s equations. It contains no gravitational waves or spatial grid and does not establish well-posedness or stability of a full numerical-relativity formulation.

Tong · Einstein equations and cosmology ↗

Measurements and notebook

Dimensionless flat-dust reduction with A(0) = V(0) = 1, zero cosmological constant and exact A(t) = (1 + 3t/2)^(2/3).

t	scale	exact	error	constraint
0	1	1	0	0
0.4	1.38	1.368	0.012019	-0.034815
0.8	1.7017	1.6915	0.010182	-0.045632
1.2	1.9893	1.9866	0.0027221	-0.050529
1.6	2.2532	2.2611	-0.0078808	-0.053199
2	2.4993	2.5198	-0.020569	-0.054829
2.4	2.7311	2.7659	-0.034798	-0.055903
2.8	2.9512	3.0015	-0.050247	-0.056651
3.2	3.1614	3.2281	-0.066714	-0.057195
3.6	3.363	3.4471	-0.084058	-0.057602
4	3.5571	3.6593	-0.10218	-0.057916

How does the computer take one time step?

The state is the pair $y=(A,V)$, and its derivative is $f(y)=(V,-1/(2A^2))$. Forward Euler follows the current slope for one step $h$: $y_{n+1}=y_n+h f(y_n)$. It treats that slope as constant across the step.

The fourth-order Runge–Kutta method, abbreviated RK4, samples a beginning slope, two trial midpoint slopes, and a trial endpoint slope:

$$k_1=f(y_n),\qquad k_2=f(y_n+hk_1/2),$$

$$k_3=f(y_n+hk_2/2),\qquad k_4=f(y_n+hk_3),$$

$$y_{n+1}=y_n+\frac h6(k_1+2k_2+2k_3+k_4).$$

The weighted slopes account for the changing derivative inside the step. For smooth solutions in the regime where truncation error dominates, halving $h$ reduces the accumulated Euler error by roughly two and the RK4 error by roughly sixteen. This is a convergence expectation to test, rather than an error bound for an arbitrary calculation.

Returning to a general spacetime. The experiment evolves a homogeneous universe with no spatial grid. A full spacetime evolution must also handle coordinate freedom, disturbances propagating across the grid, and boundaries.

Einstein’s equations contain gauge freedom, so their unreduced component form is not simply ten independent wave equations. A coordinate condition can expose the wave structure. In harmonic coordinates, for example,

and the principal, highest-derivative part of the reduced metric equations is schematically

The metric itself supplies the coefficients that determine wave propagation: the unknown metric appears in the coefficients of its own highest derivatives. Those highest derivatives enter linearly, which is the meaning of quasilinear.

A well-posed formulation has a solution, has the appropriate uniqueness, and makes that solution depend continuously on the initial data. The last requirement bounds how errors in the starting data affect the solution over a specified time interval. A formulation can be mathematically equivalent on exact constraint-satisfying solutions yet behave very differently when roundoff and discretization introduce small constraint violations. Generalized harmonic formulations prescribe the contracted connection through coordinate equations. The BSSN formulation, named for Baumgarte, Shapiro, Shibata, and Nakamura, instead separates the spatial volume factor from a unit-determinant spatial metric and separates the trace of extrinsic curvature from its trace-free part. It also evolves auxiliary connection variables. These are distinct organizations of the same physical solution, with different responses to numerical errors; deriving a full implementation goes beyond the homogeneous benchmark above.

The contracted Bianchi identity supplies constraint-propagation relations. With consistent matter evolution, exact constraints that hold initially continue to hold in a suitable exact evolution. A discretized evolution introduces errors, so constraint residuals must still be monitored and checked for convergence.

Initial conditions and boundary conditions play different roles. Initial data describe a spatial slice. At a finite simulation boundary, combinations of field disturbances propagate inward or outward at speeds set by the chosen equations; these combinations are called characteristic fields. Boundary data must handle the incoming combinations consistently, including physical radiation and changes in coordinates or constraint errors. For an isolated system one often approximates an asymptotically flat exterior; a reflecting boundary can send outgoing radiation back into the modeled region.

Some constraint formulations give spatial boundary-value equations of the same general type as Poisson’s equation, called elliptic equations. Solving them across a slice does not transmit a physical signal instantly. It constructs a mutually compatible initial state. Subsequent physical disturbances propagate according to the causal equations.

WORKED EXAMPLE

Let the constraint audit your simulation

Can a smooth-looking numerical universe solve the wrong initial-value problem?

Builds onInitial data Friedmann constraint

See the idea

A numerical solution must satisfy more than an evolution formula. Einstein’s equations also constrain the initial state and provide quantities that should remain zero as the solution evolves. A homogeneous dust universe gives a small, reproducible version of this principle.

Work it out

1. 01
   Reduce to a specified toy initial-value problem

   For spatially flat dust with zero cosmological constant, choose dimensionless time and scale factor A so the initial conditions are $A(0)=1$, $V(0)=1$. Mass conservation and the acceleration equation reduce to $\dot A=V$, $\dot V=-1/(2A^2)$. The Friedmann constraint is $V^2=1/A$.

   $$\mathcal C=V^2-\frac1A=0.$$

   Why this step works The time scale was chosen so the conserved matter coefficient is one.

2. 02
   Verify the exact benchmark

   Solving the expanding constraint gives $\dot A=A^{-1/2}$. Integrating yields $A=(1+3t/2)^{2/3}$ and $V=(1+3t/2)^{-1/3}$ for $t>-2/3$. Differentiate both expressions to verify the evolution equations.

   $$\frac{d\mathcal C}{dt}=2V\left(-\frac1{2A^2}\right)+\frac{V}{A^2}=0.$$

   Why this step works Constraint propagation follows from the exact evolution, provided the initial constraint is satisfied.

3. 03
   Expose an integration error

   Forward Euler updates $A_{n+1}=A_n+hV_n$ and $V_{n+1}=V_n-h/(2A_n^2)$ using only the old state. Starting from (1,1), one step gives $A_1=1+h$, $V_1=1-h/2$. Substitution into the constraint reveals a nonzero residual at order $h^2$.

   $$\mathcal C_1=(1-h/2)^2-\frac1{1+h}=-\frac34h^2+O(h^3).$$

   Why this step works An evolution update can accumulate truncation error even when the starting data are exact.

Go deeper

Compare the solution at the same final time for h, h/2, and h/4. In the asymptotic regime, a method of global order p has successive differences with ratio about $2^p$, provided other errors do not dominate. Convergence toward the analytic solution and shrinking constraint residual support this particular calculation. A residual alone is not a complete error estimate: even the wrong physical solution can satisfy a constraint. This is a symmetry reduction of Einstein’s equations, not a simulation of radiative three-dimensional gravity. Full numerical relativity must also handle gauge, spatial derivatives, boundaries, and a well-posed evolution formulation.

Further readingDavid Tong · General relativity

Test the idea

FIRST, PREDICT

A numerical scale factor looks smooth. Is that enough to validate the simulation?

No. Check the constraint, convergence, and an independent benchmark.Yes. Physical solutions are smooth.A zero constraint residual proves every observable is correct.

Compare the reasoning

No. Check the constraint, convergence, and an independent benchmark.

Visual smoothness does not measure truncation error or verify the equations.

Yes. Physical solutions are smooth.

Many smooth curves do not satisfy the stated equations or initial data.

A zero constraint residual proves every observable is correct.

Constraint satisfaction does not by itself establish the correct evolution or initial state.

A hint

Ask what quantity would detect a believable but inaccurate curve.

NOW CHANGE THE EXAMPLE

Take one Euler step of size h=0.2 from (A,V)=(1,1). Calculate $\mathcal C_1=V_1^2-1/A_1$.

A hint

Use $A_1=1.2$ and $V_1=0.9$.

Work through the solution

$\mathcal C_1=0.81-1/1.2=-7/300\simeq-0.023333$.

A useful simulation reports its assumptions, constraint residual, convergence, and benchmark error.

20.8 What uniqueness means when coordinates are free#

Imagine a smooth relabeling of spacetime that is exactly the identity near the initial slice but changes labels inside a later empty region—the “hole.” Apply it to the metric and all physical fields. General covariance produces a new coordinate description satisfying the same initial data. Does that destroy determinism?

Only if you assume that bare manifold points already possess observable identities independently of every field. The two descriptions preserve coincidences: where a detector meets a pulse, how much proper time its clock records, what curvature its instruments measure. They represent the same physical solution when related by the appropriate gauge diffeomorphism.

A prediction should concern “the curvature measured when this clock reads this value,” not “the curvature at a label whose attachment to any physical event I am free to change.” Boundary symmetries require care: transformations acting nontrivially on prescribed asymptotic data can carry physical charges and are not all disposable gauge.

For suitable constraint-satisfying vacuum data, and for appropriate well-posed matter systems, the relevant uniqueness statement is uniqueness of the maximal globally hyperbolic development up to diffeomorphism. “Maximal” does not promise geodesic completeness or a nonsingular future. This is the landmark result of Choquet-Bruhat and Geroch’s original Cauchy-problem paper.