Einstein Fields: A Neural Perspective To Computational General Relativity
Authors: Sandeep Suresh Cranganore
Andrei Bodnar
Arturs Berzins
Johannes Brandstetter
Venue: ICLR 2026 Poster (2026)
Keywords: neural fields (implicit neural representations)
neural compression
tensor fields
differential geometry
general relativity (GR) and numerical relativity (NR)
Sobolev training
differential geometry
finite-difference methods
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: We introduce *Einstein Fields*, a neural representation designed to compress computationally intensive *four-dimensional* numerical relativity simulations into compact implicit neural network weights. By modeling the *metric*, the core tensor field of general relativity, Einstein Fields enable the derivation of physical quantities via automatic differentiation. Unlike conventional neural fields (e.g., signed distance, occupancy, or radiance fields), Einstein Fields fall into the class of *Neural Tensor Fields* with the key difference that, when encoding the spacetime geometry into neural field representations, dynamics emerge naturally as a byproduct. Our novel implicit approach demonstrates remarkable potential, including continuum modeling of four-dimensional spacetime, mesh-agnosticity, storage efficiency, derivative accuracy, and ease of use. It achieves up to a $\mathtt{4,000}$-fold reduction in storage memory compared to discrete representations while retaining a numerical accuracy of five to seven decimal places. Moreover, in single precision, differentiation of the Einstein Fields-parameterized metric tensor is up to five orders of magnitude more accurate compared to naive finite differencing methods. We demonstrate these properties on several canonical test beds of general relativity and numerical relativity simulation data, while also releasing an open-source $\mathtt{JAX}$-based library: https://github.com/AndreiB137/EinFields, taking the first steps to studying the potential of machine learning in numerical relativity.

