1 Introduction↩︎

Finding explicit Einstein metrics is a central problem in geometry and theories of gravity. More specifically, in general relativity and its generalisations, one is particularly interested in the construction of pseudo-Riemannian Lorentzian metrics. Let \(\mathcal{M}\) be a smooth manifold equipped with a local coordinate chart \((U, x^\mu)\) on an open set \(U\subset \mathcal{M}\). Given \(U\), one may locally construct a smooth, pseudo-Riemannian metric \(g_{\mu\nu}\) on \(U\) with an associated Ricci curvature tensor \(R_{\mu\nu}\). In this notation, an Einstein metric is defined such that Ricci curvature is everywhere proportional to the metric itself, i.e.¹ \[R_{\mu\nu}=\lambda g_{\mu\nu}, \label{eq:einstein95metric}\tag{1}\] for some constant of proportionality \(\lambda\) referred to as the Einstein constant. In Riemannian geometry, 1 says that the space has the same averaged curvature in every tangent direction, whilst in general relativity, its pseudo-Riemannian version is the vacuum Einstein equation with \(\lambda\) interpreted as a cosmological constant. Although 1 is concise, it represents a non-linear second-order system for the metric components. The global problem also requires the local chart components to assemble into a single smooth tensor on the whole manifold.

In numerical relativity and stationary black-hole construction, this non-linear geometric problem is usually made tractable by turning the Einstein equations into a well-posed computational system after choosing a formulation, coordinates, gauge conditions, boundary data, and a mesh or spectral discretisation. The dynamical branch of the field is built around the Cauchy problem, where the spacetime is foliated into spatial hypersurfaces, initial data are chosen subject to the Hamiltonian and momentum constraints, and the data are evolved by a gauge-fixed form of Einstein’s equations. The ADM decomposition is the classical starting point for this viewpoint [1], while standard reviews of initial data and numerical-relativity methodology emphasise that the practical problem is not only to discretise the equations, but also to control constraints, gauges, singularities, outer boundaries, and wave extraction, see [2], [3] for instance. A major development was the conformal BSSN formulation, originating in the work of [4] and systematised by [5], which greatly improved the stability of three-dimensional evolutions relative to direct ADM evolution. The binary-black-hole breakthrough was obtained in [6] via the first stable full evolution of a binary black-hole spacetime using generalized harmonic evolution and excision, while the moving-puncture methods of [7], [8] made long-term inspiral, merger, ringdown, and gravitational-wave extraction broadly practical. Subsequent reviews describe this as the point at which binary black-hole merger simulations became routine tools for strong-field gravity and gravitational-wave astronomy [9], with high-accuracy spectral evolutions and community infrastructure such as the Einstein Toolkit further maturing the field [10], [11]. Stationary black-hole construction is different in character, instead of evolving initial data one typically solves a coupled elliptic boundary-value problem for the metric, with boundary conditions encoding asymptotic regions, axes, horizons, conformal boundaries, or compact directions. The construction of [12], for static axisymmetric vacuum solutions and non-uniform black strings, illustrates the role of adapted ansatz choices and elliptic relaxation methods. Similarly, the Einstein–DeTurck formulation of [13] gives a systematic gauge-fixed elliptic framework for static numerical relativity and Kaluza–Klein black holes. This approach was extended toward more general stationary vacuum black holes in Lorentzian signature in [14]. This stationary toolkit is reviewed in detail by [15], [16] and has been applied in less symmetric and asymptotically AdS settings, for example in numerical constructions of black resonators and related AdS instabilities. In all of these cases, the metric is represented by finite-dimensional data tied to a coordinate grid, finite-difference stencil, finite-element mesh, or spectral basis, and geometric consistency is enforced through the discretised field equations together with carefully chosen gauge, regularity, and boundary conditions.

Machine learning provides a complementary representation: the metric is treated as a differentiable function of the coordinates by representing it with a (smooth) neural network, with geometric tensors computed by automatic differentiation and inserted directly into the objective². This program of repurposing the efficient differentiation mechanisms of AI for differential geometry and geometric analysis has been recently pursued in [18], [19] for various problems concerning spheres, as well as by pioneering works in the context of Calabi-Yau metrics [20]–[24], \(G_2\) metrics [25], [26], hyperbolic spaces [27], [28], modes related to black holes [29]–[31], and for broader contexts³ such as minimal surfaces [28], [35]–[37].

AInstein⁴ [18] introduced this strategy for Riemannian Einstein metrics on spheres. In the original construction, the manifold was represented by an atlas, each chart had its own subnetwork, points in overlap regions were passed through the relevant transition maps, and an explicit overlap loss enforced the tensorial transformation law for the metric. In this way, global consistency became a learnable constraint rather than an assumption built into the architecture.

This work adapts AInstein to first find and approximate the maximally extended Schwarzschild solution, then presents a more general Lorentzian metric search. The main changes are a Lorentzian metric parametrisation, compactified Penrose-domain sampling, coordinate-invariant scalar losses, and an embedded treatment of the topologically non-trivial part of the manifold. The embedding is not only a visual convenience; it changes how global consistency is imposed. Rather than learning two independent sets of metric components and matching them with a loss term, the model learns a single ambient metric, with the two stereographic hemispheres being related analytically by pullback.

2 Theoretical Prerequisites↩︎

2.1 The Schwarzschild Problem↩︎

The Schwarzschild metric is an exact solution to Einstein’s field equations that describes the spacetime geometry outside a spherically-symmetric non-rotating mass such as a black hole. In this section, the coordinate description used throughout the work is developed.

Let \(\mathcal{P}\) denote the two-dimensional compactified Penrose domain [38], i.e. an open hexagonal subset of \(\mathbb{R}^2\) (see below), and let \(S^{2}\) denote the 2-sphere. The four-dimensional manifold is \[\mathcal{M}=\mathcal{P}\times S^{2}, \label{eq:model95manifold}\tag{2}\] where it should be noted 2 does not assume a product metric – the product notation only separates the causal coordinates from the angular coordinates. Although the Schwarzschild solution is spherically symmetric, the class of metrics considered later is a general symmetric four-dimensional pullback metric, not an ansatz restricted to the diagonal Schwarzschild form.

Let \((T,X)\) denote Penrose coordinates. A Penrose diagram is a conformally compactified representation of spacetime, infinite regions are mapped to finite boundaries while the light-cone structure is preserved. Penrose coordinates are useful because they place the event horizons inside a finite computational domain; the angular \(S^{2}\) directions are suppressed, so each point in the two-dimensional diagram represents a sphere, and radial light rays follow diagonal lines. Let \(m>0\) denote the Schwarzschild mass parameter and let \(r\) denote the Schwarzschild areal radius. In Schwarzschild coordinates, the horizon at \(r=2m\) is a coordinate singularity; in Kruskal or Penrose coordinates, it is a regular null surface. For any centre \(c\), define the two boundary functions \[\tau_{c}^{-}(X)=|X-c|-\frac{\pi}{4}, \qquad \tau_{c}^{+}(X)=-|X-c|+\frac{\pi}{4}. \label{eq:tau95boundaries}\tag{3}\]

1 shows the corresponding compactified domain. Using 3 , the four exterior/interior diamonds we consider are \[\label{eq:penrose95regions} \begin{align} B_{\rm I} &=\{0<X<\pi/2,\;\tau_{\pi/4}^{-}<T<\tau_{\pi/4}^{+}\}, \\ B_{\rm II} &=\{-\pi/4<X<\pi/4,\;|X|<T<\pi/4\}, \\ B_{\rm III} &=\{-\pi/2<X<0,\;\tau_{-\pi/4}^{-}<T<\tau_{-\pi/4}^{+}\},\\ B_{\rm IV} &=\{-\pi/4<X<\pi/4,\;-\pi/4<T<-|X|\}. \end{align}\tag{4}\] The null horizons are the diagonal lines \(T=\pm X\) in the central diamond, whilst the future and past singularities are the horizontal boundaries \(T=\pm\pi/4\) with \(|X|<\pi/4\). The remaining null sides are the conformal boundaries of the exterior regions.

