Stationary solutions in AdS

Stationary solutions in AdS [Incomplete]

by renphysics (contact: renphysics@adscft.org)

$$\boxed{R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R+\Lambda g_{\mu\nu}=8\pi GT_{\mu\nu}}$$

Aspect	Numerical relativity for astrophysics	Numerical relativity for AdS
Physical motivation	Simulate extreme gravitational events such as black holes or neutron stars merging. To understand general relativity including gravitation waves	Better understand strongly coupled systems. By simulating gravitational dynamics in AdS, one can model analogous processes like thermalization in the dual quantum system, probing strongly coupled matter far from equilibrium
Global structure	Globally hyperbolic; Cauchy data on spatial slice	Not globally hyperbolic; must supply boundary data on the timelike AdS boundary (where the CFT lives)
Stationary solutions	Limited due to uniqueness theorems - topology: spherical - rigidity - no hair theorem	Rich landscape: black droplets, black funnels
Fundations	1990s	2009-2013
Examples	Binary black holes mergers	black droplets, black funnels, black
Gauge choices	BSSN, generalized harmonic, etc.	DeTurck gauge Chesler-Yaffe gauge Bondi-Sachs gauge

Comercial coftware: Comsol. Standard equations, highly complex boundaries. Finite element method NR for AdS: non-Standard equations, simple boundaries, pseudospectral method

a strongly coupled quantum field theory can be mapped to a classical gravity system in a higher-dimensional AdS universe.

1. Exact solutions of Einstein’s equations

Symmetries

A basic and practical classification of gravitational solutions is by their symmetries, including the group and the orbit. Namely: when we solve Einstein equations, we can ask how many Killing vectors are there? A well-known case is the stationary and axisymmetric case, where there are two Killing vectors. Since Einstein’s equations have diffeomorphism invariance, which implies that different people may write the same geometry in different coordinate systems. So, to classify solutions, we need to use quantities independent of coordinates.

• the isometry group $G$ (continuous + discrete symmetries), and
• the structure of its orbits.

Besides symmetry, there are other important ways to classify spacetimes:

• Petrov classification (algebraic type of the Weyl tensor);
• classification by the Ricci tensor (e.g. Einstein manifolds, Ricci solitons);
• classification by the energy–momentum tensor (matter content, equation of state, conserved currents);
• horizon topology and global structure (especially in $D\ge 5$).

A particularly useful invariant is the cohomogeneity:

• If the metric depends nontrivially on $n$ coordinates after quotienting by symmetries, we call it cohomogeneity-$n$.
• Cohomogeneity-1 problems typically reduce to ODEs (often analytically solvable or numerically straightforward).
• Cohomogeneity-2 and higher generically require solving PDEs, and this is where the methods on this page become essential.

In asymptotically flat 4D vacuum gravity, stationary solutions are heavily constrained (e.g. Kerr uniqueness). In AdS, or with nontrivial boundary data/matter content, these uniqueness theorems typically fail, and the space of stationary solutions becomes vastly richer.

Let us start from the maximally symmetric case (10 Killing vectors): Minkowski, de Sitter, and anti–de Sitter. These are very simple solutions, and we are certainly not satisfied with only these.

When there is a black hole, we have fewer isometries. In some highly symmetric cases (for example, with enough Killing vectors), one can add a cosmological constant and still have an analytic solution. Usually the analytic solutions only depend on the AdS radial coordinate. But when there are only two Killing vectors, if we add a cosmological constant, generally there will be no analytic solution—except in special algebraically special cases, for example Petrov type D (in four dimensions). As a comparison, in the asymptotically flat vacuum case, Einstein’s equations for stationary and axisymmetric spacetimes are integrable. In the AdS case, we typically need to solve Einstein equations numerically.

境自远尘皆入咏，物含妙理总堪寻。 (Far beyond the world’s dust, the cosmic vista itself becomes poetry; within each thing lies a subtle principle, one that merits our profound inquiry.)

形骸已与流年老，诗句犹争造化功。 (My body creaks under the weight of passing years, My poems aim still to rival the perfections of nature. By Lu You, translated by C.N Yang.)

An analytic solution is like a Chinese traditional poem, and a numerical solution is like a modern poem. —Jie Ren

2. Stationary solutions in AdS spacetimes

In AdS/CFT, a standard application is that we can use a classical theory of gravity to study a strongly coupled quantum field theory. In this context, the motivation to study numerical relativity in AdS spacetime is to construct AdS spacetimes with less symmetry, so that we can study more general states in the CFT. The solutions we most commonly use are the Schwarzschild–AdS solution and the Reissner–Nordström–AdS solution. These solutions correspond to particular states in the CFT. If we want to study more general states, typically we need to find numerical solutions.

Historically, when people talked about numerical relativity (especially in the last century), they mostly meant time-evolution problems—like binary black hole mergers—because in asymptotically flat spacetimes, there are many “black hole theorems” (uniqueness/no-hair theorems, topology theorems, rigidity theorems, etc.) that strongly constrain stationary solutions. Roughly speaking, in 4D asymptotically flat Einstein–Maxwell, the stationary black hole solutions are Kerr or Kerr–Newman.

But in AdS spacetimes, black hole solutions can evade some of these theorems (or their assumptions), leading to a significantly richer variety of solutions.

An analytic solution: topological Schwarzschild–AdS

Let me first give a very simple example related to the topology theorem. In asymptotically flat spacetimes, under mild assumptions, the horizon of a black hole must be compact, and in 4D its topology must be a sphere. But in AdS we can certainly have planar and hyperbolic black holes.

For example, in $d+1$ dimensions one can write the “topological” AdS black hole in the standard form

$$ds^2=-f(r)\,dt^2+\frac{dr^2}{f(r)}+r^2 d\Sigma_{k,D-2}^2,$$

with

$$f(r)=k-\frac{\mu}{r^{D-3}}+\frac{r^2}{L^2},$$

where $k\in\{+1,0,-1\}$, and $d\Sigma_{k,D-2}^2$ is the unit metric on $S^{D-2}$ ($k=1$), $\mathbb{R}^{D-2}$ ($k=0$), or $\mathbb{H}^{D-2}$ ($k=-1$).

where $d\Sigma_{k,d-1}^2$ is the metric on a space of constant curvature $k$, and $k$ can take values $1,0,-1$ (spherical, planar, hyperbolic horizons).

If the cosmological constant is zero or positive, the “black hole” interpretation is essentially only for the spherical case. But with negative cosmological constant (AdS), we have planar and hyperbolic black holes.

For $D=4$ one often writes

$$f(r)=k-\frac{2M}{r}-\frac{\Lambda}{3}r^2,$$

which matches $\mu=2M$ and $L^2=-3/\Lambda$.

Why this matters for numerics: this family provides canonical examples of

• asymptotically AdS boundaries,
• non-extremal Killing horizons,
• different horizon topologies,
• coordinate choices that make regularity manifest.

If $\Lambda=0$, there can be black holes with different topologies in $D\ge 5$. For example, $S^3$ and $S^1\times S^2$. Even for pure gravity, higher dimensions allow richer horizon topologies (e.g. black rings with $S^1\times S^2$ in asymptotically flat 5D, and many AdS analogues). In AdS, the presence of the boundary and the cosmological constant further enriches possibilities: one can engineer boundary geometries that “force” new bulk horizon configurations.

---

More on hyperbolic black holes

Hyperbolic AdS black holes have several features that often show up in computations and in AdS/CFT applications:

• Causal structure changes with temperature.
  At high temperature, the causal structure resembles Schwarzschild–AdS. At sufficiently low temperature (even without charge), the causal structure resembles Reissner–Nordström, with two horizons. When $k=0$ and $\mu=0$, it is called the Rindler-AdS coordinates.

• Entanglement and Rényi entropies.
  Hyperbolic black holes are central in the standard mapping between Rényi entropies of a CFT and thermal physics on $\mathbb{H}^{d-1}$.

• Horizon ntersections with the boundary.
  In some coordinate choices, hyperbolic horizons can intersect the conformal boundary. This can turn a complicated domain with several curved boundaries into a rectangle in compact coordinates—ideal for tensor-product spectral grids.

---

In AdS/CFT, a stationary bulk geometry is not merely “a nice solution of Einstein’s equation”: it is often a dual description of a stationary state of a strongly coupled QFT on a prescribed background spacetime.

Two facts make AdS a natural “numerical laboratory”:

1. Boundary conditions are physical data.
   Asymptotically AdS spacetimes allow you to specify a boundary metric (a conformal class), which in AdS/CFT is the spacetime where the QFT lives. This strongly enlarges the solution space.

2. The radial direction is not the only variable.
   Once the boundary geometry or sources depend on spatial coordinates, bulk fields depend both on the holographic direction and boundary directions, turning the problem into genuinely multi-dimensional PDEs.

A convenient AdS parameterization is

$$\Lambda=-\frac{(D-1)(D-2)}{2L^2},$$

so the vacuum equation becomes

$$R_{\mu\nu}+\frac{D-1}{L^2}g_{\mu\nu}=0.$$

Now, how can uniqueness be violated more generally?

1. Add matter fields. A famous example is holographic superconductors: these are black holes that can spontaneously develop scalar hair (or vector hair). In holographic QCD, bulk scalar fields are also natural to include.

2. Even in pure gravity with $\Lambda<0$, there is still a larger variety of stationary solutions, for several reasons:

   • The horizon topology is not necessarily a sphere; it can be planar, hyperbolic, or even disconnected.
   • The rigidity theorem can be violated: the horizon need not be a Killing horizon.
   • As a result, the surface gravity is not necessarily constant: there can be heat flow along a spatial direction. In that case one expects cross terms like $dt\,dx$ in the metric.

As a comparison: in the Kerr solution there is a cross term $dt\,d\phi$, but both $t$ and $\phi$ are Killing coordinates. Here, one can have cross terms between a Killing coordinate and a non-Killing coordinate.