Published as a conference paper at ICLR 2026
EINSTEIN FIELDS: A NEURAL PERSPECTIVE TO COMPUTATIONAL GENERAL RELATIVITY
Sandeep S. Cranganore"
Andrei Bodnar* Johannes Brandstetter 1,3
Arturs Berzins"
Equal contribution
2 LIT AI Lab, Institute for Machine Learning, JKU Linz, Austria
University of Manchester, United Kingdom
3 'Emmi AI GmbH, Linz, Austria
(cranganore) berzins, brandstetter @ml jku. .tu.at andrei wwwnar@studden. aeaaarrrtttddeetmmmceette...... ac..k
manchester.ac,
ABSTRACT
We introduce Einstein Fields, a neural representation designed to compress compu- tationally intensive four-dimensionul numerical relativity simulations into compact implicit neural network weights. By modeling the metric, the core tensor field of general relativity, Einstein Fields enable the derivation of physical quantities via automatic differentiation. Unlike conventional neural fields (e.g.. signed distance, occupancy, or radiance fields), Einstein Fields fall into the class of Neural Tensor Fields with the key difference that, when encoding the spacetime geometry into neural field representations, dynamics emerge naturally as a byproduct. Our novel implicit approach demonstrates remarkable potential, including continuum model- ing of four-dimensional spacetime, mesh-agnosticity, storage efficiency, derivative accuracy, and ease of use. It achieves up a 4, 000-fold reduction in storage memory compared to discrete representations while retaining a numerical accuracy of five seven decimal places. Moreover, in single precision, differentiation of the Einstein Fields-par parameterized metric tensor is up to five orders of magnitude more accurate compared to naive finite differencing methods. We demonstrate these properties on several canonical test beds of general relativity and numerical relativity simulation data, while also releasing an open-source JAX-based library: https://github.. www.uu..oo/Annii com//nndeiB 137/EinFields, taking the first steps to studying the potential machine learning in numerical relativity.
INTRODUCTION
General relativity (GR) describes gravity as the curvature four-dimensional spacetime, encoded in the metric tensor and governed by the Einstein field equations (EFFS), system of coupled, nonlinear hyperboliceelliptic PDEs. Exact solutions are available only for idcalized cases, numerical relativity (NR) has become essential for accurate modeling of astrophysical events. Notable successes of NR include the high-precision modeling of black hole mergers (Abbott et al., 201 6a;b;c), high- precision binary neutron star merger simulations (Hayashi et al., 2025), and neutron star-black hole systems (Abbott et al., 2017). NR has also been central to the confirmation of gravitational waves (GWs) detected by LIGO and Virgo interferometers, leading to pobel-prize-winning discoveries.
Nonetheless, state-of-the-ard NR is one of the most computationally intensive domains of scientific computing, requiring massive parallelization on petascale computing infrastructures (Lovelace, 2021; Huerta et al., 2019). This is due several computational challenges in NR. including adaptive high-resolution spatial discretization and finite-difference (FD) methods on these, which, in addition, vulnerable numerical errors in sensitive regions. Moreover, NR simulations are equally storage- intensive, producing up to petabytes of data per run, prohibiting the storage and distribution of simulations HPC systems (Reed & Dongarra, 2015).
Recent progress in machine learning for scientific computing has shown the potential of hybrid neural-classic workflows (Thuerey et al., 2021; Zhang et al., 2023; Brunton al., 2020; Li et al.,

---

Published as a conference paper at ICLR 2026
2021; Brandstetter et al., 2022; Bodnar et al., 2025; Brandstetter, 2025). In addition, neural fields (NeFs) (Park et al., 2019; Müller et al., 2022; Chen & Zhang, 2019) have emerged as a powerful tool in visual computing for compact and continuous representations of traditionally discrete data, such as images, shapes, and physical fields, with ease of querying and differentiating. This raises the question of whether such hybrid approaches advance next-generation computational GR workflows, particularly in handling and compressing tensorial quantities and their derivatives defined on adaptive high-fidelity discretizations.
To this end, we propose and investigate EinFielde which provide the following contributions:
Neural compression. EinFi elds encode geometric information in compact neural representa- tions with typically fewer than two million parameters. They reproduce metric tensor components with relative precision up seven decimal digits (and up to nine in favourable coordinate charts). This yields memory-efficient approximations of complex spacetime geometries, with compression factors up to x across analytical and numerical solutions.
Discretization-free representation. EinField are trained on arbitrary point samples, includ- ing regular, irregular, and unstructured sets. They provide continuous query access to tensor fields at any resolution by learning these as continuous functions from discrete samples, which removes discretization artefacts.
Enhanced tensor differentiation. As smooth neural functions, EinF SenFielln elds support continuous evaluation of higher-order geometric quantities such as Christoffel symbols, Riemann tensors, and curvature invariants via point-wise automatic differentiation (Griewank & Walther, 2008). Initial results suggest that this approach can outperform high-order difee-difference methods on uniform grids, with accuracy gains up to 10 in FLOAT32.
High-fidelity reconstruction of tests General Relativity, We evaluate the physical fidelity of FinFie wiFielda ds analytical GR solutions and assess derived observables in addition to standard ML metrics. The models faithfully reproduce key relativistic phenomena, including orbital precession in Schwarzschild and Kсπr spacetimes, and allow accurate extraction of gravitational- wave distortions and strain. We further test a BSSN (Baumgarte & Shapiro, 202 evolution of oscillating neutron star. Early results indicate that EinFie aiields scale to realistic numerical relativity workflows despite the complexity of the solution.
BACKGROUND
Our work lies at the intersection of two domains: GR, along with its computational framework of NR, and NeFs, ubiquitous tool from computer vision. While a complete introduction to GR and its mathematical backbone, differential geometry, is beyond the scope of this Section (see detailed exposition in Appendix or succinct version Appendix D), we stress three key properties that pertain to work: (i) GR is field theory of tensor tensor-valued valued quantities, (ii) GR is intrinsically moordinaeeiindepedden and (iii) gravitational physics is entirely encapsulated in the metric and its first two derivatives.
Tensors. rank (r, s) tensor at a point I E M is the multilinear map from r covectors and s vectors to real number:
x V x....V-- -
cocops
s-copics
The r vectors and covectors pair with the respective r' covariant and s contravariant components of the tensor. As such, tensor is an element that lives in tensor product of vector and dual spaces, i.c., T E (V)00r & (V*)000. A tensor in a particular choice of basis {ca, hisnrr E V and {U hinns E V0 is given by T = Tra ary &t- B..2... Ea. .... X ear @ g X ... & gd. where gaa, CS.. , ,e,)) are the components of the tensor in this particular basis and transforms as shown in Eq. (20). A tensor field assigns to each point E M a multilinear map, i.e. a tensor, T, E 8 69 (V²*) )0 appropriate coordinates, its components T.1.... vary smoothly across the manifold.
Toias..a A..2...A. = T(001
B,..,,
General relativity extends Newtonian gravity with a geometric interpretation of gravity, where mass and energy tell spacetime how to curve, and curved spacetime tells objects how to move (Misner et al., 2017). This is formalized by the Einstein's 's field equations (EFEs)
G. Agas = G Tag.
(2)
2

---

Published as a conference paper at ICLR 2026
EFEs are set of 10 coupled non-linean tensor-valued, second partial PDEs and can be viewed as a tensorial generalization of the Newton-Poisson equation for gravity (² = pArGp(r) (Misner al., 2017: Poisson, 2004). In EFEs, Gag Ra,3 1/2 Rgg is the Einstein tensor, formed from the metric tensor field g ). which are solutions of the EFEs and tensorial generalization of the gravitational potential ((x). The Ricci curvafure tensor Rag. and the Ricci curvature scalar R are related by second derivatives 9a3. Thus, the left-hand side of EFES is entirely described by the metric and its derivatives, with A being the cosmological constant. The right-hand side depends on the stress- energy tensor Ta3 describing the matter distribution, with G being Newton's 's constant.
Metric tensor and its derivatives. The metric tensor is a rank (0, 2) symmetric bilinear form T. x T + R that generalizes the notion of an inner product on the tangent space Tg-M of a differentiable manifold M (Jost, 2008). It enables the computation of angles between vector fields and means to compute distances via the line element: ds² 9as (2*) drrd - In GR, the components Pag a particular coordinate system can be seen as 4 x symmetric matrix with ten independent components. The metric defines the causal structure and contains all geometric information of spacetime. Importantly, its partial derivatives Ô yield the Christoffel symbols I Tasy (z*), which describe the notion of parallel transport and defines covariant derivative operation Va = Ōa + ra (all detailed in Appendix A.3.3). In turn, the connection's 's derivatives (i.e., metric yield the Riemann curvature lensor Rag (I**), which encodes tidal forces of gravity. The trace part of Rep,y (index contraction w.r.t. metric: Tr, see Eq. (44)) is the Ricci tensor Ras, also a symmetric rank (0, 2) tensor. Its subsequent contraction yields the Ricci scalar R (all detailed in Appendix A.3.5). This can be summarized schematically follows (a more detailed pictorial version is available in Figure 10:
Ja 2 ag
Y R
Tre Ras
Thy R.
afty
Higer-order methods. FD methods with adaptive mesh refinement (AMR) (Berger & Oliger, 1984) have long underpinned tensor calculus in NR, discretizing space and time with high-order stencils (Appendix C). An ri-hh order stencil yields truncation errors O(h)), where the grid spacing. Widely used fourth- sixth-order schemes improve accuracy but incur larger communication costs due to broader stencil footprints in parallelized settings. In contrast, modern NR increasingly opts for (pseudo-)spectral Jspectral methods (Scheel et al., 2025), which represent fields globally through polynomial bases (Fornberg, 1996), yielding an efficiency of up to 1000- x faster on CPUs than FI) approaches on GPUs comparable accuracy (Rashti al., 2025).
Neural fields (NeFs), also known as implicit neural representations (INR) or coordinate-based neural networks, are multi-layer perceptrons (MLP) using very specific activation functions that are memory- efficient, implicit, continuous, infinitely differentiable maps, capable of capturing high-fidelity detail across complex domains (Xie et al., 2021; Essakine et al., 2024). Some well known NeFs are SIREN with sinusoidal activations (Sitzmann et al., 2020), WIRE with Gabor wavelet activations (Saragadam et al., 2023) for example. These properties have motivated their primary development and adoption in computer vision domains for representation, generation, and inversion tasks. Considering scientific computing domains, NeFs, when integrated with physics-informed losses (c.g., constraints and conservation laws), can be used for solving forward and inverse problems, including spatiotemporal dynamics governed by PDEs. In these settings, NeFs effectively act as PDE solvers, more commonly referred to physics-infommn neural networks (PINNs) (Raissi al., 2019; Karniadakis et al., 2021; Wang et al., 2025b).
METHOD PARAMETRIZING TENSOR FIELDS VIA EINFIELDS
Consider the doun-dimesiional spacetime (Al..) (a manifold equipped with metric) corresponding to an exact or numerical solution to the EFES': Gap 87 Tap. An BinField models the 10 independent components of the metric tensor field compact NeF, ultimately mapping the spacetime coordinates r = (x², _o., x², 2,23) to the symmetric rank (0, 2) metric tensor field:
9: x E M + ga (x) E Sym (TT
(4)
From now on, omit the cosmological constant term Agas-
3

---

Published as a conference paper at ICLR 2026
Soliiien to ErÈss Maiin tenser neld
Gas Alpp Se0 Tea
9ad (2
Einstan fies equations (brea)
Beef=") encroes 3pace0mm gromerry
Neural Twevor Faeide
spacotime and tur 2
Canhelds
Impl ot metri tensor Sua((") comporencs)
sugmented Biniledss
Tintwack
Jerztren Heeoderr
Improved metric tensor
connection
traectories
Curvature tonows
Figure conceptual overview of EinFielde training and downstream pipeline. (i) Premise: The Einstein field equations (EFEs) in Eq. (2) are highly non-linear PDEs defined 4D spacetime manifold, describing the geometric nature of gravitation. Their solutions define the metric tensor field (x=",, which encodes the full spacetime geometry and serves tensorial generalization of the gravitational potential. In this work, we parametrize (x*( using a neural network. (ii) Training: The training conducted such as uniform or hierarchical grids. Ei nF elds instead fit a continuous signal on these discrete representations, thus modeling 4D spacetime a 4D spacetime query coordinate
the metric tensor fields defined on 4D spacetime points,
continuum, and returning the metric tensor field for (t, x) All at arbitrary resolution. (iii) Sobolev supervision: The reconstruction quality of the metric and its derivatives is improved by augmenting Sobolev losses, i.e., metric Jacobian (neighborhood structure) and Hessian (curvatures). (iv) Validation and downstream tasks: Sobolev improved EinFields AD-based derivatives enable accurate point-wise retrieval of differential geometric quantities, such as the Levi-Civita connection (covariant derivative), geodesics, curvature tensors, and their invariants.
We deploy an MLP ğu with parameters 8, denoted ğ for simplicity. over-fit on the ground truth tensor field. Methodologically, this enables directly compressing the entire geometric information into storage-cheap NN weights, yielding continuous access (different from the training points) at arbitrary resolution of the metric and its non-trivial tensor differentiation (e.g., for Lie derivative covariant derivatives) information devoid of mesh (re)construction on curvilinear manifolds. This generalizes arbitrary rank (r, tensor field T..... (a*)). Thus, EinFiell ds posit neural alternative to address one or more of the challenges associated with traditional methods (typically utilizing higher-order finite differencing schemes) NR by not relying on costly spatiotemporal discretizations. Our framework should be considered special case of Neural Tensor Fields.
4

---

Published as a conference paper at ICLR 2026
3.1 DISTORTION
We define the distortion as the algebraic deviation of the spacetime metric from flat space,
A.s 9a.3 1a,3
(5)
with the Minkowski background fo in that particular coordinate system. From learning perspective, this decomposition acts as preprocessing step that removes the offset (fixed background metric). Flat contributions, that may even dominate numerically (e.g., 9t 1/r, Dee ~ r)) are removed, leaving only the non-trivin con-trivial geometric content, such as the curvature. Thus, the network focuses its representational capacity on meaningful deviations rather than redundantly relearning flat-space structure. As we show in 4.2, this improves scaling, accelerates convergence, and emphasizes the dynamic, physically relevant features of the metric during training.
3.2 RETRIEVING PHYSICS VIA NEURAL TENSOR FIELD DERIVATIVES
Higher derivative losses. Beyond the metric tensor itself, its first- and second-order derivatives are critical GR, as they govern geodesic motion, tidal forces, and curvature. Accurate trajectory reconstruction requires point-wise precise evaluation of Christoffel symbols and Riemann tensors. Such high-fidelit extraction from EinFielll is facilitated by Sobolev training (Czarnecki et al., 2017; Son et al., 2021), formulation that explicitly incorporates higher derivative losses (Chetan et al., 2024; Wang et al., 2025b) see Section F.4. The supervision on the metric Jacobian 0,.9. (40 independent components) and Hessian 8 (100 independent components) rectifies irregulari- ties in the metric field and its derivatives, yielding substantial accuracy gains and consequentially improving the precision of point- point-wise Christoffel symbols I l2 (z*)) and Riemann tensors R R% (x") queries by up to two orders of magnitude. This is given by the modified loss function:
200 (0) Ex (olloaa(=) (as(+)]12 + Σ N aa gap(a) 100--000 065) 065) (
(6)
with X being some coefficients and d. O₂ and a.2) o written in succinet notation. Instead of implementing higher-order FD stencils, our framework enables access to exact higher-order tensor derivatives via AD. This is illustrated in the AD workflow for differential geometry in Figure 2.
Reconstructing dynamics. Free-fall trajectories around massive objects follows a geodesic motion Eq. (56), which depends on Christoffel symbols I'(g. 0g). In our workflow (Figure 2), EinF ields reconstruct (g) ağ) with Jacobian supervision, enabling V. 0, +I. and, thus, accurate modeling of geodesic path and direct measurement of curvature via curvature tensors become possible.
Characterizing intrinsic geometry. Beyond dynamics, EinFields must reproduce the intrinsic geometry encoded curvature tensors and invariants. This constitutes Riemann Rapys tensor and the associated geodesic deviation Eq. (60), Weyl Casyá, Ricci Ras tensors, scalar curvature R, and invariants such as the Kretschmann scalar X (detailed in A.3.5.2). With Jacobian and Hessian-level supervision (see tomography plots in Appendix F.8 showcasing improvements due to higher-derivative loss inclusion), the learned fields achieve strong agreement with analytic solutions across the domain, except near singularities lim (1+) Vn N, where curvature becomes infinite.
EXPERIMENTAL VERIFICATIONS
The performance of EinFields is assessed along two axes: (i) compression, i.e., meaningful reduction in permanent storage requirements as compared to high-resolution spatiotemporal meshes utilized in NR simulations which includes the metric (10 independent components) and higher-order derivatives (20-100 components) across millions of collocation points, and (ii) reconstruction fidelity, evaluated through key GR benchmarks: geodesic dynamics around compact objects (Appendix D) and curvature diagnostics such as the curvature scalars. The evaluation criteria used either mean absolute error (MAE) or relative lz (Rel. &g) between the ground truth and NeF parametrized tensors, which is detailed in Appendix F.I.
4D training and validation data. We overfit the NeFs over synthetic data generated from exact 4D analytic solutions (see Appendix B explaining each of these use-cases) of the EFEs: (i) Schwarzschild
5

---

Published as a conference paper at ICLR 2026
Metric
g
Covariant derivative jacfwd +1'
einsum jacfud
Lie derivative jacfwd +
Christof symbol 1s
jacfwd +
Riemann tensor Riem
Bianchi II (jacfwd + r)Riem
Weyl tensor
einsum
Ricci tensor Ric
einsum
Ricci scalar
Figure 2: The directed-acyclic graph (DAG) for computing the differential geometric quantities from the metric tensor in analogy to Figure 10 and Eq. (3) The transformations include repeated differ- entiation implemented via forward mowddd-mooe Jacobian jacfwd operations and tensor index manipulation using einsum. Tensors are in depicted in teal blue, connection in light-blue, tensor derivatives in green and conservation laws (Bianchi identities (jacfud + Riem 0) in red.
(static, spherically symmetric solution), (ii) Kerr (rotating, with spin parameter a O, oblate spheroidal), and (iii) propagating gravitational waves (GW) (time-varying, linearized gravity metric). Details regarding training and validation grid resolutions and parameter ranges (e.g., mass M, spin parameter a, etc) used to generate the distortion part of the metric, as shown in Eq. 5, are described in Appendix F.2.
Training specifics. For our tasks, the most effective architectures are MLPs with SiLU activa- tions (Elfwing et al., 2018), which excel under derivative- supervision. Given the sensitivity of training dynamics to the choice of optimizer (Wang al., 2025a), employ SOAP (Vyas et al., 2025), a scalable quasi-Newton method shown to enhance gradient aligument in PINNs (Raissi et al., 2019). We adopt a GradNormmbassed scheme (Chen al., 2018), enforcing unit-norm gradients mitigating gradient imbalances induced by Sobolev supervision. The explored models span widths and depths from 64 3 to 512 totaling less than
10 parameters (~ 7 Mi The NeF training ranges between 100s (w/o Sobolev training) to 2000s (Sobolev training including metric Hessian) on NVIDIA H200 SXM GPU.
4.1 ACCURACY AND STORAGE EFFICIENCY OF EINFIELDS: METRIC AND DERIVED QUANTITIES.
Table 1: Performance evaluation (measured in Rel. l2 and MAE metries) and storage efficiency of BinFielde parametrized metric tensor fields under different representations (i.e., with and w/o Sobolev trainings). The model with the lowest MAE is selected in each row.
Representation
Rel. &
MAE
Storage
Compression
EinFields
(1.08 1 0.06)e-6
2.11e-6 0..07e-6
85 KiB
4035
EinFields (+ Jac)
(3.37 1 0.84)e-7
(9.49 1 1.51)e-7
1.1 MiB
311
EinFields (+ Jac Hess)
(1.88 = 0.16)e-7
(9.07 = 1.71)e-7
202 KiB
1698
Explicit grid
343 M13
Accurate reconstruction of higher-rank tensors from neural tensor fields is critical for recovering geodesics, tidal forces, and related physical quantities. We evaluate EinFields by comparing accuracy-memory tradeoffs for the metric (Figure 3a) and Christoffel symbols (Figure 3b), us- ing the evaluation numbers reported in Table 1. Against higher-order FD baselines in FLOAT32, EinFields ids achieves systematically lower MAE, avoiding truncation and instability issues that limit FD stencils at small h. Through AD, we compute different derived quantities listed in Table showing that Ei nF i elds outperforms FD stencils by 10 10 in accuracy across these quantities.

---

Published as a conference paper at ICLR 2026
Table Performance evaluation of EinFields reconstructed differential geometric quantities for the Schwarzschild geometry in spherical coordinates. Columns 2-3 report memory usage for full and symmetry-reduued components. Columns 4-6 report MAE relative to analytical solutions: FD stencils for h 0.01 the ground-truth (GT (FD)), EinF EinFielde ields via AD, and AD applied directly to the analytic solution (GT (AD)).
Geometric quantity
Storage [GiB]
MAE
Full
Sym.
GT (FD)
EinFields (AD)
GT (AD)
Christoffel symbol
2.6
1.6
5.37e-6
(9.98 + 2.12)e-7
5.83e-9
Riemann tensor
10.4
0.8
1.78e-2
(1.25 = e.30)e-6
2.86e-8
Weyl tensor
10.4
4.0
1.72e-2
(1.67 1.11)e5
5.89e-8
Ricci tensor
0.6
0.4
4.81e-2
(9.66 = 2.8))--6
9.02c-8
Ricci scalar
0.04
0.04
5.35e-2
(3.80 = 1.72)e-5
1.31e-8
Kretschmann invariant
0.04
0.04
1.33e-2
(1.07 0.46)e-5
3.32c-8
Bianchi identity (II)
1.68e-2
5.00e-8
4.81e-8
FinFissds EinFields (+ Jac)
fubbeeeeee TDbittherdde
10-2
Fields (+ Jac Hess)
Mean absolute error
10-5
Mean absolute error
10-3
Timiaeid i+ Tied
10
10-5
10
10-7
10-6
10
10
10a
10a
10'
10a
10a
107
10"
Memory (KiBs)
Memory (KiBs)
(a) Accuracy of K1 nFi e da for different NN sizes (measured in KiBs required to store all #1 1ACA pa- rameters) trendlines for different training schemes. Each metric tensor component
(b) Accuracy of RinFie ds' Christoffel symbols derived from the metric tensor, shown trendlines for different training schemes in 1OAT33 (KiBs). The trendlines show different training schemes (with and without Sobolev training) with MAE values ranging from le-2 to le-6. We benchmark against fourth-order and sixth-order stencils with truncation errors O(h²) and Ohh"), ), respectively. Our framework outperforms FD stencils by more than an order magnitude in accuracy.
MAE values rang- ing from le-5 to le-6. Apart from high accuracy, one additionally acquires 1000 4000 times compres- sion factors in storage memory gain, as detailed in Table
Figure 3: Trendlines of accuracy versus storage memory (KiB) requirement for the met- ric tensor and Christoffel symbols. For the explicit grid storage this is computed as num of grid collocation points 4, with 4 bytes for single precision (FLOAT32). For the NeFs, this corresponds to the storage memory of the compact implicit NN weights.
Reconstructing seminal tests of GR. As part of validation, we demonstrate high-fidelity recon- struction of seminal tests associated with general relativistic dynamics on curved manifolds: (i) geodesics curves around Schwarzschild black hole Figs.(4a, 4b, 4c) and its ray-traced rendering Fig. 6; (ii) Kerr solutions (4d, 4e, 4f); (iii) geodesic deviation describing oscillating ring of test particles due to GW distortions Fig. 5 (all detailed within Appendix F.6). Each of these use-cases shows excellent agreement with the analytic results, although they are subject to accumulated tempo- ral rollout errors (see Appendix F.6 for specifics) and are heavily affected by floating point errors requiring FLOAT64 precision.
7

---

Published as a conference paper at ICLR 2026
v/GMtc
(GMic)
(GMIc
(GMiC²)
(GMMc")
(a) Perihelion precession orbits.
(b) Circular orbit for = 3.85r..
(c) Eccentric orbits.
15
10
10
(GGMe)
(GMMt)
(GMich
-5
-5
-10
-10
-15
-15
-10
10
-10
10
-10
10
(GM/c3)
(GMic²)
(GM/cC)
(d) Retrograde orbit for 0.95.
(e) Prograde orbit for 0.90.
(f) Eccentric orbit for 0.628.
Figure 4: Row 1: Geodesics in Schwarzschild spacetime simulated in spherical coordinates -Eq. (70). Row 2: Geodesics in Kerr spacetime simulated in Boyer-Lindquist coordinates Eq. (80). Distinct regions of the geometry are indicated in solid colors. Green solid lines represent ground-truth geodesics, while the red dotted lines represent our Neñs reconstructed orbits.
0.90
t= 1.80
C 2.69
Oragedd Auabie
Drigted Asabii
Octiiial skop AakeSE
00
Asdye
an
-14
X
3.59
t= 4,49
5.39
Orizes Z2g
Gricss cog Acabie
Origtad ragg Acabts
Ociiial dkg Anaust
>,
Analyic
-3.0
+5.
X
Figure 5: Spatial deformations (stretching and squeezing) of circular ring of test particles due to "+" polarized gravitational wave Eq. (111). The NeF-reconstructed h4 cos (wit - =)) and h cos (w(t - =)) show excellent agreement with the analytic geodesic deviation for the linearized gravity use case. See Table 10 for quantitative evaluation.
Curvature associated reconstruction. Additionally, we demonstrate figh-precision reconstruction results for other seminal GR phenomena such as gravitational waves extraction via Weyl sealar Va(r, l), Kerr metric fing-singulariy structure (Kretschmann scalar) captured by BinFielde These are discussed in detail within Appendix F.6.
Oscillating neutron star NR simulation. Beyond analytical solutions, we evaluate EinF ds on widely used, dynamical test in numerical relativity: the oscillatory evolution of perturbed neutron star. Unlike the previous cases, this problem time-dependent, involves matter-spacetime coupling, has no analytical solution, and is computed using fixed mesh refinement (FMR), thus providing a more realistic NR setting for assessing model performance. The details of this setup are
8

---

Published as a conference paper at ICLR 2026
Figure 6: Neural rendering of a Schwarzschild black-hole in front celestial background. The render is constructed by tracing the geodesics of EinFie elds represented metric showcasing its compatibility with complex downstream tasks.
provided in Appendix G and summarized here. We perform simulation of a non-rotating neutron star of gravitational mass M 1.4 Mo. described by the Tolman Seeenheimer-Vookooe (TOV) equations, which is evolved under a small initial perturbation. The coupled evolution of relativistic hydrodynamics and spacetime produces the characteristic oscillation spectrum of the star. This test serves as a standard benchmark for general-relativiste simulations in the Einstein Toolkit (Löffler al., 2012).
Fixed mesh refined training data. The Einstein Toolkit software performs time-evolution using the BSSN formulation of Einstein's ' equations (Shibata & Nakamura, 1995; Baumgarte & Shapiro, 1998). The spacetime domain is discretized using fixed mesh refinement (FMR) (Schnetter et al., 2004; Hayashi et al., 2025). which hierarchy of nested grids provides higher resolution only where needed. The simulation employs five refinement levels, frl0.... r14] (cf. Fig 7b), with rl0 covering the full domain and r14 resolving the stellar interior. The evolution proceeds final time of 7 1000M, with data written every A/ = 1.00 finer levels use proportionally smaller timesteps to satisfy the CFL condition (Courant al., 1967). Training data (detailed in Table 12) at each time slice are obtained by collecting all spatial points (æ4,, 34, 2;) from the refinement levels while discarding any duplicate points lying inside finer patches. This produces a single non-uniform multi-resolution grid.
EinFields
250
Relative fzerror
5x
200
>. 150
100
10-
50
100
150
200
250
300
10"
102
X
Memory (KiBs)
ri=2
d=3
d-1
(a) Accuracy vs. compression of FinFields trained on the NR simulation data (w/ Sobolev training). The trendline indicates maximum neural compression of 2000x for Rel.6, of 3.60e-5, detailed in Table 3.
(b) FMR (excluding ghost zones) increases resolution near the center of the neutron star. The refinement lev- els (rl) are evolved independently with progressively finer Ax and At.
Figure 7: FinF SinFielde ds applied to the NR simulation of an oscillating neutron star evolved numerically using a BSSN solver fixed-mesh refinement (FMR) grid.
Training results. Table 3 and Figure 7a summarize the evaluation for the metric, comprising of compression and accuracy for the best-performing 6x256 model. The evaluation procedure

---

Published as a conference paper at ICLR 2026
explained in detail in Appendix G. For completeness, the corresponding Christoffel symbol results are also reported in Table 13.
Table 3: Performance of FinFields tested on single stable neutron star numerical relativity simulation. The table reports the relative & error, MAE, storage footprint, and resulting compression ratio obtained when training EinFields on the coalesced excd-mesh-refined (FMR) grid constructed from the simulation. EinFields achieves compression ratio of approximately 2000 while maintaining low reconstruction error, as indicated by the reported relative lz and MAE values.
Representation
Rel. l2
MAE
Storage
Compression
EinFields
3.60e-5
5.98e-5
1.4 MiB
2121
EinFields Jacobian)
6.95c-6
9.88c-6
1.4 MiB
2121
FMR coalesced grid
2.9 GiB
4.2 ABLATION STUDIES
We report ablation results relative to our best-performing baseline configuration, systematically examining the effects of matrix representations (full metric instead of distortions), activation functions, optimizers, learning rate schedulers, and Sobolev regularization. All evaluations are conducted in spherical coordinates.
Table 4: Ablation results for the Schwarzschild metric. Row 2 trains the full metric (Eq. (70)) instead of its distortion (Eq. (98)). Row and 4 ablate the learning rate schedule and the optimizer, respectively. Row 5 replaces SiLU with WIRE (Saragadam et al., 2023), well perffoerforming activation function for NeFs. Rows 6-7 ablate the derivative supervision, i.c. Sobolev training.
Ablation
Rel &
Wallclock time [s]
Metric distortion, Jac + Hes, SiL.U, SOAP, Cosine LR
Baseline:
1.40e-7
1400
Metric distortion a Metric Jas
2.13e-6
1407
Cosine + Const. LR
2.37e-5
1397
SOAP + ADAM
4.16e-6
1150
SiLU + WIRE
4.12c-6
3045
Jac Hes Jac
1.51e-7
509
Jac Hes ->
2.37e-7
364
5 CONLUSION
EinFields introduces the first implicit neural framework for compressing four four-dimensional rela- tivity simulations with differentiable modeling. By combining neural tensor fields with automatic differentiation, offers scalable, discretization-free resolution-invarian alternative grid-based methods that preserves physics and is suitable for downstream tasks. Our framework achieves accura- cies of 7 1e (see Table 8) with compression factors of 1000 4000 relative to uniform and grids. Derived quantities show up to five orders of magnitude improvement tensor derivatives (in FLOAT32) over higher-order FD schemes.
FMR
Limitations. Compression with NeFs remains lossy: even with FLOAT64 training, Rel. & errors below 1e-9 are currently unattainable. At present, the framework surpasses FD methods only in single precision. These enror propagate into Christoffel symbols, causing long-time divergence in geodesic solvers (see Figures (16a, 16b)) and curvature tensors. Moreover, the query time of these compressed NeF representations is non-trivial (a few ms for a batch of 10 queries, Figure 23), which can be prohibitive in downstream tasks that require repeated sequential evaluation. While taking the first step toward the integration of ML and NR techniques, our work does not evaluate advanced NR methods, such as adaptive mesh refinement and (pseutl-)spectral solvers.
Future work. We plan to extend the application of EinFields to other complex large-scale time- evolving NR simulations (c.g., binary black hole and binary neutron star mergers) and benchmark against advanced classical techniques such (pseudo)spectral methods.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
We sincerely thank Nils Deppe for valuable feedback on several numerical relativity-rclated aspects of the paper.
The ELLIS Unit Linz, the LIT Al Lab, the Institute for Machine Learning, are supported by the Federal State Upper Austria. We thank the projects FWF AIRI FG 9-N (10.55776/FG9), AlatnrenHeatingrrids (FFG- 899943), Stars4Waters HEORZZONCCLL222.CLLLAAEEEEEE 202--CLIMAATEE-01000 We thank NX GmbH, Audi AG, Silicon Austria Labs (SAL), Merek Healthcare KGaA, GLS (Univ. Waterloo), TÜV Holding GmbH, Software Competence Center Hagenberg GmbH, dSPACE GmbH, TRUMPF SE Co. KG.
Sandeep S. Cranganore was supported by the FWF Bilateral Artificial Intelligence initiative under Grant Agreement number 10.55776/C0E12.
REPRODUCIBILIIYY STATEMENT
We have made every effort to ensure the reproducibility of our results. The codebase, including training scripts, neural field models used, and data generation and preprocessing pipelines, will be made accessible as zip file in the supplementary material. All experiments can be reproduced using the instructions provided in the repository's 's README. md and _to__t__rre cetttririnn... EinF ields mnd, with detailed specifications of hyperparameters, optimizer settings etc.
For our synthetically generated analytic solutions data, we provide the essential configuration yaml files with appropriate parameter values in Data_generrs ion configs and the data generation scripts are within data_genera generetttttt ion. Copious example notebooks containing all the validation problems contained with the folder wxwampl__ootebooeeeeeeeeeoo visualization can be found within the bh_render folder.
The blackhole render scripts and
Software dependencies are specified in a requirement txt file, and we provide Conda virtual environments for case of setup, especially with the appropriate CUDA version. All experiments were run on (specify hardware, c.g., 1x NVIDIA A100 GPUs for prototyping and NVIDIA H200 GPUs for the main runs], with training time ranging from 100 2200 seconds depending on Jacobian, inclusion Hessian-inclusion
Hessian
losses and the specific hardware used.
REFERENCES
B. P. Abbott al. Observation of gravitational waves from binary black hole merger. Phys. Rev. Lett. 116:061102. 102, Feb 2016a. doi: N0.110//PhyReevLet.. 11.061 102.
B. Abbott et Gw150914: The advanced ligo detectors in the era of first discoveries. Phys. Rev. Lett., 116:131103, Mar 2016b. doi: 10. 10.1103/hysRevveett. 13.111103..
B. P. Abbott et al. Properties of the binary black hole merger gw 15014. Phys. Rev. Lett., 116:241102. Jun 2016c. doi: 10. 10.h 16..21100.
B. P. Abbott et al. Multi-messenger observations of a binary neutron star merger*. The Astrophysical Journal Letters, 848(2):L12, oct 2017. doi: 10.3847/2041 8211-881344 1213//a91c..
R. Arnowitt, S. Deser, and C. Misner. Dynamical structure and definition energy general Dec 1959. doi: 10. 10.110//hyy.Rev.16...222
relativity. Phys. Rev.,
Rafael Ralpolter-Ripoll, Peter Lindstrom, and Renato Pajarola. Tthresh: Tensor compression for multidimensional visual data. JEEE Transactions on Visualization and Computer Graphics, 26(9): 2891-2903, 2020. doi: 10.1 0011009TTCC..201922900
Thomas W. Baumgarte and Stuart L. Shapiro. Numerical integration of einstein's 's field equations. Phys. Rev D, 59:024007, Dec 1998. doi: 10. L.w.://yyRee...9..024.. URL https: aps. ors..or//ee doii/ 10. 1103/ 1011100///////e PhysRevD 59. 024007. www.o////oooo/oo/oooo/ooo/oooooo
Thomas W. Baumgarte and Stuart Shapiro. Numerical Relativity: Starting from Scratch. Cambridge University Press, 2021.
11

---

Published as a conference paper at ICLR 2026
Marsha J Berger and Joseph Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics, 53(3):484-555 512, 1984. ISSN 0021-9991. doi: https://doi.or//100 10.1100..111000
G.D. Birkhoff and R.E. Langer. Relativity and Modern Physics. Harvard University Press, 1923.
Cristian Bodnar, Wessel P Bruinsma, Ana Lucic, Megan Stanley, Anna Allen, Johannes Brandstetter, Patrick Garvan, Maik Riechert, Jonathan A Weyn, Haiyu Dong, Jayesh K Gupta, Kit Thambiratnam, Alexander T Archibald, Chun-Chieh Wu, Elizabeth Heider, Max Welling, Richard E Turner, and Paris Perdikaris. A foundation model for the earth system. Nature, 641(8065): 1880-11877 May 2025.
Raoul Bott and Loring W. Tu, de Rham Theory, pp. 13-88. Springer New York, New York, NY, 1982. ISBN 978-1 --557-33599---00 doi:
Robert H. Boyer and Richard W. Lindquist. Maximal analytic extension the kerr metric. Journal of Mathematical Physics, 8(2):265-28) 02 1967. ISSN 0022-2488. doi: 10. 17.56/91.7775119.
Michael Boyle al. The sxs collaboration catalog of binary black hole simulations. Classical and Quantum Graviry, 36(19): (101::991906,,, sep 2019. doi: 1.1.088/1311-663255442..
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclau- rin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumlPy programs, 2018.
Johannes Brandstetter. Envisioning better benchmarks for machine learning PDE solvers. Nat. Mac. Intell., 7(1):233, 2025. doi: 10. 10.10855425-004400962
Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022.
Michael M. Bronstein, Joan Bruna, Taco Cohen, and Petar Velickovic. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. CoRR, Abs/2104.134733 2021.
Steven L Brumton, Bernd R Noack, and Petros Koumoutsakos. Machine learning for fluid mechanics. Annual review of fluid mechanics, 52(1):4775508 2020.
S. Carroll, S.M. Carroll, and Addison Addison-Wesley. Spacetime and Geometry: An Introduction to General Relativity. Addison Wesley, 2004. ISBN 9780805387322.
S. Chandrasekhar. The Mathematical Theory of Black Holes, pp. 5-26. Springer Netherlands, Dordrecht, 1984. ISBN 978- 94-0994-40999--- doi: 00.1007/978-94 907778899--00--000 609-6469-3_2.
Zhao Chen, Vijay Badrinarayanan, Chen- Yuen-Yuloe Lee, and Andrew Rabinovich. Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks. In International conference machine learning, pp. 794 803. PMLR, 2018
Zhiqin Chen and Hao Zhang. Learning implicit fields for generative shape modeling. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 5939-5948, 2019.
Aditya Chetan, Guandao Yang, Zichen Wang, Steve Marschner, and Bharath Hariharan. Accurate differential operators for neural fields, 2024.
Alan S Cornell, Anele Neube, and Gerhard Harmsen. Using physics-informed neural networks to compute quasinormal modes. Phys. Rev. D. 124047, Dec 2022. doi: 10. 10.1103/PhsseeDD.. 124047. URL https: : / /wi.nniiiii/// aps. doi / 10 110111111101111 /et///hyyyeeiee 106. 124047.
R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathematical physics. IBM Journal of Research and Development, 11(2):215222 2344, 1967. doi:
Wojciech Marian Czarnecki, Simon Osindero, Max Jaderberg. Grzegorz Swirszez, and Razvan Pascanu. Sobolev training for neural networks. In Proceedings of the 31 International Conference on Neural Information Processing Systems, NIPS' 17, pp. 4281-425 4290, Red Hook, NY, USA, 2017.
Curran Associates Inc. ISBN 9781590860964
12

---

Published as a conference paper at ICLR 2026
Maximilian Dax, Stephen R. Green, Jonathan Gair, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf. Real-time gravitational wave science with neural posterior estimation. Phys. Rev. Lett., 127:241103 Dec 2021. doi: 10. 10.110//hyyRReeLett.2222 127.24117..
Maximilian Dax, Stephen R. Green, Jonathan Gair, Nihar Gupte, Michael Pürrer, Vivien Raymond, Jonas Wildberger, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf. Real-time inference for binary neutron star mergers using machine learning. Nature, 639(8053):49-53, Mar 2025. ISSN 1476-4687.
doi: 10. 1038//41 58-805554-5--02225-99955------
Jeffrey Dean et al. Large scale distributed deep networks. In Proceedings of the 26th International Conference on Neural Information Processing Systems- Volume 1, NIPS' 12. pp. 1223-1231, Red Hook, NY, USA, 2012. Curran Associates Inc.
Sheng Di and Franck Cappello. Fast error-bounded lossy hpe data compression with sz. 2016 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pp. 730-739, 2016. URL https: //ap:///ppppee// semant wew.chcccchieee ar. org/ Ctrg/CorpssD:::: 8296694.
Emilien Dupont, Adam Goliáski, Milad Alizadeh, Yee Whye Teh, and Arnaud Doucet. Coin: Compression with implicit neural representations, 2021. URL htttps://ara.ttt www.:.orgaaass..a 2103. 03123.
Stefan Elfwing, Eiji Uchibe, and Kenji Doya. Sigmoid-weighted linear units for neural network function approximation reinforcement learning. Neural Networks, 107:3-11, 2018. doi:
10.1006//NEUUEE2207.220001
Amer Essakine, Yanqi Cheng, Chun-Wun Cheng, Lipei Zhang, Zhongying Deng, Lei Zhu, Carola- Bibiane Schönlieb, and Angelica Aviles-Rivero.. Where do we stand with implicit neural representations? technical and performance survey. arXiv preprint arXiv:241..00 1.03688, 2024.
C. W. F. Everitt, D. B. DeBra, B. W. Parkinson, J. P. Turneaure, J. W. Conklin, M. I. Heifetz, G. M. Keiser, A. Silbergleit, Holmes. J. Kolodziejezak, M. Al-Meshari, V. G. Solomonik, K. Stahl, P. W. Worden, W. Bencze, S. Buchman, B. Clarke, A. Al-Jadaan, H. Al-Jibreen, L. J. A. Lipa, M. Lockhart, Al-Suwaidan, M. Taber, and S. Wang. Gravity probe b: Final results of space experiment to test general relativity. Phys. Rev. Lett., 106:221101, May 2011. doi: 10. 1..1100PyyseeeLctt... 106.22111......
Mester, Muhlfelder,
Bengt Fornberg. Introduction, pp. 1-3, Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 1996.
V. Frolov and I. Novikov. Black Hole Physics: Basic Concepts and New Developments. Fundamental Theories of Physics. Springer Netherlands, 1998. ISBN 97807923514500
Fuerst, V. and Wu, K. Radiation transfer of emission lines curved space-time*. A&A, 424(3): 733-746, 2004. doi: 10. 10.1001//0004666 6361:20035814,
Gianluca Galletti, Gerald Gutenbrunner, Fabian Paischer, Suresh Sandeep Cranganore, William Hornsby, Naomi Carey, Lorenzo Zanisi, Stanislas Pamela, and Johannes Brandstetter. Learning to compress plasma turbulence. In NeurlPS 2025 Al for Science Workshop, 2025. http:::: // /openreview net http///oumiiiiMAtOO..0
A. Griewank and A. Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second Edition. Other Titles in Applied Mathematics. Suciety for Industrial and Applied Mathematics, 2008 ISBN 9780898716597.
Kota Hayashi, Kenta Kiuchi, Koutarou Kyutoku, Yuichiro Sekiguchi, and Masaru Shibata. Jet from binary neutron star merger with prompt black hole formation. Phys. Rev. Lett., 134:211407, May 2025. doi: 10.1103/ 134.21 1407.
Jonathan Heek, Anselm Levskaya, Avital Oliver, Marvin Ritter, Bertrand Rondepicrre, Andreas Steiner, and Mare van Zee. Flax: A neural network library and ecosystem for JAX, 2024.
Edward Hirst, Tancredi Schettini Gherardini, and Alexander G Stapleton. Ainstein: numerical einstein metrics via machine learning. for Science, 1(2)0020001, oct 2025. doi: 10.1088/3050-2877 a 17. URL https: //doi lociioorg/110 http:/111.10///053010 088
13

---

Published as a conference paper at ICLR 2026
Langwen Huang and Torsten Hoefler. Compressing multidimensional weather and climate data into neural networks, 2023. URL https:///// www.i.owg...................oo org/ abs http/abbs/22..... /21...2211.........
Langwen Huang, Luigi Fusco, Florian Scheidl, Jan Zibell, Michael Armand Sprenger, Sebastian Schemm, and Torsten Hoefler. Error bounded compression for weather and climate applications, 2025. URL https::///// arxiv. org/ abs wtw./b/l.......... 222..
Xiaomeng Huang, Yufang Ni, Dexun Chen, Songbin Liu, Haohuan Fu, and Guangwen Yang. Czip: fast lossless compression algorithm for climate data. Int. Parallel Program., 44(6):1248-122 1267, December 2016. ISSN 0885-7458. doi: 10. (0007//00766-01--0000
URL https://das
org/10 1007 s10766-1166666 01 0403-2
E. A. Huerta, Roland Haas, Sarah Habib, Anushri Gupta, Adam Rebei, Vishnu Chavva, Daniel Johnson, Shawn Rosofsky, Erik Wessel, Bhanu Agarwal, Diyu Luo, and Wei Ren. Physics of eccentric binary black hole mergers: A numerical relativity perspective. Phys. Rev. D, 100:064003, Sep 2019, doi: 10.1
Youngsik Hwang and Young Lim. Dual cone gradient descent for training physics-informed neural networks, 2025.
C.J. Isham. Modern Differential Geometry for Physicists. World Scientific lecture notes in physics. World Scientific, 1999. ISBN 97898 10235628.
ISO Central Secretary. Information technology - jpeg 2000 image coding system part Core coding system. Standard ISO/IEC 15444-1:2024, International Organization for Standardization, Geneva, CH, 2024. URL https: /www. iso. www.:iisooo....o.... org/ httg//sandarddd. 8762. html.
Jürgen Jost. Riemannian Geometry and Geometric Analysis. Universitext. Springer Berlin, Heidelberg, Berlin, Heidelberg, 5th edition, 2008. ISBN 978-3-540-77341-2 doi: 10.1007/978-3-54--73311-22
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. informed machine learning. Nature Reviews Physics, 3(6):422 440, June 2021.
R Р Kerr and A Schild. Republication of: A new class of vacuum solutions of the einstein field equations. General Relativity and Gravitation, 41(10)24485-249, October 2009.
Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
Patrick Kidger and Cristian Garcia. Equinox: neural networks JAX via callable PyTrees and filtered transformations. Disserentiable Programming workshop at Neural Information P Processing Systems 2021, 2021.
S. Kobayashi and K. Nomizu. Foundations of Differential Geometry, Number Bd. in Foundations of Differential Geometry. Interscience Publishers, 1963. ISBN 9780470496480.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Arizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: learning maps between function spaces with applications to pdes. J. Mach. Learn. Res., 24(1), January 2023. ISSN 1532-4435.
Andrzej Krasiński. Ellipsoidal space-times, sources for the kerr metric, Annals of Physics, 112(1): 22 22-40. 40, 1978. ISSN 0003-4916 16. doi: https://doi.org/10.. 101.11166//0 0011/0003-49167788 167869999907999
V. Venkatraman Krishnan, M. Bailes, W. Straten, N. Wex, P. C. C. Freire, E. F. Keane, T. M. Tauris, P. A. Rosado, N. D. R. Bhat, C. Flynn, A. Jameson, and S. Oslowski. Lense-thirring frame dragging induced by fast-rotating white dwarf in binary pulsar system. Science, (6777)77): 577-580, 2020. doi: 10.1 10.1126/siience.axx70000
J. Lee. Introduction to Smooth Manifolds. Graduate Texts in Mathematics. Springer New York, 2012. ISBN 9781441999825.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadeneshelii Burigede Liu, Kaushik Bhat- tacharya, Andrew M. Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.
14

---

Published as a conference paper at ICLR 2026
Peter Lindstrom and Martin Isenburg. Fast and efficient compression of Soaating-pinne data. IEEE Transactions on Visualization and Computer Graphics, 12(5):1245-122 2006. doi: 10.1109/ TVCG.2006.144 143.
Bo Liu, Xingchao Liu, Xiaojie Jin, Peter Stone, and Qiang Liu. Conflict-averse gradient descent for multi-task learning, 2021.
Geoffrey Lovelace. Computational challenges in numerical relativity in the gravitational-wave era. Nature Computational Science, (()):450-452 July 2021.
Qiang Lui, Mengyu Chu, and Nils Thuerey. Config: Towards conflict-free training of physics informed neural networks, 2025.
Raimon Luna, Juan Calderón Bustillo, Juan José Seoane Martínez, Alejandro Torres-Forné, and José A. Font. Solving the teukolsky equation with physics-informed neural networks. Phys. Rev. D. 107:064025,, Mar 2023, doi: 10. 10.100/hhysevvD.1000 0740050 URL https: htapps///likkk /tt///inn....... aps. brg/doi/10/ 1103/ 3/P.///y/seaa 107 064025.
Frank Löffler, Joshua Faber, Eloisa Bentivegna, Tanja Bode, Peter Diener, Roland Haas, Ian Hinder, Bruno C Mundim, Christian D Ott, Erik Schnetter, Gabrielle Allen, Manuela Campanelli, and Pablo Laguna. The einstein toolkit: community computational infrastructure for relativistic astrophysics. Classical and Quantum Gravity, 15000) may 2012. doi:
URL https: /xxxxxx/o////// doi. org /1/
10. 1088 02264 21-4//9///////11000000
C.W. Misner, K.S. Thore,, J.A. Wheeler, and D.I. Kaiser. Gravitation. Princeton University Press, 2017. ISBN 978069
Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller. Instant neural graphics primitives with multiresolution hash encoding ACM Transactions on Graphics, 41(4):1-15, July 2022. ISSN 1557-7368.
doi: 10. 10.1145/332 10.1145/3222233530127
Ezra Newman and Roger Penrose. An approach to gravitational radiation by method of spin coefficients. Journal of Mathematical Physics, 3(3):566-578, 05 1962. ISSN 0022-2488. doi:
17244577
Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove. DeepSDF: Learning continuous signed distance functions for shape representation. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp. 165-174, 2019.
Eric Poisson. A Relativist's 's Toolkit: The Mathematics of Black-Hole Mechanics. Cambridge University Press, 2004.
M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal Computational Physics, 378:686-707, 2019. ISSN 0021-9991. doi: https://doi.or//I 10. 20006//...p.2111110000
Alireza Rashti, Rossella Gamba, Koustav Chandra, David Radice, Boris Daszuta, William Cook, and Sebastiano Bernuzzi. Binary black hole waveforms from high-resolution gr-athena++ simulations. Phys. Rev. D, 111:104078, May 2025, doi: 10. 10.1103//p2--9
Daniel A. Reed and Jack Dongarra. Exascale computing and big data. Commun. 587)5667) June 2015. ISSN 0001-0782, doi: 10. 10.1145/2699414.
Daniel Rho, Byeonghyeon Lee, Seungtae Nam, Joo Chan Lee, Jong Hwan Ko, and Eunbyung Park. Masked wavelet representation for compact neural radiance fields, 2023. URL https: //a //arxiv.org/abll /2122 2.1....00 09069.
Shawn G Rosofsky and E A Huerta. Magnetohydrodynamics with physics informed neural operators. Machine Learning: Science and Technology, 4(3):035002, July 2023. ISSN 2632-2153. doi: 10. 1088/2632-215//acc00a URL http: /ttt::///xx doi.iii.oo.........oo....o...... 100 10w..118//222222 2153 30a.
15

---

Published as a conference paper at ICLR 2026
Ian Ruchlin, Zachariah B. Etienne, and Thomas W. Baumgarte. relativity in singular curvilinear coordinate systems. Phys. Rev. D, 97:064036, Mar 2018. doi: 10.1005/hyyxecxD997.04600
Numerical
Vishwanath Saragadam, Daniel LeJeune, Jasper Tan, Guha Balakrishnan, Ashok Veeraraghavan, and Richard G. Baraniuk. Wire: Wavelet implicit neural representations, 2023.
Mark A. Scheel al. The sxs collaboration's 's third catalog of binary black hole simulations, 2025.
Erik Schnetter, Scott H Hawley, and lan Hawke. Evolutions in 3d numerical relativity using fixed mesh refinement. Classical and Quantum Gravity, 21(6): feb 2004. doi: 10. (0.108800264-381 21/6/014. URL https: //oii org/ 100. 0800088/000 04 9381/21/ 6/014/
Guangyuan Shi, Qimai Li, Wenlong Zhang, Jiaxin Chen, and Xiao-Ming Wu. Recon: Reducing conflicting gradients from the root for multi-task learning, 2023.
Masaru Shibata and Takashi Nakamura. Evolution of three-dimensional gravitational waves: Har- monic slicing case. Phys. Rev. D, 52:5428-5444 Nov 1995. doi: 10. 10//PhyyRee..52255222
Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wet- zstein. Implicit neural representations with periodic activation functions. Hugo Larochelle, Marc' Aurelio Ranzato, Raia Hadsell, Maria-Fles Aaoiiaalorina Balcan, and Hsuan- Tien Lin (eds.), Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Process- ing Systems 2020, NeurlPS 2020, December 6-12, 2020, viriual, 2020.
Hwijae Son, Jin Woo Jang, Woo Jin Han, and Hyung Hwang. Sobolev training for physics informed neural networks, 2021.
Terence Tao. What is gauge?, 2008
Edward Teo. Spherical photon orbits around a kerr black hole. General Relativity and Gravitation,
(9(1)):909-1126, November 2003.
Saul Teukolsky. The kerr metric. Classical and Quantum Gravity, 32(12):124006, jun 2015. doi: 10.1088/0264-938111 12111//64499911/2///2222240
Jonathan Thorburg. Black-hole excision with multiple grid patches. Classical and Quantum Gravity, 21(15):3665, Jul 2004. ISSN 0264-9381. doi: 10.
URL
https://dos worg/orgggoo 1088/ 2112224499922211////1
Nils Thuerey, Philipp Holl, Maximilian Mueller. Patrick Schnell, Felix Trost, and Kiwon Um. Physics-based deep learning. arXiv preprint arXiv:2109.05237, 2021.
Matt Visser. The kerr spacetime: A brief introduction, 2008
Nikhil Vyas, Depen Morwani, Rosie Zhao, Mujin Kwun, Itai Shapira, David Brandfonbrener, Lucas Janson, and Sham Kakade. Soap: Improving and stabilizing shampoo using adam, 2025.
Sifan Wang, Ananyae Kumar Bhartari, Bowen Li, and Paris Perdikaris. Gradient alignment in ihysice-informed neural networks: second-order optimization perspective, 2025a.
Yongji Wang, Mehdi Bennani, James Martens, Sébastien Racanière, Sam Blackwell, Alex Matthews, Stanislav Nikolov, Gonzalo Cao-Labora, Daniel S. Park, Martin Arjovsky, Daniel Worrall, Chongli Qin, Ferran Alet, Borislav Kozlovskii, Nenad Tomašev, Alex Davies, Pushmeet Kohli, Tristan Buckmaster, Bogdan Georgiev, Javier Gómez-Serrano, Ray Jiang, and Ching Yao Lai. Discovery of unstable singularities, 2025b. URL httphhtts/sc/ http://aaxvvvv.. org/ hta/abb//250000 /2509.1991...
Yao
Maurice Weiler, Patrick Forré, Erik Verlinde, and Max Welling. Equivariant and Coordinate Independent Convolutional Nenvorks. 2023.
Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. Neural fields in visual computing and beyond. CoRR, abs/21 1 11.11444 2021.
16

---

Published as a conference paper at ICLR 2026
Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn. Gradient surgery for multi-task learning, 2020.
Xuan Zhang, Limei Wang, Jacob Helwig, Youzhi Luo, Cong Fu, Yaochen Xie, Meng Liu, Yuchao Lin, Zhao Xu, Keqiang Yan, et al. Artificial intelligence for science in quantum, atomistic, and continuum systems. arXiv preprint arXiv:2307.0842. 2023.
17

---

Published as a conference paper at ICLR 2026
Appendices
A INTRODUCTION TO GENERAL RELATIVITY
This appendix provides the mathematical background and intuition differential geometry (covering every aspect of the paper and the library), and general relativity. We remark that the reader may appreciate several related works, such (Jost, 2008; Kobayashi & Nomizu, 1963; Isham, 1999; Lee, 2012) for mathematically rigorous coverage of differential geometry. For more physics-oriented readers, the following books extensively cover general relativity and numerical relativity (Misner et al., 2017; Carroll et al., 2004; Poisson, 2004) as alternative resources. Additionally, the Geometric Deep Learning (GDL) community can also find more ML-centric introduction to differential geometry in the following work (Bronstein et al., 2021; Weiler et al., 2023).
Table 5: Table of notations
Symbol
Description
M
Arbitrary manifold
A
4-dimensional spacetime manifold
Tu
Flat Lorentzian metric
x!
Original coordinates
Transformed coordinates
Basis
gi
Dual basis set
ey
Transformed basis
Transformed dual basis set
a
Coordinate basis (equivalent to partial derivative operator)
azn T T
Tangent space at point p Cotangent space at point
T.M
(M)
Space of one-forms
T(TM)
Smooth sections of tangent bundle (collection vector fields)
M(TM)
Smooth sections of cotangent bundle (collection of one-forms)
$
Pullback operation
4.
Pushforward operation
Riem(M)
Set (pseudo-) Riom-nniannnnn metrics on M
Diff (M)
Set of diffeomorphism maps on M
Cartesian (tensor) product
Kronecker (tensor) product
D
Directional derivative
C,
Lie derivative with respect vector field Covariant derivative Kronecker delta (identity matrix)
V
8
C
Speed of light
G
Newton's constant
A.I FUNDAMENTAL CONCEPTS OF DIFFERENTIAL GEOMETRY & AND TENSOR CALCULUS
The main concepts covered in this appendix are:
1. Fundamental concepts differential geometry and tensor calculus: We introduce contravari- ance and covariance, and further vector and dual vector spaces. This allows us to define tangent and cotangent spaces.
2. Tensors and tensor fields:
Next, define tensors and tensor fields, operations tensor fields, and the Lie derivative as generalization of the directional derivative for tensor fields.
3. Riemannian and Lorentzian geometry: This is the meat of Appendix A. We introduce 4-dimensional spacetime continuous differentiable manifold. Via the metric,
18

---

Published as a conference paper at ICLR 2026
define Riemannian manifolds, and finally Lorentzian manifolds manifold. Next, we discuss connections, covariant derivatives, and Christoffel symbols. This is all mathematical background that is required to introduce parallel transport, geodesics, the Riemann curvature tensor, the Ricci tensor, the Ricci scalar, the Weyl tensor, and finally curvature invariants and the wress-energy-momeente equations, reflecting back on the goordinate-indeeendeen of GR.
sreudo-Riemannin
tensor. We end with Einstein field
A.1.1 CONTRAVARIANT AND COVARIANT COMPONENTS
Loosely speaking, n- 7 dimensiona vector E RR can be expanded in its basis as v = vtei +u ... en. or u = n'e, if we use Einstein sum convention. In general relativity, we write v = e where Greek indices indicate 4- dimensional space-time Thus, in this non-Euclidean setting, is necessary to distinguish objects that carry an upper index (contravariant) versus objects that carry a lower index (covariant), since they satisfy different geometric properties and transformation laws.
Definition 1 (Contravariance of vector components): Let X C RR be a coordinate system (frame) that spanned by coordinate basis set { i.e., each basis vector can be expressed as ey Dzn vector E X can be expanded in its coordinate basis When transforming the vector to new coordinate system, spanned by another coordinate basis set [Sw)s i.e., e22 8g* one can express the vector components the new coordinate system 8 = (T) as
2
v (z)e,, (xx)) Uz
Or vn (x)
t"(T)
Σ Dxu
v²(x).
The ratio of change of the vector components is the inverse of the ratio of the hase components. In other words, vector components transform inversely or contravarianlly -ittn respect to basis transformations, i.e., transform in the opposite way to the change the coordinate system. Most contravariant objects represent physical quantities like displacement, velocity, and momentum, which must adjust when the coordinate basis changes.
Definition 2 (Covariance of basis set): Let X C R" be a coordinate system (frame) that is spanned by a coordinate basis set {E Isugn i.e., each basis vector can be expressed as e, = az vector v E can be expanded in its coordinate basis as v = (z()), == (*() dzn When transforming the vector u to a new coordinate system, spanned by another coordinate basis set (Evigugn i.e., Ey Der then, the basis set itself transforms as,
2_
Ox
Σ
(8)
Or"
Oz Ozy
Note that we have introduced the concept contravariant and covariant transformation by the example of vector components and the respective basis set. In general, we speak of contravariant w.r.t. their corresponding basis sets. Lee., contravariant components have covariant basis sets and covariant components have contravariant basis sets. As we introduce next, an object with covariant components an object of the dual space. These covariant vectors, or covectors, typically represent gradients, such the gradient of function. A gradient represents the change w.r.t. an infinitesimal change in direction. It is intuitive that if we make the direction larger, the change becomes larger as well. In other words, if change the basis vectors in which measure this change, the gradient transforms covariantly w.r.t. the basis vectors.
A.1.2 DUAL SPACE
While the concept of vector space well known in the machine learning community, there is closely associated concept of dual vector space (succinctly called dual space), which is an algebraic dual to the vector space itself with the same dimensions.
Definition 2 (Dual vector space): Let (V. +. :) be vector space over field F (e.g., R, C). The (algebraic) dual space (V*, +. :) is vector space of linear functionals (maps) y* {0" (v)
19

---

Published as a conference paper at ICLR 2026
Ve E V, Ve E F). satisfying:
(0* + w") (c) v²(u) + w (u) 8¹ (av Sw) = an" (n)(v) + Bwr
(9)
(10)
(cv (v) = c(v" (v))
(11)
for all ,,,w" EV", V, WVW a, 3, E F. Elements the dual space V are sometimes referred to as covectors one-forms Bott & Tu (1982).
For example, suppose, we are given basis fei...., en) of a vector space V. Then, one can introduce a dual basis set {0), ....., (²)) of the dual space V*. Let v = (aye) 02262 anC,) E V. fatlisisn E F. Thus, the action of the linear functionals O' on the vector reads: (v) = (a1e) + 0222 +... + anCn) = (e:)(ee) = a², i= L...., n
and (eee)) 6j. where we have used the orthonormality condition, and 03 is the Kronecker delta symbol. Conceptually, the covector & is (complex-conjugated if F C) row-vector, which acts column vector v produce E F. Colloquially speaking, an element of the dual space V "eats up" an element of the vector space V and returns a scalar (duality pairing).
A.1.3 TANGENT AND COTANGENT SPACES
Definition (Tangent space): Lel M be smooth (C*) manifold of dimension nt The tangent space M at point M is set of d-dimensional vectors (called langent vectors) allached at RR,, and carries the structure of real vector space.
point p. defined T.M == {(p v)
Every tangent space is spanned by an ordered basis (eulphissn {600 and vectors be expanded in this basis
8
} E M
p
s/p r" (x)
(12)
dr p
where ("(2)) (x) are the components of the vector in this basis (eplphsusn of the tangent space T.MM
It is worth noting that definition of tangent vectors.
= dim(M). With this setup, we can now formally introduce the
Definition (Tangent vector): vector dlp E Ty is called as a tangent vector if it acts as a derivation, i.e, linear map acting smooth functions C (M (M) at point E M. Specifically,
the map vlp C (M, M RR) + R satisfies 3;
v(f + g) v(f) v(g) v f.g C (M,R)) (linearity) iî) v(f) = 0, when f is a constant function, iee. v acts trivially on constants.
fii) v(fg) (g))(()) + (())x(f) Vf.g E C (M, (Leibniz product rule)
From the definition above, follows that tangent vectors should be regarded as derivation maps. This is equivalent to the notion of a "directional derivative" (Isham, 1999)
d s(v J( tw/p) 1=
0f Drt
(D.N) = ti.
(x()
Vu = -
ETAM, E
(13)
Directional derivatives are traditionally defined only for scalar-valued functions. This shall be revisited rigorously for a more generalized concept called the "Lie-derivatives", which operates on general tensors, c.f. Section A.2.3.
"A finite dimensional vector space is isomorphic to its double dual, i.e. (V*)*
For sake of ease, will drop Ip whenever possible.
20

---

Published as a conference paper at ICLR 2026
Thus, Mp.M is the space of directional derivatives. The disjoint union of all the tangent spaces at every point E M forms à structure called rangent bundles Isham (1999):
TM=
L TpM= U upp) upp E TyM}. PEM
(14)
PEM
Tangent vectors vector fields. In physics, quantities that vary spatiotemporally continuum representation are defined as fields, featuring in domains such as electrodynamics, gravity, fluid dynamics, or continuum mechanics.
vector field V is smooth assignment of tangent vector f/n to each point p E .M. Thus, a vector field is map V C (M) + C (M), and is defined as:
(V()) wln(s)-
(15)
Theorem (Cotangent space): Let M be smooth (C )-manifolld (differentiable). The colangent space M == {(p vz (p) Up)) = Vp up E T,M,K E R] at point E M the set of all linear maps v TpM + RR i.e, dual to the tangent space. The cotangent space M is Thus, any v'lp E TpM M can be expanded spanned by an ordered basis set {dz}p d ..... dx as:
v² = - = (x))dx
(16)
p
Oxtt
It follows that dx disjoint union of all the cotangent spaces at every point M are known as cotangent bundles Isham (1999):
o and dim T.M-- dim T.M = dim M. The
Ox"
p
TM- L
TpM U ((.."")): {TM) PEM
(17)
PEM
One can also construct fields of cotangent vectors (cotangent fields) by picking an element of (MM Vp E M in smooth manner. Le.. by assigning one cotangent vector smoothly each point of the manifold, one obtains a cotangent field (i.e., a smooth section of the cotangent bundle). These cotangent fields are known in mathematical literature as one-forms. The set of all smooth one-forms on M is commonly denoted as Ω (M).
A.2 TENSORS AND TENSOR FIELDS
Definition (Tensors): A rank (r, s) tensor T point p E M is described as multilinear map:
T: x..VXV
V"
x
E
(18)
P-copies
8-copies
where x denotes the Cartesian product and the resultant tensor has total rank r A tensor takes covectors and vectors, returning real number, in a multilinear way (linear each argument separately). The and input vectors and covectors pair with the and being the convariant and contravariant components, respectively. Equivalently, a tensor is an element that lives in a tensor product of vector and dual spaces, i..., T E (v) @ (V*)8s A tensor in particular basis choice (Ca,, Jisner E V and [0%a Jisnss E 1 given by
T To102...
B, Cay ②
& Ca. @ g ....
3 0
(19)
where Traaa..... basis set.
B, T(ga1 ...., ga. eg ....... eg,) the coefficients of the tensor w.r.t. the
21

---

Published as a conference paper at ICLR 2026
A.2.1 TENSOR TRANSFORMATION PROPERTIES
pivotal criterion for object to be classified tensor(field) is that transforms according to a well-defined rule under changes of coordinates. Let {ea, Jisgg E V, {0 Jisnga E V", and {ea,, hinngr E V. [f0n) Jigns E V be two coordinate systems on a smooth manifold M, related by smooth invertible map. Consider a tensor field of type (r,s) with components T S1.0,, in the original coordinate system {ea,, E V, { Jignga EV.. Under a change of coordinate systems, the components in the new coordinate system {ea, Jisngr E V, Jignss E V transform according the following tensor transformation law:
THi My..
(x) Jni ... THr "oi.,, (2) (I-1) (.-)).....-))
(20)
where Inh
F are the Jacobian and Jacobian inverse matrices in the 8atk and (J (J-1). 8r Sa
coordinate basis, respectively. The is the contravariant transformation of the contravariant compo- nents of Tar BI...," whereas (J-1) 0: the covariant transformation of the covariant components of 7 3,...,, The indices Pk, label components in the new coordinates and Oye, 3t are dummy indices summed over the old coordinates. A key feature of tensor is that, if it is zero in one coordinate system, is zero in every other coordinate system. This transformation law ensures that the tensorial nature of the object is preserved independent of the coordinate chart chosen.
Tensor fields. A tensor field is collection of tensor-valued rank quantities (r,s) such that each point E M, the multilinear function associates a value Tp V @ (V:) 8* Thus, the components Taa.ar $1...p, (p) are functions of the points of the manifold.
By definition, some known examples of tensor fields in physics and machine learning are:
Rank tensor, e.g., temperature field y : R + R (scalar field)
Rank (1, 0) tensor, e.g., (velocity, momentum, displacement) vector fields vs R" + RR" (contravariant vector field). These rank (1, 0) tensors have one component that transforms contravariantly, and "cats up" covariant component, e.g., v" to produce scalar.
Rank (0, 1) tensor, c.g., gradient vector fields V R + R" (covariant vector field). These rank (0, 1) tensors have one component that transforms covariantly, and "eat up" a contravariant component to produce scalar.
Rank (0, 2) tensor, e.g., matrix representing bilinear form that takes in two vectors and outputs a scalar. We will see the metric tensor Juv as an example. continuum and structural mechanics, known example is the strain tensor Eij representing the deformation of a crystal (body) caused by external forces such as stress.
Rank (2, 0) tensor, c.g., a matrix as multilinear map that takes in two covectors and outputs a scalar. An example for rank (2, 0) tensor is the outer product of two vectors. An example is the Cauchy stress tensor " from structural mechanics, which represents the internal forces per unit area acting inside a material body. The stress tensor takes in two vectors, i.e., the normal vector to the surface (describing orientation), and the direction vector along which the force acts (projection), and returns scalar (force per unit area in that direction).
A.2.2 OPERATIONS ON TENSOR FIELDS
For multiple tensors of the same type (r, s). the algebraic operations such as addition, subtraction or multiplication by functions are straightforward. Here, we address multiplication of tensors of different ranks.
Let T be rank (r, s) tensor and S rank (p. 9) tensor. One can construct tensor product T & S resulting in new tensor of rank (r + p, 8 + q), defined by
T 2 S(ea, ... Co.: en ..... pp 8
8 10l, so0, ...., (d0)).
8,0,, ...., (4))
(21)
TTeaa r. Ca. En......
en ((p
(22)
22

---

Published as a conference paper at ICLR 2026
and the components of this composite tensor read,
a1..arhh.... (T S)
Ta.
Sm..g BI..,,,
(23)
B BI..B.61....6
61..69"
Another useful rule is that of contracting over repeated index/indices each from the vector and dual space respectively. Consider rank (r, tensor
Tmm..ap... ap...
Tr BL...8."
(24)
I.e,, Op is summed-over in the contravariant and covariant indices, and, thus it gets contracted. The resulting tensor is of rank (r 1,3 1).
A.2.3 LIE DERIVATIVE: GENERALIZING THE NOTION OF DIRECTIONAL DERIVATIVES FOR TENSOR FIELDS
Directional derivatives are of great importance and often appear in domains such as fluid dynam- ics, where scalar field is differentiated with respect to a vector flow field, capturing infinitesimal dragging of scalar fields along flows generated by vector field. Flows can be viewed as "diffcomor- phisms" Poisson (2004) induced by these vector fields.
However, generalizing the notion of directional derivatives require defining derivatives of a set of tensor fields of arbitrary rank (r, s) w.r.t. a set of vector fields. This is often not possible on arbitrary manifolds, and requires a concept of differentiating in a tensorial setting. Geometrically, to compare tensors at infinitesimally separated points the tensor from p to (also called parallel transporting, c.f. Section A.3.3.1.
manifold V. say at points p, q M requires to "drag"
Alternatively, simpler approach to describe the dragging is via coordinate transformation from p to g. This is the idea behind the Lie derivative. The Lie derivative along a vector field vlp E TyMM measures by how much the changes in a tensor along differ from mere infinitesimal passive coordinate transformation of the tensor generated by v. In other words, the Lie derivative compares the actual rate of change of the tensor as you move along v against the change you'd 'dd get if everything were just shifted passively via a coordinate transformation. We provide a rough sketch of the derivation, but detailed explanations can be found here (Lee, 2012; Poisson, 2004).
Consider an infinitesimal coordinate transformation which maps the vector with coordinates x'I, at point p ø²I, lg point g:
3 4 x²1 + og wn (z)1p
(25)
It is to explicitly note that the original coordinates x" and the transformed coordinates ø018 are components of the same set of basis vectors. Such transformations fall under the category of active coordinate transformations that map points (or tensors at those points) at old locations to new locations in the old coordinate system in this case by "moving" a small amount og along the vector field dlp E M. In other words, an active coordinate transformation maps points (and tensors) to new locations in the old coordinate system keeping the basis set intact. Whereas, passive transformations assign new coordinates to the old points (and tensors) by transforming the basis set itself.
Assuming a coordinate basis, one can differentiate the transformation w.r.t. the original coordinates, which yields
OT" 58 + 88
(x) dr
(26)
dr
The result contains the identity matrix 5 and a small correction due to the flow field (x)).. To the
dx
Out dx"
first order, the inverse of the above Jacobian is
8
-6
Dzn
The Lie-derivative of tensor field Tу with respect to rl follows similar pattern and is defined via the limes:
TB
Ty
C.TH lim
(27)
60
23

---

Published as a conference paper at ICLR 2026
In this scheme, is important to distinguish three distinct tensor field evaluations: a) (2) (original tensor in untransformed coordinates), b) Т (ā) (transformed tensor in the transformed coordinates) and c) Ty (original tensor in transformed coordinates).
In order to compute Eq. (27), we need two important concepts from differential geometry called push- forward and pull-back operations. We direct the interested readers more advanced literature Isham (1999); Lee (2012); Kobayashi & Nomizu (1963)
These three separate tensors fields can be related in the following manner: Firstly, the tensor mapped to the new set of coordinates 73" (Z) can be obtained via Eq. (20),
Or" Ozo 12(x)
do
(J- [x + oe
T (x) + (
(28)
Ox"
Secondly, the original tensor in transformed coordinates (3) can be evaluated at q, by Taylor expansion:
ØTH
(x) Tt' (x = og 0
drd + O(EE2)
(29)
Substituting Eqs. (28, 29) into the Lie-derivative definition of Eq. (27), and og - 0 one finds the following final expression:
OTH Ox
OV dru (x)
L. 2
Dr To (x)
(30)
pullback
pushforward
The pushforward and pullback operations drag the transformed tensor field onto the original point, where differences can be computed.
Thus, tensors are being compared in the same tangent/cotangent space. Mathematically, for smooth maps* (diffeomorphismss) & .M- /N the pushforward 4 T + ( pushes vector fields forward from one tangent space of a domain T.MM to the tangent space of another tangent space To The pullback, dual linear map to pushforward, drags covectors (one-forms) living in cotangent spaces (6.)* Φ- T (p)) - T. in the reverse direction the domain. Hence, the contributions from the pushforward on the vector field components and pullback on the covector field components jointly determine the structure of the Lie derivative of mixed tensor field, as expressed in Eq. (30). These operations offer a coherent mathematical framework for transitioning between tangent and cotangent bundles mapped onto other tangent and cotangent bundles via smooth maps, acting appropriately on vector fields and one-forms, respectively.
For any arbitrary rank (r, 8) tensor, Eq. (30) can be generalized to:
0 (L,T)) vi...P v Ix
T.....p DL...,
Σ T TH ....... Delts UL....s
Σ Th
Oo
(31)
dx
vi......... Ox
Lie-derivatives do not require the notion of a connection. Connections will be introduced in detail in Section A.3.3 and intuitively stating, connects two distinct Tangent spaces different points, which is not to be confused with pullback operation. Here, is an instructive comparison table for that compares different differentiation schemes: Table 6: Comparison between actions of directional, covariant, and Lie derivatives.
Directional derivative
Covariant derivative
Lie derivative Tensor fields
Feature
Input function
Scalar fields
Tensor fields
Connection dependence
/ (Explicit)
Captures curvature
X
X
Measures
Scalar changes
Intrinsic curvature
Diffeomorphisms (flows)
for e.g., dragging of coordinates in Eq. (25) due to flows induced by vector fields.
24

---

Published as a conference paper at ICLR 2026
A.3 RIEMANNIAN AND LORENTZIAN GEOMETRY
A.3.1 FOUR DIMENSIONAL SPACETIME AS A CONTINUOUS DIFFERENTIABLE MANIFOLD
The fabric of spacetime according general relativity is combination of three-dimensional space and a strictly positively progressing time direction into a single -our-dimensional continuum. Thus, space and time mix between each other through special orthogonal transformations S0(1,3) called the Lorentz transformations. In order to rigorously define the four-dimensional spacetime, it is necessary to define the following:
Definition 4 (Manifold): A n-dimensional manifold M a space, that, locally resembles the n- dimensional Euclidean space R". However, combining these local patches together, globally, the space deviates from R".
Definition 5 (Hausdorff space): Let K. be a topological space. Then K. is said to be Hausdorff space if: For every pair of distinct points z, y E K. with + y, there exist open sets U, C succh that:
x x U, U VEV, and UnV =0
(32)
Definition 6 (Differentiable manifold): An n-dimensional differentiable manifold is Hausdorff topological space K. such that:
Locally K. is homeomorphic to RR. Thus, V K. there is an open set U such that p L and homeomorphism Φ : :u + Z with 2 an open subset RR".
ti) For two subsets Ua and Ug with Ua Nug * 0. the homeomorphisms (topologically isomorphic) da Vua Za and dbs Ug Z3 are compatible, i.e., the map dg Φa da (Ma nU4) 6 (ua nug is smooth (infinitely differentiable (²), and is its inverse map.
The 6₂ are often called charts and a collection (union) of them Ua éa is called an atlas. These charts provides a coordinate system, labeling Ua The coordinate associated p E Ua is:
(p) = (x) .
(p), (p),..., (p))
x²' (p))
Mathematically, the spacetime continuum denoted as M, differentiable manifold with the structure of Hausdorff topological space.
To summarize, a differentiable manifold a space that may be curved or complicated globally, but looks like Euclidean space up close, and allows for smooth calculus to be done on it. The Hausdorff space ensures than one can separate points nicely with open sets. This avoids weird pathological cases and makes limits and continuity well-behaved. Locally Euclidean means that one do calculus as we were on flat space even if the whole space is curved. And finally, the compatibility between overlapping charts ensures that can do calculus consistently across different charts.
A.3.2 METRIC TENSOR
Definition (Metric): metric g is rank (0, 2) tensor field that defined as a symmetric bilinear map that assigns to each p E M a positive-definite inner product g T.M x TrM + IR such that
g(llp, w| plp) g(u, w) = g(w, v) Vv, TyM (symmetric)
ii) For any E M, g(u, w) = 0 Vwlp E TrM implying = 0 (non-degenerate).
Represented in the basis set of the tangent space, the metric components at each point is given by
i
Jur Bup 9p
(33)
Ox P
and the metric can be expanded as,
9 9pd dx & dx²
(34)
25

---

Published as a conference paper at ICLR 2026
Geometrically, the metric defined in Eq. (34) generalizes the notion of distances and induces norm I1-1Ip- TpM + for generic coordinates such as curvilinear and/or manifolds possessing geometries that are intrinsically non-Euclidean in nature, for e.g., spaces of constant positive sectional curvature K ==1 (e.g.. a 2-sphere E embedded in R²',, spaces of constant sectional curvature = -1 such as hyperbolic geometry (Bolyai-Lobachevve) spaces IHT?).
The distance between two points such cases is called the line element, which is defined as,
ds² Buw (x))da dx..
(3:5)
For an n-dimensional manifold M, the metric tensor 9uv is symmetric matrix, g Φb) ( du), with a(n+1) independent components (not necessarily expanded in the coordinate basis):
(on, sho) (01 g
(0n, Φn) (11
On- (n-1)
9pv
(36)
:
(0n-1;
(0)
Definition 8 (Metric bundle): Let M be a smooth manifold and (Io, ,In) be local coordinates
on U C M. The bundle of symmetric (0, -tensos on M is the subbundle Sym²(Tr (T (M) C TO T
In fact, sections ST M) contains the symmetric bilinear forms, i.e. symmetric (0, 2)- tensor fields, and includes the neeudo-Riemannian metrics on M.
Riemannian manifolds. A metric 9 where all diagonal entries of the metric positive, i.e., Jun > 0, = 0,... (MM - 1 is called a Riemannian metric. Thus, manifold M endowed with a Riemannian metric is known as a Riemannian manifold denoted as tuple (M. g) (Jost, 2008).
For the the Euclidean space IR" with Cartesian coordinates representation
dx 2 dx +.... + dx @ dx=
(37)
the metric tensor amounts Ptj dgg and boils down to Pythagoras' theorem. A general Riemannian metric prescribes method to measure the norm of vector as Vg(u, v) ||0|| and also allows for measuring angles between any two vectors at each point cosô = Va(v.r)g( Like any
g(u.w)
other tensor, the components of the metric tensor transform under coordinate change according to Eq. (20):
Jast (7) = [[J 1)]] (x) (I
(38)
Definition 9 (Are length): Let 7 [0, - M be a piecewise smooth curve on a differentiable manifold M, with (() p and (1) q. The velocity vector along the curve denoted by (t), which lives in the tangent space M.
M. If the curve is expressed in local coordinates (t), then
the components of the tangent vector (t) are given by (t) The are length L(y) (distance) of the dt
curve is then defined by
((7) L' [1(())l)a) dt =
=f (on
dxy(t)
dx" dt
dt
(39)
dt
This are length is reparameterization invariant, i.e., does not depend on the choice parameteriza- tion of the curve (t). It is very important result that every smooth manifold admits Riemannian metric.
't is also called action in physics.
26

---

Published as a conference paper at ICLR 2026
Lorentzian manifolds. Unlike Riemannian manifolds, spacetime is actually pseudo- Riaudd--ni manifold that is, the metric is not positive definite. Thus, the underlying metric carries a signature
(-.+++. +), meaning, 9t. < 0. Consequently, spacetime is Lorentzian manifold M, and, forms the basis for electromagnetism and special relativity. The simplest example of a I Lorentzian manifold of arbitrary dimension is the Minkowksi metric, which is flat (meaning eurvature):
7 dx dx dx dx .... + dx= dr==1,
(40)
where, the components of the Minkowksi metric are Mu diag(- 11, +1,..., +1). is possible to find an orthonormal basis {eN} of Tp.M around small neighborhood of point p of Lorentzian manifold such that, "locally", the metric resembles the Minkowski metric
gurlp = Ma
(41)
In the case of Lorentzian manifolds M. the arc length L(7) in Eq. (39) is modified due to non positive-defiittness
T(o) T(o)= 1.¹ V V- do' do'
dxy dxy"(a)) (o" do',
(42)
and is sometimes is referred to as proper-time The minus sign under the square root ensures the integrand is positive for timelike paths, since for timelike intervals, the inner product of the velocity vector with itself (under the Lorentzian metric) is negative, i.e.,
ds Juu dx".
(43)
Natural isomorphism between vector spaces and dual spaces. The metric provides natural isomorphism between vector spaces and dual spaces and allows the switch between contravariant and covariant components', This is done via the following mapping g TT..- TØ M, where each point p, a one-form (covectors) obtained via contraction operation of vector field z/p with the metric g.
and the covector For vector u/p = vd Bz
U"Ip upddot the components are related by U/ = 9pm2"
Since g is non-degenerate, it is invertible. We denote the inverse metric as g such that ghd This is a rank (2, 0) tensor of the form و = grv be raised, e.g. r gd Zv-
Jav = 8 Through the inverse metric indices can oz
8_
Such index contraction rules with the metric apply to tensors of rank (r, s) or even quantities that are not tensors:
( gand, ) sh... n
S....BB Be
Ben
(44)
ai....r
A.3.3 CONNECTIONS & COVARIANT DERIVATIVE
Transporting vector and tensor fields systematically on manifolds requires mapping vector spaces at one point to vector spaces another. While this can be done trivially the Euclidean setting, for Riemannian and Lorentzian manifolds this is a non-trivial since these vector fields and tensor fields live in different vector spaces. This necessitates geometric object that behaves as "connector" between vector spaces. This is achieved via geometric entity called the affine-conecction which is a vector-valued one-form.
Definition 10 (Affine connection): Let M be smooth manifold and I(TM) be the space of vector fields on M, that the space of smooth sections of the tangent bundle (i.e., the collection of all tangent spaces). An affine connection is bilinear map
I(TM) I(TM) I(TM)
(v, w) + V,W.
In general, pscudo-Riemannian manifold has a signature
In the context of numerical relativity the switch hetween contravariant and covariant components is called "raising and lowering indices".
27

---

Published as a conference paper at ICLR 2026
The differential operator V. is the covariant derivative satisfying the following for tangent vectors v, w (short for Ulp. wlp):
v. + 2) = V.w V
ii) fow = fV.w Cf
iii) V ( (D.))w + w. C Dof=u(f) is the directional derivative.
The affine connection completely independent of the metric. However, if manifold is endowed with metric, this enables expressing the connection in terms of the metric. In GR, one looks at special subclass of affine connections called Levi-Civita connection, due to the symmetry property of the metric tensor.
Definition 11 (Levi-Civita connection):: ): An affine connection is an Levi-Civita connection for tangent vectors v, w (short for Ulp. w|p) if:
Veg VU F(TM) (metricity condition)
(45a)
V = [v, w] Vu, w E ((TM) (lorsion-free condition), VeW
(45b)
where, [u, w| (0"0 (1963).
the Lie-bracket of vector fields Kobayashi & Nomizu
Intuition behind the torsion-free condition. Imagine you are moving on smooth surface (like walking on hill), and you have two "directions" Upp and w/p at point p. Now: First move along a tiny bit, then subsequently along u. Alternatively, move along first, then along v, akin to constructing parallelogram. In flat, Fuclidean space, doing these two moves would land you at the same final point, because partial derivatives commute. On curved surface (a manifold M), they don't generally commute you end up slightly shifted. The Lie bracket [v, w] measures how far off (deficit) you are after moving in v and then w, compared to w and then v. captures the monnttmmmutttiiit"" of the vector transport along the two distinct directions, which leads to non- closure of the parallelogram. Thus, the Lie bracket is intrinsic to the manifold, and shows how the transport of u and interact. In torsion-free connections, the "commutation failure" is purely due the manifold's structure not any extra "twisting" introduced by the connection itself.
A.3.3.1 Parallel transport
Definition 12 (Parallel transport): y [0,1] M be smooth curve on the manifold, and let T be a smooth (r, -)-nnn tensor field defined along the curve 7,. The parallel transport of along the curve ((T) is defined by the condition that its directional covariant derivative along the curve's tangent vector vanishes:
V V=(+)===0 =, VT E [0. (t)) E TrM.
(46)
In local coordinates (r)), the parallel transport condition for the components of the tensor field (T) along the curve (t) is:
- dr Dinne (T) + Σ 1 (T) pT..es
d
(T)
Σ I (T) E j=1
Pl...... (T) = O. (47)
These equations a of coupled ODEs, and be solved uniquely for initial condition find a unique vector at each point along the curve ((). This ensures that as the tensor is transported along the curve, its components change in such a way that their covariant rate of change along the curve vanishes.
A.3.3.2 Christoffel symbols The Levi-Civita covariant derivative contains, apart from the partial derivative term, a correction field that calibrates the deficit between vector (tensor) fields transported along a path on the manifold. For basis (e,) } that is transported the covariant derivative is given by.
V ev (x) ( (x)
(48)
28

---

Published as a conference paper at ICLR 2026
The covariant derivative, denoted by Ve, V, - + rr. defines "modified" differentiation operator that preserves tensorial character under general coordinate transformations. The quantities I known as the Christoffel symbols, represent the components of the Levi-Civita connection, which is uniquely determined by the requirement that the connection is torsion-free and compatible with the metric. Notably, these symbols symmetric in their lower indices, i.e., - Tepe
The action of the covariant derivative general tensor field of type (r, s) ensures that derivatives of tensors transform covariantly, thereby extending the notion of differentiation from vector calculus curved manifolds.
① Taiaa (6)
Trraar PI.., (2)
Σ pa (z)7 1...0....... (x)
(49)
Or
BI..,,
81..8.
ro Taina.
mBg (x)
B........
(50)
The action of the covariant derivative on a scalar field, simply reduces to a partial derivative
0Φ(x)
V VαΦ(z)=
dr
Christoffel symbols can be solely expressed in terms of the metric and its partial derivatives:
D (x) (1) = 200h0 (Ouyov(a) + - = (x))
-
(51)
A crucial feature of any connection is that it is not a tensorial quantity. Connections don't obey the transformation law in Eq. (20) under coordinate changes. This can be easily seen through the components of the Christoffel symbols in the coordinate basis:
Oz Ox dx
0x OrP
Pe Juy
(z)
17 (2((x)
(52)
Ou" O5r Oiy
Oe" Ox
tensorial contribation
non- non-tensorial contribution
Christoffel symbols play significant role in defining most stationary trajectories (shortest longest) in the non-Euclidean setting,
Lie derivatives revisited: Levi-Civita connection included. In case of a nonzero Levi-civita connection, the partial derivatives of ordinary Lie derivative in Eq. (31) is replaced by the covariant derivatives:
(L.T) Vi..., = v V Tr
T
Tu.........l
Vouth + Σ TRi
vi...... V vo
(53)
The first term advects (drags) the tensor along the flow of u, i.e., this is the "naive" directional derivative part. The second and third terms account for how the basis vectors themselves changing, due to curvature and due to the vector field v, respectively. It is easy to show that the three terms lead to pair-wise cancellations between the Christoffel symbols present in the three different covariant derivative terms of Eq. (53). Due to this, the whole expression boils down to Eq. (31), thus corroborating the connection independence of this derivative operator. Differently put, in the covariant derivatives, the Christoffel symbols introduce extra terms. However, the extra Christoffel terms cancel out between the different contributions (first, second, and third terms). One ends up getting exactly the same final expression for the Lie derivative as if only partial derivatives had been used this is what is meant by "connection independence".
A.3.4 GEODESIC EQUATION
Geodesics are paths that correspond to the most stationary trajectories (shortest or longest distance) that connect two points p and q manifold. Often, only consider locally distance minimizing curves, and refer to them as geodesics. Geodesics are obtained by solving a calculus variations problem the distance metric L(>) in Eq. (39), i.e.,
29

---

Published as a conference paper at ICLR 2026
8C(7) := 0=
(54)
(or alternatively on T(r) ) in the Lorentzian setting of Eq. (42)). Solving the calculus of variations problem boils down to solving the Euler-Lagrange equations, which mathematically is equivalent to the condition
V (tt) =0
(55)
where (t) is the curve (path) the manifold, (t) the tangent vector (velocity vector), and V the covariant derivative. Eq. (55) intuitively says that the tangent vector is parallel transported along itself meaning, one is moving without "acceleration" relative to the curved space. Thus, parallel transporting the tangent vector along the curve preserves the tangent vector. In numerical relativity, these corresponds the equations of motion, i.e. generalization of the Newton's acceleration deg" equation dr = F# F /m. For full derivations, we direct the readers to refer to Carroll et al. (2004);
Misner al. (2017); Poisson (2004).
We shall present the final form of a very central second-order ODE describing motion (acceleration) of objects executing geodesic paths around heavy gravitating bodies, namely, the geodesic equation
dx 272
dx dx dr dr
I pa
=0
(56)
where, r is some affine paramter (typically, chosen to be the proper-time in Eq. (42)), d /dr² is the aour-accelerttion vector, dr in Eq. (51)).
Arr/dr /dr the four-velocity and T is the Christoffel symbols as seen pa
Importantly, Eq. (55) is the geometric statement of the geodesic equation. It's coordinate- i.e., it's expressed entirely in terms of geometric objects. Eq. (56) is the coordinate version of the same idea. Here, one chooses coordinate system x on the manifold, and the covariant derivative acting on vector becomes the partial derivative plus correction terms involving the Christoffel symbols".
A.3.5 CURVATURE TENSORS AND SCALARS
Curvature tensors arise naturally in differential geometry as tensorial objects that capture the intrinsic and, where appropriate, extrinsic geometric properties of a manifold. They provide a coordinate- independent way to quantify the curvature of space or spacetime by encoding how the geometry deviates from flatness through the second derivatives, constructed out of Hessians of the metric tensor. Unlike artifacts that may arise from curvilinear coordinate choices flat manifolds, curvature tensors reflect the true geometric content of a space. These generalize classical notions such as Gaussian curvature to higher dimensions and arbitrary signature. Being multilinear objects containing several tensor components, curvature tensors systematically characterize the variation of the metric across different directions. We shall introduce the key curvature related quantities, which include the Riemann relevant ones used in our paper in the following section.
A.3.5.1 Riemann curvature tensor The Riemann curvature tensor R R (2) eμ X 87 to @ 0s a rank (1,3) tensor, which quantifies the measure to which a vector that is transported along a small loop (also called holonomy) fails to return to its original orientation due to the effect of the intrinsic curvature that the vector field picks up during the transport. The Riemann curvature tensor is defined via the commutators of covariant derivatives acting on components of vector field e:
[V. V ()) (V.
VAVa) ( (a) R4,,, (z)u (a)).
(57)
The components of the Riemann curvature tensor are expressed in terms of the Christoffel symbols
0r (x) (x)
ar 2 (x)
R (x) Oz 'Ideally, concepts such connections, parallel transport and covariant derivatives are metric-indeependent formulation
dr Ga 1 Bo (x) (x) (x)
(58)
30

---

Published as a conference paper at ICLR 2026
The Riemann tensor Rabyó Iao R3, obeys the following identities:
Ras75
Ra,87
Ra375 = Rgayd
Ra375 Ryg0s
The Riemann tensor in a n-dimensional manifold has - independent components. Importantly, it satisfies two additional identities, called the Bianchi identities
1)/12
RaBy8 + Ra758 + Ra88y = 0 (Bianchi Identity I),
(59a)
VaR + V Vg Ryaad + V Ra,66 = 0 (Bianchi Identity II)
(59b)
Unlike the Christoffel symbols, which may be non-zero purely due to the choice of coordinates e.g., when imposing curvilinear coordinates such polar coordinates (r, 8) on the flat Cartesian plane- - the Riemann curvature tensor encapsulates the true geometric curvature of manifold. Since Christoffel symbols represent connection coefficients rather than tensorial objects, their non-vanishing components can give the false impression of intrinsic curvature, even on flat manifold. In contrast, coordinate-invarian statement: if the
the Riemann tensor is bona fide tensor and its vanishing Riemann tensor vanishes one coordinate system, it vanishes in all coordinate systems. Thus, provides a definitive criterion for distinguishing truly curved spaces from flat ones, independent of coordinate artifacts.
Geodesic deviation. that it encapsulates directional information about how geodesics path converge or diverge. Intuitively, it implies that Euclidean space Idd, parallel lines always remain parallel, but in the case of spherical geometry, say Sd-1 (constant positive curvature) the parallel lines converge at point, while for hyperbolic spaces HH the parallel lines continue diverging. This is captured by the geodesic deviation equation (sometimes referred to Jacobi equation) Isham (1999); Jost (2008); Poisson (2004), which shows how an infinitesimal neighborhood a given geodesics diverge or converge. Here, we shall give the equation with brief sketch.
important consequence of the existence of non-zero Riemann tensor
Theorem (Jacobi equation): Lei (T(TT) be a family of closely spaced geodesics indexed by smooth one-paramter family s and E the affine parameter. Let x* (x,T) be the coordinates of the geodesics (T),), then the tangent vector field is directional derivative expressed in these coordinates as X (s.(), Let the set of deviation vector fields S" Ox" (T.(),)) Then, the deviation vector fields that satisfy the acceleration equation are called (Jacobi fields) and read
89
5
Rd R8
X"X"X
(60)
D+2
D = X the directional covariant derivative, iee. the derivative of a vector field along where, Dr given direction manifold, while accounting for the manifold's 's curvature.
A.3.5.2 Contracted curvature tensors, scalars and invariants
Ricci tensor. From the rank (1,3) Riemann tensor, one can construct a traced (contracted) symmet- ric curvature tensor of rank (0, 2), called the Ricci tensor Rag d & 18
R7 Tr, (R()))) Rag
(61)
Mathematically, the Ricci tensor aggregates directional curvature along orthogonal planes. Thus, it can be considered as a curvature average of the Riemann tensor. is closely related to the concept of sectional curvature and reflects how volume deformations occurs evolve under geodesic flow.
Ricci scalar. The traced (contracted) Ricci tensor yields scalar field called the scalar curvature, also called the Ricci scalar. It defined as
R
Try(Rag) = R.
(62)
31

---

Published as a conference paper at ICLR 2026
Mathematically, the scalar curvature corresponds to the sum/average over all sectional curvatures, i.e., R(p) -az0 Sec(ea, eø)lp p vp E M. For point p. in an n-dimensional Riemannian manifold (M,g), it characterizes the volume of an e-radius ball in the manifold to the corresponding ball in Euclidean space, given by,
R Vol (p(() CM) -o) RR) 1 6(n+2)
+ +(02)
Weyl tensor. Another important tensor field of rank 0, 4) is the Weyl tensor, which is obtained as the "trace-free" part of the Riemann tensor. Physically, the Weyl tensor describes the tidal force experienced by body when moving along geodesics, and quantifies the shape distortion body experiences due to tidal forces (e.g., water tides caused by the gravitational pull of the moon). In an n-dimensional manifold it is defined as:
CaBy8 = Rapr5 (n Rasggy Ray988 + Rargas Readary
(63)
+ (n 1)(n -2 (gory93s Ias93r)
Mathematicallyy the Weyl tensor corresponds to the only tensor when looking at Ricci-flar manifolds, i... Rag 0. This become relevant for e.g., vaccum solutions, class of exact solutions where R., 0 of the Einstein equations in the absence of matter distribution.
components of the Riemann
Curvature invariants. Curvature invariants play central role in the analysis of spacetime ge- ometries general relativity. These scalar quantities are constructed from contractions of curvature tensors and manifestly invariant under general coordinate transformations. As such, they serve as powerful diagnostic tools for characterizing the local and global geometric and physical properties of spacetime, which includes the identification of "true" (genuine spacetime singularities) and "false" singularities (artifact of choice of coordinate charts).
Among the most prominent quadratic curvature invariants that is relevant to our simulations and features in our paper is the Kretschmann scalar, defined the full contraction of the Riemann curvature tensor with itself:
(x) Ra076 (a)))((x)))))))))) Ras)s(>)))) 9 (")) (x)g"
( (e)s((())Rass((1)))
(64)
The Kretschmann scalar provides a poordinate-indepeddeed measure of the magnitude of the cur- vature of spacetime and the singularity becomes blows-up, due infinite curvature. Examples are Kretschmann scalars X for blackholes, and Weyl scalars V for gravitational wave astrophysics. They capture the presence of intrinsic curvatures even when the Ricci tensor itself vanishes. Thus, the Kretschmann scalar encodes geometric information in frame-independent manner.
A.3.6 STRESS-ENERGY- MOMENTUU TENSOR
The stress-enes eaeeeeneergyyyooe momentum tensor (or simply called the energy morgyytummetumm tensor) is symmetric rank (2, 0) tensor
To
T ea $ eg
(65)
Physically, To.3 is a generalization of the stress tensor in continuum and fluid mechanics. It stores the information of distribution of matter fields, i.e., sources sinks a x 4 tensor, such as energy-density, energy-flux, momentum density, and momentum flux.
These matter fields satisfy the conservation laws, i.e., conservation of mass and energy via the four-dimensional continuity equation and corresponds to the divergence-free condition the energy- momentum tensor
T (x) 0.
(66)
The stress-ente
tensor features on the right hand side of the Einstein field equations (EFEs), and influences spacetime geometry by causing distortions it.
32

---

Published as a conference paper at ICLR 2026
A.3.7 EINSTEIN FIELD EQUATIONS
The EFES are set of second-order non-linear PDEs containing geometry on the left hand side and the source the right hand side. EFEs are obtained by combining all the differential geometric quantities from Eqs. (36, 61, 62, 65), together with the conservation laws for matter distribution of Eq. (66), resulting in
GoD STG To
(67)
Here, Gaß := Rass 2908 R-+ Aga, called the Einstein lensor and also satisfies the divergence-free condition V Gas = O, which, is a consequence of the IInd Bianchi identity Eq. (59b).
A.3.8 COORDINATE-INDDEEEDDENN OF GR
Fundamentally, GR posits deeper symmetry class: diffeomorphism covariance (Misner al., 2017). It asserts that the laws of physics are independent of any particular choice of coordinates or parametrization of the underlying smooth manifold. For example, the metric around the star, say sun, can be expressed in terms of the Schwarzschild metric (introduced in Section B.1). Here, the diffeomorphism acts as a gauge transformation (Tao, 2008) between two sets metries defined on Lorenizian manifold Mr in this case equivalence class of Lorentzian metrics Ri Riem(..A (..( describing the same spacetime geometry. This makes sure equations of motion, conservation laws, physical fields, etc. remain intact, hence, the term "covariance". In mathematical terms, Let, & Diff(.M(
be a smooth, invertible map between Al with smooth inverse, + : A. = M such that:
(9'9) := (v) 2.) (u))
(68)
Here, d. T.H - T.M, is the pushforward map defined on the tangent bundles. This means under diffeomorphisms the metric transforms via a pullback operation D'g = g'. T.e., go = Oz Oz"
dgo 088 (2)(z) (z) are gauge equivalent. Additionally, GR also admits changes local frames or bases (external symmetries) via the general linear group GL, IR), i.e., invertible linear transformations at each point E M. Thus, GR enjoys coordinate independence from two symmetry transformations, ie., (i) between any particular choice of coordinates or parameterization of the underlying smooth manifold M, and (ii) general linear group transformations that locally change frames of reference.
B EXACT SOLUTIONS OF THE EINSTEIN FIELD EQUATIONS
This Appendix contains detailed description of the exact solutions of EFEs corresponding class of metrics 9μv that are solutions of Eq. (67). While there exist several geometries that satisfy the EFFs, shall consider three prominent geometries: the Schwarzschild metric, the Kerr metrie, and gravitational waves. These solutions not only have a high theoretical and historical relevance, but are also of great interest in computational black hole astrophysics and gravitational wave and multi-messenger astronomy. From here on, we work in naturalized units by setting == G 1.
Our work predominantly uses the exact solutions for generating synthetic training data, which analytic expressions for (i) Schwarzschild, (ii) KerT, and (iii) linearized gravity metrics, on which we fit the NeFs.
B. SCHWARZSCHILD METRIC
The Schwarzschild metric is the simplest non-trivial solution to the EFEs. It describes the geometry around non-rotating spherical body, such as star or a black hole, constituting spherically symmetric vacuum solutions, i.e., Tr 0. A famous result of GR called the Birkhoff's 's theorem (Birkhoff & Langer, 1923) proves that any spherically symmetric vacuum solution corresponds to a static (non- rotating), iime-independent (stationary), and asymptotically flat metric (i.e., for co the metric converges to the flat Minkowski spacetime), and must essentially be equivalent to the Schwarzschild solution.
33

---

Published as a conference paper at ICLR 2026
B.1.1 COORDINATE SYSTEMS FOR SCHWARZSCHILD METRICS
Spherical polar coordinates. Schwarzschidd solution typically written in the convential spherical
polar coordinates rt, 0, Φ) where & E I,, E Rt, 0 E (0, T), and Φ E [0, 2T). The metric can be written either using the quadratic line element
ds²
2M
dt +
1
2M
-1
dr + r² (d0 sin sin?)
(69)
r
or in the equivalent matrix notation
2M
0
0
Sph
0
1
2M
0
(70)
0
22
0
0
sin
The true singularity of the Schwarzschild metric is at the origin and can be identified from the divergent Kretschmann scalar (Eq. (64)):
48M r-
X
00.
(71)
Although (fake) coordinate singularity exists at =. 2M. where the Kretschmann scalar is well defined. This special radius r, called the Schwarzschild radius. It demarcates the location of the event horizon of a non-rotating black hole and delineates a region from which no causal signal can escape to asymptotic infinity, meaning. it a point of no return for any body (including light) once it crosses this critical radius.
Cartesian Kerr-Schild coordinates. The Kerr-Schild (KS) form is a beneficial representation for finding exact solutions to the EFEs. These are perturbative corrections to a spacetime metric only upto the linear order (Kerr & Schild, 2009). The KS form enables the following simplifications to the nonlinear field equations: (i) It expresses the resultant metric lincarized perturbation to the background metric, and (ii) gets rid the coordinate singularities, which are mere artifacts of an unsuitable choice of coordinate systems. The corresponding line element expressed in the KS form reads
ds² (ā.3 + V (x) dda 'x
(72)
where fou is some background metric, la are the components of null vector & with respect to the background metric and (Ix) is scalar.
For spherically symmetric non non-rotating blackhole such as Schwarzschild, the Cartesian Kerr-Schild line element is obtained by setting Ia. 1Ja33 f.
1, z
pirrr and V 2M,
ds = - dx + dy² + d2
2M
dt + dx + y dy de
(73)
r
r
Unlike the spherical coordinate form in Eq. (69), 2M is not singular, hence removing the coordinate singularities. The metric tensor components read:
2Mx 22
2My 22
2M²
-1 2M2//r
2Mr
2M
2Mry ,3
2Mxz
,2
3
,3
g
(74)
2M
2Mry P3
y
2M
72
2M
2Mxz 3
2Myz
2.
1
,2
73
(75)
34

---

Published as a conference paper at ICLR 2026
Ingoing (advanced) BeeiingtoonFinkelttee coordinates The ingoing version of the Eddington- Finkelstein (EF) coordinates is obtained by replacing time t with an advanced time coordinate t (r)), where + M logl 2M Thus, dt in these transformed coordinates amounts to:
dt du dr. = du 2M
dr
The ingoing EF version of the Schwarzschild metric reads:
(1
2M
ds²
du + 2dl dr (d02 + sin'0 do²).
(76)
With the metric tensor being:
2M
0
1
0
(77)
0
0
0
0
sin
0
This metric is smooth (and non-degenerate), devoid coordinate singularities at the event horizon 2M, and can be continued down to the curvature singularity 0 (Carroll et al., 2004; Chandrasekhar, 1984; Frolov & Novikov, 1998).
B.2 KERR METRIC
The Kerr solution describes a massive gravitating body rotating with an angular momentum J. From the physics perspective, it is not symmetric under time-reversal symmetry, i.e., - -t, hence corresponds to a stationary but non-static solution (Teukolsky, 2015). Due to a finite angular momentum J, equivalently, rotation parameter M > 0, the Kerr metrie introduces asymmetry, and is oblate. Thus, the Kerr metric corresponds to an oblate spheroid geometry.
B.2.1 COORDINATE SYSTEMS FOR KERR METRIC
Linqqi-l coordinates. The Boyer-Lindquist (BL) coordinates are a special and convenient representation for the Kerr metric (Boyer & Lindquist, 1967; Visser, 2008; Teukolsky, 2015). The BL form (t, r, 0, Φ) is described by oblate spheroidal coordinates (Krasiński, 1978):
r2 22 + a² sind cosó
(78a)
y 22 + a² sin0 sinø
(78b)
2 = r coso
(78c)
Notice that the zenith angle + differs from the Schwarzschild case, while the azimuthal angle the same in both. As a - 0, the Kerr metric boils down to the non-rotating spherical case of the Schwarzschild metric.
d -0(1--1 2Mr
di2 4Marsin
dids+ E
(,2+02+ 2Mra2sin
sin 2 Odos
(79)
where the length scales are a = M (angular momentum per unit mass), E x2 + a cocs 0. and
= 2 2Mr + a?. The Kerr curvature singularity occurs at := ,2 + c2cos²0 = 0. implying
0 and 8 5. The metric tensor of Eq. (79) is:
35

---

Published as a conference paper at ICLR 2026
2Mr
Mar sin
Σ
Σ
Σ
0
0
A
(80)
0
Σ
0
2M or sin 0
(12 ( a²
2Ma sin"
0
sin 0
Σ
Σ
In the Boyer-Lindquis. form of the metric, there also exist coordinate singularities at A ,² 2Mr + a² 0. Thus, the roots of A = 0 are To M + VM2 a², which demarcate the outer and inner horizons.
It is easy to see the existence of curvature singularity at 0 on the equatorial plane corresponding to the zenith angle 2. Thus, unlike Schwarzschild, the singularity in Kerr geometry takes the form of ring, also known as ring singularity.
Cartesian Kerr-Schild coordinates.
The Cartesian KS form of the Kerr metric is obtained by mr setting in Eq. (72) L (1. مم py & and V +*+ in Eq. (72). The Kerr metric in Kerr coordinates are often used to write initial data for hydro simulations. The line-element in the Cartesian Kerr-Schild form reads (Teukolsky, 2015):
:
2mr ds² = - + dx dy + 7 + a² dt + r(ædx + ydy) alydr rdy) a² 42 + +-dz (81)
a² +2
Here, r(x,,,, v. is not coordinate, and given implicitly by solving the quadratic equation ty F4 i 1: The solution for the implicit function is given by the discriminant (Visser, 2008):
(x,,,,)) + d²
x2 22 + 2 2²)² 2a
The corresponding Cartesian coordinates are expressed as:
(r cospp sing) v + a² sind cos tan (a/r))
(r sinp + a cosp))sind+ = 122 + a² sind sin 4 + tan -1 (a/r))
cos
0 v /.4 translating 0 (equatorial-plang) and the ring occurring at x? y² a² In contrast, the KS representation is devoid of coordinate singularities, making it suitable to work in numerics, especially around the event-horizons.
In the BL coordinates the ring singularity for Kerr exists
Ingoing Eddington-Finkelstenn Finkelstein coordinates. In the original formulation, the Kerr metric is written in the advanced time coordinates/ingoing EF coordinates v. The line element in this representation reads (Teukolsky, 2015):
as2 1=
2Mr
(dv + sin²0 dō)²
(82)
22 + cocos²0
+ 2(dv + sin2 (dr + sin sin)
+ (r² + aaccos"0)
sin 0d0²)
where the ingoing EF coordinates are related to the Boyer-L.ndquiis coordinates Eq. (78) by the following transformation:
(f2 (12+09) dr
e ==t+/
t+
A
3 Φ + / dr
36

---

Published as a conference paper at ICLR 2026
where, A = = 2M + a², and the metric tensor components corresponding to the line element is given by
2Mr
2Mar 2Ma
0
Σ
0
a sin 2
FF Uu
(83)
0
Σ
0
2M sin" 0
asin 2 0
(22+0²
aa²
2M
sin 0
Σ
The coordinate (fake) singularities (x < oo) of the Kerr metric is given by the zeros of A = p2 2Mr + a² Solving for the zeros, one finds
T= M M a2
where, 74+ is the outer event horizon, while r- demarcates the inner event horizon.
Apart from that, rotating metrics also possess highly interesting region known the ergosphere, which fundamentally captures the non-Euclidean and non-inertial nature of general relativistic effects induced by rotation. This domain situated outside the outer to the trame-dragginn remain static within the ergosphere and is compelled to co-rotate with the black hole depending on the value of a. The location of the ergosphere is given by
horizon r4, and is created owing (einr--hiiring) effect. Consequently, no physical observer (test body) can
M + M c2cos²0
where, 1 is the outer ergosphere, while ergo demarcates the inner ergosphere. In Figure 8, the following regions of the Kerr metric are demarcated:
0.5
Z(GM/c2)
0.0
-0.5
-1.0
-1
1
X(GM/c2)
Figure 8: Kerr metric 2D slice in the x-z plane (y = 0) for a spin parameter = 0.99. The following regions plotted are: i) inner ergosphere sgo, red region, ii) inner event-horizon r.: green region, iii) outer event-horizee r+: blue region and iv) outer Ergosphere ergo, purple region.
B.3 GRAVITATIONAL WAVES
Linearized gravity models the metric as tiny fluctuations or perturbations has of the flat background metric 7n
9a nag + hag + ((has)²
where |/a,al < To describe gravitational wave propagation, it is often convenient to reduce the linearized field equations into simplified form via two gauge fixing conditions, namely, a) harmonic
37

---

Published as a conference paper at ICLR 2026
gauge, and b) transverse-tracelen (TT) gauge. Thus, the Einstein field equations for gravitational waves assume a succinct wave equation type form:
(c) Dh -66TTa
(84)
Here, Ra5 = - 2
(1() nas and D = 00 0g is the d' Alembert operator (wave operator). It can be shown that the PDEs in Eq. (84) produce gravitational wave solutions.
The transverse-traceles. (TT) perturbation (we drop the superscript (e) for the sake of ease) satisfies the following conditions:
Transverse: 0gh-Taa3 a3 0, i.e., wave propagates perpendieular perturbation direction,
Traceless: a 0,
Purely spatial: hom = 0, i.e., no time components.
Thus, gravitational wave propagating in the z-direction with frequency w is given in the TT gauge as:
0
0 h
hx
TT 0.
cos (w(( - 2))
(85)
0 hx
h+
0
0 0 0/
h4 and hx are the amplitudes of the "4" (plus) polarization and "x" (cross) polarization.
The complete metric tensor in the linearized gravity setting is given by:
0
+ h4 cos (wil - ==))
hx cos (w(t - =))
Dop Va88 + Naa
(86)
hx cos (wil - z)) 1 14 cos (wot ( - 2))
The corresponding line-clement the linearized gravity setting reads: ds² dd + [1 h4 cos(w(t =))] do² [1 h- cos(w - -t)] dy²
+ 2/x cos(w(t - z)) dx dy
(87)
Spin-weighted spherical harmonics (SWSH) metric representation. We start from the decompo- sition of the complex gravitational wave strain with the spherical harmonic basis-set expansion (New- man & Penrose, 1962). With the expansion in mode weights 1 (1, r), one can ignore (remove) the angular dependence:
M
A(1, O, 0) = (1,,,,0 ih (1, 0,
mm Yem (),
(88)
[=|+|
|m|Se
where, htt, 0, = (t,r,,, ih rt, Φ) is the complex strain. Thus, for each orbital and azimuthal indices m), one can extract the mode fm((() at a fixed radius r, one uses the orthogonality of the spig-weighted spherical harmonics (SWSHs) elements:
M fm kT
(t.r,, Ven Φ) dΩ
(89)
where, dΩ sin 9 d0 do is the spherical volume element and -2 Yem denotes the complex conjugate of the s = 2 - n = -2 where (n = 4 for V2 spin-weighted spherical harmonics. One typically carries out the integral over the 2-sphere 82 numerically finite angular grid.
38

---

Published as a conference paper at ICLR 2026
The general formula for SWSHs is:
28 + 1 (e + m)!(t m)!
1/2
0
Ytm (0, x) (-1)/tm-s
sin
28
cmm
(90)
4m
(e+ s)((t s)!
2
8-2
e ((-8) r
t+S r r m
1- 1y
2r+s-m
cot
where the parameters m are the familiar Laplace spherical harmonics (orbital-angular momentum and azimuthal indices), while is the additional spin-weight introduced by some underlying gauge
group such as (1). We especially plot for the integers -2, = m = since they are relevant for GWs and are depicted in Figure 9.
Figure Spin weighted spherical harmonics for s 2, and for |m| S The dominant contributions for the Weyl scalar 4 and the associated metric coefficients in the spherical harmonic basis (1) are shown.
B.4 MINKOWKSI METRIC
B.4.1 COORDINATE SYSTEMS FOR MINKOWSKI METRIC
The flat Minkowski metric, which a spacetime that has no curvature (M 0) can be expressed in other coordinate systems as well.
Spherical polar coordinates. In spherical coordinates (t, 1), sp), the Minkowski metric de- scribed by the quadratic line element,
ds² cd + dr² + 2262 + sin sin do²)
(91)
here, E R+ E (0, R)) and, øb [0, 2T)) the usual spherical polar coordinates. Thus, the metric tensor describing the Schwarzschild solution reads:
0
1
0
Sph Faz
0
0
2
(92)
0
0
0
2
Boyer-Lindquist coordinates. Setting M - in the Joyer-Lindquis form of the Kerr metric Eq. (79), the corresponding line element reduces to unfamiliar "oblate-spheroida" represtation:
22 2
ds² ddt
dr + (r2 co do² + (r² a²) s d (93)
22 + a²
and the components of the
0
0
BL. Muu
0
2 a coos
0
(94)
0
0
2
+ a cos2
0
0
0
(r2
+ a²) sin 20
The usual Cartesian coordinates can be related to the oblate spheroid ones via:
1 + a2 sin & cos Φ
y = vr2 + a² sin & sin Φ
cos 1
39

---

Published as a conference paper at ICLR 2026
Eddington-Finkelstein coordinates. The Minkowski metric be written in the ingoing (ad- vanced) Eddington-Finkeet Finkelstein form in two different cases, namely for the non-rotating and the rotating
i) non-rotating, a - 0 (Schwarzschild))
ds² 2 + 2 do dr + r2 sin² dø²).
and the metric tensor reads:
1
0
0
(95)
0
0
12
0
0
0
0
72
ii) rotating: (Kerr):
ds² = (dv + sin sin2 odo)² + 2(dv + sin ado) (dr + asin adin)
+ (r² + a² cos? ( sin do?),
(96)
and the metric tensor reads:
EF JAY
1
asin D
(97)
0
0
E
0
0
asin
9
(122
+ a²) sin 2
B.5 TRAINING NON-TRIVIAL. METRIC FIELDS (DISTORTIONS)
Distortion part of Schwarzschild geometry in spherical coordinates:
Spherical coordinates: obtained by subtracting Eq. (92) from Eq. (70):
Ts
0
0
9pm
(r
T
(98)
0
Kerr-Schild coordinates obtained by subtracting nu
diag( -1, +1,+1) +1,+1,+1,++1 from
Eq. (74):
2M
2Mx
2My ,2
2Mz
r
,2
y22
2Mx
2Mr²
2Mry pr3
2Mxz 733
72
p3
fjv
nl
(99)
2M
2Mzy p3
2My² 3
2Myz 3
22
2M
2Mæz n3
2Myz p3
2M
,2
3
Ingoing Eddington-Finkelen coordinates obtained by subtracting Eq. (95) from Eq. (77):
r
0
0
9
9pm
FF
Up
EF
(100)
0
0
0
0
40

---

Published as a conference paper at ICLR 2026
Distortion part of Kerr geometry:
-BoyerLLindquis coordinates: obtained by subtracting Eq. (94) from Eq. (80):
2Mr
or sin
0
0
Σ
2ME
0
9pz
BL. f/
rA
(101)
aM sin?
2M sin 0
0
Σ
Σ
Kerr-Schild coordinates: obtained by subtracting
diag (1,,+1,+1,+1) +1) from
Eq. (81):
Eddington-Finkeen Finkelstein coordinates: obtained by subtracting Eq. (97) from Eq. (83):
2Mr
2Mar ar sin?
0
0
Σ
0
0
(102)
0
0
0
2M ar sin
2M sin
0
C FINITE-DIFFERENCE METHOD (FDM) FOR TENSOR DIFFERENTIATION
The main concept in this appendix section details numerical differentiation methods for tensor-valued quantities, focusing on the practical use higher-order difee-diffeeence schemes (in particular, sixth-order stencils). We outline the treatment of discretization errors and the use of neighboring grid collocation points as part of a numerical tensor calculus toolbox.
To compare the performance against automatic differentiation tensor fields defined the four dimensional spacetime, throughout the paper we opt for the highly accurate sixth-order order forward difference (n = 6 accuracy). This scheme queries six neighboring points per evaluation and the general formula of the differential operators are given by:
49
6
15
20
[0x ~
20h
fx)
T.
fx h
2h
f(x + 2h) +
3h'
f( + 3h)
15
6
f(x + 4h) + 5h
f(x 5h)
6h f(x 6h) + O(h))
(103)
4h
i-th index h
Here, = (x1, E Rdd and the h = he, s.t. = (0,
0), depending with respect
to the variable the partial derivative is performed over. Thus, this is accurate upto O(h)), and the truncation error occurs at seventh-order.
In general, for an n-th order niffe-difference approximation, the stencil is constructed by querying neighboring collocation points on the voxel grid. The resulting truncation error on function evaluated on the grid (gridfunctions) scales as follows (Ruchlin et al., 2018):
G = f1
Thus, higher order stencils enable larger step size h choices since the error scales exponential to the stencil order, i.e, & x h., This results in not only better accuracy, but also lesser memory consumption due lower grid resolution.
41

---

Published as a conference paper at ICLR 2026
Finite-difference method bottlenecks in NR
Higher-order finite-difference stencils require a collection of padded grid points exterior to the cube as boundary handling. These are often called ghost cells (zones). For e.g., if ni ghost cells are required for the n-th order forward difference stencil, the grid size increases of a 3D spatial voxel grid (Ruchlin et al., 2018), N₂ Ny N - (N, + +gg) (Ny ngg) X [N= + nz).
Sensitive to numerical noise especially for tensor- tensor-valued functions defined on multidimen- sional voxel grids.
D SUCCINCT INTRODUCTION TO GENERAL RELATIVITY, EQUATIONS OF MOTION AND EXACT SOLUTIONS
Derivative operators. symbols
The metric and its partial derivatives can be used to construct the Christoffel
T (x) : (Qugov(x)
Ougou + (00
- Dogur (x))
The Christoffel symbols denote how the metric varies across spacetime and define a parallel transport machinery to translate tensor fields around the manifold. With these, it is possible to construct two pivotal modified tensor differentiation operators, namely: (i) The covariant derivative (also called the Levi- Civita connection), which can be seen as "calibration" of the partial derivative operator for parallel transportation in the curvilinear setting:
VoToi
Talman
Er T no
-r T
BI....
Oz! PI...e
91...A
PI........."
2m1
jm1
(ii) The Lie derivative, which generalizes the notion of a directional derivative that is connection independent (cf. Appendix A.3.3.2). The Lie derivative captures infinitesimal dragging of the tensor field along the flow generated by the vector field E:
() Orr SI...9, = To mir Br...e,
ETamp Br... 0p5 + +
T. B.... 0g, 0z
Differential geometric objects. Using the modified derivatives, we can construct a hierarchy of higher-rank differential geometric quantities, such the Riemannian curvature tensor 13, or the Ricci tensor Rag, via a series of derivatives 0 covariant derivatives V = Ô + r. and tensor index contractions 8; Vr+1 - Vrr (typically, Try).
V
C Casy8
9a
r ag
R c,9y
RaB
R
Rias3jo =
Lx
V
a Ra
Figure 10: Differential geometry workflow in general relativity (only lhs of the EFEs Eq. (2)): We describe each quantity starting left: The metric tensor 9a3 defines the spacetime geometry. Its partial derivatives i yield the Christoffel symbols 1 ag. which describe the notion of parallel transport and defines a covariant derivative operation V. θ, I... The connection also defines the Lie derivative C. along vector fields v. The connection's 's derivatives give the Riemann curvature tensor Recy. which encodes tidal forces. The contraction operator & Tr, contracts with the metric, producing the trace part, i.e., the Ricci tensor R.,9. Its subsequent contraction yields the Ricci scalar R. The Riemann tensor further splits into the Weyl tensor C. (trace-free curvature part) and satisfies the algebraic and differential Bianchi identities R|a,>)0 0 and VioRasha O. Together, these geometric objects form the backbone of general relativity, ultimately entering the Einstein field equations through the Einstein tensor Gas-
42

---

Published as a conference paper at ICLR 2026
Conservation laws. It follows from the contracted Bianchi identities, i.e., cyclic sum of Riemann curvature tensor covariant derivatives (II Bianchi identity see Eq. (59b)) vanishes identically:
V RaRoy + V VaRaaōo + V 0
This is geometric identity that holds for any (torsion-free ffre))) connection compatible with the metric. The identity above consequently leads to the covariant derivative of the stress-energy tensor vanishing, that is. VaTa3 := 0 (see Eq. (66)), which corresponds to the general relativity. If required, conservation laws typically feature soft constraints in PDEs, and relevant especially when matter distribution/fields are considered.
conservation in
Equations of motion. The geodesic equation is a central second-order ODE that describes the motion of objects following geodesic paths in the curved spacetime background
d²
dx dx² pa
+
dr
dr dr
Here, T E. IR is some affine parameter (often the proper time; not be confused with the coordinate time t). The geodesic equation is the general relativistic analogue of Newton's 's second law and generalizes the concept of a "straight line" curved spacetime by describing the trajectories of bodies under the influence of gravitational field. Related, the geodesic deviation equation describes how nearby geodesics diverge or converge due to the intrinsic curvature of manifold, quantified by the separation vector (T(T):
S²
Rd a
X
Dr2
where, X is vector field and Dr X denotes the directional covariant derivative (see Definition 2). Thus, it encodes information about the tidal forces of gravitation.
D.0.1 ANALYTICAL (EXACT) SOLUTIONS
Exact solutions of the EFEs are metric tensor fields Jap that satisfy Eq. (2). Many exact solutions are known, which can be classified into exterior (vacuum) solutions and interior solutions (Misner et al., 2017). While there exist several geometries that satisfy the EFEs, we shall focus on three prominent vacuum solutions: the Schwarzschild metric, the Kerr metrie, and gravitational waves. These geometries not only have a high theoretical and historical relevance, but are also of great interest in computational black hole astrophysics and gravitational wave and multi-messen messenger astronomy. Appendix B discusses these solutions in more detail, including other prominent coordinates, as well real and fake (coordinate) singularities. From here on, work in naturalized units by setting = 1.
Schwarzschild metric It is the simplest non-trivial solution to the EFEs and describes static spherically symmetric geometry around spherical non-rotating gravitating bodies, such as stars or black-holes. Although simple, the Schwarzschild metric predicts many phenomena beyond Newtonian gravity, most notably the precession of elliptical orbits and the bending of light rays. Both of these predictions have been experimentally verified the Solar system, using the motion of Mercury perihelion and in the Eddington experiment during the 1919 Solar eclipse, respectively. The metric typically written in spherical polar coordinates (t, 0, Φ))) where E R, r E IRt, E (0,z). x), and (0,
2 = -- 2M
2t + (1-
2M (122M)
2M
-1
dr + 22002 + sin?
Kerr metric generalizes the Schwarzschild solution to rotating bodies with the angular momentum or, equivalently, the rotation parameter a J. The solution forms rotating, stationary (but not static) geometry, which is oblate around the rotation axis that breaks spherical symmetry. This geometry again permits new phenomena, notably the geodetic effect and frame dragging, both of which have been experimentally verified in the Earth's 's' orbit by the Gravity Probe B.
The metric can be described in the corresponding oblate spheroidal coordinates also known as the Royer-Lidqute (BI.) coordinates (1, r, 0, Eq. (78) (Boyer & Lindqiist, 1967):
43

---

solving vaccum Published as a conference paper at ICLR 2026
ds² = -11 2Mr
(1-2M) dt² 4Marsin
Σ dtdo+ dr + (1+ 2Mra?sin's sin??
E
Linearized gravity models the metric as tiny fluctuations or perturbations hags Ihal < flat background metric na33
of the
gap R Mo8 + has + O(hag)²
Famously, this model can describe GW propagation using a periodic dieedependent perturbation, which has served as the theoretical basis for the Nobel-prize winning detection of GWs generated by binary black hole mergers (Abbott et al., 2016c). As detailed in Appendix B.3, the choice of
a certain gauge essentially transforms the vacuum EFEs into the wave equation DMa
where (c)
C = 1) 0g0g is the d' d'Alembert or wave operator and a = h h 1 m a family of GW solutions: we will use the plane wave propagating in the z-direction with the angular frequency w expressed in the Cartesian coordinates:
Na,3. This equation admits
0
0
0
h4
hx
TT 0B
cos (w(t 2))
0
hx
h+
0
0
0
0
Here, h4 and hx are the amplitudes of the "4" (plus) polarization and "x" (cross) polarization
E RELATED WORK CONTINUED
Compression techniques scientific computing. classical compression strategies have been a versatile tool in reducing data sizes of large-scale numerical simulation data, which constitute storage memory-intensive domain. Simulation runs can range between petabytes to exabytes of data, thus requiring compression strategies be integrated within the simulation pipeline, for efficient data volume reduction. These range from lossless and lossy compression techiques (Lindstrom & Isenburg, 2006; Di & Cappello, 2016; ISO Central Secretary, 2024) for multidimensional weather modeling tasks (Huang et al., 2016; 2025), discrete wavelet tranform (bae)- compression (Rho et al., 2023) or multidimensional data via tensor decomposition (Ballester- (Ballester-Rpooll al., 2020).
Recently increasing works have leveraged implicit neural representations for lossy neural com- pression (Dupont et al., 2021) by embedding nigh-dimensional (explicit grids). time-dependent physical simulations into compact, differentiable network weights. Such representations achieve compression factors of several orders of magnitude while retaining physical accuracy and offering efficient gradient access for downstream analysis control. Applications include multidimensional weather and climate modeling (Huang Hoefler, 2023) or even hybrid compression techniques using turbulent plasma plasma-based based simulations (Galletti et al., 2025).
MI applied to gravitational physics. Gravitational wave modeling and numerical relativity problems have been tackled by state-of-the-are deep learning methods. These include DINGO rapid gravitational wave parameter estimation toolkit using NNs as surrogates for Bayesian posterior distributions (Dax al., 2021). They show orders of magnitude reduction in inference time, bringing it down from O(day) to 20s. Similar lines of work include DINGO-BNS that performs real-tin val-time inference for binary neutron star (BNS) mergers and applicable for multi-messenger astronomy (Dax et al., 2025). On the other hand, physics informed neural networks (PINNs) (Raissi al., 2019), have found applications in general relativistic phenomena and NR, such as solving the Tekoulsky equation inorder to compute the first quasinormal modes of the blackholes. for e.g.. Kerr geometry (Luna et al., 2023; Cornell et al., 2022). Some other works explore physics informed neural operators (PINO) (Rosofsky & Huerta, 2023) for magnetohydrodynnmic Einstein equations (Hirst et al., 2025)
Physics-informed neural networks, neural operators and neural fields. Physics-Informt Neural Networks (PINNs) augment neural network training with physical constraints derived from governing differential equations (Karniadakis et al., 2021). This is typically achieved by adding loss terms that
44

---

Published as a conference paper at ICLR 2026
penalize violations of the residuals of the underlying PDEs ODEs (Raissi et al., 2019), together with constraints from boundary conditions and conservation laws. PINNs are primarily data-free approaches and have been especially effective for solving forward and inverse problems. While PINNs share with neural fields the use of coordinate-based neural networks, their purpose is fundamentally different: PINNs solve physical equation by enforcing its residual during optimization, whereas neural fields represent given physical field directly from data without solving the governing equations. In this sense PINNS can be considered as special case of neural fields (Xie et al., 2021), that satisfy the govering physical equations at each step.
Neural operators, contrast, aim learn mappings between unfnite-dimensional function (Banach) spaces (Kovachki et al., 2023), providing data-driven approach to approximate solution operators for entire families of PDEs. They prioritize generalization across different input functions, geometries, or forcing conditions, typically requiring large training sets covering many PDE instances.
Implicit neural representations (Müller et al., 2022) occupy a distinct position relative to both paradigms. Rather than enforcing PDE residual (as in PINNs) learning operator over families of solutions (as neural operators), neural fields provide a compact, continuous, and fully differentiable representation of single high-dimensional physical field. This makes them especially well suited for encoding scientific data with high fidelity, enabling continuous spatial (and temporal) query access, implicit compression, and differentiable downstream analysis. In this sense, neural fields are not a competing method for PDE solution or operator learning, but a complementary representation framework for capturing and reconstructing complex physical domains. nntructing the dynamics with high fidelity, setting them apart from PINNs and neural operators in general.
F EXPERIMENTAL DETAILS
This appendix provides detailed experimental specifics, including: (i) AD as a superior differentiation framework for tensor fields compared to higher-order
methods; (ii) Effectivity of higher derivative losses for retrieving high-precision dynamics and higher-order curvature tensors & invariants; (iii) setup data preparation (iv) gradient alignment aspects relevant SOAP optimizers: (v) Component-wii wise tomography of compression eтror on the implicit metric and its derived higher- rank differential geometric quantities; (vi) training across varied coordinate systems illustrate the pardinate-choiie flexibility of NeFs; (vii) hyperparameter configurations employed; and (vii) the hardware and software environments used for these experiments.
F.1 EVALUATION CRITERIA
We flatten the ground truth tensor at point p E M with its components indexed by k be denoted by (p) E R², with < < n and the corresponding EinFiel Felds parametrized tensors are denoted by (p). The dimensionality depends on the tensor under consideration. For instance, for a symmetric metric tensor n 10, corresponding to its independent components, while for the Riemann curvature tensor, = 256 when considering all components explicitly, or n 20 when accounting only for the independent components under the symmetries inherent to the tensor, respectively.
We evaluated these quantities over a set of m s 125, 000 validation collocation points D {p:) (pi)isinn and use standard error criteria in discretized form, which includes double sums: one over the total number of tensor components {fr]isken while the other for the total number collocation points {p.)isimm:
Mean-absolute error (MAF) =
(fa(p)) - fx(P))
(104a)
mm
Σn (a( fa (p:)) mM
Relative &g error (Rel. E2)
(104b)
Ei= (p:))
These are applied to the metric tenors and their derived quantities, illustrated in Figures 10 and 2. Recall that the tensor components cooddinate-eependent (and even more so, the metric Jacobian,
45

---

Published as a conference paper at ICLR 2026
metric Hessian, and Christoffel symbols not even tensors), and, hence, these enror lack immediate physical meaning. This is improved with the consideration of scalar quantities such as the Ricci scalar, Kretschmann invariants, and Weyl scalars, which by definition are coordinate- independent quantities.
F.2 DATA GENERATION
Our use cases are exact analytic solutions the EFEs, i.e., the set of metrics 9a3 that satisfy the Eq. (2). These solutions describe the exterior (vacuum) solutions around massive gravitating objects. For our main Section D.O.I, each having different features and spatio-temporal symmetries:
experiments, we a NeF against the analytic solutions introduced in
Schwarzschild metric in spherical coordinates - Eq. (98),
Kerr metric in Boyer-Lindquist and Kerr-Schild coordinates Eqs. (101, 81),
gravitational waves metric (TT gauge) Cartesian coordinates - Eq. (110).
For each, compute the distortion after subtracting the flat background metric in that particular coordinate chart. Detailed information on data specifications provided in Table 7.
Additionally, train each geometry in different coordinate systems investigate how the choice of coordinates impacts NeFs (recall: the physical laws do not depend on the coordinate system).
Table Training data generation specifications: spacetime metric, coordinate system, domain extent, grid resolution, and physical parameters.
Metric
Coordinates
Domain
Resolution
Parameters
Schwarzschild
Spherical (t,r,0,))
0 [2.5,150]
= 1
128
0 (0, T) Φ E [0, 2T)
128
128
Kerr
Boyer- Lindquist (f, Φ)
0 [3,14]
M 1 [0.628,0.95]
128
0 E (0, m) Φ E [0, 2m) 0 x E [-3, 3] y E [] 2 E [0.1, 3]
128
128
Kerr-Schild (t, v. y,
M 1
128
a = 0.7
128
128
Linearized gravity
Cartesian (t, v, y,
I E [0, 10]
140
w
I [0, 10]
10
10
E [0, 10]
10
2 [0, 10]
140
46

---

Published as a conference paper at ICLR 2026
F.3 COMPARING AD FD BASED METHODS.
We quantify the performance of automatieediiferenttattoo operations on the ground truth metric against the 6-th order finite difference stencils. We test it against the Kretschmann scalar X By R which is prone to errors, especially due to floating point errors accumulated in the Riemann curvature tensor:
0
10-14
50
10-15
100
10-16
0
100
X
Figure 11: Absolute error aaaytyi MADl profile plotted for = 0.3 between the analytic
Kretschmann scalar and the ground truth Kretschmann scalar obtained via implemented on the ground truth (analytic) metric.
47

---

Published as a conference paper at ICLR 2026
F.4 HIGHER TENSOR DERIVATIVE LOSSES SOBOLEV TRAINING
Sobolev training (Czarnecki et al., 2017) refers to a class of learning paradigms where NNs are trained not only match target function values but also additionally its derivatives. Formally, given target function f X R, and NN approximation fe. Sobolev training minimizes a joint loss involving the functional and its derivatives:
N
Lsab = Ez )ollf(z) fe())) + ΣM f(0) - DO0 fecor) j=1
(105)
where D denotes the 7th derivative operator, which in our case could be the partial derivatives @ or covariant derivatives V and A are weighting coefficients. This loss promotes alignment not only in function space but also in the Sobolev space W which encodes both value and derivative information. Sobolev training enhances generalization, stability, and accuracy of NeF derivatives (Chetan et al., 2024).
Algorithm I: "'Einniell ield training scheme
1: Input: Training dataset {(x', g(x),, a g(x"), )(2)((')))
g(x'() Je=1 number epochs Nepochs:
learning rate 7/, optimizer O, Sobolev order N (0, 1.2)
2: Initialize neural field parameters 0 device D (c.g., GPU) in single (FLOAT32) precision
3: for epoch to N Nepochs do
4: for each mini-batch (Featch- Poatch: 0.21) Pbatch: 2.21 Doatch) in dataset do
5: Move (Zbatch Pbatch to device D
6: Abatch EinFieill. ds (Tbatch; 0)
7: loss + (bbaach- Ibauch)
8:
if N 2 then [Jacobian supervision)
Compute 0=2t Joarch through AD loss +- loss X L(a) Iranch 0% ghatch) end if
9:
10:
11:
12:
if N 2 2 then (Hessian supervision)
Compute a!²²) Ooarch through AD
13:
14:
loss +- loss λ2 L(D.2) (2)
Ibatch 0=22) gbatch)
15:
end
16:
Compute gradients: Og Ög loss
17:
Update parameters: + ( 0g. 17)
18:
Optionally: synchronize gradients across devices if using distributed training
19:
end for
20:
Optionally: evaluate validation set, log MAE and memory usage for monitoring
21:
Optionally: checkpoint for fault tolerance and reprodueibility
end for
23:
return optimized parameters 0
The expected losses c(ol") g. a04 g) put-forth in Algorithm is short-hand notation for Eg ao4)) (a)(x) a.2) (x)(x))
48

---

Published as a conference paper at ICLR 2026
F.5 GRADIENT ALIGNMENT
Competing tasks a well-known problem in multi-objective learning (Yu et al. (2020), Liu et al. (2021), Shi et al. (2023)), the gradients of the loss functions pull the weights in different directions. In Scientific Machine Learning (SciML), a lot of work emerged analyzing and mitigating gradient conflicts in the context of PINNS (Wang et al. (2025a), Lui et al. (2025), Hwang & Lim (2025)). Although Sobolev training differs from PINNs, particularly from supervision perspective, it exhibits the same problem where some loss terms dominate others. In PINNs, the training is highly dependent on first satisfying the initial/boundary conditions, which provide uniqueness to the solution. The different levels complexity between these and the residual loss create different optimization priori- ties, but both losses are equally important. Similarly, Sobolev training faces analogous challenges with competing loss components. The Jacobian data serve to constrain the model's derivatives, while the target function outputs determine the integration "constant", both components being equally valuable. However, the sources complexity differ between these approaches. In PINNs, the primary challenge stems from determining solution through unsupervised learning on PDE losses, whereas in our Sobolev training specifically, the complexity arises from managing optimization stability in high-dimensional spaces: 16 16-dimensional output space, 64-dimensional Jacobian, and 256-dimensional Hessian. Moreover, this complexity accompanied by the challenge of handling gradient imbalances. Depending on the point spacetime, the metrie its derivatives dominate the loss. Generally speaking, an analogy is to think 9ur x _,0p9 p and popgp & A.. making clear how gradient magnitudes differ depending on the radius.
Mitigating gradient conflicts does not necessarily result in better accuracy, but explains a possible reason why the loss does not improve further, the optimization being stuck in the local minimum of one of the objectives. The intra-step gradient alignment scores presented in Wang et al. (2025a) demonstrate SOAP as a far superior alternative compared to other well-established optimizers, or at least for the experiments considered in that study. To provide direct comparison using the same methodology as in Wang et al. (2025a), we evaluated both ADAM and SOAP on the Cartesian KS representation of the Kerr metric, chosen the most complex metric investigated in this work. The Sobolev training contains only two objectives: metric and Jacobian supervision. For our experimental setup, employed an MLP architecture with hidden layers, 190 hidden units per layer, and SILU activation function to compare gradient conflicts between the two optimizers. The training utilized cosine decay learning rate schedule, starting from initial learning rate of 1E-2 and decaying to final rate of 1e-8 over 200 epochs. For weighting the losses, gradient normalization was used with and without an exponential moving average. As shown in Figure 15, even though Adam is providing twice as much gradient alignment score in almost all epochs, SOAP's second-orde and preconditioning capabilities allow for 100x training loss improvement.
1.0
0.8
10 3
Algnment
0.6
10-5
0.4
10-7
0.2
10-9
0
50
100
150
200
0
50
100 Epochs
150
200
Epochs
Adam
SOAP
Figure 12: Average gradient alignment per epoch (left) and MSE loss during training for Adam and Soap optimizers (right). The shaded light color in the alignment plot represents a minimum and maximum deviation compared to using an exponential moving average or not for the weights multiplying the gradients.
49

---

Published as a conference paper at ICLR 2026
F.6 RECONSTRUCTING GENERAL RELATIVISTIC DYNAMICS AND CURVATURE SCALARS
We use synthetic data generated from analytical solutions to validate and characterize Einnieldee We primarily focus two interesting aspects of general relativity: (i) general relativistic dynamies, particularly geodesic motion around massive gravitating objects (see Section D) and (ii) global curvature structures encoded in tensorial invariants.
F. FLOATING POINT PRECISION CAVEAT.
NR simulations are inherently high-precision endeavors, with the accurate modeling of complex gravitational phenomena critically reliant on high-fidelity numerical computations. In contrast to traditional machine learning domains, such as large language models (LLMs), where reduced- precision arithmetic (FLOAT16 6 or BFLOATIC 6) yields strong results in both training accuracy and memory efficiency (Dean et al., 2012), this paradigm does not extend to NR workflows, where floating-point precision is a dominant factor influencing the fidelity of the results.
EinFiele Medeeddbsaed geodesics. To compute geodesic motion, we numerically integrate the trajecto- ries using a fifth-order explicit singly diagonally implicit Runge- Kutta (ESDIRK) solver. Specifically, we evolve Eq. (56) with respect to the affine parameter (proper time)", generating ground-truth geodesics from the analytic Christoffel symbols for the Schwarzschild and Kerr spacetimes. These are then compared against the rollouts obtained using the EinF Leld-reconstructee Christoffel symbols.
To accurately retrace geodesic orbits, it is essential to incorporate Jacobian supervision within the Sobolev training framework. In contrast, additional Hessian supervision results in only marginal improvements for geodesic simulations and not required in practice. Following Section F.6. 1, all geodesic solvers executed in double precision ensure numerical stability and high-fidelity trajectory reconstruction.
Following the strategy presented in Section 3.2, we leverage EinFields' ability to yield high- precision derivatives the spacetime metric, which includes Christoffel symbols, Riemann curvature tensors, and scalar invariants. In this section, demonstrate how accurately model particle trajectories derived from geodesic integration with our implicit parameterizations.
F.6.2 SCHWARZSCHILD METRIC GEODESIC SETUP
Geodesics in the Schwarzschild spacetime are of fundamental interest, as they underlie phenomena such as gravitational lensing and the perihelion precession of Mercury, as well as the motion of planets in the solar system more generally.
The initial conditions of the trajectories chosen in the experiments are fully specified by the initial position
(t,r, ()) = 0) = (0, agT म 7/2
(106)
and the initial velo-iel
0,0, do
(v',"", v", 2 (t (V- ((1-verro(C 1/))
(107)
Where to = bo 11(ro r.) r.) and ag, by E R can be chosen freely to select the desired orbit in the 0 2 plane. The geodesics in Figure 4 demonstrate good qualitative agreement several orbits. The error quantified and discussed further Section F.7.2.
F.6.3 SCHWARZSCHILD BLACKOOLE RENDER
Being able to compute geodesics is sufficient to perform rendering. We use the Schwarzschild EinField ds metric render black-hole on celestial background. This requires propagating pawpaireerpppppaaaaa geodesics from the camera observer via the spacetime terminating at the distant background. The re- sulting ray-traced image as shown in the main text provides visual evidence for the global consistency and quality of the metric and the derived Christoffel symbols.
Not to be confused with the coordinate time
50

---

Published as a conference paper at ICLR 2026
F.6.4 KERR METRIC GEODESIC SETUP
Geodesies in a Kerr spacetime around rolating body (see details Appendix B.2) play a central role in several key astrophysical observations and experimental tests of GR. Notably, photon geodesics determine the black hole shadow images captured by the Event Horizon Telescope (Fuerst, S. V. Wu, K., 2004), and trame-dragging Kerr geometry. These have been measured experimentally by the Gravity Probe mission (Everitt et al., 201 I) and recently via radio pulses arriving from pulsars (Krishnan et al., 2020).
Thirrin) effects (Misner et al., 2017) are a hallmark of the
While, the mathematical description of geodesics around Kerr metric beyond scope for this work (see Teo (2003) for detailed exposition), consider solving three different cases numerically, These include Zackiger orbits (retrograde geodesics stable geodesics with larger radii), prograde orbits (stable geodesics with smaller radii), and arbitrary eccentric orbits, which depend on the initial conditions, including choice of energy E and angular momentum Lg of the test partiele.
F.6.5 KERR KRETSCHMANN INVARIANT
The Kretschmann invariant (scalar), X 1a (as)) Rapys is key curvature invariant distinguishing true (curvature) singularities from coordinate (apparent) singularities. The nontrivial part of the Kretschmann invariant for the Kerr metric reads:
Rosy R Ra876 -Capra 2r2 a²cos² (r² + 00) 16r coo=2/
(r2 + c))
(108)
This guarantees that the curvature singularity (.V + oc) occurs at the ring E = r² + coco0²0 = with zeros at:
r 0, and, 0 2
Thus, the rotation induces ring singularity radius a on the equatorial plane 0 = 7/2, where the curvature diverges. Accurately capturing this geometric structure requires isolating true singularities from coordinate artifacts, which can otherwise lead incorrect classification of singularities.
We perform training in Cartesian KS coordinates (see Eq. (81) eliminate coordinate singularities that would otherwise impede convergence. We first train EinF i elds (+Jac + Hess) on Cartesian KS coordinates, subsequently constructing the Riemann tensor see Section A.3.5.1 via successive automatic differentiation steps and raising indices using the parametrized metric 9 (see Eq. (108)).
The NeF reconstructed X Figure 13b captures the ring singularity structure and agrees well with the analytical solution, shown in Figures 13a. However, the reconstruction remains sensitive to floating-poing errors and requires high NeF accuracy for stability as seen from Figure 13c.
7788
0.00
0,00
99009
-0.89
2
0099
2
- 10.32
-178
1198
2.67
1
11.55
-3.56
33996
-12.7
-4.55
0
0
-14.01
-1
-5.34
-1
-6.93
-5.4
-1
-16447
-2
-7.12
-2
7792
11.70
-8.01
-891
-18.93
-3
-2
-3
22
(a) Analytic form X
(b) EinFields reconstructed X
(c) Absolute enror I.X
X1
Figure 13: The Kretschmann sealar X of the Kerr metrie computed in Cartesian Kerr-Schild form
(Eq. (81)) in the z-y plane for 0.3.
51

---

Published as a conference paper at ICLR 2026
F.7 TRAINING ON VARIED COORDINATE SYSTEMS
Table 8: Relative L. enror considered on a grid of validation collocation points EinFields, (ii)
EinFields (+Jac) and, (iii) EinFields (+Jac + Hess) supervision. As described above in the text, we quantify the effect of inputs queries varied coordinate charts and how EinFie inFiiels training generalizes over these different metric (geometry) representations.
Metric
Representation
Coordinate
Rel. L2
Spherical
2.26e-7
EinFields
Cartesian Kerr-Schild
1.37e-5
Mudington-FFnkelstein
9.21e-9
Spherical
1.37e-7
Schwarzschild
EinFields (+ Jac)
Cartesian Kerr-Schild
3.00e-6
Finkelstein
6.47c-9
Spherical
1.20c-7
EinFields (+ Jac + Hess)
Cartesian Kerr-Schild
1.53e-6
Fidinngton-Finellsteinn
9.08e-9
Boyer-Lindquis.
6.95e-8
EinFields
Cartesian Kerr-Schild
4.47e-6
Eddington- Finkelstein
6.44e-8
4.72e-8
Kerr
EinFields (+ Jac)
Cartesian Kerr-Schile
8.83e-7
Faddington-Fineelstenn
4.95e-8
Moyer-Lindqure
4.69c-8
EinFields (+ Jac + Hess)
Cartesian Schr-S
4.95c-7
Fddinntoon
4.72e-8
Results reported in Table suggest that the choice of coordinates has strong impact on the metric up to three orders of magnitude. This aspect should be investigated further in future work.
NeFs take physical coordinates inputs and map them directly to field values. Unlike traditional machine learning architectures that ingest abstract learned feature spaces (such as token embeddings or extracted features). INRs operate directly on the physical coordinate space. enabling them to represent continuous signals in domain-agnostic manner.
In the context of GR, this implies that lour-dimensional representation of metric tensor fields by an INR explicitly depends on the input coordinate system, more generally, the chosen frame reference. Despite this apparent dependency, GR possesses the fundamental property of diffeomorphism covariance (see Section A.3.8), which asserts that the laws of gravitation remain invariant under smooth coordinate transformations, However, the choice of coordinate system remains an essential practical tool for simplifying the form of the metric tensor. For example, while the Schwarzschild metric is diagonal spherical coordinates (albeit with coordinate singularity at the event horizon), transforming Cartesian KS coordinates produces dense, off-diagonal metric representation, or, for that matter, moving to Sädinntonn-innellee coordinates, which both remove coordinate-related artifacts (see Paragraph A.3.5.2).
Understanding this behavior is essential for developing robust INR-based frameworks for represent- ing geometric quantities in numerical relativity, while respecting the underlying diffeomorphism invariance of general relativity. For the Schwarzschild case, we initiate the training by sampling query spacetime coordinates spherical representation (t. r, 6, Φ). These sampled collocation points are then transformed into their corresponding collocation points in Cartesian coordinates (1, V, =) and ingoing Beddington-Finkesteee coordinates (v, 0, Ø) (see Section B. for explicit transfor- mation details). Subsequently, FinF weldd outputs the metric tensors corresponding to these coordinate systems, yielding Eqs. (70, 74, 77). For the Kerr metric, which is characterized by its oblate spheroidal geometry, we sample query collocation points the Boyer-Lindquis coordinates (t, r, 0, σ). followed by the collocation points transformed into Cartesian (1, y, and ingoing
52

---

Published as a conference paper at ICLR 2026
Eddington Fäddington-Finkelstenn Finkelstein coordinates (t, 0. 0). EinFields then outputs the Kerr metric tensors in these respective coordinate systems, resulting in Eqs. (80, 81, 83).
Table 9: Col. 1 lists the spacetime metrics (Schwarzschild and Kerr). Cols. 2-4 indicate the coordinate charts used for NeF training: spherical-like, Cartesian-like, and lightcone-like. For Schwarzschild, these correspond to spherical coordinates (t, r. B, Ø), Cartesian Kerr-Schild (KS) coordinates (, V/ =). and ingoing Fadington-Finkelsteenn (EF) coordinates (v, 0, Φ). trained on the metrics described Eqs. (98-100) respectively. For Kerr, these correspond to Boyer-Lindquist (BL) coordinates (t, r, 0, Φ). Cartesian KS coordinates (t, z, y, 2), and ingoing EF coordinates (v, 7, 0, Ø). trained on the metries in described in Eqs. (101-102) respectively.
Metric
Spherical-like
Cartesian-
Lightcone-like
Schwarzschild
Kerr
This multi coordinate training strategy ensures that the neural tensor field learns consistent repre- sentations across coordinate systems while maintaining geometric and physical consistency under diffeomorphisms, facilitating generalization and stability in downstream geometric learning tasks.
53

---

Published as a conference paper at ICLR 2026
F.7.1 LINEARIZED GRAVITY: GEODESIC DEVIATION, GRAVITATIONAL WAVE STRAINS AND WEYL SCALARS
Linearized gravity models the solution of the EFEs via periodic perturbations on fixed background metric. These linearized solutions are highly relevant in numerical relativity, as they describe the groundbreaking, experimentally verified discovery of gravitational waves generated by binary black hole mergers (Abbott et al., 2016c). The metric tensor can written as
ga MaB + hag + (haa))²
(109)
where |haal < is the perturbation term. As detailed in Section B.3, plane gravitational wave propagating in the 2-direction with angular frequency w can be described in the rranverse-tracelerr (TT) gauge as
0
0
0
0
0
h+
hx
0
cos (w(l -2]] z))
(110)
0
hx
h+
0
0
0
0
Here, h4 and hx are the amplitudes of the "4" (plus) polarization and "x" (cross) polarization.
Validation problems for GW metric and derivatives quality. Compared to Schwarzschild and Kerr metrics, a key distinction of the linearized gravity setting describing gravitational waves is its time dependence - see Eq. (1 (110). Although it does not depend on æ and y. the temporal dependence motivates us to consider our model trained on full spacetime grid of size N. N₂ Ny Nz.
Distortion of ring of test-particles. When the described gravitational wave interacts with a ring of freely falling test particles initially at rest in the x-y plane, induces periodic deformations of the ring. For a purely + polarized wave, the resulting motion causes the ring to stretch and squeeze along the and y-axes, leading characteristic "plus" deformation pattern.
The motion of the test particles under the influence of this gravitational wave is obtained by solving the geodesic deviation equation, up to leading order in the strain amplitude h4. As result, the particle trajectories in the TT gauge are
x(t) 2 cos (w(t 10---111 2)) =))) =0() s(t) 10:00(111 1 nh
1 th
(1- cos (wit z)) 3(0).
(111)
Here, x0) and y(0) denote the initial coordinates of test particle, and the time-dependent pertur- bations reflect the tidal nature of gravitational waves. The cosine dependence captures the periodic stretching and squeezing of spacetime caused by the wave as it traverses the partiele ring. Figure and Table 10 show how the famous ring oscillation experiment can be reproduced with EinF deldss This is done by parametrizing the perturbation hTT and captures the famous stretching and squeezing effect.
Weyl scalars of gravitational radiation field. The Weyl scalars are five complex quantities Vo, 11, V2 3 Va that arise in the Newman-Penrase formalism of GR (Newman & Penrose, 1962). They encode all the independent components of the Weyl tensor Capys (see Eq. (63)), representing the "free" gravitational field the part of spacetime curvature that can propagate gravitational waves, distinct from the curvature directly caused by matter. In NR and GW modeling, V1(()) the primary scalar quantity used to extract observable GW signals from simulations. It is defined as
(112) with n, A being particular choice of Newman-Penross tetrads and Ie its complex conjugate The central relation in an asymptotically flat spacetime (cf. Boyle et al. (2019) for details) is that v (t is cquivalent to the second coordinate-ting derivative of the strain h(t) h4 (t) ih (t):
V T
V = h + ihx
(113)
We compute V from the NeF-parameterrized strain hag in two distinct ways:
Note that the Weyl scalars not invariant and depend particular choice of the tetrad fields.
54

---

Published as a conference paper at ICLR 2026
1. indirectly via the Weyl tensor obtained with the differential geometric chain (see also Figures (112),
10 and 2): aB
+1aa
-
r
-
R
487
1/
2. directly via the second time-derivatives
(113)
44
Spin-weighted spherical harmonic representation for GW extraction. A quantity of central interest in gravitational waveform construction the mode decomposition of the GW strain into its angular components. The complex strain h(t, O, Φ) = h4 (1, r,0, Φ) expanded in terms of spin-weighted spherical hurmonics (SWSHs) as
i i (, 0, Φ) can be
M
hit, r. O, 0)
h -2m)m2 (0,0).
(114)
m=-É
where -2Y 6) are the SWSHs (see Eq. (90)) with spin-weight - reflecting the helicity of GWs in the TT gauge (see Eqs. (88 and 89). In practice, the dominant contributions to the strain arise from the quadrupole (E = 2, 2 modes, denoted by h 2,2=2(t)) (t), which capture the leading-order gravitational radiation (detailed in Section B.3).
1e-8
| h222 hi I/M
3.0
|R(h@iniees - han \/
2.5
Amplitude error
2.0
1.5
0.5
0.0
t/M
Figure 14: The absolute error of the amplitude between the EinFields and analytic values 22. (see Eq. (89)) at fixed radial distance R plotted against I/M. The amplitudes agree to 1E-8, indicating that EinF i elds can capture the complex strain h and subsequently 22) GW signals.
55

---

Published as a conference paper at ICLR 2026
Radiated power of GWs. Another important physical observable for GWs is the radiated power loss given by the famous quadrupole formula (Carroll et al., 2004). The ume-averaged power or luminosity radiated by GWs is given
dE p2 32
/ an
5) - - 122)
(115)
de
The particular perturbation metric in the above experiments (see Eq. (110) has equal amplitude A h4 hx for both and polarizations. As consequence, the radiated power loss simplifies
to
de
w²
(116)
dt
4
10
1e-10
10-66
10-11
1.00
Analytac: 30/M))
4.105
EinPaalia: (K/M))
4.0
8
3.650
0.75
3.5
3.195
0.50
6
2.740
2.5
2.285
0.00
2.0
4
1.830
-0.25
1.5
1.375
1.0
2
-0.70
0.5
0.465
00
2
8
10
0.010
-1.00
t
(a) Absolute error in the (z, t) plane, averaged over w-y (b) Temporal evolution of the Weyl scalar computed from the analytic and Ri nii e da perturhations (left axis) and their absolute error (right axis) at a fixed position.
slices.
Figure 15: Comparison of the real part of the Weyl scalar R(L,) (Eq. (112)) computed from the EinFields and the analytic metric. The eπrrrs are on the order of E-10 and E-11, respectively, indicating highly accurate gravitational waveform reconstruction capacities of finniilll
Table 10: Rel. for key quantities in the linearized gravity case with two different NeF architectures: (i) perturbation metric. (ii) perturbation metric trained with Sobolev loss and gradient normalization. (iii) reconstructed Ricci scalar, and (iv) reconstructed real part of the Weyl scalar 1 where V x h4 (t, x). The final column reports the absolute difference in the predicted gravitational radiation energy loss, integrated over the unit sphere.
h TT (+Jac Hess) 0g (GradNorm)
Ricci scalar R
Weyl scalar R(V,))
Luminosity dE/dt
Model
SiLU
8.56e-4
5.90e-13 1.08e-12
2.53e-5
2.7le44
SIREN
3.78e-2
9.56e-5
3.34e-4
WIRE
1.68e-2
1.55e-13 13
1.81e-5
3.69e-4
56