In the Penrose coordinates of 4 , the target Schwarzschild line element may be written as \[\mathrm{d}s^{2} =F(T,X)(-\mathrm{d}T^{2}+\mathrm{d}X^{2}) +r(T,X)^{2}\mathrm{d}\Omega_{2}^{2}, \label{eq:schw95metric}\tag{5}\] where \(\mathrm{d}\Omega_{2}^{2}\) is the round metric on the unit two-sphere, \(r(T,X)\) is the Schwarzschild areal radius and \(F(T,X)\) is the Penrose conformal factor. This conformal factor is \[F(T,X)= \frac{32m^{3}e^{-r(T,X)/(2m)}}{r(T,X)\left(\cos^{2}T-\sin^{2}X\right)^{2}} \, , \label{eq:F95factor}\tag{6}\] with \(r(T,X)\) defined below. For stereographic angular coordinates \((q_{1},q_{2})\) on \(S^{2}\), \[\mathrm{d}\Omega_{2}^{2} =\frac{4}{(1+q_{1}^{2}+q_{2}^{2})^{2}} (\mathrm{d}q_{1}^{2}+\mathrm{d}q_{2}^{2}). \label{eq:stereo95s2}\tag{7}\] The conformal factor \(F\) is finite and non-zero at the horizons; the singular behaviour of the Schwarzschild solution is instead carried by \(r(T,X)\rightarrow0\) on the spacelike singularity (in addition to the conformal boundary divergence). Thus, the singular behaviour is restricted to the boundary of the Penrose domain. Let \(W\) denote the real Lambert-\(W\) function, and let \(\mathcal{B}\) label one of the four regions in 4 . The radius \(r\) is computed over \(B_{\rm I},\ldots,B_{\rm IV}\) by \[r=2m(1+W_{\mathcal{B}}), \label{eq:radius95lambert}\tag{8}\] where \(\chi=\cos 2T/\cos 2X\) and \[\begin{align} W_{\mathcal{B}} &= W\!\left[\exp\{-2\,\operatorname{arccoth}\chi-1\}\right], \qquad\qquad\qquad\quad \mathcal{B}=B_{\rm I},B_{\rm III},\\ W_{\mathcal{B}} &= W\!\left[-\exp\{-2\,\operatorname{arctanh}\chi-1\}\right], \qquad\qquad\qquad \mathcal{B}=B_{\rm II},B_{\rm IV}. \label{eq:piecewise95W} \end{align}\tag{9}\] For a detailed account on these coordinates the reader is referred to [39], for instance. The tensor equation to be solved remains the vacuum condition \(R_{\mu\nu}=0\).

2.2 Birkhoff’s Theorem and Scalar Invariants↩︎

This section collects the invariant statements that characterise the Schwarzschild and Type-I geometries. The key point here is not that any single scalar uniquely determines the metric in general, but rather the specification of increasing numbers of invariants constrains the space of possible solutions.

2.2.1 Petrov Classes, Weyl Invariants, and Curvature Scales↩︎

Any smooth region of a 4D classical black-hole spacetime has a local Petrov type. A classifier of the type may be defined through the complex speciality index ([40], [41]). Let \(C_{abcd}\) denote the Weyl tensor and define the self-dual Weyl tensor \[\mathcal{C}_{abcd} := C_{abcd} - i\,{}^{\star}C_{abcd} \qquad \text{where} \qquad {}^{\star}C_{abcd} := \frac{1}{2}\epsilon_{ab}{}^{ef}C_{efcd}. \label{eq:self95dual95weyl95petrov}\tag{10}\] With our normalisation, the complex Weyl invariants are \[I := \frac{1}{32}\mathcal{C}_{abcd}\mathcal{C}^{abcd} \qquad \text{and} \qquad J := \frac{1}{384} \mathcal{C}_{ab}{}^{cd} \mathcal{C}_{cd}{}^{ef} \mathcal{C}_{ef}{}^{ab}. \label{eq:weyl95petrov95invariants}\tag{11}\] The corresponding speciality index is \[\mathcal{S} := \frac{27J^{2}}{I^{3}}, \label{eq:speciality95index}\tag{12}\] following the invariant normalisation used by [40], and is defined away from the conformally flat locus \(I=0\). In this convention, an algebraically special Kerr background has \(\mathcal{S}=1\) ([40]). Schwarzschild is the \(a=0\) Kerr limit; in the Schwarzschild principal tetrad the only non-zero Newman–Penrose Weyl scalar is real, \(\Psi_{2}=-m/r^{3}\) (see [42]). Therefore \[I=3\Psi_{2}^{2}=\frac{3m^2}{r^6}, \qquad J=-\Psi_{2}^{3}=\frac{m^{3}}{r^{9}}, \qquad \mathcal{S}_{\rm Schw}=1+0i . \label{eq:schwarzschild95weyl95speciality}\tag{13}\] Equivalently, the algebraically special Petrov condition is measured by the discriminant \[\Delta_{\rm Petrov} := I^{3}-27J^{2}. \label{eq:petrov95discriminant}\tag{14}\] Schwarzschild is Petrov type \(\mathrm{D}\), so the invariant target is \(\mathcal{S}=1\), or equivalently \(\Delta_{\rm Petrov}=0\), at every non-flat point. The implication is not reversible: \(\Delta_{\rm Petrov}=0\) is the algebraically special condition and also includes type \(\mathrm{II}\), while the more degenerate types \(\mathrm{III}\), \(\mathrm{N}\), and \(\mathrm{O}\) have \(I=J=0\) and hence a discriminant with an undefined speciality index [41]. Petrov type I is associated with a non-zero discriminant. Thus the discriminant is a test for repeated principal null directions, not a complete type-\(\mathrm{D}\) classifier, yet corroborates type D and can be used to determine the solution is not type D. The same invariant, \(\mathcal{S}=27J^{2}/I^{3}\), is also used in later local Petrov-classification diagnostics ([43]).

The Kretschmann invariant fixes the curvature scale of Schwarzschild giving partial information on the coordinate representation. We write \[K := R_{abcd}R^{abcd} \label{eq:kretschmann95definition}\tag{15}\] for the quadratic curvature scalar. For the Schwarzschild solution, \[K_{\rm Schw}=\frac{48m^{2}}{r^{6}}, \label{eq:schwarzschild95kretschmann95scale}\tag{16}\] where \(r\) is the areal radius [42]. Thus the dimensionless product \[\frac{Kr^{6}}{m^2}=48 \label{eq:schwarzschild95kr6}\tag{17}\] is constant on the non-singular Schwarzschild region. This provides an invariant radial scale against which the learnt metric can be checked.

One may also employ a useful Schwarzschild-specific cubic curvature ratio. Let \[J_{\rm cub} := R_{ab}{}^{cd}R_{cd}{}^{ef}R_{ef}{}^{ab}. \label{eq:cubic95riemann}\tag{18}\] With the Riemann-sign convention fixed throughout this work, the Schwarzschild solution has \(J_{\rm cub}=96m^{3}/r^{9}\), using [44]. In the Weyl-invariant normalisation of 11 , this is equivalent to \(J_{\rm cub}=96\,\mathrm{Re}(J)\) on Schwarzschild, allowing definition of the dimensionless reference scalar \[\rho_{\rm Schw} := -\frac{96\,\mathrm{Re}(J)}{|K|^{3/2}} = -\frac{J_{\rm cub}}{|K|^{3/2}} = -\frac{1}{\sqrt{12}}, \label{eq:schwarzschild95cubic95ratio}\tag{19}\] everywhere the metric is defined, with \(|\rho_{\rm Schw}|=1/\sqrt{12}\); further details are provided in 6. The value follows from the chosen normalisation of the diagnostic \(\rho\); the literature input is the standard Schwarzschild Weyl and Kretschmann invariants [42], not the diagnostic name itself. This ratio is simply an additional scalar check of the Schwarzschild curvature pattern, independent of the (areal) radius.

2.2.2 Killing Vectors and Spherical Symmetry↩︎

The analytic Schwarzschild solution is spherically symmetric. The relevant infinitesimal symmetries are the three rotations of the round two-sphere. Equivalently, the Killing algebra of the round \(S^{2}\) is \(\mathfrak{so}(3)\), generated by the infinitesimal rotations inherited from the ambient Euclidean space ([45]). In the Cartesian embedding \(S^{2}\subset\mathbb{R}^{3}\) these are the familiar angular-momentum vector fields \[J_x=-z\partial_y+y\partial_z,\qquad J_y=z\partial_x-x\partial_z,\qquad J_z=-y\partial_x+x\partial_y . \label{eq:cartesian95rotation95generators}\tag{20}\] To obtain the coordinate formulae used here, one may write the northern inverse stereographic map as \[x=\frac{2q_{1}}{1+|q|^{2}}, \qquad y=\frac{2q_{2}}{1+|q|^{2}}, \qquad z=\frac{1-|q|^{2}}{1+|q|^{2}} \, , \label{eq:north95stereo95rotation95pullback}\tag{21}\] which yields the following generators in coordinates: \[\begin{align} \xi_{z}&=-q_{2}\,\partial_{q_{1}}+q_{1}\,\partial_{q_{2}}, \nonumber\\ \xi_{x}&=-q_{1}q_{2}\,\partial_{q_{1}} -\frac{1}{2}(1-q_{1}^{2}+q_{2}^{2})\,\partial_{q_{2}}, \nonumber\\ \xi_{y}&=\frac{1}{2}(1+q_{1}^{2}-q_{2}^{2})\,\partial_{q_{1}} +q_{1}q_{2}\,\partial_{q_{2}}. \label{eq:s295killing95vectors} \end{align}\tag{22}\] The corresponding formulae on the southern chart differ only by the expected stereographic sign changes for two of the generators. A metric is invariant under these rotations precisely when \[(\mathcal{L}_{\xi_A}g)_{\mu\nu} = \xi_A^{\rho}\partial_{\rho}g_{\mu\nu} +g_{\rho\nu}\partial_{\mu}\xi_A^{\rho} +g_{\mu\rho}\partial_{\nu}\xi_A^{\rho} =0, \qquad A=x,y,z. \label{eq:killing95equation95s2}\tag{23}\] The Killing equations are the infinitesimal statement of spherical symmetry. In the numerical workflow below, the same operators are used both as a training loss and as validation residuals. This choice makes the Birkhoff–Jebsen hypothesis an explicit part of the Schwarzschild targeted search rather than a purely posterior observation. Equivalently, one may construct and implement any other symmetry property in this manner to perform other targeted searches.