3. The AdS boundary can be conformal to (almost) any spacetime. For example, if we want to study $\mathcal{N}=4$ super-Yang–Mills theory on a black hole background, we can impose boundary conditions such that the induced metric on the AdS boundary is conformal to (say) a Schwarzschild black hole metric. In this case one calls the bulk geometry asymptotically locally AdS.

Also, due to the violation of uniqueness theorems in AdS, we have a larger variety of stationary solutions. Actually, even when the cosmological constant is zero, in higher dimensions one can have other horizon topologies like black rings. In AdS there are many possibilities: black holes, black strings, black rings, and even black bottles and black globules. I will give some examples.

Beyond analytic solutions: boundary black holes, funnels, droplets

A key AdS/CFT application is to place the boundary CFT on a black hole spacetime. The boundary horizon then extends into the bulk in ways that allow qualitatively different bulk horizon geometries, commonly described as:

• black funnels: a connected horizon reaching the boundary;
• black droplets: a horizon hanging from the boundary but disconnected from a separate “planar” bulk horizon.

These are cohomogeneity-2 (or higher) problems, and are typically accessible only numerically.

From the numerical point of view, the moral is:

Choose coordinates so that the computational domain is as simple as possible (ideally a rectangle or a product of intervals), and factor out known singular/vanishing behavior analytically so that the remaining unknown functions are smooth.

3. Ellipticity and DeTurck gauge fixing

3.1 Character of Einstein’s equations

When we try to solve PDEs, we first want to ask: what type of PDE is it? Elliptic, hyperbolic, or parabolic?

These different classes have quite different mathematical and numerical features.

• For elliptic systems, after discretization the differential operator becomes a matrix, and solving the PDE is essentially inverting this operator (with appropriate boundary conditions).
• For hyperbolic systems, the operator is not invertible in that sense; one typically has a time-evolution problem and must use evolution schemes.
• Parabolic problems (diffusion-type) again have different features.

Now, Einstein’s equation is difficult in part because it does not have a definite character “by itself”. Let me explain why, and how to deal with it.

People often say Einstein’s equation is “highly nonlinear”. In mathematics language, it is a quasi-linear second-order system: the highest-derivative part is linear, so in some sense it is controllable. But it does not have a definite character because of gauge invariance.

To determine the character, we linearize Einstein’s equation around a background, $g_{\mu\nu}\to g_{\mu\nu}+h_{\mu\nu}$, and look only at the highest-derivative terms (the second-derivative part), which is called the principal symbol.

Einstein’s equation is a quasi-linear second-order PDE system for $g_{\mu\nu}$. A robust way to understand whether a boundary value problem is well-posed is to study the principal symbol of the operator.

Take a perturbation

$$g_{\mu\nu}\to g_{\mu\nu}+h_{\mu\nu}.$$

The principal symbol of the Ricci tensor acting on $h_{\mu\nu}$ can be written (schematically) as

$$(P_g h)_{\mu\nu} =\frac12 g^{\alpha\beta} \bigl( -2\,\partial_\alpha\partial_\beta h_{\mu\nu} +2\,\partial_\alpha\partial_{(\mu}h_{\nu)\beta} -\partial_\mu\partial_\nu h_{\alpha\beta} \bigr).$$

In Fourier space, replace $\partial_\alpha\to \xi_\alpha$ to get

$$(P_g(\xi)h)_{\mu\nu} =\frac12\bigl( -\xi^2 h_{\mu\nu} +2\,\xi^\alpha\xi_{(\mu}h_{\nu)\alpha} -\xi_\mu\xi_\nu h \bigr),$$

where $h\equiv g^{\alpha\beta}h_{\alpha\beta}$, $\xi^2\equiv g^{\alpha\beta}\xi_\alpha\xi_\beta$. For an elliptic system, one wants $P_g(\xi)$ to be invertible for all nonzero covector $\xi$, i.e.,

$$P_g(\xi)h=0 \quad \Rightarrow\quad h=0.$$

But Einstein’s equation fails this test because of diffeomorphism invariance. Infinitesimally, a diffeomorphism generated by a vector field $v_\mu$ changes the metric by a pure gauge perturbation

$$h_{\mu\nu}=\nabla_{(\mu}v_{\nu)},$$

which always satisfies $(P_g(\xi)h)_{\mu\nu}=0$. In other words, pure gauge perturbations lie in the kernel of the principal symbol, so the system is not elliptic as written.

This is the stationary analog of a familiar counting statement: the metric has $D(D+1)/2$ components, but Einstein’s equations do not give that many independent equations without gauge fixing, because of the Bianchi identity. In 4D one often summarizes it as “10 components = 6 independent Einstein equations + 4 gauge conditions.”

We can decompose

$$h_{\mu\nu}=\hat h_{\mu\nu}+\nabla_{(\mu}v_{\nu)},$$

where $\hat h_{\mu\nu}$ is the physical (transverse) part and $\nabla_{(\mu}v_{\nu)}$ is pure gauge (longitudinal).

So: ellipticity requires gauge fixing: it is what makes the PDE problem well-posed.

3.2 Harmonic gauge and the DeTurck vector

At the linearized level, harmonic (de Donder) gauge is

$$\partial_\nu \hat h^\nu_\mu-\frac12\partial_\mu \hat h=0.$$

A geometrically clean way to implement harmonic gauge for the full nonlinear problem is to introduce a covector field $\xi_\mu$ built from the connection. One convenient form is

$$\xi^\mu = g^{\alpha\beta}\Gamma^{\mu}_{\alpha\beta} = -\nabla^2 x^\mu,$$

which is the standard harmonic-coordinate condition. It can be written as

$$C_\mu \equiv \nabla^\nu h_{\mu\nu}-\frac{1}{2}\nabla_\mu h = 0, \qquad h\equiv g^{\mu\nu}h_{\mu\nu}.$$

Equivalently, in terms of coordinates one can view harmonic gauge as a condition that makes coordinates satisfy wave/Laplace-type equations.

The most flexible and globally well-defined formulation used in stationary numerical relativity is the Einstein–DeTurck method, which introduces the DeTurck vector (“generalized harmonic” vector)

$$\xi^\mu = g^{\alpha\beta} \big( \Gamma^{\mu}_{\alpha\beta}(g) - \bar\Gamma^{\mu}_{\alpha\beta}(\bar g) \big),$$

where $\bar g_{\mu\nu}$ is a fixed reference metric with the same boundary structure (and typically the same “qualitative” causal/regularity structure) as the desired solution. Because the difference of two Levi-Civita connections is a tensor, $\xi^\mu$ is a globally defined vector field. Equivalently, one can write $\xi_\mu=g_{\mu\nu}\xi^\nu$; in coordinates it reduces to a familiar “harmonic coordinate” expression in the special case $\bar g$ is flat.

Then define the DeTurck-modified Ricci tensor

$$R^H_{\mu\nu}=R_{\mu\nu}-\nabla_{(\mu}\xi_{\nu)}.$$

The Einstein–DeTurck equation is

$$R^H_{\mu\nu}-\frac{2\Lambda}{D-2}g_{\mu\nu}=0, \qquad\text{(vacuum)}$$

or in AdS notation,

$$R^H_{\mu\nu}+\frac{D-1}{L^2}g_{\mu\nu}=0.$$

A key fact: this system has the same number of independent equations as unknown metric components (schematically $\frac12 D(D+1)$), because the DeTurck term supplies the missing gauge conditions. This is summarized by counting equations vs unknowns

• $g_{\mu\nu}$: $\frac12 D(D+1)$ unknowns,
• $R_{\mu\nu}+\frac{D-1}{L^2}g_{\mu\nu}=0$: $\frac12 D(D-1)$ independent equations (Bianchi reduces by $D$),
• $R^H_{\mu\nu}+\frac{D-1}{L^2}g_{\mu\nu}=0$: $\frac12 D(D+1)$ independent equations.

In $D=4$, for example, $g_{\mu\nu}$ has 10 components, while the Einstein equation provides only 6 independent equations until you fix 4 coordinate/gauge conditions. The DeTurck formulation packages this gauge fixing into the PDE system itself.

---

Counting independent equations

The statement $10=6+4$ corresponds to the most general situation where all components of the metric are nonzero and the metric depends on all coordinates. Now consider special cases. For instance, suppose the metric depends only on one coordinate — say a single “holographic” coordinate — or perhaps it depends on a holographic coordinate plus one spatial coordinate. Then the counting changes: $3=2+1$.

If we substitute the ansatz

$$ds^2=-A(r)dt^2+B(r)dr^2+C(r)d\Omega^2$$

into Einstein’s equations, we will not get a fully determined set unless you fix one function by a gauge choice, because there is a coordinate freedom $r\to f(r)$. We can choose $C(r)=r^2$ or $B(r)=1/A(r)$ to fix the gauge. By the DeTurck method, we must start with the above ansatz.

The number of unknown functions equals “number of independent Einstein equations” plus “number of gauge conditions,” but the numbers depend on the cohomogeneity of the ansatz.

How do we know it is two? One way is brute force: actually try it yourself. Plug different ansätze into the Einstein equations and see what equations you get. You will likely find three equations that look different. But seeing that only two of them are independent is not so obvious. Sometimes you need a trick — using the (contracted) Bianchi identity:

$$\nabla_\mu G^{\mu\nu} = 0,$$

to understand the dependence between equations.

---

$$R_{\mu\nu}\to R^H_{\mu\nu}=R_{\mu\nu}-\nabla_{(\mu}\xi_{\nu)}.$$

The payoff is visible at the level of the principal symbol:

$$(P_g^H h)_{\mu\nu}=-\frac12 g^{\alpha\beta}\partial_\alpha\partial_\beta h_{\mu\nu},$$