---

Published as a conference paper at ICLR 2026
F.7.2 ACCUMULATION OF ROLLOUT ERRORS FOR GEODESICS
Minute floating-point inaccuracies (around 1E-5 IE-6) arising from Christoffel symbols retrieved via EinF ields autoregressively accumulate when evolving the equations of motion for test particles along geodesics.
To quantify the inaccuracies between the ground truth and NeF-evolved geodesics, we compute the deviation between the position vectors r(r) E R as a function of the affine parameter (proper time) r in Cartesian coordinates. Specifically, for the ground truth trajectory, the spatial coordinates corresponding the position vector are given by r(r) (a(r), U(T), (()), while for the NeF- evolved trajectory, we denote f(r) = (()))) The deviation at each proper time T is then computed as the Euclidean norm,
5r(r) = ||((r) - f(r)ll = (x(r)
(u(r) - ((T))
(117)
In practice, the geodesic trajectories computed in (r, 8, σ) (e.g.. Boyer-Lindquist) coordinates and subsequently transformed into Cartesian coordinates before evaluating the deviation using the above expression.
While single-precisin (FLOAT32) arithmetic is sufficient for training EinFie sinFields in most ex- periments and downstream tasks presented, geodesic simulations indicate the need for FLOAT64 precision results in MAE and relative &g error for the reconstructed metric and its derivatives. Only FLCAT64 ensures the mitigation of error accumulation during temporal rollout, preserving the accuracy necessary for reliable scientific inference in gravitational physics.
Given the high sensitivity of time-stepped trajectories to such numerical inaccuracies, we quantify this error accumulation by explicitly presenting the deviation т. especially for eccentric orbits for both Schwarzschild and Kerr metric for a = 0.628. These reported for the Schwarzschild use case in Figure 16a, and Figure 16b for the Kerr metric use case, respectively. For Schwarzschild, the error accumulates stably, while for Kerr it is erratic. We hypothesize this is likely due to the stable versus chaotic nature of orbits in the respective spacetimes. Eventually, orbits diverge significantly, especially when leaving the NeF training domain.
function of the affine parameter
0.10
12.5
0.08
10.0
6r(m)
0.06
6r(t)
7.5
0.04
5.0
A
0.02
2.5
0.00
0.0
2000
4000
250
500
750
Unmolds
Biniheldd (+ Jac)
Innihelds ++ Jac + Hlsss)
Uniields
Limthhllds (+ Jac)
Iinitedds Jac Hesss
(a) Schwarzschild eccentric orbits.
(b) Kerr eccentric orbits for 0.628.
Figure 16: Geodesic rollout deviation ör over proper time T.
The results suggest that incorporating the Hessian supervision into training may introduce noise that can hinder convergence, performing worse than using metric Jacobian supervision or, for that matter, metric alone. For geodesic equations, supervising second derivatives is often unnecessary, and Jacobians alone provide significant improvements in trajectory reconstruction. However, Hessians become essential when computing Riemann tensors and curvature-related quantitics, and are required in applications such numerically solving the geodesic deviation equation (see Eq. (60)), which are typically encountered for solving for the test ring oscillation in linearized gravity use cases.
57

---

Published as a conference paper at ICLR 2026
F.8 TOMOGRAPHY: METRIC, METRIC JACOBIAN AND METRIC HESSIAN COMPONENTS
Here, we demonstrate the quality of EinFie lds parametrized metric tensor fields for the Kerr metric with spin parameter 0.7" report the mean absolute (MAE) between the ground truth and the NeF-fitted metric tensors in Figure 17. The evaluation is performed on validation grid with collocation points sampled arbitrarily within the training range but distinct from the training collocation points. Using a model configured with SiLU activations, SOAP optimizer, GradNorm, and without Sobolev regularization, we observe agreement with the ground truth up to six decimal places, achieving an MAE on the order of 1E 6. The effect of introducing losses pertaining to metric
1e-5
Jty
1.0
Oxt
Qxx
Gay
0.8
0.6
gvt
9yx
9/
9vz
0.4
Oxt
O2x
Gay
O2z
0.2
X
Figure 17: Kerr metric absolute error between ground truth (analytic) metric and the EinFields parametrized metric. The metrics depicted in the Cartesian Kerr-Schild (KS) representation as presented in Eq. (81). The 2D slice of all the metric components captured in the x-y plane at fixed 0.7.
1.4 for spin parameter value
Jacobian and Hessian supervision, apart from the metric loss that EinFields predominantly uses, can be quantified and visualized with the following plots below, Here, for the sake of visualization, we do tomography (2D cuts) of different metric components along particular axis for the Kerr metric in Cartesian KS coordinates (Eq. (81)).
The first, second, and third columns in each figure correspond to EinFields training without Sobolev supervision, FinFis e ds (+Jac), and EinFis el (+Jac + Hess) trained, respectively, for randomly sampled components of differential geometric quantities.
"aCattssian Kerr-Schild coordinates chosen to avoid coordinate singularities, enabling tomography over larger coordinate ranges.
58

---

Published as a conference paper at ICLR 2026
12,3, component
[2.3,11 component
[2.3.1] component
1.3
1.5
1.0
1.0
0.3
0.5
[0.2.2] pomponent
compotent
[0.2.2] compunent
0.0
Jacobian value
-50
-100
-1.0
component
componest
X coordinate
Ground truth
DinFields
Errur
EinPields (+joc)
EinCields (+jac I[ees)
Figure 18: 2D Tomography of Kerr metric Jacobian components in Cartesian KS representation.
1,,3,100] companent
[2,3,110] component
12,3, component
10
2.22
1.0
Hossian value
compunent
component
component
X coordinate
Ground truth
Finliss
Error
(+Jac)
EinFiolds (+jac Hees)
Figure 19: 2D Tomography of Kerr metric Hessian components in Cartesian KS representation.
59

---

Published as a conference paper at ICLR 2026
[3.2.] companent
3.2.1] component
[3.2]] compunent
0.0
-6-
Christoffel symbols value
12, 23 component
12,12 component
12.1.2 companent
-10
compunent
[0.3.2] component
component
x coordinate
Ground truth
FinFields
10.0x Firor
FinFields (+Joc)
FinFields (ajac Hess)
Figure 20: 2D Tomography of Kerr Christoffel symbols components in Cartesian KS representation.
component
[10.2,1] component
component
20
-20
Riemann tensor value
[1.10.3 31 component
111.0.33 0.331 compunent
18.1.0.3 component
40
20
-20
10
component
11.133.11 component
component
-2.5
+5.0
coordinate
Grouzdtruth
EinPields
10.0x Hrror
BinPields (+fec)
BinFields (+Jor Hess)
Figure 21: 2D Tomography of Kerr Riemann tensor components in Cartesian KS representation.
Kretschmamm hennnn value
1000
2000
-1000 1000
10000
-2000
coordinate 10.0k Error
Ground truth
EinFlelds
BinFields (+Jac)
EinFields (+Jac + Hess)
Figure 22: 2D Tomography of Kerr Kretschmann invariant in Cartesian KS representation.
60

---

Published as a conference paper at ICLR 2026
F.9 TRAINING HYPERPARAMETER
Table 11: Training configurations for Schwarzschild, Kem, and GWs used in the geodesies, Kretschmann plots, Table 8, and linearized gravity section.
Parameter
Schwarzschild
Kerr
GWs
Architecture
MLP
Depth
3/3/5/7
Width
64 / 128 256 512
190
128 / 128 /90
Activation
SiLU
SiLU
SiLU / SIREN / WIRE
Input dimension
4
4
4
Output dimension
10 / 16
16
16 85K
# Parameters
13.5K 50K 332K/ 1.5M
185K
Optimizer
SOAP
31
0.95
82
0.95
Precondition frequency
Learning rate schedule
Initial learning rate
E-2 E-3
Decay steps
10 E-5 / B-6
6 x 10 22//22,4 4, 6 x 10 E-7/B-8
x 103 /4 x 10* E-9
Final learning rate
Training
Epochs
100
200 30 / 100 200 / 300
200 20 200
Number of batches
100
Gradient weighting scheme
None / GradNorm
F.10
HARDWARE & LICENSES
For primary computational work, utilize high- pegh-performanc CPU system equipped with 2x32-core Intel Xeon@ Platinum 8452Y+ processors, each operating at 4.1 GHz, and 2048 GiB of RAM. All NeF-related training is performed single NVIDIA H200 SXM GPU with 144 GiB of HBM3e memory. For prototyping and preliminary experiments, we employ single NVIDIA Tesla A100 GPU with 40 GiB of memory.
This work would not have been possible without the open-source software ecosystem. Our implemen- tation is built upon multiple community-maintaiinn libraries, and gratefully acknowledge their licenses below. The core computations were performed using JAX [cuda12] (Bradbury et al., 2018) 8) with CUDA support, licensed under the Apache 2.0 License. For model definition and training, we relied on Equinox (Kidger& & Garcia, 2021) and Flax. .. Linen (Heek et al., 2024), both also under Apache 2.0. For solving differential equations, we employed Di ffra (Kidger, 2021), distributed under the Apache 2.0 License. All libraries used are permissively licensed, enabling free academic and non-commercial research.
G ADDITIONAL EXPERIMENTS
G.1 OSCILLATING NEUTRON STAR NR SIMULATION SPECIFICATIONS
Simulator setup and solver. The system (BSSN) formulation of Einstein's equations (Shibata & Nakamura, 1995; Baumgarte & Shapiro, 1998) as implemented in the McLachlan code, together with reneral--elatiiistti hydrodynamics provided by GRHydro of Einstein Toolkit (Löffler et al., 2012). Reconstructing the full 4D spacetime metric 9.æ around the neutron star from the numerical BSSN solver output variables is done by collecting the ADM variables (Arnowitt et al., 1959) (via
evolved using the Baum-
ShhapiroSSiibtaaNaaaa
61

---

Published as a conference paper at ICLR 2026
the ADMBase Thorn in EinsteinToolki) in the following manner.
9pw -
0² B,8
B1
(118)
B
76t
where a corresponds to the lapse function, 3; :- [.3z. By, 32) is the shift-vector and the 3-metric Tij == {Taz, Taas ".. J. J. The simulation employs fixed mesh refinement (FMR) (Schnetter et al., 2004), i.e., grid patches of con-uniform resolution. The coarsest grid usually encloses the whole simulation domain. Successively finer grids overlay the coarse grid at those locations where a higher resolutions is needed. This done by the Carpet infrastructure (Löffler et al., 2012). Four distinct, static grid configurations are used, namely low, medium, high, and highest resolutions, which differ only by their number of collocation points per refinement level. FMR requires far less resources than globally increasing the resolution. No dynamic or adaptive regridding is performed during the evolution; instead, each run is performed independently at fixed resolution to assess convergence and scalability. into unigrid application with minimal changes its structure. Instead of only one grid, there are several grids or grid patches with different resolutions
Training data generation. The training data are obtained by aggregating all spatial collocation points (mi, Vi. 2;) from the refinement levels (r10, rl1..... r14] at single time slice t. The NeFs are trained on this coalesced hierarchy, resembling an effective multiresolution representation, while discarding any duplicate points, yielding single non-uniform grid.
Table 12: Grid refinement hierarchy and resolution parameters for the fixed mesh refinement (FMR) configuration used in the TOV benchmark. The spatial resolution Ar and number of time slices N₂ increase with refinement level, while the physical domain shrinks correspondingly. The last column lists the full 4D grid shape (N,, Nı, Ny N₂) for each refinement level. Apart from the medium resolution simulation grid, all other grids don' 't contain the ghost zones.
4D grid shape (N,, Nx. N.)
NR Simulation Dataset
Refinement level
Spatial domain
Ax
r10
38.4, 345.6]
12.8
(313, 31, 31, 31)
r11
19.2, 217.6
6.4
(626, 38, 38)
Medium resolution (simulation grid)
r12
9966, 108.8
3.2
(1251, 38, 38, 38)
r13
4.8, 52.8)
1.6
(2501, 37, 37, 37)
r14
2.4, 24.8]
0.8
(5001,
35, 35, 35)
r10
[0.0, 307.2
12.8
(313, 25, 25, 25)
rl1
[0.0, 198.4
6.4
(626, 32, 32, 32)
Medium resolution (training grid)
r12
[0.0, 99.2
3.2
(626, 32, 32,
32)
r13
[0.0, 48.0
(626, 31, 31, 31)
r14
[0.0, 22.4
0.8
(626,
29, 29,
29)
r10
[0.0, 307.2|
8.0
(501,
40, 40, 40)
High resolution (evaluation grid)
rl1
[0.0, 198.4
(1001,
40, 10, 40)
r12
[0.0, 99.2
2.0
(2001,
10,
40,
10)
r13
[0.0, 4
1.0
(4001,
40,
40, 40)
r14
[0.0, 22.4]
0.5
(8001,
40,
40,
40)
Tensor differentiation on NR grids. Numerical relativity solvers augment the primary computational domain with ghost zones (Thornburg, 2004), which provide additional layers of points surrounding the grid, These zones support high-order finite difference stencils, ensure that boundary conditions are applied consistently, and help manage coordinate singularities. The simulations used here employ three ghost zones in each spatial direction. Tensor derivatives are evaluated using a sixth-order central finite difference stencil applied the metric tensor grid functions:
((()) f(z + 3h) 9f + 2h) (z + h) + 45f(z h) -9((x -2h) f(x 3h)
60h
1eSee http://ehheinnooliiioolllleeieie for details
62

---

Published as a conference paper at ICLR 2026
The three ghost zones required because this stencil accesses three neighbouring points on both sides of each evaluated location.
Training specifies. The full four-dimensional numerical relativity dataset contains approximately 160 million collocation points, including contributions from the refinement hierarchy. This scale necessitates hyperparameters distinct from those used for the analytic benchmarks in EinFields. We retain SiLU activations and the SOAP optimizer, and introduce the following modifications: (i) Fourier embedding layer of size 256 for input normalization with embedding frequency scale of 0.01, (ii) a cosine learning-rate schedule with initial learning rate 10-3 annealed upto 10-5, (iii) a batch size of 10 and (iv) 200 training epochs with 71666 decay steps.
Evaluation procedure. EinFields is trained on the Table 12). As no analytical solution exists for the oscillating TOV configuration, we adopt the highest-resolution BSSN simulation (Table 12) our ground truth metric field values. Although this evaluation grid features finer spatial (and temporal) resolution at every refinement level rlx, its spatial domain restricted the same coordinate ranges corresponding the medium-resolution training grid data. We therefore evaluate EinFields on this high-resolutie grid within the identical domain range, enabling a direct comparison against the ground truth solution. This procedure yields the results reported in Table 3 and assesses the ability of the implicit representations to provide continuous query access of the metric tensor field values at coordinates not present in the training grid.
coalesced FMR grid (see
Run-time tradeoffs and query speed. We plot the query speeds of FinFields and its AD-derived
Jacobian values over the collocation points.
Metric
Jacobian
Width
Wiath 512
Query Time [s]
10
10-3
10-3
NN Depth
NN Depth
Figure 23: The time to query EinFields 10 batch of points. The model is trained the neutron star numerical relativity simulation containing approximately 71 million collocation points. There is clear trend of increasing query time the MLP model size increases. Shaded regions indicate small uncertainty (only visible when zoomed in). Timings are performed on an NVIDIA H200 GPU with jAX. JIT.
Retrieving higher derivative quantities. To extend our analysis to the differentiation of tensorial quantities, we additionally evaluate the Christoffel symbols. Specifically, we compare the Christoffel symbols obtained by our automatic differentiation based FinFields pipeline with those computed using a fourth-order central finite-difference stencil applied to the reconstructed four-metric on this multi-resolution coalesced mesh. The finite-difference operator acts on metric data that include three padded ghost zones and are generated from the BSSN formulation. We use a fourth-order stencil rather than sixth-order one to maintain numerical stability. All comparisons are carried out at high resolution, using the simulation configurations listed in Table 12, and as done for the metric tensor field. The accuracy results for the NR simulation Christoffel symbols reported in Table 13. Additionally we report the compression ratio obtained for storing such quantities in explicit form as compared to implicit machinery.
63

---

Published as a conference paper at ICLR 2026
Table 13: Performance and compression for oscillating neutron star: EinFields retrieved Christoffel symbols as compared to FD methods applied on the multi-resolution coalesced FMR grids.
Geometric quantity
Storage [GiB]
MAE
Compression factor (Sym.)
Full
Sym.
EinFields (AD)
Christoffel symbol
17.0
10.6
8.61e-4
7753
Static mesh refinement evaluation. the previous set of numerical relativity experiments, both the simulation and training grids preserved fixed block structure over the full BSSN time evolution. Each refinement level was pre-defined and spatially static, with no patch motion or regridding. However, technical inconsistency arises from the fact that coalesced FMR grids, as commonly used, contain coarse-level cells embedded inside regions where finer patches are available. As result, multiple grid resolutions represent the same physical region simultaneously. Retaining all levels produces redundant sampling density, and therefore does not reflect a strictly hierarchical representation of the spacetime domain.
A more consistent procedure is to remove coarse-level grid points that fall within the bounding region of any finer-level patch. practice, we retain only the resolon- data available at each spatial location. This construction yields what is effectively static mesh refinement (SMR) hierarchy, analogous to a static octree, in which the resolution is high only near the neutron star and decreases smoothly with distance. Unlike adaptive mesh refinement (AMR). SMR does not dynamically move patches, perform tagging, or trigger online mesh adaptation. Nonetheless, this setting forms strong and meaningful test of implicit field models, since the resulting sampling distribution is highly non-ume uniform and reflects true octree-like scaling.
Formally, the refinement bounding boxes satisfy:
BBox(r14) C BBox(r13) C... BBooxr10),
(119)
where BBox(rlk) denotes the spatial extent of refinement level rlk, with r14 the finest and r10 the coarsest. Coarse-level points inside the domain of any finer box are removed, eliminating multi-resolution overlap. The result for this strategy are also detailed in Table 14.
Table 14: Evaluation of EinFields (best model 6 x 256) on a single stable neutron star simulation under static octree/static mesh refinement (SMR). We remove coarse-level samples that lie inside the regions covered by finer patches (Eq. 119,, retaining only the highest available resolution at each spatial location. The table reports the relative &2 error, MAE, storage footprint, and the resulting compression ratio. Under SMR reduction, EinFields attains a compression ratio of approximately 1900x while preserving low reconstruction error, suggesting that climinating overlapping coarse- resolution points improves metric fidelity in high-resolution regions.
Representation
Rel. 62
MAE
Storage
Compression
FinFields
1.09-55
4.4e-5
1.4 MiB
1901
EinFields (+ Jacobian)
6.82e-6
9.24e-6
1.4
1901
SMR grid
2.6 GiB
This benchmark provides a representative NR workflow in which matter and spacetime are tightly coupled. It therefore offers ideal context for assessing the robustness and efficiency of EINFIELDS when applied to high-resolution simulations of relativistic astrophysical systems.
64


Reviews:
[{"review_id": "2Tw74AD6K1", "reviewer": "Reviewer_kkPc", "summary": "This paper presents a neural tensor‑field (tensors in the General Relativity sense) parametrization of GR metrics with a JAX implementation for differential geometry operators wired through automatic differentiation. The demos (Schwarzschild/Kerr orbits, Weyl scalars, ring deformation under a linearized GW) look correct and generate polished images, and the ablations around Jacobian/Hessian supervision are good.\n\nThat said, I have major concerns with the claims and evaluations. I believe these are not aligned with numerical relativity (NR) practice, and so this is not yet at the point of being useful for actual science. Thus, the current narrative is a bit misleading. The headline is \"compressing 4D NR simulations by 1000-4000x with better derivatives than FD in FLOAT32,\" but most experiments are static analytic 3D cases (t = 0) plus one simple time‑varying GW; storage comparisons are made against a dense \"explicit grid\" strawman (modern NR code uses AMR or pseudo-spectral methods); baselines omit spectral/ROM methods; coordinate chart sensitivity is large; and long‑horizon dynamics need FLOAT64 to avoid divergence. I think it's a promising direction for ML to help, but I think this paper needs to be honest about the current status of such an approach. Doing so would not only be better for the paper, but I think also benefit the authors in that it would point out to the ML community where more work is needed. Basically I would like to see the narrative of this paper modified to be honest about the practicality of this, and about the toy baselines, before I consider acceptance, especially as general ML audiences will not know how to evaluate this.", "questions": "- Please clarify the precise data format used for training (\"4D spacetime grid\" vs \"arbitrary samples\"), since the paper simultaneously describes the approach as discretization-free yet refers to training on regular grids.\n- Could you provide any comparison, even small-scale, against a spectral or AMR baseline (ideally an actual code used by the GR community) to contextualize the claimed compression?\n- Please explain how derivative accuracy scales in FLOAT64 and whether the observed long-time geodesic divergence persists under higher precision.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- I like the clear pipeline and library, and a JAX code for this seems useful for the NR community. The graph from metric to derived quantities is explicit and leverages forward‑mode Jacobians/Hessians with einsum operations which is nice. This is a useful contribution of reusable tooling for GR in ML, and I think it is a good contribution by itself.\n- I liked the ablations accounting done in 4.2, it is nice and I think quite useful to see the effect of every modification to the training process. It is interesting to see that Jacobian/Hessian supervision indeed helps.\n- Neat canonical tests: precession, circular/eccentric Kerr orbits, ring response under a linearized GW, and Weyl/Kretschmann diagnostics are shown and largely match analytics on short horizons.\n- Multi‑chart training attempt: training/evaluating in multiple coordinate charts acknowledges a real pain point in GR ML (though this is near the end of the appendix, I think it should come earlier)", "weaknesses": "First, my main concerns:\n\n1. First, I think the scope and storage comparisons are misaligned with NR. The paper advertises 4D compression of NR simulations, but the primary experiments are analytic snapshots at t = 0 (Schwarzschild/Kerr); only the linearized‑gravity toy has time evolution. The compression factors compare MLP weights to explicit dense grids counted as \"#points x 4 bytes\" in FLOAT32, which is not how NR codes actually store data (they would use adaptive mesh refinement stored with an octree). So the 1000–4000x headline is basically comparing against a strawman and misleading to the ML community about the state of this domain.\n2. The paper itself notes that modern NR \"increasingly opts for (pseudo‑)spectral methods ... up to 1000–5000x faster on CPUs than FD on GPUs at comparable accuracy.\" But all quantitative baselines in this paper are finite difference stencils (on a uniform grid - but the only actual finite difference codes used in NR are based on adaptive meshes) and an analytic AD. There's no head‑to‑head vs actual NR codes used in GW modeling. This makes it hard to situate their method, even for someone who knows NR, let alone the ML community.\n3. The paper claims up to five orders of magnitude derivative gains over FD in FLOAT32, but geodesic integration requires FLOAT64 and long‑time rollouts still diverge (the authors show this in their own figures and explicitly state that they only surpass FD in single precision). This puts the method far from NR‑readiness, where double precision (or even higher) is standard.\n4. There is large coordinate‑chart sensitivity which is a bit worrisome. Table 8 (deep in the appendix) reports up to three orders of magnitude variation in \"Rel‑L2\" error across charts for the same spacetime. That undermines generality claims unless the representation or training explicitly handles diffeomorphisms or evaluates with chart‑invariant metrics.\n5. The physics is not enforced or audited. The pipeline mentions Bianchi identities, but experiments focus on pointwise tensor errors, scalar invariants, and geodesic tests. While the pointwise tensor errors are no doubt useful in clarifying (to a NR person) that these methods aren't yet ready, there's no reporting of physics checks, like conservation laws, etc., which are exactly the diagnostics one needs to trust a compressed metric in downstream NR workflows.\n6. I am a bit confused about the \"discretization‑free\" claims, since it seems the method is ultimately trained on a grid. Several places describe training on 4D spacetime grids or \"4D training and validation grid data,\" which undercuts the claim of being discretization‑free. Even an INR is ultimately a finite parametrization/basis.\n7. The throughput trade‑off is not discussed. Even if the file is tiny, post‑processing requires many MLP queries to reconstruct fields, whereas decompressing spectral grids is reading coefficients + evaluating polynomials. You still likely win on storage, but the compute/runtime trade‑off for analysis & viz should be stated.\n\nSecond, some other suggestions/comments:\n\n- I make the following comment purely to help the authors improve their work, and do not include this comment as part of my score for the paper, so feel free to not address this in your rebuttal. It is simply a suggestion/idea. So I think the branding of the method as \"Einstein Fields\" is not an optimal choice, because it would likely conflict with \"Einstein Field Equations\" in any search. Also, to the physicist, who I assume you would like to include in the audience of the paper, it does not give them any idea that this is related to machine learning.\n- Obviously GR is quite complicated to someone with no background. I am not sure it is possible to give it much of an introduction here, and I worry the current introduction might give the wrong ideas. I think it is better to simply direct the reader to an online resource, rather than give an inevitably incomplete description of the mathematics in the appendix. Try to think think about what purpose it serves: (1) for those who don't know GR, this is not going to be nearly enough to introduce them even to the basics; (2) for those who do know GR, they will not need this anyways. So, why include it at all? Consider, for audience (1) I think you need to simply target the intuition for each variable you model. That is the fundamentally useful thing to write about. And for audience (2) (curious physicists, maybe) I think you simply need to _translate_ GR concepts to machine learning for them. So, when you consider these two audiences, the appendix seems to serve little purpose. I recommend trying the split approach above: focus on intuition of the key target variables for the non‑physicists (and point them to other resources), and focus on translation of the machine learning concepts for the physicists. This would be much more effective in my view. \n  - Also, I did not check through all the math in the appendix. Thus, there could be errors.\n- It might be worth calling a GR tensor exactly that: a \"GR tensor,\" to differentiate from the ML meaning.\n- Use proper scientific notation instead of \"5.37E‑6\" in tables. \n- Even if the file is small, you still must run tens of thousands of MLP queries for post‑analysis/viz; be transparent about that runtime trade‑off versus reading spectral coefficients. \n- Figure 1 is noisy. Consider simplifying or splitting. (The caption itself also states training on a 4D spacetime grid, which conflicts with \"discretization‑free\" messaging.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8foyZsKxH5", "reviewer": "Reviewer_Lt1j", "summary": "This paper derives a novel neural representation method to compress the relativity simulations.", "questions": "See above.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "- Outstanding efficiency and accuracy are shown when representing the symmetry of the simulations.\n- This paper is especially well written and well presented.\n- The problems addressed by the new tool is of interest to a wide community.", "weaknesses": "- I am not entirely sure that how interesting this paper will be for the readership of ICLR, of whom so few are well versed in this area.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "GF5GU3kxny", "reviewer": "Reviewer_DSBV", "summary": "This paper introduces Einstein Fields (EinFields), a neural field approach to compress and represent 4D spacetime metrics from general relativity simulations. The method achieves up to 4,000× compression of metric tensor fields with 7-9 decimal digit accuracy (1E-7 to 1E-9 relative precision) while providing discretization-free continuous representations that can be trained on arbitrary point samples and queried at arbitrary resolutions. A key contribution is improved differentiation accuracy through automatic differentiation (AD), computing Christoffel symbols, Riemann tensors, and other derived quantities with up to 10^5 better accuracy than finite difference methods in FLOAT32. The approach parametrizes the metric distortion (deviation from flat space) using MLPs and employs Sobolev supervision. The validation focuses on three analytical solutions to Einstein's field equations: Schwarzschild, Kerr, and linearized gravitational waves, successfully reconstructing key relativistic phenomena.", "questions": "## Questions\n\n\n* Actual NR simulation data: Your validation strategy using analytical solutions (Schwarzschild, Kerr, linearized GW) with known ground truth is appropriate for demonstrating the method's capabilities. As a natural next step, have you tested EinFields on any actual numerical relativity simulation outputs, even at small scale? What additional challenges arise with real NR data?\n\n* Table 1 mentions selecting \"the model with the lowest MAE\" - how many training runs were performed? Can you report mean ± standard deviation over multiple random seeds for the key results in Tables 1-3? This is important for assessing reproducibility and typical vs. best-case performance.\n\n* Parameter generalization: Do you train a separate network for each physical configuration (M, a, etc.), or can one network generalize across parameter ranges?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "## Strengths \n\n**Originality**: This paper presents a novel application of neural fields to general relativity, introducing the first implicit neural representation for tensor-valued spacetime geometries. The approach creatively adapts neural field techniques from computer vision to computational physics, with several original contributions.\n\n**Quality**: The paper demonstrates strong technical rigor with comprehensive validation across multiple canonical GR test cases (Schwarzschild, Kerr, gravitational waves). The evaluation methodology is sound. Ablation studies (Table 3) properly isolate the contributions of different design choices. The authors are transparent about limitations.\n\n**Clarity**: The paper is well-written and accessible. The background section (Section 2) effectively introduces both GR concepts and neural fields. Figure 1 provides an excellent conceptual overview of the pipeline. The mathematical notation is consistent and properly defined (though dense in places).\n\n**Significance**: This work addresses genuine computational bottlenecks in numerical relativity—storage (petabytes per simulation) and accurate tensor differentiation. The 4,000× compression factor and 10^5 improvement in derivative accuracy (FLOAT32) represent substantial practical gains.", "weaknesses": "## Weaknesses\n\n**Limited experimental scope**: The validation is restricted to three analytical solutions to Einstein's field equations (Schwarzschild, Kerr, linearized gravitational waves). While these are canonical test cases, they represent idealized scenarios far simpler than realistic numerical relativity (NR) simulations. \n\n**Limited contextualization within scientific computing:** While the introduction mentions neural fields and ML for scientific computing, it lacks: (1) discussion of prior ML work specifically targeting numerical relativity or gravitational physics, (2) comparison with traditional compression methods used in scientific computing, and (3) detailed positioning relative to neural operators and PINNs. A dedicated related work section would help readers better understand the landscape and the paper's specific contributions.\n\n**Missing Error Quantification:** Tables 1-3 report single-valued metrics without error bars or confidence intervals. Table 1 mentions selecting 'the model with the lowest MAE,' suggesting multiple runs were performed but statistics are not reported.\n\n### Minor issues:\n\nPage 10, line 490: \"supplimentary\" → \"supplementary\"\n\nFigure 4 a caption: \"Perihilion precession\" → \"Perihelion precession\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)