BSSN Formulation

For two decades after Arnowitt, Deser and Misner recast general relativity as an initial-value problem in 1962, numerical relativists faced a maddening obstacle: the ADM evolution equations, though mathematically exact, are only weakly hyperbolic. Feed them clean, constraint-satisfying black-hole data and the simulation still dies — small truncation and gauge errors, instead of propagating away at the speed of light, sit and amplify under the nonlinear feedback, so the code crashes within roughly a single light-crossing time. The physics was right; the formulation was numerically unstable. Every attempt to simulate the one problem the field existed to solve — two black holes spiralling together — ran into the same wall. The BSSN formulation dismantles that wall without touching the physics. Following Nakamura, Oohara and Kojima, then Shibata and Nakamura, and finally Baumgarte and Shapiro, one rewrites the same ADM equations in new variables: a conformal rescaling $\gamma_{ij}=e^{4\phi}\tilde\gamma_{ij}$ that separates the metric's "size" from its "shape", a trace split of the extrinsic curvature into $K$ and $\tilde A_{ij}$ , and — the decisive move — the promotion of the contracted conformal connection $\tilde\Gamma^i=-\partial_j\tilde\gamma^{ij}$ to an independently evolved field with its own equation. That single step pulls the dangerous mixed second derivatives out of the Ricci tensor and turns the principal part into a clean set of coupled wave operators. The reshuffled system is strongly hyperbolic; errors now propagate and leave the grid. Married to the moving-puncture gauge — $1+\log$ slicing and a Gamma-driver shift — BSSN let black holes advect freely across the grid without excision, and in 2005–2006 produced the first stable binary-black-hole inspiral–merger–ringdown waveforms. It remains, with the generalized-harmonic formulation, one of the two engines of modern numerical relativity and the source of the templates LIGO–Virgo–KAGRA match against data.

Key takeaways

BSSN rewrites the ADM 3+1 evolution equations in conformal-traceless variables — the same physics, a numerically stable form — converting a weakly hyperbolic system into a strongly hyperbolic one.
The evolved fields are the conformal factor $\phi$ , the unit-determinant conformal metric $\tilde\gamma_{ij}$ , the trace $K$ , the conformal traceless curvature $\tilde A_{ij}$ , and the conformal connection functions $\tilde\Gamma^i$ — 17 fields where ADM evolved 12.
The pivotal trick is promoting $\tilde\Gamma^i=-\partial_j\tilde\gamma^{ij}$ to an independent variable, which removes the mixed second metric derivatives from the Ricci tensor and restores a wave-like principal part.
Two algebraic conditions, $\det\tilde\gamma=1$ and $\tilde\gamma^{ij}\tilde A_{ij}=0$ , must be actively re-imposed, and the momentum constraint substituted into the $\tilde\Gamma^i$ equation, or stability is lost — BSSN is not automatically constraint-damped.
With the moving-puncture gauge ( $1+\log$ slicing plus a Gamma-driver shift) BSSN evolves black holes without excision and powered the 2005–2006 binary-black-hole breakthroughs.
Key idea: promote the contracted conformal connection $\tilde\Gamma^i=-\partial_j\tilde\gamma^{ij}$ to an independently evolved variable — this eliminates the mixed second derivatives of the metric in the conformal Ricci tensor, replacing them with first derivatives of $\tilde\Gamma^i$ and restoring strong hyperbolicity to the otherwise weakly hyperbolic ADM system.

Where this lives today

BSSN is running right now on supercomputers around the world. The open-source Einstein Toolkit ships BSSN-based thorns (McLachlan and its successors) that hundreds of groups use to evolve black-hole and neutron-star spacetimes, and production codes such as BAM, LazEv and the newer Z4c and CCZ4 constraint-damped descendants of BSSN generate the numerical-relativity waveform catalogs — the SXS, RIT and Georgia Tech collections — that calibrate the effective-one-body and phenomenological templates matched against LIGO–Virgo–KAGRA data. Every confident statement about the masses and spins of GW150914 or the hundred-plus mergers since traces back through those templates to a BSSN-family evolution. Coupled to relativistic magnetohydrodynamics, BSSN also drives binary-neutron-star simulations of the kind that interpreted GW170817 and its kilonova. The frontier now is accuracy and reach: constraint-damped formulations (CCZ4) that actively suppress rather than merely propagate constraint violations, extreme mass ratios and high spins where the moving-puncture gauge strains, and the eccentric and precessing corners of parameter space that next-generation detectors — LISA, Einstein Telescope, Cosmic Explorer — will demand. The reformulation that first made the problem tractable is still the workhorse of the theory it made possible.

The physical picture