which is Laplace-like and therefore elliptic in a Riemannian setting (and becomes elliptic in stationary Lorentzian reductions, as explained next).

Finally, after solving the Einstein–DeTurck equation, one must check the gauge condition

$$\xi^\mu=0.$$

If $\xi^\mu=0$, the solution is an actual Einstein solution. If $\xi^\mu\neq 0$, the solution is a Ricci soliton (a solution of the DeTurck-modified system but not of Einstein’s equations). The possibility of solitons is real and must be controlled.

These are exactly the conceptual points emphasized in the lectures: you solve the modified system, then verify it is truly an Einstein solution by checking $\xi^\mu=0$.

3.3 Why stationary Lorentzian problems become elliptic

A subtle but central point: we are typically solving Lorentzian stationary problems. Why do we still get an elliptic boundary value problem?

The answer is that stationary symmetry reduction effectively removes the time evolution direction from the principal symbol. A standard decomposition (used explicitly in the review) is:

$$ds^2=G_{AB}(x)\big(dz^A+A^A_i(x),dx^i\big)\big(dz^B+A^B_j(x),dx^j\big) + h_{ij}(x),dx^i dx^j,$$

where $z^A$ are coordinates along Killing directions (e.g. $t,\phi_1,\dots$) and $x^i$ are coordinates on a base manifold. If the base metric $h_{ij}$ is Riemannian (positive definite), then the principal symbol of the Einstein–DeTurck operator is controlled by $h^{ij}\partial_i\partial_j$ and is elliptic.

This is the mathematical reason stationary problems can be posed as elliptic boundary value problems.

Caveat (important): there are stationary setups where the reduced equations can become mixed type (or fail to be elliptic), such as certain “heat-flow” or “flowing” configurations. The lectures emphasize that the DeTurck method is best viewed as a flexible gauge-fixing method, not a magic “ellipticizer” in all circumstances.

When is a stationary Einstein problem elliptic?

A crucial (and sometimes underappreciated) point is that “stationary” does not automatically mean “elliptic.”
Ellipticity is a property of the reduced PDE system after you impose symmetries and gauge fixing.

A very general way to organize stationary solutions is to split coordinates into:

• Killing coordinates $y^A$ (for example $t$ and azimuthal angles $\phi_a$),
• base coordinates $x^i$ on the orbit space $\mathcal{B}=\mathcal{M}/G$.

Locally, many stationary metrics can be written in a Kaluza–Klein-like form

$$ds^2 = G_{AB}(x)\big(dy^A+A^A{}_i(x)\,dx^i\big)\big(dy^B+A^B{}_j(x)\,dx^j\big) + h_{ij}(x)\,dx^i dx^j,$$

where all unknown functions depend only on the base coordinates $x^i$.

For the Einstein–DeTurck equations, the principal symbol on perturbations is controlled by the inverse of the base metric $h^{ij}$ in the directions where the fields vary. A sufficient practical condition for ellipticity is that the base metric $h_{ij}$ is Riemannian (positive definite) on the computational domain.

• For static solutions, one can often Wick rotate $t\to i\tau$ and interpret the problem as a Riemannian boundary value problem, where ellipticity is manifest.
• For rotating but stationary solutions with Killing horizons, the base metric can still be Riemannian outside the horizon, and the DeTurck system is elliptic on that domain (with the horizon treated as a regular “fictitious boundary”).

However, there are important situations where ellipticity can fail:

• Flowing/nonequilibrium steady states (e.g. heat transport setups) can involve horizons that are not generated by a single Killing field; in such cases the reduced equations can become mixed hyperbolic–elliptic.
• Certain coordinate/gauge choices can accidentally reintroduce residual gauge modes, leading to a singular Jacobian in Newton’s method even if the continuum equations are formally elliptic.

For this reason, in serious computations one should always:

1. choose an ansatz that is closed under residual coordinate transformations of the base,
2. monitor the conditioning of the linear solves in Newton steps,
3. verify a posteriori that the intended elliptic character is realized (typically by observing robust spectral convergence and well-conditioned Jacobians).

3.4 Einstein–DeTurck with matter fields (and gauge fields)

In holography and many GR applications, we solve Einstein coupled to matter. The DeTurck idea extends cleanly.

A convenient “trace-reversed” form is:

$$R_{\mu\nu} = \frac{2\Lambda}{d-2}g_{\mu\nu} + T_{\mu\nu} -\frac{1}{d-2}T,g_{\mu\nu},$$

and the Einstein–DeTurck version is simply

$$R_{\mu\nu}-\nabla_{(\mu}\xi_{\nu)} = \frac{2\Lambda}{d-2}g_{\mu\nu} + T_{\mu\nu} -\frac{1}{d-2}T,g_{\mu\nu}.$$

This form is explicitly discussed in the review.

For scalar fields, the principal symbol is already well-behaved. For gauge fields, one must also fix gauge. For example, a Maxwell equation $\nabla_\mu F^{\mu\nu}=0$ is gauge-degenerate; the review describes adding a gauge-fixing term so that the principal symbol becomes Laplace-like:

$$\nabla_\mu F^{\mu\nu}+\nabla^\nu\chi=0, \qquad \chi=\nabla_\mu A^\mu-\nabla_\mu \bar A^\mu,$$

which yields a well-posed elliptic operator for stationary problems.

3.5 Ricci solitons and how to rule them out

The Einstein–DeTurck equation is not identical to Einstein’s equation. It admits additional solutions with $\xi\neq 0$, called Ricci solitons (in this context). So we must ensure $\xi=0$. One generally deals with this in two ways:

1. A posteriori check: verify numerically that $|\xi|^2=\xi^\mu\xi_\mu$ goes to zero up to numerical tolerance (improving with resolution).

2. Analytic nonexistence arguments: in some classes (especially static problems with $\Lambda\le 0$ and appropriate boundary conditions), one can prove solitons do not exist using maximum principle arguments for a scalar built from $\xi^\mu$.

3.6 Choosing the reference metric $\bar g_{\mu\nu}$

By choosing different reference metrics $\bar g_{\mu\nu}$, we effectively implement different gauge choices—so gauge fixing is encoded in $\bar g_{\mu\nu}$ rather than being imposed by explicit coordinate conditions.

Practical guidelines:

• $\bar g$ should share the same conformal boundary, the same “locations” of fictitious boundaries (axes/horizons in your coordinate chart), and the same general topology as the desired solution.

• Often $\bar g$ is built by patching together known analytic geometries near each boundary (e.g. near horizon looks like a Rindler-type product; near axis looks like polar coordinates; near infinity looks like AdS). This “patching” viewpoint and its importance are stressed in the lectures.

• When writing the metric ansatz, do not prematurely gauge-fix by hand; let the DeTurck term handle gauge fixing. This is one of the major practical advantages of the method.

DeTurck method

Comparison with other gauge choices

To use the DeTurck method, we typically write a metric ansatz that is closed under the coordinate transformations we want to allow (for example, for a cohomogeneity-two problem we write the metric with the most general form compatible with the symmetries and with dependence on two coordinates).

$$ds^2=-A(r,x)dt^2+B(r,x)(dr^2+dx^2)+C(r,x)dy^2$$ $$ds^2=-Tdt^2+Adr^2+Bdx^2+Sdy^2+2Fdrdx$$

One can compare this to more direct coordinate gauge fixings.

For example, for a static solution depending on $(r,x)$ one could try a gauge where one reduces the number of unknown functions (for instance by adopting a “conformal gauge” form of the metric). Then the ansatz may look much simpler. However, in that approach there is often residual gauge freedom that must still be fixed; otherwise the elliptic operator is not invertible numerically. In many cases it is difficult to identify and fix all residual gauge freedom, and also the chosen coordinate gauge may not be the most suitable for complicated geometries.

Besides conformal gauge, another option that has been explored in some contexts is Bondi–Sachs gauge, which can be convenient for certain problems, but I have not spent enough time exploring it in detail.

So the conclusion is: the Einstein–DeTurck method (solving Einstein–DeTurck instead of Einstein directly) is a very flexible way to choose coordinates and solve stationary Einstein equations, because gauge fixing is built into the PDE system through the reference metric.

---

Why not fix the gauge directly in the metric ansatz?

A tempting approach is to reduce the number of unknown functions in the metric ansatz by using coordinate freedom “by hand.”
For example, in a cohomogeneity-2 static problem one might try to choose a conformal gauge on the base so that the $(x,y)$ part of the metric is explicitly conformally flat, thereby eliminating some functions at the level of the ansatz.

This can work, but it often comes with two practical hazards:

1. Residual gauge freedom.
   Even after simplifying the ansatz, there can remain coordinate transformations of the base $(x,y)$ that preserve the form of the ansatz. These residual diffeomorphisms reintroduce non-uniqueness and typically show up numerically as a non-invertible (or badly conditioned) Jacobian in Newton’s method.

2. The “right” gauge is problem dependent.
   A gauge choice that is convenient for one geometry can be a poor choice for another (especially when the geometry develops multiple scales or complicated boundary structure). Over-fixing the gauge too early can paint you into a corner.

The Einstein–DeTurck method avoids these issues by:

• keeping an ansatz that is closed under reparameterizations of the base coordinates,
• implementing gauge fixing through the DeTurck vector built from a reference metric,
• yielding a determined elliptic system whose solutions can be checked by verifying $\xi=0$.

In GR, the “best coordinates for solving the PDE” are often not the same as the “best coordinates for interpreting physics.”

---

In summary, the DeTurck method has two major advantages at once:

• It helps make the system elliptic (in the Euclidean/Riemannian stationary setting), because the troublesome gauge-degenerate terms in the principal symbol are canceled by the DeTurck term.
• It fixes the gauge automatically, in a flexible way: changing the reference metric $\bar g_{\mu\nu}$ corresponds to changing how you fix the gauge, but you retain a lot of coordinate flexibility because you do not hard-fix coordinates directly at the ansatz stage.