The Birkhoff–Jebsen theorem states that a four-dimensional spherically symmetric vacuum solution is locally Schwarzschild; for \(\lambda\neq0\) the corresponding solution is Kottler/Schwarzschild–de Sitter ([46], [47]). Israel’s theorem gives a complementary static black-hole result, under the standard assumptions of asymptotic flatness, regular horizon, vacuum equations, and horizon topology, the regular static vacuum black hole is Schwarzschild, as discussed in [48]. These criteria are intentionally redundant. Ricci flatness and spherical symmetry already place the metric on the local Schwarzschild branch by the Birkhoff–Jebsen theorem, while non-degeneracy excludes collapsed metric limits and the Kretschmann scale fixes the mass parameter. Since each condition is stated tensorially or through scalar curvature invariants, the combined identification is insensitive to the coordinates used to represent the metric – it is a diffeomorphism-invariant identification of the geometry as Schwarzschild.

2.3 Spacetime Embedding↩︎

The \(S^2\) factor is represented by two stereographic hemispherical charts. This two-chart representation is only a coordinate device. Let \[w=(T,X,q_{1},q_{2},p), \label{eq:intrinsic95input}\tag{24}\] be an intrinsic coordinate together with a chart label, where \(p=0\) denotes the northern chart and \(p=1\) denotes the southern chart. We write \(q=(q_{1},q_{2})\), \(|q|^{2}=q_{1}^{2}+q_{2}^{2}\), and define \(s_{0}=+1\), \(s_{1}=-1\), allowing the embedding to be written as \[W(w)= \left(T,X,\frac{2q_{1}}{1+|q|^{2}}, \frac{2q_{2}}{1+|q|^{2}}, s_{p}\frac{1-|q|^{2}}{1+|q|^{2}}\right) \, . \label{eq:embedding}\tag{25}\] 25 maps the four intrinsic coordinates in 24 to \(\mathbb{R}^{2}\times S^{2}\subset\mathbb{R}^{5}\). The first two coordinates are unchanged, while the final three are the inverse stereographic projection of the angular point.

The embedding Jacobian \(J^{A}{}_{\mu}=\partial W^{A}/\partial w^{\mu}\) reads \[J= \begin{pmatrix} 1&0\\ 0&1 \end{pmatrix} \oplus \frac{1}{(1+q_{1}^{2}+q_{2}^{2})^{2}} \begin{pmatrix} 2(1+q_{2}^{2}-q_{1}^{2}) & -4q_{1}q_{2}\\ -4q_{1}q_{2} & 2(1+q_{1}^{2}-q_{2}^{2})\\ -4s_{p}q_{1} & -4s_{p}q_{2} \end{pmatrix}. \label{eq:sphere95embedding95block}\tag{26}\] The pullback of the flat Euclidean metric on the last three coordinates through 25 gives 7 . The sign \(s_{p}\) changes which pole is represented by \(q=0\), while the induced round metric is the same in both charts. This equality of the induced metric is the mathematical reason the two hemispheres can be used as one embedded angular manifold.

All tensorial calculations on \(\mathcal{P}\times S^{2}\) are taken with respect to \((T,X,q_{1},q_{2}) = w\). The two sets of coordinates therefore define a single embedded sphere, rather than two a-priori independent charts (which was the case in [18]).

2.4 Trapped Surfaces↩︎

The invariants of 2.2 characterise the Schwarzschild solution among vacuum geometries, but the defining property of a black hole is causal rather than curvature-based: the presence of a horizon, and locally of trapped surfaces. The event horizon is a global, teleological (physical) object, the boundary of the causal past of future null infinity, and cannot be evaluated pointwise on a bounded domain. The quasi-local notion of a trapped surface [49]–[51] is therefore used, which is local and requires no symmetry.

Let \(\Sigma\) be a closed spacelike two-surface with the two future-directed null normals \(\ell,n\) (normalised by \(g(\ell,n)=-1\)), and let \(q_{\mu\nu}=g_{\mu\nu}+\ell_\mu n_\nu+n_\mu\ell_\nu\) project onto \(\Sigma\). The null expansions \[\theta_\ell=q^{\mu\nu}\nabla_\mu\ell_\nu, \qquad \theta_n=q^{\mu\nu}\nabla_\mu n_\nu, \label{eq:null95expansions}\tag{27}\] give the fractional rate of change of the area element of \(\Sigma\) along the two light fronts emanating from it. In a normal region the outgoing family diverges and the ingoing family converges, \(\theta_\ell>0>\theta_n\). The surface is (future) trapped when both expansions are negative, \(\theta_\ell,\theta_n<0\), so that even the nominally outgoing light is dragged inward; \(\theta_\ell=0\) is the marginal (apparent-horizon) case, and \(\theta_\ell,\theta_n>0\) is the time-reversed, anti-trapped (white-hole) case. By the Penrose singularity theorem ([49]) a closed trapped surface, together with the null energy condition (saturated in vacuum), forces geodesic incompleteness and lies within the black-hole region; trapped surfaces are the standard local certificate of a black-hole interior.

The manifold 2 supplies a natural family of such surfaces: the \(S^2\) fibres at fixed \((T,X)\). Their induced metric is exactly the angular block of the pullback 32 , with no symmetry assumed, so each fibre has a well-defined area \(A=4\pi R^{2}\) and areal radius \(R\). It is convenient to encode the trapping in the single scalar \[\Xi:=g^{\mu\nu}\,\partial_\mu R\,\partial_\nu R, \qquad \operatorname{sign}\Xi=-\operatorname{sign}(\theta_\ell\theta_n), \label{eq:trapped95scalar}\tag{28}\] as the squared norm of the areal-radius gradient, where \(\Xi>0\) (spacelike \(\nabla R\)) is untrapped, \(\Xi=0\) is marginal, and \(\Xi<0\) (timelike \(\nabla R\)) is trapped or anti-trapped, the two being distinguished by the common sign of the expansions. In spherical symmetry \(\Xi\) reduces to the Misner–Sharp combination \(1-2M_{\rm MS}/R\) of 54, which for Schwarzschild equals \(1-2m/r\) and is negative precisely on the interior \(r<2m\).

Two remarks matter for the non-spherically-symmetric metrics learnt in 4.4. First, the Misner–Sharp mass is only invariantly defined under spherical symmetry and is unavailable in general, whereas the trapped condition \(\theta_\ell,\theta_n<0\) is symmetry-free and remains the operative definition. Second, \(\Xi\) is a smooth scalar built from the metric and its first derivatives, hence a natural differentiable training signal, while the individual expansions \(\theta_\ell,\theta_n\) – which alone distinguish a black hole from its time reverse – are reserved for the posterior certification. The identification 28 is exact for round fibres and is corrected only by their departure from roundness, which we find to be small.

3 AInstein Embedding Neural Networks↩︎

3.1 Warm up – Local Lorentzian Schwarzschild↩︎

Before considering the embedded Schwarzschild problem, we study a local two-dimensional Lorentzian simplification. In this mode the manifold is a single coordinate ball, the dimension is set to \(d=2\), and no \(S^{2}\) embedding is used. The network therefore represents the metric directly on the local coordinates \(u=(u^{0},u^{1})\). Its outputs are three independent entries of a lower-triangular \(2\times2\) matrix \(L_{\theta}(u)\), and the Lorentzian metric is obtained from an indefinite-Cholesky parametrisation, \[g^{(\theta)}_{\mu\nu}(u) = \left(L_{\theta}(u)\eta_{(2)}L_{\theta}(u)^{\mathsf{T}}\right)_{\mu\nu}, \qquad \eta_{(2)}=\operatorname{diag}(-1,1), \qquad \mu,\nu=0,1. \label{eq:local95lorentz95cholesky}\tag{29}\] This local architecture isolates the Lorentzian signature and curvature calculation from the additional global structure of \(\mathcal{P}\times S^{2}\) later considered.

The network has a sets of four 128 neuron-wide linear layers, and employs GELU activations; GELU allows derivatives to be \(C^2\) smooth, as needed for metric derivatives. To mitigate confusion, it is emphasised that the network is trained such that all Schwarzschild losses later defined in 3.4.3 are disabled. Only this Lorentzian Einstein condition is used for each of \(\lambda \in \{-1,0,1\}\).

3.2 Network Architecture↩︎

For the central, global, architecture, let \(\theta\) denote the trainable weights and biases of the numerical model. The neural architecture is built using an internal multilayer perceptron \[\mathcal{N}_{\theta}:\mathbb{R}^{5}\longrightarrow\mathbb{R}^{15}. \label{eq:network95map}\tag{30}\] This map in 30 takes the embedded point \(W(w)=(T,X,x(q_1, q_2),y(q_1, q_2),z(q_1, q_2))\), where \((x,y,z)\) are the Cartesian sphere coordinates implicitly defined in 25 , and outputs the independent entries of a lower-triangular \(5\times5\) matrix \(L_{\theta}\). The use of this lower-triangular matrix enforces the symmetric property of the metric under the below combination meaning there are fewer components to learn. The ambient Lorentzian metric is \[G_{\theta}=L_{\theta}\eta L_{\theta}^{\mathsf{T}}, \qquad \eta=\operatorname{diag}(-1,1,1,1,1). \label{eq:lorentz95cholesky}\tag{31}\] The four-dimensional metric used in the loss is the pullback \[g^{(\theta)}_{\mu\nu}(w) = (G_{\theta})_{AB}(W(w)) J^{A}{}_{\mu}(w)J^{B}{}_{\nu}(w) \, . \label{eq:pullback95metric}\tag{32}\] Here \(\mu,\nu=0,\ldots,3\) and \(A,B=0,\ldots,4\). The biases of the architecture’s output layer are initialised so that the model is initialised near the flat ambient metric \(\eta\). In the limit that the standard deviation of the sample distribution tends to 0, one recovers the flat Minkowski space metric.