Picture the ADM Ricci tensor as a wave operator for the metric that has been sabotaged. Its principal symbol is degenerate: alongside the physical modes that propagate at light speed sit zero-frequency gauge and constraint modes that neither travel nor decay. Any error you seed into those directions has nowhere to go — it stays on the grid and the nonlinear terms feed it. BSSN's cure is structural, not dissipative: give the offending second-derivative combination its own name, $\tilde\Gamma^i$ , and its own first-order evolution equation. Once that combination is a genuine field rather than a hidden derivative of the metric, the remaining metric second derivatives collapse into a clean Laplacian, and the system reads as coupled wave equations with honest characteristic speeds. The other half of the construction is the conformal split, borrowed from the way one solves the initial-data constraints and recycled here for evolution. Separating the volume degree of freedom $\phi$ from the shape $\tilde\gamma_{ij}$ (frozen to unit determinant) lets a black hole's diverging curvature be packed into $\phi$ alone — handled through the smooth variable $\chi=e^{-4\phi}$ that vanishes gently at a puncture — while the shape variables stay finite. The parallel trace split of $K_{ij}$ into $K$ and $\tilde A_{ij}$ sends the slice's expansion into a single scalar that the $1+\log$ lapse can tame, and leaves the radiative content in the traceless $\tilde A_{ij}$ .

Figure

The conformal-traceless split at the heart of BSSN. On a spacelike slice the physical metric factorizes into a volume part carried by e to the 4φ and a unit-determinant shape part, γ_ij = e^{4φ} γ̃_ij with det γ̃ = 1, while the extrinsic curvature splits into its trace K and a traceless remainder Ã_ij. The lapse α and shift βⁱ carry the slice forward in coordinate time

The construction is deliberately gauge-aware, which is why BSSN never travels alone. Hyperbolicity is a property of the coupled metric-plus-gauge system: a badly chosen shift can spoil it even with the conformal machinery in place. So the lapse and shift come with their own evolution equations — the $1+\log$ slicing that feeds on $K$ and the Gamma-driver shift that chases $\tilde\Gamma^i$ — and those gauge equations are as much a part of the formulation as the metric equations. The payoff is best seen not in an equation but in a stability curve: where ADM's constraint violation runs away and the code dies, BSSN's stays bounded and convergent for as long as you care to run it.

Figure

Why BSSN replaced ADM: schematic growth of the constraint-violation norm during a black-hole evolution. The weakly hyperbolic ADM system amplifies seeded error exponentially and crashes near one light-crossing time, t ≈ 40 M; the strongly hyperbolic BSSN system with 17 evolved fields keeps the violation bounded and convergent out to hundreds of M, letting the moving punctures orbit indefinitely

The core equations

BSSN is defined by three algebraic changes of variable and the evolution equations they induce. The first is the conformal rescaling that separates volume from shape. The unit-determinant condition on $\tilde\gamma_{ij}$ makes $\phi$ the only variable that carries the local volume, and forces $\phi=\tfrac{1}{12}\ln\det\gamma$ because a conformal factor $e^{4\phi}$ scales a three-metric's determinant by $e^{12\phi}$ .

\gamma_{ij}=e^{4\phi}\,\tilde\gamma_{ij},\qquad \det\tilde\gamma_{ij}=1,\qquad \phi=\tfrac{1}{12}\ln\det\gamma_{ij} \qquad (1)

where: $\gamma_{ij}$ — physical spatial 3-metric on the slice; $\tilde\gamma_{ij}$ — conformal metric, normalized to unit determinant; $\phi$ — logarithmic conformal factor, a BSSN variable; $e^{4\phi}$ — conformal scaling carrying the local volume and curvature. Reading it: the size lives in $\phi$ , the shape in $\tilde\gamma_{ij}$ . The second change of variable splits the extrinsic curvature the same way, peeling off its trace $K$ (the local expansion rate of the foliation) from a conformal traceless remainder $\tilde A_{ij}$ that carries the gravitational-wave content. One checks the split is consistent by tracing it with $\gamma^{ij}$ , which returns $K$ exactly because $\tilde\gamma^{ij}\tilde\gamma_{ij}=3$ .

K_{ij}=e^{4\phi}\!\left(\tilde A_{ij}+\tfrac{1}{3}\tilde\gamma_{ij}K\right),\qquad K=\gamma^{ij}K_{ij},\qquad \tilde\gamma^{ij}\tilde A_{ij}=0 \qquad (2)

where: $K_{ij}$ — physical extrinsic curvature of the slice; $K$ — its trace, the expansion of the foliation; $\tilde A_{ij}$ — conformal traceless extrinsic curvature, a BSSN variable carrying the radiative shear. The third and decisive change introduces the contracted conformal connection. For a unit-determinant metric the trace of the Christoffel symbols vanishes, so this contraction collapses to a pure divergence of the inverse conformal metric — and it is this object, not any derivative of it, that BSSN evolves as a field.

\tilde\Gamma^i \equiv \tilde\gamma^{jk}\tilde\Gamma^i{}_{jk}=-\partial_j\tilde\gamma^{ij} \qquad (3)

where: $\tilde\Gamma^i{}_{jk}$ — Christoffel symbols of the conformal metric $\tilde\gamma_{ij}$ ; $\tilde\Gamma^i$ — contracted conformal connection, promoted to an independent BSSN field. Its independence from $\tilde\gamma_{ij}$ during evolution is the whole trick. Taking the trace of the ADM metric-evolution equation and dividing by twelve gives the evolution of $\phi$ : the volume of a fluid element shrinks with the expansion $K$ and is advected and sheared by the shift.