A remark on history: the “DeTurck trick” was originally introduced in the early 1980s in the study of Ricci flow, where adding such a term makes the flow strictly parabolic; later it was adapted to stationary Einstein problems in AdS numerics.

4. Boundary conditions

An operational way of defining stationarity

A stationary solution is a solution to a well-posed boundary value problem. (Dias-Santos-Way 1510.02804)

The traditional definition of stationary solutions in GR (existence of a timelike Killing vector), solutions of well-posed boundary value problems, and solutions of elliptic problems are not identical sets, but they overlap strongly in the applications we care about.

4.1 Manifolds with a boundary

Let $\partial\mathcal{M}$ be a (timelike or spacelike) boundary with outward unit normal $n^\mu$. The boundary geometry is described by the induced metric $\gamma_{ij}$ and the extrinsic curvature $K_{ij}$.

Common boundary conditions include:

• Dirichlet: fix $\gamma_{ij}$ (i.e. fix the boundary metric).
• Neumann: fix $K_{ij}$ (or some part of it).
• Mixed / Robin: fix a combination like $a\,\gamma_{ij}+b\,K_{ij}$. For example, in the entanglement island / brane setups, one imposes mixed (Neumann-type) boundary conditions.
• Modified Dirichlet (very common in gravity): fix the conformal class $[\gamma]$ and the trace $K\equiv \gamma^{ij}K_{ij}$ (or another scalar combination). This is often a natural choice when a conformal factor diverges (as in AdS).

Which one is appropriate depends on the variational principle and on the physical interpretation (e.g. which quantities are “sources” vs “responses”).

There are also asymptotic boundaries, such as asymptotic infinity of Minkowski spacetime and asymptotic boundary of AdS. These are not physical boundaries, but we still impose boundary conditions there. In practice, prescribing asymptotic boundary conditions is often similar in spirit to prescribing boundary conditions on a physical boundary.

4.2 Asymptotic boundary in AdS

For asymptotically AdS spacetimes, the boundary is at infinite proper distance, but one can perform a conformal compactification. Roughly, if $z$ is a defining function with $z=0$ at the boundary, one writes

$$g_{\mu\nu} \sim \frac{L^2}{z^2}\,\tilde g_{\mu\nu}, \qquad z\to 0,$$

where $\tilde g_{\mu\nu}$ is regular at $z=0$. The boundary metric is the restriction $\tilde g|_{z=0}$, defined only up to Weyl transformations.

From the numerical point of view, the practical rule is:

• factor out the known divergent conformal factor;
• impose Dirichlet-type conditions on the remaining smooth fields so that the induced boundary metric is the desired one.

In AdS/CFT, this boundary metric (plus any boundary values of matter fields) is the physical input specifying the QFT background and sources.

Fictitious boundaries: axes and Killing horizons

Many “boundaries” in numerical coordinates are not true boundaries of the manifold — They are coordinate boundaries. They are loci where a cycle shrinks smoothly (axis) or where a Killing vector degenerates (Killing horizon). The correct boundary conditions are regularity conditions.

Axis of symmetry

Near an axis where an angular coordinate $\phi$ degenerates, regularity typically demands that the metric in the $(\rho,\phi)$ plane looks like flat polar coordinates:

$$d\rho^2+\rho^2 d\phi^2$$

up to smooth warping by functions of the transverse coordinates. The angular coordinate $\phi$ must have the correct period: $2\pi$. Otherwise we get a conical defect, but the spacetime we want does not have such a defect. This implies conditions such as:

• no cross terms that would be singular,
• evenness of functions in $\rho$,
• specific equalities between metric functions multiplying $d\rho^2$ and $\rho^2 d\phi^2$ (to avoid conical defects).

In spectral implementations, these often translate into Neumann conditions $\partial_\rho(\cdots)=0$ at $\rho=0$ plus one algebraic condition fixing the conical deficit.

Non-extremal Killing horizon

In Euclidean signature, a non-extremal Killing horizon behaves like the origin in polar coordinates (regularity fixes the relation between the Euclidean time periodicity and surface gravity).

For a static black hole, after Wick rotation to Euclidean signature, the horizon behaves just like an axis of rotation. So the horizon regularity condition becomes the analog of “no conical singularity,” and it fixes the periodicity of Euclidean time:

$$\tau \sim \tau + \beta,\qquad \beta = \frac{2\pi}{\kappa} = \frac{1}{T},$$

where $\kappa$ is surface gravity and $T$ is temperature.

A non-extremal Killing horizon generated by

$$K=\partial_t+\sum_a \Omega_H^{(a)}\,\partial_{\phi_a}$$

is locally analogous to an axis after Wick rotation. In suitable coordinates, one expects near the horizon ($\rho=0$)

$$ds^2\approx -\kappa^2\rho^2 dt^2 + d\rho^2 + \gamma_{ab}(x)dx^a dx^b + \cdots,$$

where $\kappa$ is the surface gravity and $T=\kappa/(2\pi)$ is the Hawking temperature.

Regularity imposes:

• relations among metric functions at the horizon (often $g_{tt}$ and the inverse of $g_{\rho\rho}$ must vanish at the same rate);
• absence of divergent cross terms in a corotating frame;
• constancy of $\kappa$ over the horizon for a genuine Killing horizon.

A practical implementation strategy is:

• choose a reference metric $\bar g$ that already has the desired horizon location and temperature;
• factor out the expected $(1-y)$ or $\rho^2$ behavior in the ansatz so that unknown functions remain smooth and finite;
• impose the remaining regularity as Neumann/Dirichlet conditions at the horizon boundary in compact coordinates.

Saying in words “$\beta$ has period $2\pi$” is not yet an equation. How do you translate that physical statement into the actual boundary conditions you implement in the PDE system? This “translation from physics language to mathematical equations” depends on the coordinate system you are using. Different coordinates lead to different explicit forms of the regularity conditions.

Extremal horizon

Non-extremal horizons are relatively straightforward: there are universal regularity conditions. Extremal horizons are more complicated: in the zero-temperature limit, the IR geometry can be of different types. For example, for the extremal Reissner–Nordström black hole, the near-horizon region is $AdS_2\times S^2$ and the extremal horizon has finite entropy. A standard extremal black hole is the case where an inner horizon and outer horizon merge. In terms of metric functions, you can think of it as two simple zeros merging into a double zero, e.g. schematically

$$g_{tt}\sim -(r-r_+)^2 \quad \text{(double zero, extremal)}$$

rather than

$$g_{tt}\sim -(r-r_+) \quad \text{(simple zero, non-extremal)}.$$

Extremal horizons are qualitatively different:

• the surface gravity vanishes, $\kappa=0$;
• the near-horizon geometry typically develops an $\mathrm{AdS}_2$ throat (or another scaling region);
• the horizon can sit at infinite proper distance in suitable slices.

Numerically, extremal horizons often behave more like an additional asymptotic region than like a smooth “axis-like” boundary. Spectral methods can still work, but one must:

• pick coordinates adapted to the scaling region;
• impose boundary conditions through near-horizon expansions rather than naive “regular at $\rho=0$” rules;
• expect that continuation from finite temperature down to $T\to 0$ may become increasingly stiff.

In physics, extremal (zero-temperature) horizons are very important, especially in condensed matter applications, because zero temperature corresponds to the ground state. In holography language, you often think of the AdS boundary as the UV and the deep interior as the IR. At zero temperature, the IR geometry can take different forms: it might be an extremal horizon, it might be a “pinching” geometry, it might be singular, and so on. So extremal horizons are difficult and diverse.

Also, at zero temperature, whether a spectral method works becomes more delicate. It’s not that spectral methods can never be used, but they often run into difficulties. Spectral methods tend to work better for finite-temperature problems.

Non-Killing horizons

For a Killing horizon, the temperature (surface gravity) is well-defined and constant. There are also cases where surface gravity is not well-defined globally, or the “effective temperature” varies along the horizon — that corresponds to a non-Killing horizon. For example, when there is heat flow along the horizon, the surface gravity is not constant. Both non-Killing horizons and extremal horizons are harder than the standard finite-temperature Killing horizon case.

---

It is conceptually helpful to distinguish:

• Defining boundary conditions: the data that defines the physical problem (e.g. the boundary metric for an AdS/CFT setup, the temperature and angular velocities of a horizon, the values of external sources).
• Derived boundary conditions: conditions that are forced by smoothness and/or by the equations of motion once the defining data is chosen (e.g. certain Neumann relations at an axis, or relations among metric functions at a regular horizon).

In practice, most “fictitious boundary” conditions (axes, horizons) are of the derived/regularity type, while asymptotic AdS data is typically defining.

---

What happens if boundary conditions are missing or wrong?

Consider a schematic equation

$$\mathcal{D}u = f,$$

and try to invert the operator. If $\mathcal{D}$ is not invertible, you can get:

• non-uniqueness (many solutions),
• or even no solution (if constraints are inconsistent).

A simple example: take

$$\frac{d}{dx}u(x) = f(x).$$

Without boundary conditions, the solution is not unique: $u(x)+C$ is also a solution. In linear algebra language, the matrix representation of $\frac{d}{dx}$ is singular unless you fix the boundary conditions appropriately. I have a slide specifically about this linear algebra point.

If the boundary value problem is truly well-defined, then after discretization the resulting matrix should be invertible and you should get a unique solution. If boundary conditions are not correct, you might get non-uniqueness; if they are inconsistent, you may get no solution.

---

4.5 Boundary conditions and the DeTurck vector

Boundary conditions must also be compatible with the requirement $\xi=0$ for Einstein solutions.

A useful identity is that the contracted Bianchi identity implies that $\xi$ satisfies a linear, elliptic, second-order equation on a solution of the Einstein–DeTurck system:

$$\nabla^2\xi_\nu + R_{\nu\mu}\xi^\mu = 0$$