As mentioned, this is where the embedded sphere changes the original AInstein construction. In [18], separate chart metrics were matched by an overlap loss. Here both stereographic hemispheres are obtained by the push-forward of the same ambient metric through 25 , so their agreement is built into the parametrisation and no angular overlap loss is needed; any learnt metric is automatically globally consistent on the manifold.

3.3 Data Generation↩︎

To sample a general point on the \(\mathbb{R}^2 \times S^2\) manifold, the Penrose sampler first draws a point in a 2-dimensional disc and then maps it to the compactified Penrose domain. This approach is more efficient than rejection sampling in the polygonal Penrose diagram, and gives direct control over the radial density of a sample.

Let \(R_{\mathcal{P}}\) be the maximum radius used in the auxiliary unit disc. We draw a disc point \(z=\rho u\), where \(0\leq\rho\leq R_{\mathcal{P}}\leq1\), \(u=(u_{T},u_{X})\) is a unit vector, and \(|u|=1\). In the implementation \[\begin{align} \rho&=R_{\mathcal{P}}\,b, \nonumber\\ b&\sim{\rm Beta}(\alpha_{\mathcal{P}}^{-1},\alpha_{\mathcal{P}}), \qquad \arg u\sim{\rm Unif}[0,2\pi), \label{eq:disc95sampling} \end{align}\tag{33}\] where \(\alpha_{\mathcal{P}}\) is the density-power hyperparameter. Define the offset distance \[\delta(\rho)=\frac{\pi}{4}(1-\rho) \label{eq:penrose95offset}\tag{34}\] and the three radial boundary scales \[\begin{align} A(u)&=\frac{\pi/4-\delta}{|u_{T}|},\\ B(u)&=\frac{\pi/2-\sqrt2\,\delta}{|u_{T}+u_{X}|},\\ C(u)&=\frac{\pi/2-\sqrt2\,\delta}{|u_{X}-u_{T}|}. \label{eq:boundary95scales} \end{align}\tag{35}\] The sampled Penrose point is obtained using the directional scale \(\ell(u)\), \[(T,X)=\ell(u)\,u, \qquad \ell(u)=\min\{A(u),B(u),C(u)\}, \label{eq:disc95to95penrose}\tag{36}\] with the limiting prescription implied by 35 when a denominator vanishes. A representative sample of 2000 points for \(R_\rho=0.85\) is shown in Figure 2, where the ill-defined nature of the singularity at the boundaries motivates this use of a \(R_\rho < 1\). At \(\rho=1\), 36 reaches the full diagram boundary; at \(\rho=0\), it collapses to the centre.

3.4 Training Objectives↩︎

The first step of training a neural architecture is the specification of a training objective in the form of a loss function. In this section, the training objectives for supervised training against the known target solutions, and then four related problems are outlined. First, local two-dimensional Lorentzian Einstein runs are used as a controlled signature test. We then consider the task of constructing a simulacrum of the classic Schwarzschild black hole, before introducing two general searches for Petrov type I solutions: the first a search agnostic to solution style, and the second specifically searching for black hole trapped surfaces.

3.4.1 Supervised Network↩︎

The supervised model, where training instead approximates a known construction directly, is used for initialisation and benchmarking. The base network still outputs the ambient Cholesky vector, while a wrapper applies 32 and returns the flattened intrinsic metric. For training samples \(w_{a}\) and target metric \(g^{\rm tar}_{\mu\nu}\), the default supervised loss is the componentwise relative mean-squared error \[\mathcal{L}_{\rm sup} = \frac{1}{16N}\sum_{a=1}^{N} \sum_{\mu,\nu} \frac{ \left(g_{\mu\nu}^{(\theta)}(w_{a})-g_{\mu\nu}^{\rm tar}(w_{a})\right)^{2}}{\left(g_{\mu\nu}^{\rm tar}(w_{a})\right)^{2}+1}. \label{eq:supervised95loss}\tag{37}\] In these investigations, the target metric is either the Lorentzian identity metric on \(\mathbb{R}^{2}\times \mathbb{R}^{3}\) (useful for a neutral initialisation), or the analytic Schwarzschild metric in 5 (useful for benchmarking the general unsupervised search performance).

3.4.2 Unsupervised Local Lorentzian Network↩︎

For the local two-dimensional calibration runs, for each training sample \(w_{a}\), define an Einstein error tensor \[E_{\mu\nu}(w_{a}) := R_{\mu\nu}[g^{(\theta)}](w_{a}) -\lambda g^{(\theta)}_{\mu\nu}(w_{a}), \label{eq:einstein95residual}\tag{38}\] which is used to define the Einstein loss \[\label{eq:einstein95loss95local} \mathcal{L}_{E, \text{Local}} := \frac{1}{N} \sum_{a=1}^N \bigg|E_{\mu\nu}\big(g_E^{(\theta)}\big)^{\mu\alpha}\big(g_E^{(\theta)} \big)^{\nu\beta}E_{\alpha\beta} \bigg|_{w_a}\tag{39}\] where \(g_E\) is the Euclideanised variant of the metric⁵. The network described in 3.1 is then trained with this \(\mathcal{L}_{E, \text{Local}}\) loss only, for \(\lambda\in\{+1,0,-1\}\).

3.4.3 Unsupervised Schwarzschild Network↩︎

For the four-dimensional Schwarzschild and general Petrov type searches, the vacuum equation is \(R_{\mu\nu}=0\), so \(\lambda=0\) is set⁶. The default Schwarzschild Lorentzian Einstein loss is then defined as the inverse-metric contraction \[\mathcal{L}_{E,\text{Schw}} := \frac{1}{N}\sum_{a=1}^{N} \frac{1}{\text{min}(|K|, K_\text{cap}) + \epsilon}\left| E_{\mu\nu}\big(g^{(\theta)}_E\big)^{\mu\alpha}\big(g^{(\theta)}_E \big)^{\nu\beta}E_{\alpha\beta} \right|_{w_{a}}, \label{eq:contracted95einstein95loss95sch}\tag{40}\] where once more \(g_E\) is the Euclideanised version of \(g\), producing a well-defined scalar, since all objects involved are tensorial (as opposed to using MSE). The numerical prefactor⁷ involves a local scalar weighting which uses \(K_\text{cap} = \kappa K_\text{hor}\) to cap this factor to the horizon value \(K_\text{hor}=0.75\). The numerator of 40 is simply the most basic statement of the Einstein constraint which is consistent between all problems (and indeed [18]).

Secondly, the curvature-scale loss is imposed through the quadratic Weyl invariant \(I\), rather than directly through the full Kretschmann scalar. For a Ricci-flat metric the Riemann and Weyl tensors agree, but away from the vacuum locus the Kretschmann scalar also contains Ricci contributions. Using \(I\) therefore decouples the scale-fixing term from the Einstein residual: \(\mathcal{L}_{\rm E}\) drives the Ricci tensor to zero, while the Weyl term fixes the non-Ricci curvature scale of the Schwarzschild branch.

With the normalisation in 11 , the exact Schwarzschild solution obeys \[|I_{\rm Schw}|=\frac{3m^{2}}{r^{6}} \quad \implies \quad \sqrt{|I_{\rm Schw}|}\,r^{3}=\sqrt{3}\,m . \label{eq:schwarzschild95weyl95scale}\tag{41}\] For each sampled point, the implemented Weyl loss is therefore \[\mathcal{L}_{W} := \frac{1}{N}\sum_{a=1}^{N} \left( r(w_a)^{3} \sqrt{|I_\theta(w_a)|+\epsilon} - \sqrt{3}\,m \right)^{2}. \label{eq:weyl95invariant95loss}\tag{42}\] Here \(r(w_a)\) is the analytic areal radius of the sampled Penrose point, and \(I_\theta\) is computed from the self-dual Weyl tensor of the learnt metric. The square-root form removes the six-power radial dynamic range of the raw invariant, while still excluding the flat \(I=0\) branch.

The Birkhoff-motivated symmetry term is \[\mathcal{L}_{SO(3)} := \frac{1}{3N} \sum_{A=x,y,z}\sum_{a=1}^{N} \frac{ \left| (\mathcal{L}_{\xi_A}g)_{\mu\nu} (g_E)^{\mu\alpha}(g_E)^{\nu\beta} (\mathcal{L}_{\xi_A}g)_{\alpha\beta} \right|_{w_a}}{ \left| g_{\mu\nu}(g_E)^{\mu\alpha}(g_E)^{\nu\beta}g_{\alpha\beta} \right|_{w_a} +\epsilon}, \label{eq:symmetry95operator95loss}\tag{43}\] where the three operators \(\xi_A\) are the rotational fields in 22 , \(\mathcal{L}_{\xi_A}\) are the usual symmetry generators, and \(g_E\) is the same Euclideanised metric used in 40 . In the implementation this term is evaluated by a small finite-difference pullback along each rotational flow. Its role is to impose the spherical-symmetry assumption needed for the Birkhoff–Jebsen identification, which it does by measuring the metric-contracted Killing-vector action along each basis direction relative to the metric-contracted size of \(g\). This action should be zero in each case according to 23 . Together with the vacuum Einstein loss, this loss selects the Schwarzschild branch locally, while the Weyl curvature-scale target fixes scaling over the Penrose diagram to ensure it is fully covered and not a more trivial diffeomorphism.