\partial_t\phi=-\tfrac{1}{6}\alpha K+\beta^i\partial_i\phi+\tfrac{1}{6}\partial_i\beta^i \qquad (4)

where: $\alpha$ — lapse function; $\beta^i$ — shift vector; $\partial_t$ — coordinate-time derivative along the foliation. The trace-free part of the same ADM equation drives $\tilde\gamma_{ij}$ via $\tilde A_{ij}$ .

\partial_t\tilde\gamma_{ij}=-2\alpha\tilde A_{ij}+\beta^k\partial_k\tilde\gamma_{ij}+\tilde\gamma_{ik}\partial_j\beta^k+\tilde\gamma_{kj}\partial_i\beta^k-\tfrac{2}{3}\tilde\gamma_{ij}\partial_k\beta^k \qquad (5)

where: $\tilde A_{ij}$ — conformal traceless curvature, the shape driver; the $-\tfrac23\tilde\gamma_{ij}\partial_k\beta^k$ term keeps the right-hand side trace-free so $\det\tilde\gamma=1$ is preserved to leading order. The extrinsic-curvature sector splits the same way. The traceless part $\tilde A_{ij}$ is driven by the trace-free part of the lapse Hessian and the 3-Ricci tensor — this is the equation that carries the gravitational-wave content. The trace $K$ obeys an equation in which the 3-Ricci scalar has been eliminated using the Hamiltonian constraint — the substitution that removes second spatial derivatives of the metric from the trace sector and leaves only the Laplacian of the lapse and a quadratic curvature source.

\partial_t\tilde A_{ij}=e^{-4\phi}\big[-D_iD_j\alpha+\alpha R_{ij}\big]^{\mathrm{TF}}+\alpha\!\left(K\tilde A_{ij}-2\tilde A_{ik}\tilde A^k{}_j\right)+\beta^k\partial_k\tilde A_{ij}+\tilde A_{ik}\partial_j\beta^k+\tilde A_{kj}\partial_i\beta^k-\tfrac{2}{3}\tilde A_{ij}\partial_k\beta^k,\qquad \partial_t K=-\gamma^{ij}D_iD_j\alpha+\alpha\!\left(\tilde A_{ij}\tilde A^{ij}+\tfrac{1}{3}K^2\right)+\beta^i\partial_i K \qquad (6)

where: $D_i$ — covariant derivative of the physical 3-metric $\gamma_{ij}$ ; $R_{ij}$ — its 3-Ricci tensor; $[\,\cdot\,]^{\mathrm{TF}}$ — the trace-free part, taken with $\tilde\gamma_{ij}$ ; $\tilde A^k{}_j=\tilde\gamma^{kl}\tilde A_{lj}$ ; $\tilde A_{ij}\tilde A^{ij}$ — quadratic curvature source, indices raised with $\tilde\gamma$ . The first equation is where the radiative degrees of freedom live (in vacuum; matter adds a trace-free stress term); the second is the equation the $1+\log$ lapse feeds on. The conformal connection functions evolve by an equation obtained from differentiating their definition and substituting the momentum constraint for the divergence $\partial_j\tilde A^{ij}$ — the step that buys strong hyperbolicity.

\partial_t\tilde\Gamma^i=\tilde\gamma^{jk}\partial_j\partial_k\beta^i+\tfrac{1}{3}\tilde\gamma^{ij}\partial_j\partial_k\beta^k-2\tilde A^{ij}\partial_j\alpha+2\alpha\!\left(\tilde\Gamma^i{}_{jk}\tilde A^{jk}-\tfrac{2}{3}\tilde\gamma^{ij}\partial_jK+6\tilde A^{ij}\partial_j\phi\right)+\cdots \qquad (7)

where: $\tilde A^{ij}=\tilde\gamma^{ik}\tilde\gamma^{jl}\tilde A_{kl}$ ; $\cdots$ — advective and shift-derivative terms $\beta^j\partial_j\tilde\Gamma^i-\tilde\Gamma^j\partial_j\beta^i+\tfrac23\tilde\Gamma^i\partial_j\beta^j$ . Note the double spatial derivatives of the shift, which couple the gauge into the principal part. Finally the constraints, rewritten in conformal variables, that a BSSN evolution must monitor even though it does not solve them at every step.

\mathcal H=\tilde R+R^\phi+\tfrac{2}{3}K^2-\tilde A_{ij}\tilde A^{ij}=0,\qquad \mathcal M^i=\tilde D_j\!\left(e^{6\phi}\tilde A^{ij}\right)-\tfrac{2}{3}e^{6\phi}\tilde\gamma^{ij}\partial_jK=0 \qquad (8)

where: $\mathcal H$ — Hamiltonian constraint in conformal form; $\mathcal M^i$ — momentum constraint; $\tilde R,\,R^\phi$ — the conformal-metric and conformal-factor pieces of the 3-Ricci scalar; $\tilde D_j$ — covariant derivative of the conformal metric $\tilde\gamma_{ij}$ , not the physical $D_j$ of equation (6). A rising $\mathcal H$ is the first sign a run is going bad.