(up to conventional overall factors depending on definitions).

Therefore, if your boundary conditions enforce $\xi|_{\partial\mathcal{M}}=0$ (or otherwise rule out nontrivial solutions of this elliptic equation), then $\xi$ must vanish everywhere.

In many practical coordinate systems where a fictitious boundary sits at (say) $\rho=0$ with $\rho$ the normal coordinate, regularity conditions for the metric typically imply boundary conditions of the schematic form

$$\xi^{\rho}\big|_{\rho=0}=0, \qquad \partial_{\rho}\xi^{\tilde i}\big|_{\rho=0}=0,$$

where $\tilde i$ labels directions tangent to the boundary (including Killing directions). These are compatible with $\xi=0$ and ensure that the elliptic problem for $\xi$ is well-posed. In many static AdS problems, one can turn this into a rigorous nonexistence theorem for solitons; in general stationary problems, treat $\xi$ as a stringent numerical diagnostic.

5. Numerical implementation: pseudospectral method and Newton–Raphson

Methods

So, from the numerical relativity perspective for stationary solutions, a good “toolbox” is the combination of three methods:

1. Einstein–DeTurck method (for gauge fixing and getting a good boundary value formulation),
2. Spectral method (for discretization/approximation),
3. Newton’s method (for solving nonlinear equations).

Of course each has alternatives:

• DeTurck vs fixing gauge directly in the ansatz (but residual gauge problems become hard).
• Spectral vs finite difference vs finite element:
  • finite difference approximates locally by low-degree polynomials,
  • spectral methods approximate globally using a complete basis on the whole domain.
• Newton vs flow/relaxation methods (e.g. Ricci(-DeTurck) flow–type relaxation).

The package GRSpectral implements the above method.

5.1 Elliptic vs hyperbolic vs parabolic (prototype equations)

Elliptic

$$(\nabla^2+m^2)u=f$$

Hyperbolic

$$(\partial_t^2-\nabla^2)u=f$$

Parabolic

$$(\partial_t-\nabla^2)u=f$$

Stationary Einstein problems, once correctly gauge-fixed and symmetry-reduced, are designed to fall into the elliptic category so they can be solved as boundary value problems.

The core philosophy of spectral methods is:

Expand the unknown functions in a complete basis on the domain (a Hilbert-space viewpoint), and aim to achieve high accuracy with as few basis functions as possible.

This contrasts with:

• Finite differences: approximate derivatives by local polynomial fits (stencils). They are conceptually simple and flexible for nonsmooth problems, but achieving high accuracy can require many grid points.
• Finite elements: approximate the solution by piecewise polynomials on a mesh; very powerful for complicated geometries and nonsmoothness, but typically more infrastructure-heavy.

For smooth stationary geometries on simple domains, spectral collocation often wins because it can reach very high accuracy with comparatively few degrees of freedom.

---

Functions in Hilber space

Think of a 3-dimensional Euclidean space. A vector can be expanded in a basis. Now think of a Hilbert space: a function can be expanded in a basis. When a differential equation does not have a closed-form solution, one way to treat it is to regard the unknown function as an element of a function space — a “vector” in an infinite-dimensional space — and then solve for its coordinates in a chosen basis.

Think about Descartes’ contribution. If you study a circle, you can draw a circle; if you study an ellipse, it’s harder to draw precisely. Descartes introduced coordinates: once you have a coordinate system, you describe points by numbers. You don’t need to draw the picture.

Similarly, how do we treat a function? We also build a “coordinate system” — by choosing a basis. In principle this coordinate system has infinitely many axes. A function corresponds to a “vector” in this infinite-dimensional coordinate system. That space is a Hilbert space. Euclidean space is really just a special case of a Hilbert space viewpoint.

But choosing a basis is not trivial. You want an orthogonal and complete basis. Once you choose it, you can expand the function in that basis, and solving the equation becomes solving for the expansion coefficients.

In practice, we hope that we can use as few basis elements as possible, because that is what makes computation efficient. For periodic functions, Fourier modes are natural: a finite number of modes can already approximate the function well.

This is not merely an analogy; mathematically, it is the same structure:

• Orthogonality corresponds to the inner product.
• Normalization corresponds to choosing basis functions of unit norm.
• Completeness corresponds to the ability to represent the function arbitrarily well.
• In finite-dimensional linear algebra, you get components by dot products with basis vectors.
• In function spaces, you get Fourier coefficients by inner products.

---

For a finite interval, the most common and effective choice is the Chebyshev basis. A standard set of collocation points on $[x_-,x_+]$ is the Chebyshev–Gauss–Lobatto grid

$$x_j = \frac{x_+ + x_-}{2} + \frac{x_+ - x_-}{2}\cos\!\left(\frac{j\pi}{N}\right), \qquad j=0,1,\dots,N.$$

(These points cluster near the endpoints, avoiding the Runge phenomenon and enabling exponential convergence for smooth functions.)

At collocation points, derivatives are approximated by differentiation matrices:

$$f^{(n)}(x_i)\approx \sum_j D^{(n)}_{ij} f(x_j).$$

Consider a 1D example. The Chebyshev grid (8 points) is

$$x=(-1.,\, -0.9009,\, -0.6234,\, -0.2225,\, 0.2225,\, 0.6234,\, 0.9009,\, 1.)^T$$

The corresponding first-derivative differentiation matrix is

$$D\equiv\frac{d}{dx}= \left( \begin{array}{cccccccc} -16.5 & 20.19 & -5.311 & 2.572 & -1.635 & 1.231 & -1.052 & 0.5 \\ -5.048 & 2.392 & 3.603 & -1.473 & 0.8900 & -0.6559 & 0.5549 & -0.2630 \\ 1.327 & -3.603 & 0.5100 & 2.493 & -1.182 & 0.8019 & -0.6559 & 0.3079 \\ -0.6431 & 1.473 & -2.493 & 0.1170 & 2.246 & -1.182 & 0.8900 & -0.4089 \\ 0.4089 & -0.8900 & 1.182 & -2.246 & -0.1170 & 2.493 & -1.473 & 0.6431 \\ -0.3079 & 0.6559 & -0.8019 & 1.182 & -2.493 & -0.5100 & 3.603 & -1.327 \\ 0.2630 & -0.5549 & 0.6559 & -0.8900 & 1.473 & -3.603 & -2.392 & 5.048 \\ -0.5 & 1.052 & -1.231 & 1.635 & -2.572 & 5.311 & -20.19 & 16.5 \\ \end{array} \right)$$

Once derivatives become matrix multiplications, a linear differential equation schematically becomes a linear algebra problem:

$$Lu=f \quad\Rightarrow\quad u=L^{-1}f,$$

and nonlinear equations become nonlinear algebraic systems.

Multi-dimensional tensor-product grids

If your domain is a rectangle (or a product of intervals), a major advantage is that you can use tensor-product grids and Kronecker products to build derivative operators efficiently.

For example, on a 2D grid with $N_x+1$ points in $x$ and $N_y+1$ points in $y$, you can represent $f(x_i,y_j)$ as a vector $f_I$ using a fixed index map $I\leftrightarrow (i,j)$, and write

$$D_x = I_{N_y+1}\otimes D_x^{(1)}, \qquad D_y = D_y^{(1)}\otimes I_{N_x+1}.$$

This is one reason why choosing coordinates that make the domain a rectangle is so valuable.

Enforcing boundary conditions in collocation form

A very practical technique is:

• keep the PDE residual equations at interior grid points;
• replace the PDE residual at boundary collocation points by the boundary condition residual.

This row-replacement viewpoint is helpful when thinking about invertibility: boundary conditions are not “extra”; they are the equations that make the discrete operator invertible.

---

Eigenvalue problems vs non-eigenvalue problems

In linear algebra, a system $A u = b$:

• If $b\neq 0$ (at least one component nonzero), and $A$ is invertible, you solve for $u$. This corresponds to a typical boundary value problem with sources.
• If $b=0$, then nontrivial solutions exist only when $A$ is not invertible: $$\det A = 0.$$ That is an eigenvalue-type situation: only for special parameter values can you have nontrivial solutions.

This is the same as for differential equations. Many problems we want to solve are eigenvalue problems (normal modes, instabilities, etc.), and many are non-eigenvalue problems (driven boundary value problems).

---

5.3 Newton–Raphson for nonlinear PDEs

After discretization, the unknown fields become a vector $U\in\mathbb{R}^N$ (often with $N$ in the thousands to millions), and the PDEs become a nonlinear system

$$F(U)=0.$$

Newton’s method linearizes around a current iterate $U^{(k)}$:

$$J(U^{(k)})\,\delta U^{(k)} = -F(U^{(k)}), \qquad U^{(k+1)}=U^{(k)}+\delta U^{(k)},$$

where $J=\partial F/\partial U$ is the Jacobian.

A compact pseudocode summary is:

Given: reference metric and fields → initial guess U₀
for k = 0,1,2,...
    compute residual  F(U_k)
    if ||F|| and ||ξ|| below tolerance: stop (converged)
    assemble / apply Jacobian J(U_k)
    solve linear system   J δU = -F
    choose damping λ ∈ (0,1] (optional line search)
    update U_{k+1} = U_k + λ δU
end

In large systems one often avoids forming $J$ explicitly and instead implements Jacobian–vector products for iterative Krylov solvers; in moderate-size problems direct sparse factorization can be faster and simpler.

A robust practical implementation usually includes:

• a line search / damping parameter $\lambda\in(0,1]$ so $U\to U+\lambda\delta U$;
• a good linear solver for $J\delta U=-F$ (direct for moderate size, iterative for large size);
• a careful stopping criterion based on residual norms and $\|\xi\|$.

The unavoidable issue: an initial guess