Finally, in order to aid numerical stability, it is necessary to add a ‘determinant loss’, \(\mathcal{L}_{\det}\), which prevents the \(R^2\) block becoming degenerate during training. Let \(g^{(\theta)}_{(T,X)}\) denote the upper-left \(2\times2\) block of the learnt metric in the Penrose directions. The determinant loss is defined by \[\mathcal{L}_{\det} := \frac{1}{N}\sum_{a=1}^{N} \frac{1}{\big|\det g^{(\theta)}_{(T,X)}(w_{a})\big|+\epsilon}, \label{eq:r295det95loss}\tag{44}\] and prevents the degenerate attractor in which the two-dimensional Lorentzian block components flow to 0 whilst reducing the curvature terms.

The Schwarzschild-search objective becomes \[\mathcal{L}_{\rm tot}^{\rm Schw}(t) := \frac{\alpha_{\rm E}(t) \mathcal{L}_{E,\text{Schw}} + \alpha_{W}(t)\mathcal{L}_{W}+ \alpha_{SO(3)}(t)\mathcal{L}_{SO(3)} + \alpha_{\det}(t)\mathcal{L}_{\det}}{\alpha_{\rm E}(t)+\alpha_{W}(t)+\alpha_{SO(3)}(t)+\alpha_{\det}(t)}, \label{eq:total95loss}\tag{45}\] where \(t\) denotes the training epoch, and the \(\alpha_i(t)\) are the scheduled multipliers of the active loss terms which change the weighting of various loss components throughout training. The training samples may optionally be regenerated at every epoch using the current scheduled density powers, while the validation set is held fixed. Batches are shuffled through a dataset pipeline, and the implementation facilitates Ricci, Kretschmann, and Killing-diagnostic quantity evaluation either by the standard two-tape kernel or by the optimised forward-mode kernel; validation is computed in batches. Post training, on an independent test data sample the training losses are computed to evaluate performance. In addition to the losses, the 3 invariant scalars introduced in 2.2.1 are computed and compared to the expected values for the Schwarzschild solution as further validation of the found solution the architecture is approximating. These include matching the speciality index \(\mathcal{S} = 1 + 0i\), the Kretschmann scalar \(\frac{Kr^6}{m^2} = 48\), and the cubic curvature ratio \(\rho_{\text{Schw}} = -\frac{1}{\sqrt{12}}\).

The Schwarzschild-targeted experiment trains Ricci flatness, the Weyl curvature scale, spherical symmetry, and non-degeneracy, and then asks whether these independent invariant diagnostics agree with the Schwarzschild/type-\(\mathrm{D}\) values. By Birkhoff’s theorem, if the resulting solution is spherically symmetric and vacuum, it is guaranteed to be locally Schwarzschild; these scalar checks provide additional confidence that the resulting metric is indeed a Schwarzschild solution.

3.4.4 Unsupervised Petrov Type-I Network↩︎

The Schwarzschild experiment above uses the quadratic Weyl loss to fix the curvature scale and then checks the Petrov speciality index and Schwarzschild cubic invariant a posteriori. A complementary search can instead use the speciality index in the loss function as a repeller from the Schwarzschild/type-\(\mathrm{D}\) value \(\mathcal{S}=1\). This is motivated by the role of \(\mathcal{S}\) as an invariant measure of deviation from algebraically special Kerr/Schwarzschild geometry [40] and by its later use in local Petrov-classification diagnostics [43]. These runs do not include the analytic Schwarzschild Kretschmann loss or the \(SO(3)\) symmetry loss, so the network is not supplied with the Schwarzschild radial curvature profile.

Beginning by stating a variant of the Einstein loss, 40 , which is homothety⁸ invariant, the type I Einstein loss is defined as \[\mathcal{L}_{E,I} := \frac{1}{N} \sum_{a=1}^{N} \frac{1}{|I_\theta(w_a)|+\epsilon} \left. \left| E_{\mu\nu} \left(g_E^{(\theta)}\right)^{\mu\alpha} \left(g_E^{(\theta)}\right)^{\nu\beta} E_{\alpha\beta} \right| \right|_{w_a}. \label{eq:type-i-einstein-loss}\tag{46}\]

For each sample with Weyl invariants \(I(w_a)\) and \(J(w_a)\), define the speciality index as in 12 . The implemented Petrov-profile term fits a prescribed nonconstant profile \(\mathcal{S}_{\rm targ}(w_a)\) centred away from the algebraically special value \(1\): \[\mathcal{L}_{\rm P}^{\rm profile} := \frac{1}{N}\sum_{a=1}^{N} \left| \frac{27J(w_a)^{2}}{I(w_a)^{3}+\epsilon_I} - \mathcal{S}_{\rm targ}(w_a) \right|^{2}, \qquad \mathcal{S}_{\rm targ}(w_a)\not\equiv1 . \label{eq:petrov95repeller95loss}\tag{47}\] In the runs reported below the profile is quadrupolar in nature⁹ and centred near the real value \(2\), so agreement with \(\mathcal{S}_{\rm targ}\) explicitly penalises collapse back to the Schwarzschild/type-\(\mathrm{D}\) value. The target profile is an arbitrary numerical training choice rather than a canonical closed-form Ricci-flat solution from the literature, and any complex value would work equivalently well. Since \(\mathcal{S}=1\) characterises the algebraically special type-\(\mathrm{D}\) value for Kerr/Schwarzschild, this term is best understood as a conservative device for searching away from that basin, not on its own a positive classifier for a particular alternative Petrov type. The retained \(I\) and \(J\) diagnostics are therefore still needed to check that the model has not collapsed to a conformally flat or otherwise degenerate case where the speciality index is ill-conditioned.

In practice, we also include a bounded curvature repeller \[\mathcal{L}_{K{\rm rep}} := \frac{1}{N}\sum_{a=1}^{N} \frac{\epsilon_K}{|K_\theta(w_a)|+\epsilon}, \label{eq:k95repeller95loss}\tag{48}\] to discourage collapse to a conformally flat or nearly flat state. This term does not supply the Schwarzschild Kretschmann target; it only acts as a soft barrier against the \(K_\theta\simeq0\) failure mode.

The corresponding objective is thus \[\mathcal{L}_{\rm tot}^{\rm I}(t) := \frac{ \alpha_{\rm E}(t) \mathcal{L}_{E,I} +\alpha_{\rm P}(t)\mathcal{L}_{\rm P}^{\rm profile} +\alpha_{\det}(t)\mathcal{L}_{\det} +\alpha_{K{\rm rep}}(t)\mathcal{L}_{K{\rm rep}}}{\alpha_{\rm E}(t)+\alpha_{\rm P}(t)+\alpha_{\det}(t)+\alpha_{K{\rm rep}}(t)}. \label{eq:non95petrov95d95loss}\tag{49}\]

3.4.5 Unsupervised Petrov Type-I Black Hole Network↩︎

The type-I search of 3.4.4 drives the geometry away from the algebraically special locus while keeping it numerically Ricci-flat, but it fixes neither an absolute curvature scale nor any causal structure. The homothety-invariant loss \(\mathcal{L}_{E,I}\) is insensitive to constant rescalings \(g\mapsto\Lambda^{2}g\), so the resulting type-I vacua carry no distinguished horizon and are not, in themselves, black holes. To promote them to black holes we add two terms: a horizon curvature anchor that fixes the scale and enforces a regular finite-curvature horizon, and a trapped-surface term that installs the causal structure of 2.4.

The anchor uses the same degree-two Weyl scalar as the Schwarzschild curvature-scale loss 42 , but Gaussian-localised (sample-wise weighted by \(\phi_a^{-1}\)) to the horizon band \(r\simeq2m\) rather than imposed everywhere, \[\mathcal{L}_{\rm hor} := \frac{\sum_{a}w_a\big(r(w_a)^{3}\sqrt{|I_\theta(w_a)|+\epsilon}-\sqrt3\,m\big)^{2}}{\sum_{a}\phi_a}, \qquad \phi_a=\exp\!\left[-\left(\frac{r(w_a)-2m}{\sigma\,2m}\right)^{2}\right]. \label{eq:horizon95anchor95loss}\tag{50}\] As the target \(\sqrt{|I_{\rm Schw}|}\,r^{3}=\sqrt3\,m\) of 41 is a finite constant, anchoring it at \(r=2m\) simultaneously fixes the overall scale – so that the Kretschmann invariant attains its physical value \(K_{\rm hor}\) at the horizon, rather than being sent to zero by the homothety freedom – and forbids a divergent, naked curvature there. Imposing the constraint at a single radius, instead of throughout the domain as in 42 , leaves the bulk profile free to remain type-I rather than collapsing back to Schwarzschild. This term subsumes and replaces the curvature repeller 48 .