Derivation AWhy the contracted conformal connection must be an independent field

The goal here is to expose, term by term, exactly where ADM's weak hyperbolicity hides and how promoting $\tilde\Gamma^i$ removes it. We decompose the physical 3-Ricci tensor under the conformal rescaling into a conformal-metric part $\tilde R_{ij}$ and a conformal-factor part $R^\phi_{ij}$ , isolate the second-derivative structure of $\tilde R_{ij}$ , and show that the mixed metric derivatives which spoil the principal symbol survive only if $\tilde\Gamma^i$ is treated as a derivative of the metric — and vanish the moment it is an independent variable with its own evolution equation.

Step 1. Under $\gamma_{ij}=e^{4\phi}\tilde\gamma_{ij}$ the Ricci tensor splits cleanly, $R_{ij}=\tilde R_{ij}+R^\phi_{ij}$ , with the conformal-factor piece $R^\phi_{ij}=-2\tilde D_i\tilde D_j\phi-2\tilde\gamma_{ij}\tilde D^k\tilde D_k\phi+4\tilde D_i\phi\,\tilde D_j\phi-4\tilde\gamma_{ij}\tilde D^k\phi\,\tilde D_k\phi$ , where $\tilde D$ is the connection of $\tilde\gamma_{ij}$ . Why: the conformal-factor part contains only derivatives of the scalar $\phi$ , which are harmless — all the hyperbolicity trouble must therefore live in the tensor part $\tilde R_{ij}$ .

Step 2. Write the conformal-metric Ricci in the standard form $\tilde R_{ij}=\partial_k\tilde\Gamma^k{}_{ij}-\partial_{(i}\tilde\Gamma^k{}_{j)k}+\tilde\Gamma^k{}_{ij}\tilde\Gamma^l{}_{kl}-\tilde\Gamma^k{}_{il}\tilde\Gamma^l{}_{kj}$ . The two derivative-of-Christoffel terms are the principal part; the products of Christoffels are lower order. Why: hyperbolicity is decided entirely by the highest-derivative (principal) terms, so we only need to track $\partial_k\tilde\Gamma^k{}_{ij}$ and $\partial_{(i}\tilde\Gamma^k{}_{j)k}$ .

Step 3. Expand the first term. Using the Christoffel formula, $\partial_k\tilde\Gamma^k{}_{ij}=-\tfrac12\tilde\gamma^{kl}\partial_k\partial_l\tilde\gamma_{ij}+\tilde\gamma_{k(i}\partial_{j)}\tilde\Gamma^k+\dots$ , where the ellipsis is quadratic in first derivatives. The second spatial derivative of $\tilde\gamma_{ij}$ appears — but only as a Laplacian $\tilde\gamma^{kl}\partial_k\partial_l\tilde\gamma_{ij}$ , together with a term carrying $\partial\tilde\Gamma^k$ rather than a mixed second derivative of the metric. Why: a pure Laplacian is exactly the principal symbol of a wave operator; the mixed derivatives $\partial_i\partial_k\tilde\gamma_{jl}$ that would degrade the symbol have been folded into the single object $\partial_j\tilde\Gamma^k$ .

Step 4. Now make the choice. If $\tilde\Gamma^k$ is left as shorthand for $-\partial_l\tilde\gamma^{kl}$ , then $\partial_j\tilde\Gamma^k$ reintroduces exactly the mixed metric second derivatives, and the principal symbol degenerates back to weakly hyperbolic conformal ADM. If instead $\tilde\Gamma^k$ is an independent field evolved by equation (7), the principal part of $\tilde R_{ij}$ is just $-\tfrac12\tilde\gamma^{kl}\partial_k\partial_l\tilde\gamma_{ij}+\tilde\gamma_{k(i}\partial_{j)}\tilde\Gamma^k$ — a wave operator coupled to a first-order source. Why: this is the entire content of BSSN's stability. The reformulation is not cosmetic; the independence of $\tilde\Gamma^i$ is what makes the difference between a code that crashes at one light-crossing time and one that runs for hundreds of $M$ .

Go deeper · Derivation B · 1 minFrom the ADM metric equation to the conformal pair

Derivation BFrom the ADM metric equation to the conformal pair

The aim is to derive both the $\phi$ and $\tilde\gamma_{ij}$ evolution equations from a single ADM equation, showing that the trace and trace-free parts of $\partial_t\gamma_{ij}$ split cleanly into a volume equation and a shape equation and that the split is self-consistent with $\det\tilde\gamma=1$ . The seed is the ADM metric evolution $\partial_t\gamma_{ij}=-2\alpha K_{ij}+\mathcal L_\beta\gamma_{ij}$ , which says the intrinsic geometry of a slice changes because the slice bends into the next one (extrinsic curvature) and because the coordinates slide (the shift Lie derivative).