But regardless of which nonlinear solver you use, one thing is basically unavoidable: you need an initial guess. For example, Mathematica’s FindRoot asks for an initial guess; this is a universal feature of nonlinear solving. So finding a good initial guess is an important technique.

Two standard strategies:

1. Reference-metric guess: start from the reference metric $\bar g$ (and reference matter fields).
   Works when the solution is close to $\bar g$ or the nonlinearity is mild.

2. Parameter continuation (highly recommended):
   start from a parameter regime where the solution is known (often high temperature, weak deformation, small amplitude), solve there, then slowly change parameters and use the previous solution as the next initial guess.

For example, if a family of solutions depends on a parameter $\alpha$, solve first at $\alpha=\alpha_0$, then use that converged solution as a seed for $\alpha=\alpha_0+\Delta\alpha$, and so on. This continuation strategy becomes essential at low temperatures or near critical points, where Newton’s basin of attraction shrinks.

5.4 Diagnostics: what to check

A credible stationary solution should pass (at least) the following checks:

• DeTurck norm: verify $\|\xi\|$ is small and decreases with increasing resolution (ideally exponentially in $N$ for spectral methods).
• Residual norm: verify $\|F(U)\|$ is small and convergent.
• Convergence test: repeat at multiple resolutions and confirm exponential convergence for smooth problems.
• Constraint/identity checks: verify independent identities (e.g. Smarr-like relations, first law, Ward identities for holographic stress tensor) when applicable.
• Regularity: inspect curvature invariants and confirm no unphysical singularities.

The package GRSpectral implements this overall strategy (DeTurck + pseudospectral + Newton) in a way tailored to stationary AdS problems.

---

---

$$Lu=f\quad \Rightarrow\quad u=L^{-1}f$$

Now let me briefly talk about how to implement the numerical techniques to construct these solutions. The basic idea is:

• The unknown functions form a vector in some function space.
• The differential operator acts like a (nonlinear) operator.
• After discretization, the operator becomes a matrix acting on a finite-dimensional vector.

To solve an elliptic system, we discretize the coordinates and the differential operator, impose boundary conditions, and then invert the resulting operator (or solve a linear system).

Concretely:

1. Construct a grid (discretize the coordinates and functions).
2. Write the discretized differential operator as a matrix.
3. Impose boundary conditions (often by replacing appropriate rows corresponding to boundary points).
4. Solve the resulting linear system.

But Einstein–DeTurck is a nonlinear equation, so we need an iterative method. Typically we use the Newton–Raphson method.

In one dimension, Newton’s method for solving $f(x)=0$ is

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.$$

For a discretized nonlinear system $F(u)=0$ (where $u$ is a vector collecting all grid values of all unknown functions), Newton’s method takes the form

$$J(u_n)\,\delta u = -F(u_n),\qquad u_{n+1}=u_n+\delta u,$$

where $J$ is the Jacobian matrix of $F$.

---

6. Examples

The general workflow for solving physics problems are

1. Pose the physical problem.
2. Translate the physical problem into a mathematical problem.
3. Solve the mathematical problem.
4. Interpret the physical meaning of the solution.

Note that step 2 is crucial: you must translate the physical question into a well-posed boundary value problem—write down the equations and the boundary conditions explicitly—before you can feed anything into a solver.

Examples of stationary AdS solutions: funnels, droplets, tunnels, hammocks, binaries, islands

If we put $\mathcal{N}=4$ super-Yang–Mills theory on a Schwarzschild black hole background (i.e. the boundary metric has a black hole), then there are two phases:

• The higher-temperature phase is a deconfined phase, and the dual bulk solution is called a black funnel.
• The lower-temperature phase is a confined phase, and the bulk solution is called a black droplet. Black droplets involve two disconnected horizons. Typically the domain is more complicated: for example one may have boundaries corresponding to the boundary black hole, an axis, a bulk horizon, and also a planar horizon. This is significantly more difficult than the simpler planar black hole case, and often requires clever coordinate choices in order to implement boundary conditions cleanly.

In the typical schematic figure, the middle plot corresponds to the black droplet, and the right plot corresponds to the black funnel.

If we put the CFT on a “cosmological” black hole background (for example a Schwarzschild–de Sitter boundary metric), similarly there are also two phases. The deconfined phase is associated with a black hammock, and the confined phase is associated with a black tunnel. The AdS boundary has two horizons: one is the black hole horizon, and the other is the cosmological horizon.

There are also black binaries. Usually when people talk about numerical relativity, they mean binaries orbiting each other and a time-evolution problem. But in AdS, there can be static (or stationary) configurations of binaries. Those are very difficult solutions to construct.

Another interesting direction is entanglement islands in higher dimensions. This is related to important progress on the black hole information paradox, and is a successful application of numerical methods for stationary solutions.

In that setup, there is a brane in AdS, and the brane typically comes with mixed (Neumann-type) boundary conditions. When we put a brane in AdS spacetime, the nontrivial boundary condition backreacts on the geometry. So the geometry depends on a spatial coordinate in addition to the radial coordinate.

And of course there are many other applications in holographic quantum matter: holographic superconductors, charge density waves, pair density waves, and so on. One can attempt to construct phase diagrams of strongly coupled systems using these gravitational solutions. Holographic QCD also uses numerical techniques to solve Einstein’s equations.

A useful strategy (especially for funnel/droplet-like geometries) is to introduce compact coordinates $(x,y)$ so that the computational domain is a rectangle, with different boundaries corresponding to:

• an asymptotic AdS boundary (after factoring the conformal divergence),
• a boundary black hole horizon (or another “boundary feature”),
• a bulk horizon,
• an axis or another asymptotic end (e.g. the hyperbolic black hole infinity).

A concrete metric ansatz of this type is:

$$ds^2=\frac{L^2}{x\,y(1+x)^2}\bigg\{ -x(1-y)(1+x\,y)\,T\,dt^2 +\frac{x(1+x)^2A\,dy^2}{4y(1-y)(1+x\,y)} +\frac{r_0^2B\,[dx+x(1-x)^2F\,dy]^2}{x(1-x)^4} +\frac{r_0^2S}{(1-x)^2}\,d\Omega_{2}^2 \bigg\},$$

where $T,A,B,S,F$ are functions of $(x,y)$ to be solved for, and the prefactors have been chosen so that:

• known vanishing/divergent behavior near boundaries is built in;
• the remaining unknown functions are expected to be smooth on the closed rectangle.

A simple reference metric is obtained by setting

$$T=A=B=S=1,\qquad F=0.$$

What boundary conditions look like in practice

For an ansatz of this form, boundary conditions typically include:

• Asymptotic boundary (AdS): Dirichlet conditions fixing the induced boundary metric (equivalently, fixing the leading coefficients in the asymptotic expansion).
• Bulk horizon: regularity conditions enforcing a smooth non-extremal Killing horizon (often Neumann conditions plus an algebraic equality among some functions, and matching the horizon temperature to that of the reference metric).
• Axis: regularity conditions preventing conical singularities.
• Corner consistency: when two boundaries meet, the expansions implied by each boundary must be compatible. Choosing $\bar g$ to satisfy all boundary conditions to sufficient order dramatically simplifies this.

After imposing these, one solves the Einstein–DeTurck equations for $T,A,B,S,F$ with Newton–Raphson, often using parameter continuation in temperature, boundary deformation strength, or rotation.

---

Other common stationary targets (where the same toolkit applies)

Once you are comfortable with the DeTurck + spectral + Newton pipeline, you can tackle a wide range of stationary problems, for example:

• Rotating AdS black holes with deformed horizons (“lumpy” solutions and branches connected to superradiant instabilities).
• Horizonless stationary solutions such as AdS geons and boson stars (Einstein–scalar), obtained by imposing regularity at the origin instead of a horizon.
• Black resonators and other solutions with only a single Killing field (often requiring careful treatment of symmetry assumptions).
• Black rings/strings and other higher-dimensional horizon topologies, where cohomogeneity-$\ge 2$ is generic and numerics is essential.
• CFT steady states with transport, where one must pay close attention to whether the reduced equations remain elliptic or become mixed type.

The details differ (ansatz and boundary conditions), but the numerical backbone is the same.

Metric ansatz:

$$ds^2=\frac{L^2}{x\,y(1+x)^2}\bigg\{-x(1-y)(1+x\,y)Tdt^2+\frac{x(1+x)^2Ady^2}{4y(1-y)(1+x\,y)}+\frac{r_0^2B[dx+x(1-x)^2Fdy]^2}{x(1-x)^4}+\frac{r_0^2S}{(1-x)^2}d\Omega_{2}^2\bigg\}$$

The reference metric given by $T=A=B=S=1$ and $F=0$.

{\hfill [J.E. Santos, B. Way, 1208.6291]}

{\bf Planar black hole} at $x=1$. Boundary conditions are $T=A=B=S=1$ and $F=0$.

$$ds^2=\frac{L^2}{y}\left[-\frac{1}{4}(1-y^2)dt^2+\frac{dy^2}{4y(1-y^2)}+\frac{r_0^2dx^2}{4 (1-x)^4}+\frac{r_0^2}{4(1-x)^2}d\Omega_2^2\right]$$

After coordinate transformation $t = 2 \tau$, $x = 1-r_0/2R$, $y= z^2$:

$$ds^2=\frac{L^2}{z^2}\left[-(1-z^4)d\tau^2+\frac{dz^2}{1-z^4}+dR^2+R^2d\Omega_2^2\right]$$

{\bf Conformal boundary} at $y=0$.

$$ds^2=\frac{L^2}{y}\left[\frac{dy^2}{4 y}+\frac{1}{x(1+x)^2}\;ds_\partial^2\right],$$ $$ds_\partial^2=-x\,dt^2+\frac{r_0^2dx^2}{x(1-x)^4}+\frac{r_0^2}{(1-x)^2}d\Omega_2^2$$

After $x= 1-r_0/r$, it becomes the Schwarzschild metric:

$$ds_\partial^2=-\left(1-\frac{r_0}{r}\right)dt^2+\frac{dr^2}{1-\frac{r_0}{r}}+r^2d\Omega_2^2\;.$$

{\bf Horizon} at $y=1$.

Expanding the equations of motion about the horizon will give the condition $T=A$ and other conditions on $\partial_yA|_{y=1}$, $\partial_yB|_{y=1}$, $\partial_yF|_{y=1}$, and $\partial_yS|_{y=1}$.

{\bf Hyperbolic black hole} at $x=0$.

$$ds^2=\frac{L^2}{y}\left[-(1-y)dt^2+\frac{dy^2}{4y(1-y)}+\frac{dx^2}{4x^2}+\frac{1}{4x}d\Omega_2^2\right]$$

After coordinate transformation $y=z^2$ and $x=\exp(-2\eta)$,

$$ds^2=\frac{L^2}{z^2}\left[-(1-z^2)dt^2+\frac{dz^2}{1-z^2}+d\eta^2+\frac{e^{2\eta}}{4}d\Omega_2^2\right]$$

Zero energy hyperbolic black hole:

$$ds^2_{\mathbb{H}}=\frac{L^2}{z^2}\left[-(1-z^2)dt^2+\frac{dz^2}{1-z^2}+d\eta^2+\sinh^2\eta\, d\Omega_{d-3}^2\right]$$

7. A Mathematica package: GRSpectral

Mathematica’s NDSolve handles Dirichlet boundary conditions relatively well. But once you go beyond Dirichlet — e.g. to Neumann or mixed boundary conditions — things quickly become extremely difficult.

There are big difference between industrial problems and theoretical physics problems.

Industrial problems are often:

• standard equations (e.g. elasticity equations, Maxwell equations, etc.),
• but with very complicated boundaries (complex shapes of parts).

To solve complicated boundaries, industrial software relies heavily on the finite element method. But FEM usually requires writing the PDE in a weak form (integration by parts, test functions, etc.). For relatively simple standard PDEs, one can do this systematically.

However, for equations as complicated as the Einstein–DeTurck equations, writing the weak form by hand is basically impossible, and it’s not clear how to automate it robustly.

In theoretical physics research, the situation is often the opposite:

• the equations are complicated and nonstandard (Einstein–DeTurck is a prime example),
• but the boundaries are often relatively simple, and complicated boundaries can often be simplified by coordinate transformations.

So industrial FEM software is not naturally suited to our type of problems.

This is why I felt it was necessary to develop a package that, given:

• equations (in “strong form”),
• and boundary conditions (also in a direct form), can solve the boundary value problem using spectral methods, without requiring a weak form.

Mathematica technically allows you to input equations without writing weak form, but internally it treats Dirichlet conditions as “the boundary condition,” while Neumann/mixed conditions are handled in a more complicated way (often treated as part of the equation system). For the Einstein–DeTurck system and general boundary conditions, this becomes almost impossible to handle directly in NDSolve.

---

So I wrote a package with a flow-chart style design: you input equations and boundary conditions and it solves.

It is available on GitHub. You can download the .m file and install it according to the instructions (save it into an appropriate system folder).

The package has only two core functions:

• EqProcess: specify the equations and boundary conditions.
• SP... (the solver function): construct the collocation grid and compute the numerical solution (this is the “spectral version” of a solve routine).

Concretely, you must input the essential things (nothing redundant):

1. Unknown functions (the fields to be solved). For example, in a complicated Einstein system you might have many functions $G_1,G_2,\dots$. You don’t need to understand the particular example I showed on the slide; the point is it can handle many coupled unknown functions.

2. Equations. Often you derive them once and store them in a file (since they can be extremely long), and then load them.

3. Boundary conditions. Boundary conditions must be written as equations. For instance:

   • at $z=0$, impose $G_1=1$,
   • or impose $\partial_x G_1 = 0$, etc.

4. An initial guess (initial trial solution). This is required for nonlinear solving.

5. Optionally, specify what outputs you want.

If you don’t specify outputs, it outputs the solved functions by default. But anyway: these are the necessary inputs. This is what it means to define a well-posed boundary value problem.

Then you call the solver, which I named something like spndsolve. This name has nothing to do with Mathematica’s NDSolve; I just chose a similar naming because it plays the same “role,” but it solves using spectral methods. Here sp refers to “spectral.”

After you call it, it constructs the collocation grid and the problem is solved. Ideally, you do not need to learn special numerical tricks — you just call these two functions.

Now about the grid input:

You specify the coordinates, say $z$ and $x$, because the metric depends on these two coordinates.

You typically map the coordinate range to the Chebyshev interval $[-1,1]$. For example, in Poincaré coordinates for a black hole geometry, the radial coordinate might be $z\in[0,1]$, where $z=0$ is the AdS boundary and $z=1$ is the horizon. Then you map $[0,1]\to[-1,1]$ with a simple affine map.

You also specify the grid sizes, e.g. $20\times 20$. In some problems, one direction might require a periodic grid (for a periodic coordinate). In that case you specify it accordingly.

Finally, if the problem has parameters (e.g. $a$, $c$, etc.), you pass their values. Writing parameters this way has a practical advantage: you don’t have to constantly Clear symbols, and you can sweep parameters systematically.

So, in summary: input functions, equations, boundary conditions, initial guess, and grid/coordinate mapping. Then the package automates the rest.

---

Now, about multiple solutions and initial guesses.

Sometimes if the problem has non-uniqueness (for example in an eigenvalue-type setting), you might get different solutions depending on the initial guess. This is an important issue. One key technique is to try two different initial guesses.

For example, if your solution is expected to be close to a Schwarzschild-like black hole, you can use the Schwarzschild solution as the initial guess. Also, when you vary a parameter like temperature: at high temperature you might have a solution close to a known one; as you lower temperature, the solution may flow toward a different IR geometry. This is often how you follow branches.

Check precision and stability.

After you obtain a solution, you should compute some quantities that should be zero (or should satisfy some identity). For instance, you can compute:

• constraints that should vanish,
• or other consistency checks (depending on the physical problem).

They will never be exactly zero numerically. The key is: as you increase the number of collocation points $N$, the values of these “should be zero” quantities should trend toward zero.

This is a very important way to test both accuracy and stability.

A more strict approach is: increase $N$ significantly, say from $20$ to $60$, and check the convergence trend. This trend is extremely important. If it doesn’t behave as expected, something may be wrong. Sometimes you can find a solution that “converges,” but the convergence rate is not as expected — then there might be a subtle bug. Such subtle issues occur frequently in designing iterative methods.

As $N$ increases:

• Accuracy should improve (and ideally for spectral methods you often expect very fast convergence).
• Computation becomes slower (bigger matrices).

You might see quantities go from $10^{-2}$ to $10^{-3}$ to $10^{-4}$, etc. But it can plateau: for example, it may get stuck around $10^{-5}$ and then stop improving even as you increase $N$.

One possible reason is that if $N$ becomes too large, the matrix becomes more ill-conditioned and you hit limitations of machine precision. Then you must increase working precision. Roughly speaking, the number of digits of working precision should exceed what the condition number demands; otherwise, increasing $N$ can even make results worse.

The boundary conditions are written as equations, and treated on the same footing as the equations. Boundary conditions can be just as complicated and nonlinear as the equations — that makes no essential difference.

Finally, let me say something about the applicability of spectral methods.

Any spectral method has a range of applicability. Typically you need the solution to be analytic inside the domain. At the boundary, you can sometimes tolerate mild non-analyticity — for example, expansions like

$$a_0 + a_1 x + a_2 x^2 + \cdots + a_4 x^4 + \cdots + (\text{polynomial})\times \log x,$$

or fractional powers with sufficiently “mild” behavior.

But if the two independent asymptotic behaviors are truly two different fractional powers and neither can be made analytic in a suitable way, then spectral methods become very hard to apply. Sometimes you can fix this by factoring out a fractional power so that one branch becomes analytic, and then imposing that the coefficient of the non-analytic branch is zero as the boundary condition. That can work in some cases.

In any case, after you solve, you should check whether the solution matches the expected asymptotics and physical behavior.

I call this package spectralNDSolve (sometimes I say “spndsolve”), because I was disappointed with Mathematica’s built-in NDSolve for this kind of problem. NDSolve uses finite element methods. Finite element methods are suitable for problems with complicated boundary shapes but standard equations; and there are many expensive commercial packages that do this very well.

But for our problems, it is the opposite: Einstein equations are highly non-standard PDEs, while the boundary shapes are usually simple (or we can make them simple by choosing appropriate coordinates). So I found that many standard packages are not well suited to stationary solutions of Einstein equations in practice, and I wrote this package.

The basic usage works in two stages:

1. Symbolic stage: specify equations and boundary conditions using a function like EQProcess.
2. Numerical stage: solve the system numerically by calling spectralNDSolve with grid information.

I have not explained what a “holographic pair density wave” is here; it is a gravitational dual of a kind of holographic superconductor. What I want to emphasize is that the equations are highly nontrivial PDEs—very complicated. But this package automatically handles many internal complexities: from the user side, the input is only the equations and boundary conditions (of course, after you derive the equations).

Workflow

1. The equation we solve numerically is not the raw Einstein equation, but the Einstein–DeTurck equation.

2. After solving, we must verify that the DeTurck vector vanishes $\xi=0$.

3. When writing the metric ansatz, we should not fix gauge prematurely. The gauge-fixing is handled automatically by the DeTurck formulation.

8. References for NR

Stationary

[1] M. Headrick, S. Kitchen and T. Wiseman, A new approach to static numerical relativity, and its application to Kaluza-Klein black holes, Class. Quant. Grav. 27 (2010) 035002 [arXiv:0905.1822].