The trapped-surface term acts on the scalar \(\Xi\) of 28 , evaluated for the learnt metric by automatic differentiation of the fibre areal radius \(R\). A black-hole interior requires \(\Xi<0\) inside \(r<2m\) and \(\Xi\to0\) at the horizon; we impose this one-sidedly, \[\mathcal{L}_{\rm trap} := \frac{1}{N}\sum_{a=1}^{N} \Big[\max\!\big(0,\;\beta(r_a)-\operatorname{sign}(r_a-2m)\,\Xi_\theta(w_a)\big)\Big]^{2}, \qquad \beta(r)=\Big(\beta_0+\beta_1\frac{\max(0,2m-r)}{2m}\Big)\mathbf{1}_{r<2m}. \label{eq:trapped95loss}\tag{51}\] The sign factor flips the target across the horizon, requiring \(\Xi\le-\beta\) on the interior and \(\Xi\ge0\) on the exterior. The margin \(\beta(r)\) is the Schwarzschild-linearised trapped depth \(\beta_1(2m-r)/2m\), which vanishes at the horizon and grows inward, plus a small floor \(\beta_0>0\) that maintains a non-zero gradient just inside \(r=2m\), where the linear part is weak. The one-sided \(\max(0,\cdot)\) penalises only insufficient trapping and leaves \(\Xi\) free to be more negative, so that the type-I radial profile is not pinned to the Schwarzschild one. Since \(\Xi\) is invariant under time reversal, \(\mathcal{L}_{\rm trap}\) enforces \(\Xi<0\) but not the black-hole orientation \(\theta_\ell,\theta_n<0\) specifically; the orientation is checked a posteriori.

The black-hole search objective is therefore \[\mathcal{L}_{\rm tot}^{\rm BH}(t) := \frac{\alpha_{\rm E}(t)\mathcal{L}_{E,I} +\alpha_{\rm P}(t)\mathcal{L}_{\rm P}^{\rm profile} +\alpha_{\rm hor}(t)\mathcal{L}_{\rm hor} +\alpha_{\rm trap}(t)\mathcal{L}_{\rm trap} +\alpha_{\det}(t)\mathcal{L}_{\det}}{\alpha_{\rm E}(t)+\alpha_{\rm P}(t)+\alpha_{\rm hor}(t)+\alpha_{\rm trap}(t)+\alpha_{\det}(t)}. \label{eq:blackhole95loss}\tag{52}\] Relative to the type-I objective 49 , this drops the curvature repeller and adds \(\mathcal{L}_{\rm hor}\) and \(\mathcal{L}_{\rm trap}\); the trapping weight \(\alpha_{\rm trap}\) is switched on only after a warm-up during which the vacuum, type-I and horizon terms are first established.

4 Results↩︎

This section reports the results from each numerical investigation in the search of Einstein metrics using the architectures described in 3.2. All runs hyperparameters are stated in 5.

4.1 Local Lorentzian Searches↩︎

As a local calibration before the four-dimensional searches, we train two-dimensional Lorentzian Einstein metrics on a single ball patch (contractible, topologically trivial) for Einstein constants \(\lambda\in\{+1,0,-1\}\). These runs use the same Lorentzian metric parametrisation and the Einstein residual in 38 , with no other terms in the loss function.

1 reports the global test losses for the three Einstein constants, together with the corresponding mean values of \(\det(g)\) over the sampled points. The losses are small in all three cases, with the Ricci-flat case giving the lowest Einstein loss. The determinant values remain negative and away from zero, indicating that the learnt metrics retain Lorentzian signature rather than collapsing towards a degenerate solution.

The component plots in 3 show representative learnt metric and Ricci components for each sign of \(\lambda\). In two dimensions the Einstein condition fixes the Ricci tensor component-wise once the metric and \(\lambda\) are specified, so the agreement between the small global losses, nonzero negative determinants, and the plotted components provides a compact sanity check for the Lorentzian local training setup used by the higher-dimensional experiments.

4.2 Schwarzschild Search↩︎

Having corroborated the most important term among the PINN’s loss components, the next natural step is to test our method for the prototypical Ricci-flat solution: the Schwarzschild black hole. As discussed in 2.1, the choice of Penrose coordinates for this experiment provides a finite computational domain and a unified description of the whole spacetime (except the centre point of the black hole singularity).

In order to have a reference benchmark for results obtained through unsupervised training, we first performed a series of runs in supervised mode. As described at the beginning of 3.4, the supervised objective simply consists of approximating the target metric, which in this case is the Schwarzschild black hole written in Penrose coordinates. Since infinities cannot be represented computationally¹⁰, points very close to the boundary of the Penrose diagram, which importantly include the both spatial infinity and singularity (where curvature blows up), are excluded from the experiments. The results presented in rest of the paper were obtained by sampling according to the strategy described in 3.3 with \(\rho = 0.85\). This value provides a substantial covering of the Penrose diagram, spanning a region from well beyond the event horizon up to multiple event horizon radii outside of the black hole.

For this benchmark \(10\) independent supervised training runs were performed, which all saturated at a supervised training losses of \(\sim 10^{-3}\). For each supervised trained model an independent test dataset was used to evaluate the performance, as presented in the respective columns of 2. The test evaluation metrics were the same physical losses used in the unsupervised runs (none of which are used in supervised training), and the relevant invariants for the Schwarzschild solution, as introduced in 2.2.1. These provide the benchmarks to assess the testing metrics of the unsupervised runs, whereas the supervised runs have exact knowledge of the analytic Schwarzschild solution to train against, and represents the limits of approximation with the allocated computational resources. The unsupervised runs to compare against will be performing a general search with the Birkhoff-Jebsen aligned training objectives of 3.4.3.

Proceeding to the unsupervised experiment. These runs were trained subject to the losses presented in 3.4.3, i.e. the optimiser is not given the analytic Schwarzschild metric components as pointwise targets, but instead learns a Lorentzian metric by minimising geometric residuals computed from its own output by automatic differentiation. The active constraints are Ricci flatness, the Weyl curvature-scale target, the \(SO(3)\) symmetry residual, and the determinant barrier preventing degeneration of the Penrose \((T,X)\) block. These losses encode the invariant conditions needed to select the Schwarzschild branch, rather than an exact fitting of the metric component-by-component.

Each of the 10 Schwarzschild-targeted runs were evaluated on an independent test sample of 2000 points. The average test losses (with standard deviation), and invariant values, are given in the respective columns of 2. Each of the test losses are orders of magnitude below the benchmark supervised training, validating that this training process has indeed found a Schwarzschild solution in each case, in alignment with the Birkhoff-Jebsen conditions. As independent validation, each of the invariants computed for the trained models on the independent test data are also orders of magnitude closer to the target values for the Schwarzschild solution, despite the training process providing no information about them. Additionally, the low value of the maximum determinant across the manifold validates the solutions are consistently Lorentzian. With both the strong loss results, and strong independent invariant corroboration, there is solid confidence in concluding the AInstein architecture can find and model the Schwarzschild solution.

For further visual corroboration, 4 shows a learned \(g_{00}\) component for a representative run together with the corresponding Ricci component \(R_{00}\), as well as the analytic values of each from the known Schwarzschild solution for comparison. This metric component has the expected large variation across the compactified domain, while the Ricci component remains small on the same samples, validating the Ricci-flat behaviour of the found Schwarzschild solution for this representative component. The remaining components of the trained metric and respective Ricci tensor for a representative Schwarzschild run are shown in 7 of 7, which further corroborates the Ricci-flat nature of the found solution; additional plots on \(S^2\) domain are included in the associated GitHub repository.

The independent scalar checks in 5 visualise the deviation of each invariant of 13 for a representative unsupervised trained Schwarzschild solution from the expected values of Schwarzschild. Corroborating the low deviation from expectation shown in 2, all residuals in this figure show low values, demonstrated by the colour scales. This again provides independent corroboration that this training process and AI architecture can indeed reproduce the Schwarzschild solution via a targeted search.

4.3 Petrov-Type I Vacuum Search↩︎

The second experiment uses the same embedded Lorentzian architecture and compactified sampling domain, but changes the target from Schwarzschild targeted search to a general type-I vacuum search. As described in 3.4.4, the Einstein residual is replaced by the homothety-invariant \(\mathcal{L}_{E,I}\), so the search is insensitive to overall constant metric rescalings, whilst the Schwarzschild Weyl-scale loss is removed and replaced with the speciality-index profile loss in 47 . The vacuum type-I objective of 49 therefore combines this Ricci-flatness condition with an arbitrary target speciality index profile which is representative of type I solutions, and includes two non-physical artefact losses to deter degeneracy (determinant barrier and the boundary curvature repeller). The \(SO(3)\) symmetry term is not used as a training constraint in this run, so the model is free to find non-symmetric solutions; the Killing residual is retained only as a diagnostic of departure from this spherically-symmetric basin.

3 reports the test diagnostics for the 10 seed runs of 2000 points, showing test losses and relevant invariant values. The low average Einstein loss (still lower than supervised Schwarzschild benchmark) shows that these runs do find numerically Ricci-flat solutions; and that these solutions do match the arbitrary speciality index profile. This strongly supports the idea that the AInstein architecture can be used for more general Ricci-flat vacuum solutions, with freedom to change the profile to match desired conditions.