Step 1. Take the determinant. Using $\partial_t\ln\det\gamma=\gamma^{ij}\partial_t\gamma_{ij}$ and $\gamma^{ij}K_{ij}=K$ , the trace of the seed gives $\partial_t\ln\det\gamma=-2\alpha K+2\partial_i\beta^i+\beta^i\partial_i\ln\det\gamma$ . Why: contracting with $\gamma^{ij}$ projects out the pure volume response, and the only place the extrinsic curvature enters is through its trace $K$ — precisely the volume-changing part.

Step 2. Insert $\phi=\tfrac{1}{12}\ln\det\gamma$ and divide by twelve: $\partial_t\phi=-\tfrac16\alpha K+\tfrac16\partial_i\beta^i+\beta^i\partial_i\phi$ , which is equation (4). Why: because $\det\tilde\gamma=1$ is fixed once and for all, $\phi$ is the sole carrier of volume, so its evolution is nothing but the rescaled trace equation.

Step 3. Differentiate $\tilde\gamma_{ij}=e^{-4\phi}\gamma_{ij}$ to get $\partial_t\tilde\gamma_{ij}=e^{-4\phi}\partial_t\gamma_{ij}-4\tilde\gamma_{ij}\partial_t\phi$ . Why: this subtracts the volume evolution we just found from the full metric evolution, leaving only the shape-changing remainder.

Step 4. Substitute the seed equation and the trace split $K_{ij}=e^{4\phi}(\tilde A_{ij}+\tfrac13\tilde\gamma_{ij}K)$ . The trace parts cancel by construction, leaving $\partial_t\tilde\gamma_{ij}=-2\alpha\tilde A_{ij}+\beta^k\partial_k\tilde\gamma_{ij}+2\tilde\gamma_{k(i}\partial_{j)}\beta^k-\tfrac23\tilde\gamma_{ij}\partial_k\beta^k$ , which is equation (5). Why: the result is manifestly trace-free with respect to $\tilde\gamma_{ij}$ , so it preserves the unit determinant to leading order — the conformal and traceless splits are mutually consistent, and only the physical shear $\tilde A_{ij}$ drives the shape. *The trace-free structure of equation (5) guarantees $\partial_t\det\tilde\gamma=0$ analytically.*

Go deeper · Derivation C · 2 minMoving-puncture gauge — 1+log slicing and the Gamma-driver shift

Derivation CMoving-puncture gauge — 1+log slicing and the Gamma-driver shift

This derivation motivates the live gauge conditions that make BSSN usable for black holes without excision, and explains why their specific form is dictated by hyperbolicity rather than convenience. Because hyperbolicity is a property of the coupled metric-plus-gauge system, the lapse and shift cannot be prescribed freely; they must be evolved by equations built from the very BSSN fields ( $K$ and $\tilde\Gamma^i$ ) so that the combined principal symbol stays well posed and the coordinate singularity at a puncture is handled automatically.

Step 1. Impose $1+\log$ slicing, $\partial_t\alpha=-2\alpha K+\beta^i\partial_i\alpha$ , a member of the Bona–Masso family with $f(\alpha)=2/\alpha$ . Near a puncture $K$ grows, driving $\alpha\to0$ , so the lapse collapses and the evolution freezes exactly where the singularity forms. Why: this is singularity avoidance without excision — the slice simply refuses to advance into the strong-field region, and it feeds on $K$ , a field BSSN already evolves.

Step 2. Drive the shift by the conformal connection, $\partial_t\beta^i=\tfrac34 B^i+\beta^j\partial_j\beta^i$ with $\partial_t B^i=\partial_t\tilde\Gamma^i-\eta B^i+\dots$ . The auxiliary field $B^i$ makes the shift a damped chase of $\tilde\Gamma^i$ , counteracting the coordinate distortion (slice-stretching) that would otherwise pile up grid points behind an orbiting hole. Why: $\tilde\Gamma^i$ measures how non-uniform the conformal coordinates are, so steering the shift by it keeps the grid sensibly distributed; the damping $\eta$ , of order $1/M$ to $2/M$ , prevents the chase from overshooting into oscillation.

Step 3. Check that the coupling preserves hyperbolicity. The double shift derivatives in equation (7) and the gauge equations share the same principal structure as the metric sector, so the full metric-plus-gauge system is strongly hyperbolic with real characteristic speeds. In the far field $\alpha\to1$ and the $1+\log$ gauge mode propagates at $\sqrt{f}=\sqrt{2}\approx1.41$ times the speed of light — a harmless coordinate speed, not a physical one, exactly analogous to a superluminal phase velocity. Why: the gauge is not an afterthought bolted onto BSSN but part of the well-posed system; both gauge variables are BSSN fields, so the punctures advect across the grid as the black holes physically orbit. *The coefficients $2$ ( $1+\log$ ) and $3/4$ (Gamma-driver) and the damping $\eta\sim1/M$ are the standard production values. The far-field gauge speed $\sqrt2$ follows directly from $f(\alpha\!\to\!1)=2$ .*

Limiting cases