[2] P. Figueras, J. Lucietti and T. Wiseman, Ricci solitons, Ricci flow, and strongly coupled CFT in the Schwarzschild Unruh or Boulware vacua, Class. Quant. Grav. 28 (2011) 215018 [arXiv:1104.4489].

[3] A. Adam, S. Kitchen and T. Wiseman, A numerical approach to finding general stationary vacuum black holes, Class. Quant. Grav. 29 (2012) 165002 [arXiv:1105.6347].

[4] O. J. C. Dias, J. E. Santos and B. Way, Numerical Methods for Finding Stationary Gravitational Solutions, Class. Quant. Grav. 33 (2016) 133001 [arXiv:1510.02804]. ---Review

physics-oriented review: holographic thermal field theory

Pseudospectral method: two books

Time evolution

[5] P. M. Chesler and L. G. Yaffe, Numerical solution of gravitational dynamics in asymptotically anti-de Sitter spacetimes, JHEP 07 (2014) 086 [arXiv:1309.1439].

[6] K. Balasubramanian and C. P. Herzog, Losing Forward Momentum Holographically, Class. Quant. Grav. 31 (2014) 125010 [arXiv:1312.4953].

[7] Z. Ning, Q. Chen, Y. Tian, X. Wu and H. Zhang, Spontaneous Deformation of an AdS Spherical Black Hole, Phys. Rev. D 109 (2024) 064082 [arXiv:2307.14156].

Talk

6. The Einstein–DeTurck method and Ricci solitons

If we look at the principal symbol of this modified equation, then in Euclidean signature it becomes elliptic. For static solutions in Euclidean signature, this is a clean elliptic boundary value problem. For stationary but non-static Lorentzian problems, it is still elliptic under certain conditions and coordinate choices, but it is not always elliptic in general. Still, one can often proceed as long as the discretized operator is invertible.

One advantage of the method is that $\xi$ depends on the reference metric $\bar g$. Choosing different reference metrics corresponds to choosing different gauges / coordinate conditions, so this is a flexible way to choose coordinates suited to the geometry.

Example 1: Black funnels and hyperbolic/planar coordinates

A useful technical point is that one wants to choose coordinates so that the computational domain becomes as simple as possible—ideally a rectangle in two coordinates $(x,y)$ (a direct product domain). In some constructions, using hyperbolic black hole coordinates is helpful because the horizon intersects the AdS boundary in a way that makes the domain rectangular.

Then the metric ansatz should be written in a form that is closed under the coordinate transformations in $x$ and $y$, and one introduces appropriate factors so that the unknown functions satisfy boundary conditions in a manifestly consistent way. A reference metric is obtained by taking the unknown functions to simple values consistent with the same boundary conditions.

One then examines each boundary:

• the AdS conformal boundary,
• axes (if present),
• horizons (bulk and/or planar),
• and other coordinate boundaries.

For instance, at one boundary the geometry approaches a planar AdS black hole; on the conformal boundary the induced metric is essentially a Schwarzschild metric (after an appropriate coordinate transformation).

Lecture

This is also a good place to compare two approaches:

• Approach A: Fix a gauge first (i.e. fix coordinates) at the level of the ansatz, then solve the resulting reduced system.
• Approach B: Use the Einstein–DeTurck method, where you do not gauge-fix at the ansatz level, and instead add the DeTurck term so the system becomes suitable for a boundary value problem.

If you are interested in actually solving Einstein’s equations, I think many of you have probably done the following before: write an ansatz, plug it in, then look at the resulting equations. Typically some equations look simpler and some look more complicated. Then you try to form combinations that look simpler, and solve the simpler part first. That is a common style when one aims for analytic solutions.

But this method is not fully satisfactory. The key issue is that “when can you combine equations into something simple?” is not predictable. The process of combining equations itself does not have a systematic, guaranteed method — it relies on experience and tricks (the “art” of finding analytic solutions).

A feature of the DeTurck method is that at the very beginning, when you write an ansatz, you do not impose a gauge condition there. Then the Einstein equations you get are in a sense “on equal footing.” You do not need to hunt for clever combinations. In many situations, if you put an ansatz into the Einstein–DeTurck equations, you simply get the right number of equations (matching the number of unknown functions), and they are treated uniformly.

The cons:

• You have more unknown functions (because you didn’t gauge-fix early).
• The DeTurck term introduces extra contributions, so the equations become more complicated (more terms).

The pros:

• You gain a very large amount of gauge flexibility (coordinate freedom).
• You do not need to “massage” the equations into simpler forms. In fact, when solving with the DeTurck method, you can generate the equations, store them in a file, and not even look at them.

Even though the equations become long and complicated, the object you really care about — both at the start and in the end — is the metric $g_{\mu\nu}$, not the printed form of the equations. The complexity of the printed equations is not so important if you have a robust numerical pipeline.

---

For example, for a scalar field you impose boundary conditions by analyzing the asymptotics at the AdS boundary: for a second-order equation you have two linearly independent solutions, and a boundary condition means restricting the two degrees of freedom down to one.

---

After discretization, a function becomes a vector. A differential operator becomes a matrix.

For example, suppose we solve an ODE on an interval. Discretize the interval into grid points:

• the first grid point corresponds to $x=0$,
• the next grid point is the next point,
• and so on,
• the last grid point corresponds to the other boundary.

Then the values of the function at these grid points form a vector.

What about the differential operator? It becomes a matrix. How do you write the matrix? This is basic numerical analysis / computational physics knowledge: use finite difference formulas — two-point, three-point, etc. The more points you use in the stencil, the more off-diagonal components the matrix has. If you go to very high order, the matrix may stop being sparse.

But I don’t want you to get stuck in these details right now. There are standard formulas, and in many contexts you don’t need to implement them manually. You only need to understand the principle: a differential operator is represented by a matrix acting on the vector of function values.

At this point you might think: but the Einstein equations are nonlinear, not linear.

Yes. Nonlinear problems are solved by reducing them to many linear problems and iterating. This is Newton’s method. So to solve nonlinear problems, you must first understand how to solve linear problems.

In practice, in the process of solving stationary Einstein equations numerically, the most time-consuming step is solving the linear algebra problem in each Newton step. It is not that you explicitly compute the inverse matrix; rather, you solve a linear system efficiently. But still, solving that system is the bottleneck.

Now one more important analogy:

In linear algebra there is just “a matrix equation.” For differential equations, you also have boundary conditions. How do boundary conditions enter the linear algebra picture?

Look at the vector: its first component corresponds to one endpoint, and its last component corresponds to the other endpoint. Correspondingly, in the matrix equation, the first row and last row are naturally associated with boundary constraints. When you discretize a boundary value problem, you must replace the appropriate rows of the matrix equation with the discretized boundary conditions.

Once you do that, a well-defined boundary value problem becomes a well-defined linear algebra problem.

---

There is also a subtle point: an $n\times n$ matrix has $n$ eigenvalues. But the continuum differential equation typically has infinitely many eigenvalues (think of quantum mechanics with infinitely many energy levels). After discretization you only get $n$ eigenvalues, and most of them are not reliable. Only some small subset — typically the low-lying part — approximates the true continuum spectrum. You must check convergence by increasing resolution and selecting the stable eigenvalues.

The “essence” of numerical methods: discretize the differential equation, turn it into linear algebra, solve it, and check convergence.

In principle we could use Legendre polynomials, Hermite polynomials, and so on. But in practice, Chebyshev polynomials have a major advantage: the collocation grid points can be written down explicitly in closed form. This makes them extremely convenient in actual spectral computations. (We won’t go into the details right now.)

---

Let me return to the big picture again.

We said the DeTurck method is a gauge-fixing method. It is flexible, and in many cases it plays a role in turning Einstein’s equations into a system suitable for a boundary value problem — often an elliptic system, or at least close to elliptic.

But here there are subtleties:

Mathematically, on a Euclidean-signature (Riemannian) manifold, the Einstein–DeTurck equation is elliptic.

Physically, our spacetimes are Lorentzian. Then the Einstein–DeTurck equation is not automatically elliptic. You need to consider cases.

• If the problem is static, you can analytically continue $t\to i\tau$. Then the problem becomes Riemannian, and you get an elliptic system.
• If the problem is stationary but not static, then in many cases (because we are not doing time evolution) the Einstein–DeTurck system is still elliptic in the relevant sense.
• But in some situations, such as geometries with heat flow or non-Killing horizons, the system is not elliptic; it can become mixed-type.

Even if it is not elliptic in a strict mathematical sense, if the discretized linear algebra problem (the Jacobian in Newton’s method) is still well-posed/invertible, then the method can still be usable.

---

Now you can compare the “numerical relativity” we are discussing now with numerical relativity in the last century.

In the last century, when people said “numerical relativity,” they usually meant time evolution, because stationary problems had little to solve (due to uniqueness theorems and fewer known families). But now, in AdS (for example in AdS$_4$), there are extremely rich families of stationary solutions — many different black objects, many branches, many phases.

Traditional numerical relativity emphasized the $3+1$ decomposition: splitting Einstein equations into:

• evolution equations (hyperbolic),
• and constraint equations (elliptic).

Here, with the Einstein–DeTurck approach for stationary problems, the constraints are essentially “promoted” into the same footing as the dynamical equations. You solve as many independent second-order equations as there are unknown metric components, and you no longer have the same constraint/evolution split.

Also, time evolution problems and stationary boundary value problems differ significantly in techniques. Next week perhaps another lecturer will talk about time evolution methods; those have their own specialized tools.

They are complementary. Stationary problems often require more flexible coordinate choices; time evolution typically requires you to commit to a particular coordinate/gauge setup early, which limits flexibility. So time evolution cannot solve all stationary problems. Conversely, stationary methods cannot solve time evolution.

With the DeTurck formulation, stationary problems can often be turned into a relatively systematic pipeline:

• write the ansatz,
• generate the equations,
• impose boundary conditions,
• solve.

---