In addition to these test losses, the Killing loss was also computed, and the reported value is much larger than in the Schwarzschild search, indicating the found solutions are not spherically symmetric (which is expected because spherical symmetry is no longer imposed as a loss in training). The determinant check remains strictly negative, confirming the consistent Lorentzian nature, and the speciality index is centred near the intended type-I profile value.

The full component-level metric and Ricci diagnostics for the Petrov type-I vacuum search are collected in 8 of 7. These plots provide the component-level counterpart to the invariant diagnostics in 3: the learnt metric components are non-Schwarzschild in scale and structure, while the Ricci components remain small relative to the metric components. Taken together, the loss diagnostics, determinant check, and speciality-index profile indicate that the model has found non-degenerate, numerically Ricci-flat geometries outside the Schwarzschild/type-\(\mathrm{D}\) basin.

4.4 Petrov-Type I Black Hole Search↩︎

The final experiment applies the black-hole objective 52 on the same embedded domain, adding the horizon anchor 50 and the trapped-surface term 51 to the type-I recipe of 4.3. Two variants are reported here, differing only in the trapped-surface margin and weight, \((\beta_0,\alpha_{\rm trap})=(0.05,15)\) and \((0.08,25)\) (both with \(\beta_1=1\)), each over \(10\) random seeds, with \(\alpha_{\rm trap}\) warmed up over the first \(30\%\) of training. 4 collects the diagnostics, evaluated on the independent test samples of 2000 points; the interior/exterior of each fibre and its expansions \(\theta_\ell,\theta_n\) are computed by finite differences of the learnt areal radius.

Across all \(20\) runs the black-hole criteria hold simultaneously and robustly. The geometry is numerically Ricci-flat, shown by the average Einstein loss matching the supervised Schwarzschild scale. The full component-level metric and Ricci diagnostics for the black-hole search are collected in 9 of 7. The components of the trained metric and respective Ricci tensor for a representative Schwarzschild run in the \((\beta_0,\alpha_{\rm trap})=(0.08,25)\) case are shown in 9 of 7, which visually validate the Ricci-flat nature of the found solutions.

4 also shows the solutions are algebraically general, \(\mathcal{S}\simeq 1.99\), well away from the type-D value. The horizon curvatures are regular and at the Schwarzschild scale, \(\sqrt{|I|}\,r^{3}\to1.75\approx\sqrt3\) and \(K_{\rm hor}\to0.77\approx0.75\), so the solutions are not naked singularities; the negative max(det(\(g\))) confirms the solutions are consistently Lorentzian; and the metric is manifestly non-Schwarzschild, with relative Frobenius distance \(\simeq0.6\) to the analytic solution. Crucially in each run, the interior is genuinely trapped: \(\Xi<0\) on \(97\)–\(99\%\) of the interior samples, with interior median \(\Xi\simeq-0.1\).

Because \(\Xi\) is invariant under time reversal, the loss fixes \(\Xi<0\) but not the black-hole orientation. We resolve it a posteriori through the two null expansions \(\theta_\ell,\theta_n\), with the convention calibrated so that the future interior \(B_{\rm II}\) of the analytic Schwarzschild solution comes out future-trapped.¹¹

6 shows where the trapped region sits for this model, and reproduces the expected causal structure of 1: the trapped region coincides with the interior \(r<2m\), bounded by the null horizons \(T=\pm X\), with the exteriors untrapped.

Finally, as these metrics are not spherically symmetric, we confirm that \(\Xi\) faithfully reports the trapping. Recomputing the expansions with the exact fibre normals – which, owing to the small residual \(\mathbb{R}^{2}\)–\(S^{2}\) cross terms of the learnt metric, do not lie exactly in the \((T,X)\) plane – shifts \(\theta_\ell,\theta_n\) by only about \(5\%\), far below their interior magnitude \(\simeq0.1\), and leaves the whole interior future-trapped. The classification is therefore robust against the departure from spherical symmetry.

These final results therefore validate that the AInstein architecture can be specialised to trapped surface solutions too, setting up the codebase and this methodology for general black hole search.

5 Conclusions and Outlook↩︎

This work upgrades AInstein from Riemannian sphere searches to Lorentzian black-hole geometries. The essential ingredients are an indefinite metric parametrisation, compactified Penrose-domain sampling, a single ambient network whose pullback handles the \(S^{2}\) topology globally, and differentiable losses built directly from curvature, invariants, non-degeneracy, and trapping diagnostics. This allows for general vacuum solution search while retaining tensorial control over the geometry being learnt.

The results show this machinery working in increasing generality. Local Lorentzian tests recover two-dimensional Einstein metrics with all three signs of \(\lambda\). The Schwarzschild experiment then recovers the maximally extended vacuum black hole from geometric constraints rather than supervised metric fitting: Ricci flatness, the Weyl scale, and spherical symmetry select the Schwarzschild branch, while independent invariant checks confirm the expected type-\(\mathrm{D}\) curvature structure. The Petrov-profile search generalises from Schwarzschild-specific desired curvatures and uses the same framework to find non-degenerate, numerically Ricci-flat metrics in a different Petrov class, type I. Finally, adding a horizon curvature anchor and a trapped-surface loss leads the architecture to find and approximate algebraically general vacuum black-hole candidates with regular Schwarzschild-scale horizons and trapped interiors.

This work develops the AInstein methodology into Lorentzian solution reproduction, and then further into a flexible variational search engine for vacuum geometry. One may specify invariant profiles and causal conditions, and the network searches the corresponding metric space directly. This is especially promising for black-hole problems where analytic ansätze are restrictive or unavailable, yet this AI methodology can provide numerical approximations. Natural next targets include Kaluza–Klein black strings and localized black holes, non-uniform or distorted horizons, AdS black holes with prescribed boundary data, and stationary vacuum searches where symmetry assumptions should be relaxed rather than built in from the start.

Data Availability↩︎

The respective codebase was written in Python, and scripts used to generate all the results are available at the AInstein repository:
https://github.com/xand-stapleton/ainstein.

Acknowledgements↩︎

The authors thank Michael Douglas, Fabian Ruehle, and Tomás Silva for helpful comments during the “Mathematics and Machine Learning Program” at Harvard University; Gianfranco Cortes and Yueqing Feng for discussions on the embedding architecture; Kymani Armstrong-Williams for comments on the Kretschmann invariant of the Schwarzschild solution; and David Berman for helpful comments on well-defined tensorial integration measures, as well as bringing us all together.

EH is supported by São Paulo Research Foundation (FAPESP) grant 2024/18994-7. TSG is supported by the 2024 Max Planck-Humboldt Research Award bestowed on Geordie Williamson and Catharina Stroppel by the Max Planck Society and the Alexander von Humboldt Foundation. AGS and TSG acknowledge support from Pierre Andurand over the course of this research.

Computations were performed on the HPC system Viper at the Max Planck Computing and Data Facility. This research utilised computational resources of the Centro Nacional de Processamento de Alto Desempenho em São Paulo (CENAPAD-SP). This research used Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT.

6 Invariant Derivations for the Schwarzschild Solution↩︎

Spherical symmetry gives a preferred geometrical radius, the areal radius \(r\), defined by the area \(A\) of each two-sphere, \[A=4\pi r^{2}. \label{eq:areal95radius95def}\tag{53}\] In any spherically symmetric spacetime there is a corresponding Misner–Sharp mass \(M_{\rm MS}\), defined invariantly by \[1-\frac{2M_{\rm MS}}{r} = \nabla^{a}r\nabla_{a}r. \label{eq:misner95sharp95def}\tag{54}\] This definition depends only on the areal radius and the spacetime metric, not on a particular Schwarzschild coordinate chart. The Einstein equations imply the local mass-balance relation \[\nabla_a M_{\rm MS}=4\pi r^2\left(T_a{}^b\nabla_b r-T^{A}{}_{A}\nabla_a r\right) \, , \label{eq:misner95sharp95balance}\tag{55}\] with \(T^A{}_A=T^t{}_t+T^r{}_r\) being the trace over the 2D space normal to the symmetry spheres. Therefore, in vacuum, \(T_{ab}=0\), the Misner–Sharp mass is constant: \[M_{\rm MS}=m. \label{eq:misner95sharp95constant}\tag{56}\]