Case	Result	Meaning
Conformally flat, momentarily static	$\tilde\gamma_{ij}=\delta_{ij},\ K_{ij}=0\ \Rightarrow\ \tilde A_{ij}=0,\ K=0,\ \tilde\Gamma^i=0$	Brill–Lindquist and Bowen–York puncture data reduce BSSN at the initial instant to evolving $\phi$ alone; the entire traceless apparatus starts trivially.
Drop the $\tilde\Gamma^i$ promotion	$\tilde\Gamma^i\to-\partial_j\tilde\gamma^{ij}$ (not evolved)	Substituting the definition back everywhere returns conformal ADM, which is only weakly hyperbolic — the stability gain is lost, proving the independence is essential, not cosmetic.
Linearized weak field	$\gamma_{ij}=\delta_{ij}+h_{ij},\ \ \Box\,h^{\rm TT}_{ij}=0$	Around flat space BSSN reduces to wave equations for the transverse-traceless perturbation, reproducing gravitational-wave propagation at the speed of light.
Maximal slicing limit	$K=0,\ \partial_tK=0$	Freezing the trace sector leaves $\phi,\ \tilde\gamma_{ij},\ \tilde A_{ij},\ \tilde\Gamma^i$ evolving with $K$ held to zero — a classic singularity-avoiding gauge, decoupled from the $1+\log$ lapse.

Common trapBSSN is automatically stable once you rewrite the variables

It is not. Rewriting ADM in conformal-traceless variables makes the system strongly hyperbolic, which is necessary for well-posedness, but it does not make it constraint-damped. Two things must be done by hand every few steps or the evolution still drifts and dies. First, the algebraic conditions $\det\tilde\gamma=1$ and $\tilde\gamma^{ij}\tilde A_{ij}=0$ must be actively re-imposed — codes rescale $\tilde\gamma_{ij}$ back to unit determinant and subtract the trace from $\tilde A_{ij}$ — because truncation error makes both drift, and once they drift the conformal split loses meaning. Second, the momentum constraint must be substituted into the $\partial_t\tilde\Gamma^i$ equation rather than using the raw divergence $\partial_j\tilde A^{ij}$ ; skip it and the very mode you promoted goes unstable. And hyperbolicity is gauge-dependent: an arbitrary shift can spoil it, which is why the specific Gamma-driver condition is part of the recipe, not an optional garnish.

Worked examples

1 · Counting the BSSN state vector (fully worked)

How many fields does a moving-puncture BSSN code evolve at each grid point, how does that compare with ADM, and how many physical degrees of freedom survive?

Step 1. Count the ADM variables. ADM evolves the symmetric spatial metric $\gamma_{ij}$ (6 components) and the symmetric extrinsic curvature $K_{ij}$ (6 components), for 12 evolved fields. All 12 are stored and stepped.

Step 2. Count the BSSN variables. The conformal factor $\phi$ (1), the conformal metric $\tilde\gamma_{ij}$ (6), the trace $K$ (1), the conformal traceless curvature $\tilde A_{ij}$ (6), and the conformal connection functions $\tilde\Gamma^i$ (3) give $1+6+1+6+3=\mathbf{17}$ evolved fields. BSSN deliberately stores more numbers than ADM — the extra bookkeeping is the price of strong hyperbolicity.

Step 3. Add the live gauge. Moving-puncture runs also evolve the lapse $\alpha$ (1), the shift $\beta^i$ (3), and the Gamma-driver auxiliary $B^i$ (3), i.e. 7 gauge fields, for a total of $17+7=\mathbf{24}$ evolved quantities per grid point.

Step 4. Recover the physical degrees of freedom. Two algebraic conditions, $\det\tilde\gamma=1$ and $\tilde\gamma^{ij}\tilde A_{ij}=0$ , are enforced (removing 2), and the four Hamiltonian-plus-momentum constraints together with the four-fold coordinate freedom remove the rest, leaving the same 2 physical polarizations of the gravitational field that ADM describes. The extra fields are redundancy in the service of stability, not new physics. ADM 12 → BSSN 17 (+7 gauge = 24), yet still 2 physical polarizations.

2 · The gauge speed of 1+log slicing

The $1+\log$ condition is $\partial_t\alpha=-2\alpha K$ . At what coordinate speed do gauge pulses propagate in the asymptotically flat region, and does a value above $c$ signal a problem?

Step 1. Identify the slicing family. Bona–Masso slicing is $\partial_t\alpha=-\alpha^2 f(\alpha)K$ ; matching to $\partial_t\alpha=-2\alpha K$ gives $\alpha^2 f=2\alpha$ , hence $f(\alpha)=2/\alpha$ .

Step 2. Recall the characteristic speed. For Bona–Masso slicing the gauge (lapse) modes propagate, relative to the normal observers, at $c_{\rm gauge}=\sqrt{f}$ in units of the physical light speed.

Step 3. Evaluate in the far field. As $r\to\infty$ , $\alpha\to1$ , so $f\to2$ and $c_{\rm gauge}=\sqrt{2}\approx\mathbf{1.41}\,c$ . The $1+\log$ gauge pulse outruns light.