In an orthonormal frame adapted to the static Schwarzschild geometry, with \(e_{0}\) timelike, \(e_{1}\) radial, and \(e_{2},e_{3}\) tangent to the symmetry two-sphere, the independent non-zero Riemann components are \[\begin{align} R_{0101}&=-\frac{2m}{r^{3}}, & R_{0202}=R_{0303}&=\frac{m}{r^{3}}, \nonumber\\ R_{1212}=R_{1313}&=-\frac{m}{r^{3}}, & R_{2323}&=\frac{2m}{r^{3}}, \label{eq:schwarzschild95orthonormal95riemann} \end{align}\tag{57}\] up to the overall sign convention for the Riemann tensor. The Kretschmann scalar is insensitive to this overall sign, and the final Schwarzschild value agrees with the standard curvature-invariant expression [42]. Using the usual Riemann symmetries, the contraction counts each independent sectional component four times, giving \[\begin{align} K_{\rm Schw} &= R_{abcd}R^{abcd} \nonumber\\ &= 4\left[ \left(\frac{2m}{r^{3}}\right)^{2} +\left(\frac{m}{r^{3}}\right)^{2} +\left(\frac{m}{r^{3}}\right)^{2} +\left(\frac{m}{r^{3}}\right)^{2} +\left(\frac{m}{r^{3}}\right)^{2} +\left(\frac{2m}{r^{3}}\right)^{2} \right] \nonumber\\ &= \frac{48m^{2}}{r^{6}}. \label{eq:schw95kretschmann95geometric} \end{align}\tag{58}\] Equivalently, this gives the invariant check \(Kr^{6}=48m^{2}\). Substituting the Penrose-domain radius \(r=2m(1+W_{\mathcal{B}})\) from 8 then gives the analytic target used in the implementation, \[K_{\rm Schw} = \frac{48m^{2}}{r^{6}} = \frac{3}{4m^{4}(1+W_{\mathcal{B}})^{6}}. \label{eq:analytic95kretschmann}\tag{59}\] The same orthonormal components give the cubic contraction. To see the numerical factor explicitly, raise the second antisymmetric index pair. In the ordered bivector basis \[(01),(02),(03),(12),(13),(23),\] the non-zero diagonal entries of \(R_{ab}{}^{cd}\) are \[R_{01}{}^{01}=\frac{2m}{r^3},\quad R_{02}{}^{02}=R_{03}{}^{03} =R_{12}{}^{12}=R_{13}{}^{13}=-\frac{m}{r^3},\quad R_{23}{}^{23}=\frac{2m}{r^3} . \label{eq:riemann95bivector95eigenvalues}\tag{60}\] The signs differ from the lowered components whenever a time index is raised. Since the full Einstein summation is over ordered antisymmetric pairs, each independent bivector eigenvalue appears eight times: for example the \(01\) bivector contributes through the choices \(ab=01\) or \(10\), \(cd=01\) or \(10\), and \(ef=01\) or \(10\). The three antisymmetry signs multiply to give the same contribution in each case. Equivalently, the contraction is eight times the trace of the cube of the curvature operator on \(\Lambda^{2}\): \[\begin{align} R_{ab}{}^{cd}R_{cd}{}^{ef}R_{ef}{}^{ab} &= 8\left[ \bigg(\frac{2m}{r^3}\bigg)^3+\bigg(-\frac{m}{r^3}\bigg)^3+\bigg(-\frac{m}{r^3}\bigg)^3+\bigg(-\frac{m}{r^3}\bigg)^3+\bigg(-\frac{m}{r^3}\bigg)^3+\bigg(\frac{2m}{r^3}\bigg)^3 \right] \nonumber\\ &= 96\bigg(\frac{m}{r^3}\bigg)^{3}. \label{eq:schwarzschild95cubic95counting} \end{align}\tag{61}\] Thus \[J_{\rm cub} = R_{ab}{}^{cd}R_{cd}{}^{ef}R_{ef}{}^{ab} = 96\frac{m^{3}}{r^{9}}. \label{eq:schwarzschild95cubic95invariant}\tag{62}\] Combining 58 and 62 gives \[\frac{J_{\rm cub}}{K_{\rm Schw}^{3/2}} = \frac{96m^{3}r^{-9}}{(48m^{2}r^{-6})^{3/2}} = \frac{1}{\sqrt{12}}. \label{eq:schwarzschild95cubic95ratio95derivation}\tag{63}\] The sign of this cubic contraction changes with the opposite Riemann-tensor convention, whereas the Kretschmann scalar does not.

The speciality-index calculation uses the complex Weyl invariants \(I\) and \(J\), not the real Riemann contraction \(J_{\rm cub}\). The relevant invariant construction is the one used by Baker and Campanelli for the speciality index, and by later Petrov-classification diagnostics [40], [43]. In Newman–Penrose notation [53] the five Weyl scalars \(\Psi_{0},\ldots,\Psi_{4}\) are the components of the Weyl tensor in a complex null tetrad, and the scalar polynomial invariants may be written as \[\begin{align} I &= \Psi_{0}\Psi_{4} -4\Psi_{1}\Psi_{3} +3\Psi_{2}^{2}, \nonumber\\ J &= \begin{vmatrix} \Psi_{4} & \Psi_{3} & \Psi_{2}\\ \Psi_{3} & \Psi_{2} & \Psi_{1}\\ \Psi_{2} & \Psi_{1} & \Psi_{0} \end{vmatrix}. \label{eq:np95weyl95invariants95appendix} \end{align}\tag{64}\] These are the Newman–Penrose form of the self-dual Weyl contractions in 11 , with the same normalisation. The speciality index is then \(\mathcal{S}=27J^{2}/I^{3}\), which equals unity for the algebraically special Kerr/type-\(\mathrm{D}\) limit in this convention [40]. In vacuum \(R_{abcd}=C_{abcd}\), and in the Schwarzschild principal tetrad the only non-zero Weyl scalar is \[\Psi_2=-\frac{m}{r^3}, \label{eq:schwarzschild95psi2}\tag{65}\] which is real [42]. All other Weyl scalars vanish in this tetrad. Substituting \[\Psi_0=\Psi_1=\Psi_3=\Psi_4=0\] into 64 gives \[I=3\Psi_2^2=\frac{3m^2}{r^6}, \qquad J=-\Psi_2^3=\frac{m^3}{r^9}, \qquad \mathcal{S} = \frac{27J^2}{I^3} = 1+0i . \label{eq:schwarzschild95speciality95appendix}\tag{66}\] Thus the constants in the Schwarzschild speciality-index target come directly from the standard Newman–Penrose invariant polynomials and the single Coulomb Weyl scalar of the type-\(\mathrm{D}\) Schwarzschild geometry. Comparing 62 and 66 gives \(J_{\rm cub}=96\,\mathrm{Re}(J)\) for the convention used here. Therefore the dimensionless scalar reported in the diagnostics is \[\rho := -\frac{96\,\mathrm{Re}(J)}{|K|^{3/2}} = -\frac{J_{\rm cub}}{|K|^{3/2}}, \qquad \rho_{\rm Schw} = -\frac{1}{\sqrt{12}}, \qquad |\rho_{\rm Schw}| = \frac{1}{\sqrt{12}}. \label{eq:rho95schwarzschild95appendix}\tag{67}\] This last normalisation defines the diagnostic used in the numerical plots; the cited invariant results are the Schwarzschild Weyl scalar, Weyl invariants, and Kretschmann scalar. Hence the invariant targets are not additional coordinate ansatz data: they follow from the vacuum equation, spherical symmetry, and the areal-radius function \(r(T,X)\).

7 Component-Level Metric and Ricci Visualisations↩︎

This appendix collects the full component-level validation plots for the three main four-dimensional runs. In each figure, the left block shows the learnt metric components \(g_{\mu\nu}\) and the right block shows the corresponding Ricci tensor components \(R_{\mu\nu}\), both arranged by tensor index over the validation samples. The shared colour bar beneath each pair fixes the component scale used across the metric and Ricci panels, while the captions report the mean absolute component sizes used in the diagnostics.

8 Model Hyperparameters↩︎

Summaries of the optimal training hyperparameters for the Schwarzschild experiment are presented in 5, whilst 6 highlights the key changes for the Type-I vacuum search, and 7 the key changes for the black hole search. An exhaustive list of parameters is presented in the GitHub repository.

References↩︎

---

1. We work with the ‘mostly-positive’, or ‘East-Coast’ signature \((-,+,+,+)\), and use the conventional definition \(R^a{}_{bcd}=\partial_c\Gamma^a_{db}-\partial_d\Gamma^a_{cb}+\Gamma^i_{bd}\Gamma^a_{ic}-\Gamma^i_{bc}\Gamma^a_{id}\) for the Riemann tensor.↩︎

2. This is the Physics Informed Neural Network (PINN) formalism [17].↩︎

3. Here we note complementary applications of pure supervised learning to the black hole context [32]–[34].↩︎

4. AInstein is available at: https://github.com/xand-stapleton/ainstein.↩︎

5. By this we mean \(g_E\) is obtained by the same vielbein as \(g\), but using \(\delta\) as the flat metric, instead of \(\eta\).↩︎

6. Note the code also has functionality for general \(\lambda\) search.↩︎

7. This and other losses have \(\epsilon\) terms in denominators as numerical guards. Unless otherwise specified, \(\epsilon=10^{-12}\).↩︎

8. A homothety transform is simply a transformation that scales distances from a fixed point, known as the centre, by a given ratio.↩︎

9. The bounded quadrupolar target profile used was inspired by [52], with motivation from the Legendre-polynomial multipolar structure of Weyl/Zipoy-Voorhees-type static axisymmetric solutions. This choice can easily be generalised in the codebase for other searches.↩︎

10. For a fixed architecture there is a fixed limit of expressibility within the function class the architecture approximates. As infinities are approached this leads to numerical instabilities and overflows which need to be mitigated.↩︎

11. Interestingly, but a-posteriori unsurprisingly, the seeds divide almost evenly: half realise a genuine black hole, with \(\theta_\ell,\theta_n<0\) throughout \(B_{\rm II}\), and half the time-reversed anti-trapped (white-hole) branch.↩︎