Step 4. Interpret. This is a coordinate speed, not a signal speed — the same status as a phase velocity above $c$ in optics. It carries no information and violates no causality; the physical light cones are untouched, and the metric and matter modes still propagate at $\le c$ . In fact a gauge speed slightly above the light cone is desirable, because it lets coordinate transients escape the domain slightly ahead of the physical fields. The far-field 1+log gauge speed is √2 ≈ 1.41 c, a harmless coordinate effect.

Physics of the portrait

The BSSN formulation is where three decades of frustrated numerical relativity finally cohered — a Japanese conformal reformulation, a pair of American authors who proved its stability and gave it its name, and a gauge invented on two continents at once that let black holes move.

Origins · 1 minBSSN formulation and how this was discovered

BSSN formulation

the conformal-traceless 3+1 evolution system γ_ij = e^{4φ} γ̃_ij · numerical relativity of strong-field spacetimes

BSSN is one of the rare physics results named for four people who never all sat in a room together. The conformal-traceless idea first appeared in 1987, when Takashi Nakamura, Ken-ichi Oohara and Yasufumi Kojima, working in the Japanese numerical-relativity tradition, rewrote the ADM equations in conformal variables to tame gravitational collapse. Masaru Shibata and Nakamura sharpened it in 1995, introducing the contracted conformal connection functions as evolved fields — the move that quietly fixed the hyperbolicity, though its structural reason was not yet understood. The formulation might have stayed a specialist's tool had Thomas Baumgarte and Stuart Shapiro not, in 1999, run ADM and the Nakamura–Shibata system head to head on the same data and shown the latter stable where the former exploded; their paper crystallized why, and the community soldered the four surnames together as "BSSN". The last piece came in 2005–2006. Frans Pretorius had just cracked the binary-black-hole problem with a rival generalized-harmonic code, and within months two groups — Manuela Campanelli, Carlos Lousto and collaborators, and John Baker, Joan Centrella and collaborators — independently found that BSSN plus a moving-puncture gauge could do the same, letting the punctures simply advect across the grid. Together those results ended forty years of failure and gave gravitational-wave astronomy its waveform templates: the reformulation that changed no physics changed everything about whether we could compute it.

History

1962

— Arnowitt, Deser and Misner cast general relativity as a 3+1 initial-value problem, later found only weakly hyperbolic

1987

— Nakamura, Oohara and Kojima introduce a conformal-traceless reformulation for collapse

1995

— Shibata and Nakamura add the conformal connection functions as evolved variables

1999

— Baumgarte and Shapiro demonstrate the system's stability and fix its name as BSSN

2006

— Campanelli et al. and Baker et al. add the moving-puncture gauge and simulate binary-black-hole mergers

Going further

Texts. Baumgarte and Shapiro, Numerical Relativity: Solving Einstein's Equations on the Computer (2010), Chapter 11 · Alcubierre, Introduction to 3+1 Numerical Relativity (2008), Chapter 5 · Gourgoulhon, 3+1 Formalism in General Relativity (2012) · Wald, General Relativity (1984), Chapter 10 (initial-value formulation) · Misner, Thorne and Wheeler, Gravitation (1973), Chapter 21 (ADM).

Papers. Nakamura, Oohara and Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987) · Shibata and Nakamura, Phys. Rev. D 52, 5428 (1995) · Baumgarte and Shapiro, Phys. Rev. D 59, 024007 (1999) · Campanelli, Lousto, Marronetti and Zlochower, Phys. Rev. Lett. 96, 111101 (2006) · Baker, Centrella, Choi, Koppitz and van Meter, Phys. Rev. Lett. 96, 111102 (2006). Key self-consistency checks: $\tilde\Gamma^i=-\partial_j\tilde\gamma^{ij}$ follows from the vanishing trace of unit-determinant Christoffels, and $\phi=\tfrac1{12}\ln\det\gamma$ from $\det\gamma=e^{12\phi}$ ; equations (4) and (5) follow from the ADM $\partial_t\gamma_{ij}$ by trace and trace-free splitting; the variable counts run 12 → 17 (+7 gauge); the far-field $1+\log$ gauge speed is $\sqrt f=\sqrt2$ with $f=2/\alpha$ ; and the gauge coefficients are $2$ ( $1+\log$ ), $3/4$ (Gamma-driver), and $\eta\sim1/M$ .

Prerequisites. ADM 3+1 Decomposition of Spacetime · Gauss–Codazzi–Mainardi Equations · Hamiltonian and Momentum Constraints · Initial-Value (Cauchy) Problem in General Relativity. Leads to. Gauge Conditions in Numerical Relativity · Black-Hole Excision and Puncture Methods · Numerical Relativity and Binary Black-Hole Mergers · Apparent Horizons and Horizon Finders.

Ask the tutor: Ask me to walk through why the promoted $\tilde\Gamma^i$ turns the conformal Ricci tensor into a wave operator, or how the moving-puncture gauge lets a black hole slide across the grid without excision.

Why was ADM unstable if its equations are exactly Einstein's?

Because stability is a property of the formulation, not just the physics. ADM writes the correct evolution equations, but as a partial-differential system they are only weakly hyperbolic: the principal symbol is not diagonalizable, so it cannot bound the growth of every mode. Some gauge and constraint combinations neither propagate away nor decay, and the nonlinear terms then amplify any truncation error seeded into them. The result is that even perfect initial data develops exponentially growing violations and the code crashes within about one light-crossing time. BSSN keeps the identical physical content but reshuffles it into a strongly hyperbolic system whose modes all propagate with real speeds, so errors leave the grid instead of blowing up.

What exactly does promoting Γ̃ⁱ to an independent variable mean in code?

It means $\tilde\Gamma^i$ gets its own storage array and its own time-integration, updated by equation (7), instead of being recomputed from spatial derivatives of $\tilde\gamma^{ij}$ at every step. Analytically the two are equal, $\tilde\Gamma^i=-\partial_j\tilde\gamma^{ij}$ , but numerically they are different objects: the independently evolved field appears in the Ricci tensor as a first derivative $\partial_j\tilde\Gamma^k$ , whereas the recomputed version would inject mixed second derivatives of the metric. Keeping them separate is what makes the principal part wave-like. The definition is not discarded, though — it becomes a constraint, monitored during the run, and codes sometimes gently drive $\tilde\Gamma^i$ back toward $-\partial_j\tilde\gamma^{ij}$ to control drift.

Why the conformal factor χ = e^{-4φ} instead of φ near a puncture?

At a puncture the physical metric diverges, so $\phi\to+\infty$ and $e^{4\phi}$ blows up — a bad thing to store and difference numerically. The variable $\chi=e^{-4\phi}$ (or the related $W=e^{-2\phi}$ ) instead goes smoothly to zero at the puncture while staying finite and differentiable everywhere, so no grid point holds an infinite number and no region has to be carved out (excised). Evolving $\chi$ rather than $\phi$ was one of the practical keys to the 2006 moving-puncture breakthrough: it let the coordinate singularity sit on the grid harmlessly. The trade-off is that $\chi$ appears in denominators in a few places, so codes floor it away from exactly zero.

How is BSSN different from the generalized-harmonic formulation?

They are the two independent solutions to the same well-posedness problem, and both work. The generalized-harmonic approach, used by Pretorius in the 2005 breakthrough, imposes harmonic-type coordinate conditions that make Einstein's equations manifestly wave-like from the start and typically pairs with black-hole excision. BSSN keeps the 3+1 split explicit, achieves hyperbolicity through the conformal-traceless variables and the promoted $\tilde\Gamma^i$ , and pairs with the moving-puncture gauge instead of excision. In practice generalized-harmonic codes (like SpEC) often use spectral methods and excise, while BSSN codes use finite differencing with adaptive mesh refinement and punctures. That the two very different formulations agree on the resulting waveforms to high precision is a strong cross-check on numerical relativity as a whole.

What is slice-stretching and how does the Gamma-driver shift fight it?

Slice-stretching is a coordinate pathology: with singularity-avoiding slicing the lapse collapses near a black hole, so coordinate time nearly stops there while it runs normally far away, and the spatial coordinates get dragged and piled up steeply near the horizon. Gradients sharpen until the finite-difference grid can no longer resolve them and the run fails. The Gamma-driver shift counteracts this by steering the coordinates outward: it chases $\tilde\Gamma^i$ , which measures exactly how non-uniform the conformal coordinates have become, redistributing grid points so they stay usable. Without a good shift condition, even a strongly hyperbolic BSSN evolution of a black hole eventually chokes on its own coordinates.

Is BSSN still state of the art, or has it been superseded?

Both, depending on what you count. Plain BSSN is superseded for the most demanding runs by its constraint-damped descendants Z4c and CCZ4, which add terms that actively drive constraint violations toward zero rather than merely letting them propagate — a real robustness gain for long neutron-star inspirals and high-spin black holes. But those formulations are direct evolutions of the BSSN idea, sharing its conformal-traceless variables and its gauge, so in the sense that matters BSSN's framework is very much alive. It remains the default in the widely used Einstein Toolkit and the conceptual entry point everyone learns first. For the majority of production binary-black-hole simulations, moving-puncture BSSN is still entirely adequate.

Where do matter sources enter the BSSN equations?

Through the right-hand sides, as projections of the stress-energy tensor onto the slice. The energy density $\rho$ measured by normal observers sources the trace equation for $K$ and the Hamiltonian constraint; the momentum density $S^i$ sources the momentum constraint and the $\tilde\Gamma^i$ equation; and the spatial stress $S_{ij}$ (with its trace $S$ ) sources the traceless curvature $\tilde A_{ij}$ and, through $S-3\rho$ -type combinations, the trace $K$ . In vacuum all of these vanish and one recovers the equations quoted above. Coupling BSSN to relativistic hydrodynamics or magnetohydrodynamics just means computing these projections from the fluid variables each step, which is exactly how binary-neutron-star and stellar-collapse simulations are built.