Paper deep dive
Gauge-Equivariant Intrinsic Neural Operators for Geometry-Consistent Learning of Elliptic PDE Maps
Pengcheng Cheng
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%
Last extracted: 3/22/2026, 5:13:49 AM
Summary
The paper introduces Gauge-Equivariant Intrinsic Neural Operators (GINO), a novel architecture for learning solution operators of elliptic PDEs on manifolds. GINO leverages intrinsic spectral multipliers of the Laplace-Beltrami operator and gauge-equivariant nonlinearities to ensure representation-independent, discretization-robust, and geometry-consistent performance for form-valued fields.
Entities (5)
Relation Signals (3)
GINO → parameterizes → Elliptic solution maps
confidence 95% · GINO, a class of neural operators that parameterize elliptic solution maps primarily through intrinsic spectral multipliers
GINO → utilizes → Laplace-Beltrami operator
confidence 92% · parameterize elliptic solution maps primarily through intrinsic spectral multipliers acting on geometry-dependent spectra
Gauge Equivariance → improves → Discretization robustness
confidence 90% · enforcing intrinsic structure and gauge equivariance yields operator surrogates that are more geometry-consistent and discretization-robust
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Learning solution operators of partial differential equations (PDEs) from data has emerged as a promising route to fast surrogate models in multi-query scientific workflows. However, for geometric PDEs whose inputs and outputs transform under changes of local frame (gauge), many existing operator-learning architectures remain representation-dependent, brittle under metric perturbations, and sensitive to discretization changes. We propose Gauge-Equivariant Intrinsic Neural Operators (GINO), a class of neural operators that parameterize elliptic solution maps primarily through intrinsic spectral multipliers acting on geometry-dependent spectra, coupled with gauge-equivariant nonlinearities. This design decouples geometry from learnable functional dependence and enforces consistency under frame transformations. We validate GINO on controlled problems on the flat torus ($\mathbb{T}^2$), where ground-truth resolvent operators and regularized Helmholtz--Hodge decompositions admit closed-form Fourier representations, enabling theory-aligned diagnostics. Across experiments E1--E6, GINO achieves low operator-approximation error, near machine-precision gauge equivariance, robustness to structured metric perturbations, strong cross-resolution generalization with small commutation error under restriction/prolongation, and structure-preserving performance on a regularized exact/coexact decomposition task. Ablations further link the smoothness of the learned spectral multiplier to stability under geometric perturbations. These results suggest that enforcing intrinsic structure and gauge equivariance yields operator surrogates that are more geometry-consistent and discretization-robust for elliptic PDEs on form-valued fields.
Tags
Links
- Source: https://arxiv.org/abs/2603.14734v1
- Canonical: https://arxiv.org/abs/2603.14734v1
Trouble viewing inline? Open PDF directly →
Full Text
159,680 characters extracted from source content.
Expand or collapse full text
Gauge-Equivariant Intrinsic Neural Operators for Geometry-Consistent Learning of Elliptic PDE Maps Cheng 1022@mails.jlu.edu.cn of Mathematics Jilin University Changchun, 130012, China Abstract Learning solution operators of partial differential equations (PDEs) from data has emerged as a promising route to fast surrogate models in multi-query scientific workflows. However, for geometric PDEs whose inputs and outputs transform under changes of local frame (gauge), many existing operator-learning architectures remain representation-dependent, brittle under metric perturbations, and sensitive to discretization changes. We propose Gauge-Equivariant Intrinsic Neural Operators (GINO), a class of neural operators that parameterize elliptic solution maps primarily through intrinsic spectral multipliers acting on geometry-dependent spectra, coupled with gauge-equivariant nonlinearities. This design decouples geometry from learnable functional dependence and enforces consistency under frame transformations. We validate GINO on controlled problems on the flat torus (2T^2), where ground-truth resolvent operators and regularized Helmholtz–Hodge decompositions admit closed-form Fourier representations, enabling theory-aligned diagnostics. Across experiments E1–E6, GINO achieves low operator-approximation error, near machine-precision gauge equivariance, robustness to structured metric perturbations, strong cross-resolution generalization with small commutation error under restriction/prolongation, and structure-preserving performance on a regularized exact/coexact decomposition task. Ablations further link the smoothness of the learned spectral multiplier to stability under geometric perturbations. These results suggest that enforcing intrinsic structure and gauge equivariance yields operator surrogates that are more geometry-consistent and discretization-robust for elliptic PDEs on form-valued fields. Keywords: Neural operators, geometric deep learning, gauge equivariance, elliptic PDEs, Laplace–Beltrami operator, spectral multipliers, Helmholtz–Hodge decomposition, discretization consistency, cross-resolution generalization, scientific machine learning 1 Introduction Partial differential equations (PDEs) govern a wide range of phenomena in physics, engineering, and the natural sciences. In many applications, the central computational task is not the solution of a single PDE instance, but the repeated evaluation of a solution operator that maps an input field (e.g., a forcing term, boundary data, coefficients) to an output field (e.g., a state, potential, or flux). Classical numerical solvers—finite differences, finite volumes, finite elements, and spectral methods—offer systematic accuracy and stability guarantees, but can be computationally expensive when embedded in design loops, inverse problems, uncertainty quantification, or multi-query regimes Hughes (2003); Brenner and Scott (2008); Temam (2024). This has motivated substantial interest in scientific machine learning approaches that learn surrogate models for PDE solution maps from data, including physics-informed methods Raissi et al. (2019); Karniadakis et al. (2021); Brunton et al. (2020) and operator-learning methods Lu et al. (2021); Li et al. (2020); Kovachki et al. (2023). Among data-driven approaches, neural operators aim to learn maps between infinite-dimensional function spaces, enabling discretization-independent inference in principle Lu et al. (2021); Li et al. (2020); Kovachki et al. (2023). Representative examples include DeepONet Lu et al. (2021), the Fourier Neural Operator (FNO) Li et al. (2020), and related integral-operator architectures Kovachki et al. (2023). These models have demonstrated impressive performance on benchmark PDE families, including fluid dynamics and diffusion-type problems Brunton et al. (2020); Stachenfeld et al. (2022); Takamoto et al. (2022). However, many existing operator-learning architectures are trained and evaluated on fixed grids with coordinate-dependent representations. As a result, their behavior can degrade under changes of discretization (cross-resolution evaluation) and under geometric changes (e.g., varying metric, mesh, or coordinate charts), and they often lack built-in guarantees of equivariance with respect to representation choices. A core challenge is that many physically meaningful PDEs are geometric: their inputs/outputs may be vector fields or differential forms, and the governing operators depend on the underlying Riemannian metric through objects such as the Laplace–Beltrami operator and the Hodge star. In such settings, the same geometric field can be represented by different component vectors under different choices of local frame or gauge. Standard convolutional networks, including coordinate-augmented CNNs, are typically not designed to respect these symmetries; their outputs may vary substantially under benign changes of representation, even when the underlying geometric field is unchanged. This issue is closely related to the broader theme of geometric deep learning, which emphasizes architectures that respect symmetry and structure—on groups, manifolds, and meshes—through equivariance principles Bronstein et al. (2021); Cohen and Welling (2016b); Cohen et al. (2019). In Euclidean settings, group-equivariant CNNs and steerable CNNs have yielded strong gains in sample efficiency and generalization by hard-coding symmetry Cohen and Welling (2016a); Worrall et al. (2017); Weiler and Cesa (2019); Bekkers (2019). More recently, gauge-equivariant CNNs extend equivariance to settings where features transform under spatially varying frames, providing a principled way to handle fields on manifolds and meshes Cohen et al. (2019); Köhler et al. (2020); Kondor and Trivedi (2018); de Haan et al. (2020). In parallel, equivariant neural networks for 3D geometry have matured rapidly, including tensor field networks and SE(3)-equivariant architectures Thomas et al. (2018); Fuchs et al. (2020); Geiger and Smidt (2022); Satorras et al. (2021), reinforcing the practical value of enforcing symmetry at the architectural level. In this work we focus on a complementary question: how to design an operator-learning architecture that is intrinsic and gauge-consistent, and that remains stable across geometric and discretization changes. Our approach is guided by classical geometric analysis. For elliptic operators, the spectrum of the Laplace-type operator encodes both geometry and scale, and solution operators can be expressed via spectral multipliers (e.g., resolvents and Green’s functions). Meanwhile, the Helmholtz–Hodge decomposition provides a canonical splitting of vector fields (or 1-forms) into exact, coexact, and harmonic components, underpinning stable formulations in both continuous and discrete settings Jost (2005); Warner (1983); Desbrun et al. (2006); Hirani (2003); Arnold et al. (2006, 2010). Discrete exterior calculus (DEC) and finite element exterior calculus (FEEC) make these ideas computationally concrete, yielding discretizations that preserve topological and geometric structure Desbrun et al. (2006); Hirani (2003); Arnold et al. (2006, 2010). Building on these insights, we introduce a Gauge-Equivariant Intrinsic Neural Operator (GINO) that parameterizes the learned operator primarily through intrinsic spectral multipliers of the form mθ(λ)m_θ(λ), where λ denotes a geometry-dependent spectral quantity (e.g., λg(k)=k⊤Ak _g(k)=k Ak for constant-metric problems). This design explicitly separates (i) geometry-dependent spectral structure from (i) learnable functional dependence, and we couple it with a gauge-equivariant nonlinearity to ensure consistent transformation under frame changes. The resulting model is intended to behave like a geometry-aware, discretization-consistent surrogate for elliptic solution operators on fields with nontrivial transformation laws. We validate the approach in a controlled setting where the ground-truth operator is available in closed form in the Fourier domain, enabling diagnostics that directly probe the theory. Across experiments E1–E6, we show that: (i) GINO matches the target operator accurately; (i) it achieves gauge equivariance up to numerical precision, while coordinate-based CNN baselines exhibit large gauge sensitivity; (i) it remains stable under structured metric perturbations; (iv) it generalizes across resolutions and approximately commutes with restriction/prolongation operators, indicating reduced discretization dependence; (v) it supports a structure-preserving (regularized) Hodge decomposition task; and (vi) ablations link multiplier smoothness to stability, aligning with theoretical expectations about spectral regularity and perturbation amplification. We propose an intrinsic operator-learning architecture that parameterizes elliptic solution maps via learned spectral multipliers coupled with gauge-equivariant nonlinearities, aimed at representation- and discretization-consistent behavior. We design a suite of theory-aligned diagnostics for gauge equivariance, metric-perturbation stability, and discretization commutation, enabled by exact Fourier-domain ground truth. We demonstrate strong empirical advantages over coordinate-based CNN baselines on accuracy, gauge consistency, metric stability, and cross-resolution generalization, and we provide ablations linking multiplier roughness to perturbation amplification. 2 Problem Setup and Target Operator 2.1 Geometric and Functional-Analytic Setting Let (ℳ,g)(M,g) be a connected, compact, oriented, d -dimensional smooth Riemannian manifold without boundary. We denote by Ω1(ℳ) ^1(M) the space of smooth differential 11 -forms on ℳM . The Riemannian metric g induces a pointwise inner product ⟨⋅,⋅⟩g ·,· _g on 11 -forms and the associated L2L^2 inner product ⟨ω,η⟩L2:=∫ℳ⟨ω(x),η(x)⟩gvolg(x),ω,η∈Ω1(ℳ), ω,η _L^2:= _M ω(x),η(x) _g\,dvol_g(x), ω,η∈ ^1(M), with norm ‖ω‖L22=⟨ω,ω⟩L2\|ω\|_L^2^2= ω,ω _L^2 . We write HsΩ1(ℳ)H^s ^1(M) for the Sobolev space of 11 -forms of order s∈ℝs , defined via local charts and a partition of unity (equivalently, through the functional calculus of a fixed elliptic operator on 11 -forms). When no ambiguity arises, we abbreviate HsΩ1(ℳ)H^s ^1(M) as HsΩ1H^s ^1 . We denote by d:Ωk(ℳ)→Ωk+1(ℳ)d: ^k(M)→ ^k+1(M) the exterior derivative and by δ:Ωk+1(ℳ)→Ωk(ℳ)δ: ^k+1(M)→ ^k(M) its L2L^2 -adjoint (the codifferential). The (Hodge) Laplacian on 11 -forms is ΔH:=dδ+δd:Ω1(ℳ)→Ω1(ℳ). _H:=dδ+δ d: ^1(M)→ ^1(M). On a compact manifold without boundary, ΔH _H is essentially self-adjoint, nonnegative, and has discrete spectrum. Throughout, we use the same notation for its continuous extension to Sobolev spaces. To make the inverse problem well-posed without explicitly quotienting out harmonic 11 -forms, we introduce a positive shift parameter α>0α>0 . Define the shifted operator ℒg:=ΔH+αI,L_g:= _H+α I, where I denotes the identity on 11 -forms. The shift ensures ℒgL_g is strictly positive and invertible on L2Ω1L^2 ^1 , and improves conditioning both analytically and numerically. Finally, since vector fields are often the physically meaningful objects, we note that g provides an isomorphism between 11 -forms and vector fields via the musical isomorphisms ♭,♯ , . All statements below may be equivalently read in the language of tangent vector fields, but we adopt the 11 -form formulation because it is naturally compatible with Hodge theory and with discretizations based on discrete exterior calculus. 2.2 Elliptic PDE on 1-Forms and the Solution Operator We consider the following family of elliptic problems on 11 -forms: for a given forcing f∈Hs−1Ω1f∈ H^s-1 ^1 (with s≥0s≥ 0 typically), find ω∈Hs+1Ω1ω∈ H^s+1 ^1 such that ℒgω=f,i.e.,(ΔH+αI)ω=f.L_gω=f, .e., ( _H+α I)ω=f. (2.1) Equation(2.1) covers, as special cases or close relatives, the vector Laplacian and Hodge–Poisson problems arising in geometry processing, computational electromagnetism, and fluid mechanics on surfaces (e.g., via vorticity–stream formulations or Hodge decompositions). The Hodge Laplacian is the canonical second-order elliptic operator acting on differential forms, and thus provides a clean testbed for developing intrinsic, coordinate-free operator learning theory. Standard elliptic regularity on compact manifolds implies that ℒgL_g defines a continuous isomorphism between Sobolev spaces, ℒg:Hs+1Ω1→Hs−1Ω1,L_g:H^s+1 ^1→ H^s-1 ^1, and therefore admits a bounded inverse (the resolvent at −α-α ), g:=ℒg−1:Hs−1Ω1→Hs+1Ω1.S_g:=L_g^-1:H^s-1 ^1→ H^s+1 ^1. (2.2) We refer to gS_g as the (shifted) Hodge–Poisson solution operator. In this work, gS_g is the primary target operator to be learned. Spectrally, since ΔH _H is self-adjoint and nonnegative on L2Ω1L^2 ^1 , there exists an L2L^2 -orthonormal eigenbasis ψkk≥0⊂Ω1(ℳ)\ _k\_k≥ 0⊂ ^1(M) with eigenvalues 0≤λ0≤λ1≤⋯0≤ _0≤ _1≤·s , λk→∞ _k→∞ , such that ΔHψk=λkψk _H _k= _k _k . Consequently, for sufficiently regular f , gf=∑k≥01λk+α⟨f,ψk⟩L2ψk,S_gf= _k≥ 0 1 _k+α f, _k _L^2\, _k, (2.3) which exhibits gS_g as a spectral multiplier (a functional calculus) of ΔH _H . This representation motivates the model class introduced later: intrinsic neural operators parameterizing spectral multipliers while preserving geometric equivariances. 2.3 Learning Objective and Data Model We assume access to samples (fi,ωi)i=1n\(f_i, _i)\_i=1^n drawn from a distribution D over forcings f and corresponding solutions ω=gfω=S_gf to (2.1), where the metric g is fixed unless stated otherwise. Concretely, the training pairs may be generated by a classical numerical solver on a discretization of ℳM , or obtained from observations when the underlying elliptic model is appropriate. Our goal is to learn a parameterized operator family ^θ∈Θ\ S_θ\_θ∈ mapping 11 -form forcings to 11 -form solutions, ^θ:Hs−1Ω1→Hs+1Ω1, S_θ:H^s-1 ^1→ H^s+1 ^1, that approximates gS_g in a strong, geometry-aware sense. Specifically, we target guarantees of the following form: Operator approximation in Sobolev topology: ‖^θ−g‖Hs−1→Hs+1\| S_θ-S_g\|_H^s-1→ H^s+1 is small, with explicit dependence on model capacity (e.g., spectral truncation level) and parameterization error. Gauge equivariance: the learned operator is equivariant with respect to changes of local orthonormal frames, reflecting the intrinsic nature of 11 -forms (formally defined in Section3). Discretization consistency: when ^θ S_θ is implemented on a mesh or point cloud via a convergent discrete Hodge Laplacian, the discrete learned operator converges to its continuum counterpart as resolution increases. In Sections 4–6 we formalize and prove approximation, stability, and continuum–discrete consistency results under the above setup. 3 Gauge-Equivariant Intrinsic Neural Operators This section introduces a class of intrinsic neural operators acting on 11 -forms that (i) are defined without reference to coordinates, (i) are equivariant to local changes of orthonormal frames (gauge transformations), and (i) admit discretizations that commute with refinement under standard convergent discrete Hodge Laplacians. Our construction is motivated by the observation that the target solution operator g=(ΔH+αI)−1S_g=( _H+α I)^-1 is a spectral multiplier of ΔH _H (cf.(2.3)), and hence can be approximated by learning appropriate functional calculi of ΔH _H . 3.1 Gauge Structure on 11 -Forms and an Equivariance Principle Although a differential 11 -form is an intrinsic geometric object, any numerical representation necessarily uses local frames. Let Uℓ=1L\U_ \_ =1^L be an open cover of ℳM and, on each UℓU_ , let Eℓ(x)=(eℓ,1(x),…,eℓ,d(x))E_ (x)=(e_ ,1(x),…,e_ ,d(x)) be a smooth local orthonormal frame of TxℳT_xM . Any ω∈Ω1(ℳ)ω∈ ^1(M) admits local coordinates [ω]Eℓ(x)∈ℝd[ω]_E_ (x) ^d defined by [ω]Eℓ(x)j:=ω(x)(eℓ,j(x)).[ω]_E_ (x)_j:=ω(x) (e_ ,j(x) ). If the frame is changed to Eℓ′(x)=Eℓ(x)Rℓ(x)E _ (x)=E_ (x)R_ (x) with Rℓ(x)∈O(d)R_ (x)∈ O(d) , then the coordinate representation transforms as [ω]Eℓ′(x)=Rℓ(x)⊤[ω]Eℓ(x).[ω]_E _ (x)=R_ (x) [ω]_E_ (x). (3.1) This gauge action (pointwise orthogonal transformations of local frames) is unavoidable in computations on tangent/cotangent bundles. A learning architecture intended to represent intrinsic operators on 11 -forms should respect this transformation law. We therefore adopt the following axiom. Definition 3.1 (Gauge equivariance). An operator :Ω1(ℳ)→Ω1(ℳ)T: ^1(M)→ ^1(M) is gauge-equivariant if for any choice of local orthonormal frames Eℓ\E_ \ and any smooth Rℓ:Uℓ→O(d)R_ :U_ → O(d) defining new frames Eℓ′=EℓRℓE _ =E_ R_ , the local representations satisfy [(ω)]Eℓ′(x)=Rℓ(x)⊤[(ω)]Eℓ(x)whenever[ω]Eℓ′(x)=Rℓ(x)⊤[ω]Eℓ(x),[T(ω)]_E _ (x)=R_ (x) [T(ω)]_E_ (x) [ω]_E _ (x)=R_ (x) [ω]_E_ (x), (3.2) for all x∈Uℓx∈ U_ . Equivalently, T commutes with the pointwise action of O(d)O(d) on the fiber coordinates. In words, changing the computational gauge (local orthonormal frame) should only rotate the coordinate representation of the output accordingly; the operator itself is independent of this choice. 3.2 Intrinsic Linear Layers via Functional Calculus of ΔH _H A central design choice is to build linear operator layers from intrinsic differential operators. Since ΔH _H is self-adjoint and nonnegative on L2Ω1L^2 ^1 , it admits a spectral decomposition ΔHψk=λkψk _H _k= _k _k with an L2L^2 -orthonormal eigenbasis ψkk≥0\ _k\_k≥ 0 . For any bounded measurable function m:[0,∞)→ℝm:[0,∞) , the functional calculus defines a linear operator m(ΔH)m( _H) by m(ΔH)ω:=∑k≥0m(λk)⟨ω,ψk⟩L2ψk,m( _H)ω:= _k≥ 0m( _k) ω, _k _L^2\, _k, (3.3) with convergence in L2Ω1L^2 ^1 (and in HsΩ1H^s ^1 under appropriate growth/regularity of m ). This construction is coordinate-free and depends only on (ℳ,g)(M,g) . We now define the basic learnable intrinsic linear layer as a truncated multiplier. Definition 3.2 (Truncated spectral multiplier layer). Fix a truncation index K∈ℕK and let mθ:[0,∞)→ℝm_θ:[0,∞) be a parameterized multiplier (e.g., an MLP in the scalar variable λ ). The spectral multiplier layer is (θ,Kω):=∑k=0Kmθ(λk)⟨ω,ψk⟩L2ψk. (T_θ,Kω ):= _k=0^Km_θ( _k) ω, _k _L^2\, _k. (3.4) The truncation serves two purposes: (i) it yields a finite-dimensional parameterization compatible with discretizations, and (i) it isolates a controllable approximation error that can be quantified in Sobolev operator norms. A complementary and often numerically convenient representation uses heat kernels. If m is completely monotone or admits a Laplace transform representation m(λ)=∫0∞a(t)e−tλtm(λ)= _0^∞a(t)e^-tλ\,dt , then m(ΔH)ω(x)=∫0∞a(t)(e−tΔHω)(x)t,m( _H)ω(x)= _0^∞a(t)\,(e^-t _Hω)(x)\,dt, (3.5) where e−tΔHe^-t _H is the heat semigroup on 11 -forms. This emphasizes that our convolution is intrinsic: it is induced by the geometry through the heat flow rather than by Euclidean translation. 3.3 Gauge-Equivariant Pointwise Nonlinearities on 11 -Forms To construct expressive nonlinear operator architectures, we interleave intrinsic linear layers with pointwise nonlinearities. However, generic coordinate-wise nonlinearities (e.g., applying a scalar activation to each fiber coordinate in an arbitrary gauge) break gauge equivariance. We therefore restrict to nonlinearities that are equivariant under the O(d)O(d) action in (3.1). A standard characterization is that continuous O(d)O(d) -equivariant maps φ:ℝd→ℝd :R^d ^d are radial: φ(v)=ρ(‖v‖)vfor some scalar functionρ:[0,∞)→ℝ. (v)=ρ(\|v\|)\,v some scalar function ρ:[0,∞) . (3.6) Accordingly, we define the fiberwise nonlinearity using intrinsic norms. Definition 3.3 (Fiberwise gauge-equivariant nonlinearity). Let ρθ:[0,∞)→ℝ _θ:[0,∞) be a learnable scalar function. Define σθ:Ω1(ℳ)→Ω1(ℳ) _θ: ^1(M)→ ^1(M) by (σθ(ω))(x):=ρθ(‖ω(x)‖g)ω(x),( _θ(ω))(x):= _θ (\|ω(x)\|_g )\,ω(x), (3.7) where ‖ω(x)‖g2=⟨ω(x),ω(x)⟩g\|ω(x)\|_g^2= ω(x),ω(x) _g . Then σθ _θ is intrinsic and gauge-equivariant. In practice, ρθ _θ can be parameterized by a small MLP, a spline, or a polynomial basis with bounded range (useful for stability). Importantly, (3.7) depends only on the norm induced by g , hence it is invariant to local orthonormal frame rotations. We also allow an intrinsic gated pointwise linear map that remains gauge-equivariant: (ℬθω)(x):=bθ(x)ω(x),(B_θω)(x):=b_θ(x)\,ω(x), (3.8) where bθ:ℳ→ℝb_θ:M is a learned scalar field (potentially produced by another intrinsic operator acting on invariants such as ‖ω‖g\|ω\|_g , curvature scalars, or known metadata). Since multiplication by a scalar commutes with O(d)O(d) , (3.8) preserves gauge equivariance. 3.4 The Gauge-Equivariant Intrinsic Neural Operator Architecture We now assemble the above ingredients into a deep operator network. Let ω0:=f _0:=f denote the input forcing 11 -form. For ℓ=0,…,L−1 =0,…,L-1 , define the hidden state recursion ωℓ+1:=σθℓ(θℓ,Kωℓ+ℬθℓωℓ), _ +1:= _ _ (T_ _ ,K _ +B_ _ _ ), (3.9) and output ^θ,K(f):=ωL S_θ,K(f):= _L . Here: • θℓ,KT_ _ ,K is the truncated intrinsic multiplier layer (3.4), • ℬθℓB_ _ is a gauge-compatible pointwise linear map (3.8), and • σθℓ _ _ is the gauge-equivariant nonlinearity (3.7). This architecture can be interpreted as learning a nonlinear perturbation of functional calculi of ΔH _H while remaining intrinsic. In the elliptic setting, one may also use purely linear versions (setting σ to identity) to directly approximate resolvents m(λ)≈(λ+α)−1m(λ)≈(λ+α)^-1 . We keep the nonlinear form (3.9) because it allows modeling operator families with parameter dependence, mild nonlinearity in constitutive relations, or data/model mismatch, while preserving the geometric symmetries. 3.5 Intrinsicness and Gauge Equivariance: Formal Guarantees We record the two basic structural properties required later for stability and discretization analysis. Proposition 3.1 (Intrinsicness). Let (ℳ,g)(M,g) be a compact smooth Riemannian manifold without boundary. Let φ:ℳ→ℳ :M be a diffeomorphism and set g′:=φ∗g := ^*g . Consider a network operator ^θ,K(g) S^(g)_θ,K built by composing layers of the form ω↦σθ(θ,K(g)ω+ℬθ(g)ω),ω _θ (T^(g)_θ,Kω+B^(g)_θω ), where: (1) Spectral layer θ,K(g)T^(g)_θ,K is a truncated functional calculus of the Hodge Laplacian ΔH,g _H,g on 11 -forms: θ,K(g):=mθ(ΔH,g)Π≤Λ(g),T^(g)_θ,K:=m_θ( _H,g)\, ^(g)_≤ , (3.10) for some cutoff Λ=Λ(K,g) = (K,g) , and where Π≤Λ(g) ^(g)_≤ is the L2L^2 -orthogonal spectral projector of ΔH,g _H,g onto the direct sum of eigenspaces with eigenvalues ≤Λ≤ . (This covers top- K modes truncation up to multiplicity; it is equivalent to truncating by the first K eigenmodes when eigenvalues are ordered.) (2) Pointwise linear term ℬθ(g)B^(g)_θ is multiplication by a scalar field bθ(g)∈C∞(ℳ)b^(g)_θ∈ C^∞(M) : (ℬθ(g)ω)(x)=bθ(g)(x)ω(x),(B^(g)_θω)(x)=b^(g)_θ(x)\,ω(x), (3.11) and the scalar field is natural under pullback: bθ(g′):=φ∗bθ(g).b^(g )_θ:= ^*b^(g)_θ. (3.12) (3) Nonlinearity σθ _θ is fiberwise radial with respect to the metric: (σθ(g)(ω))(x)=ρθ(‖ω(x)‖g)ω(x),(σ^(g)_θ(ω))(x)= _θ(\|ω(x)\|_g)\,ω(x), (3.13) where ρθ:[0,∞)→ℝ _θ:[0,∞) is fixed and |⋅|g|·|_g is the pointwise norm induced by g on 11 -forms. Then the operator is intrinsic in the sense that ^θ,K(g′)(φ∗f)=φ∗(^θ,K(g)(f))for all f∈Ω1(ℳ). S^(g )_θ,K( ^*f)= ^* ( S^(g)_θ,K(f) ) all f∈ ^1(M). (3.14) Proof We first establish the key intertwining identity: φ∗(ΔH,gω)=ΔH,g′(φ∗ω)∀ω∈Ω1(ℳ). ^*( _H,g\,ω)= _H,g ( ^*ω) ∀\,ω∈ ^1(M). (3.15) Recall ΔH,g=dδg+δgd _H,g=d _g+ _gd on 11 -forms, where d is the exterior derivative (metric-independent) and δg _g is the codifferential (metric-dependent). Two standard identities are: Pullback commutes with d : φ∗(dη)=d(φ∗η)∀η∈Ωk(ℳ). ^*(dη)=d( ^*η) ∀\,η∈ ^k(M). (3.16) Pullback intertwines Hodge star with pulled-back metric: φ∗(∗gη)=∗g′(φ∗η)∀η∈Ωk(ℳ). ^*(*_gη)=*_g ( ^*η) ∀\,η∈ ^k(M). (3.17) This follows from the defining property of ∗* : for any α,β∈Ωkα,β∈ ^k , α∧∗gβ=⟨α,β⟩gdvolg,α *_gβ= α,β _g\,dvol_g, and the facts that φ∗(α∧γ)=φ∗α∧φ∗γ ^*(α γ)= ^*α ^*γ , and that φ∗⟨α,β⟩g=⟨φ∗α,φ∗β⟩g′ ^* α,β _g= ^*α, ^*β _g , φ∗(dvolg)=dvolg′ ^*(dvol_g)=dvol_g . Using δg=(−1)d(k+1)+1∗gd∗g _g=(-1)^d(k+1)+1*_gd*_g on (k+1)(k+1) -forms, combine 3.16–3.17 to obtain, for any η∈Ωk+1η∈ ^k+1 , φ∗(δgη) ^*( _gη) =cφ∗(∗gd∗gη) =c\, ^* (*_gd*_gη ) =c∗g′φ∗(d∗gη) =c\,*_g ^* (d*_gη ) =c∗g′dφ∗(∗gη) =c\,*_g d\, ^* (*_gη ) =c∗g′d∗g′(φ∗η)=δg′(φ∗η), =c\,*_g d\,*_g ( ^*η)= _g ( ^*η), (3.18) where c=(−1)d(k+1)+1c=(-1)^d(k+1)+1 . Thus φ∗δg=δg′φ∗ ^* _g= _g ^* . Finally, φ∗(ΔH,gω) ^*( _H,gω) =φ∗(dδgω+δgdω) = ^*(d _gω+ _gdω) =dδg′(φ∗ω)+δg′d(φ∗ω)=ΔH,g′(φ∗ω), =d\, _g ( ^*ω)+ _g \,d( ^*ω)= _H,g ( ^*ω), (3.19) which proves 3.15. Let Ag:=ΔH,gA_g:= _H,g viewed as a self-adjoint nonnegative operator on L2Ω1L^2 ^1 . Define the pullback operator Uφ:L2Ω1(ℳ,g)→L2Ω1(ℳ,g′)U_ :L^2 ^1(M,g)→ L^2 ^1(M,g ) by Uφω:=φ∗ωU_ ω:= ^*ω . One checks that UφU_ is unitary because the L2L^2 inner product on forms is preserved when both the metric and volume form are pulled back. From 3.15 we have the operator intertwining relation (on smooth forms, hence by density on the operator domain): UφAg=Ag′Uφ.U_ A_g=A_g U_ . (3.20) A standard consequence of 3.20 for self-adjoint operators is that the entire spectral calculus intertwines: UφF(Ag)=F(Ag′)Uφfor all bounded Borel functions F.U_ \,F(A_g)=F(A_g )\,U_ all bounded Borel functions F. (3.21) In particular, taking F=mθF=m_θ yields φ∗(mθ(ΔH,g)ω)=mθ(ΔH,g′)(φ∗ω). ^* (m_θ( _H,g)\,ω )=m_θ( _H,g )\,( ^*ω). (3.22) Likewise, taking F=[0,Λ]F=1_[0, ] gives the spectral projector identity φ∗(Π≤Λ(g)ω)=Π≤Λ(g′)(φ∗ω). ^* ( ^(g)_≤ ω )= ^(g )_≤ ( ^*ω). (3.23) Combining 3.22 and 3.23, and using that mθ(ΔH,g)m_θ( _H,g) commutes with Π≤Λ(g) ^(g)_≤ (both are functions of the same self-adjoint operator), we obtain the key commutation for the truncated spectral layer 3.10: φ∗(θ,K(g)ω)=θ,K(g′)(φ∗ω). ^* (T^(g)_θ,Kω )=T^(g )_θ,K( ^*ω). (3.24) Pointwise linear term. Using 3.11 and the naturality assumption 3.12, φ∗(ℬθ(g)ω) ^* (B^(g)_θω ) =φ∗(bθ(g)ω)=(φ∗bθ(g))(φ∗ω)=bθ(g′)(φ∗ω)=ℬθ(g′)(φ∗ω). = ^* (b^(g)_θ\,ω )=( ^*b^(g)_θ)\,( ^*ω)=b^(g )_θ\,( ^*ω)=B^(g )_θ( ^*ω). (3.25) Radial nonlinearity. We need the pointwise norm compatibility: ‖φ∗ω‖g′=φ∗(‖ω‖g).\| ^*ω\|_g = ^*(\|ω\|_g). (3.26) This is a direct computation from the definition of the pullback metric g′=φ∗g = ^*g and the pullback of a 11 -form. In particular, the pointwise inner product on 11 -forms induced by g is preserved under pullback when the metric is pulled back, hence the norm transforms as 3.26. Therefore, φ∗(σθ(g)(ω)) ^* (σ^(g)_θ(ω) ) =φ∗(ρθ(‖ω‖g)ω) = ^* ( _θ(\|ω\|_g)\,ω ) =ρθ(φ∗‖ω‖g)(φ∗ω) = _θ( ^*\|ω\|_g)\,( ^*ω) =ρθ(‖φ∗ω‖g′)(φ∗ω)=σθ(g′)(φ∗ω). = _θ(\| ^*ω\|_g )\,( ^*ω)=σ^(g )_θ( ^*ω). (3.27) Define the layer map (for fixed ℓ ) as ℱθℓ,K(g)(ω):=σθℓ(g)(θℓ,K(g)ω+ℬθℓ(g)ω).F^(g)_ _ ,K(ω):=σ^(g)_ _ (T^(g)_ _ ,Kω+B^(g)_ _ ω ). Using 3.24, 3.25, and 3.5, φ∗(ℱθℓ,K(g)(ω))=ℱθℓ,K(g′)(φ∗ω). ^* (F^(g)_ _ ,K(ω) )=F^(g )_ _ ,K( ^*ω). (3.28) Let ω0=f _0=f and ωℓ+1=ℱθℓ,K(g)(ωℓ) _ +1=F^(g)_ _ ,K( _ ) , and similarly ω0′=φ∗fω _0= ^*f , ωℓ+1′=ℱθℓ,K(g′)(ωℓ′)ω _ +1=F^(g )_ _ ,K(ω _ ) . We prove by induction that ωℓ′=φ∗ωℓω _ = ^* _ for all ℓ . The base case ℓ=0 =0 holds by definition. If ωℓ′=φ∗ωℓω _ = ^* _ , then ωℓ+1′=ℱθℓ,K(g′)(ωℓ′)=ℱθℓ,K(g′)(φ∗ωℓ)=φ∗(ℱθℓ,K(g)(ωℓ))=φ∗ωℓ+1,ω _ +1=F^(g )_ _ ,K(ω _ )=F^(g )_ _ ,K( ^* _ )= ^* (F^(g)_ _ ,K( _ ) )= ^* _ +1, where we used 3.28. Thus ωL′=φ∗ωLω _L= ^* _L , i.e., ^θ,K(g′)(φ∗f)=φ∗(^θ,K(g)(f)), S^(g )_θ,K( ^*f)= ^* ( S^(g)_θ,K(f) ), Proposition 3.2 (Gauge equivariance). Let (ℳ,g)(M,g) be a compact Riemannian manifold. Fix an open set U⊂ℳU and two smooth local orthonormal frames E=(e1,…,ed)E=(e_1,…,e_d) and E′=(e1′,…,ed′)E =(e _1,…,e _d) on U . Assume they are related by a smooth map R:U→O(d)R:U→ O(d) via ej′(x)=∑k=1dek(x)Rkj(x)(x∈U).e _j(x)= _k=1^de_k(x)\,R_kj(x) (x∈ U). (3.29) Let [ω]E(x)∈ℝd[ω]_E(x) ^d denote the coordinate representation of a 11 -form ω in frame E , i.e., [ω]E(x)j:=ω(x)(ej(x)),[ω]E′(x)j:=ω(x)(ej′(x)).[ω]_E(x)_j:=ω(x) (e_j(x) ), [ω]_E (x)_j:=ω(x) (e _j(x) ). (3.30) Consider the layer map ℱ(ω):=σθ(θ,Kω+ℬθω),F(ω):= _θ (T_θ,Kω+B_θω ), (3.31) where: • θ,K:Ω1(ℳ)→Ω1(ℳ)T_θ,K: ^1(M)→ ^1(M) is a (truncated) spectral multiplier of ΔH _H as in (3.4), hence a well-defined operator on 11 -forms; • ℬθω=bθωB_θω=b_θ\,ω for a scalar field bθb_θ ; • σθ(ω)(x)=ρθ(‖ω(x)‖g)ω(x) _θ(ω)(x)= _θ(\|ω(x)\|_g)\,ω(x) for a scalar function ρθ _θ . Then for every ω∈Ω1(ℳ)ω∈ ^1(M) and every x∈Ux∈ U , [ℱ(ω)]E′(x)=R(x)⊤[ℱ(ω)]E(x).[F(ω)]_E (x)=R(x) \,[F(ω)]_E(x). (3.32) Consequently, the depth- L network ^θ,K S_θ,K obtained by composing such layers satisfies the gauge-equivariance condition in Definition 3.1. Proof Fix ω∈Ω1(ℳ)ω∈ ^1(M) and x∈Ux∈ U . Using 3.29 and linearity of ω(x):Txℳ→ℝω(x):T_xM , [ω]E′(x)j [ω]_E (x)_j =ω(x)(ej′(x))=ω(x)(∑k=1dek(x)Rkj(x)) =ω(x)(e _j(x))=ω(x)\! ( _k=1^de_k(x)R_kj(x) ) =∑k=1dRkj(x)ω(x)(ek(x))=∑k=1dRkj(x)[ω]E(x)k. = _k=1^dR_kj(x)\,ω(x)(e_k(x))= _k=1^dR_kj(x)\,[ω]_E(x)_k. (3.33) In vector form, [ω]E′(x)=R(x)⊤[ω]E(x).[ω]_E (x)=R(x) [ω]_E(x). (3.1.1) This identity holds for any 11 -form, independently of how the form is produced. Spectral multiplier θ,KT_θ,K . Since θ,Kω∈Ω1(ℳ)T_θ,Kω∈ ^1(M) is itself a 11 -form (by definition of the spectral expansion), we can apply(3.1) directly with η=θ,Kωη=T_θ,Kω : [θ,Kω]E′(x)=R(x)⊤[θ,Kω]E(x).[T_θ,Kω]_E (x)=R(x) [T_θ,Kω]_E(x). (3.34) No additional structure is needed: this follows purely because θ,KT_θ,K maps forms to forms in an intrinsic manner (i.e., it does not depend on the chosen frame). Scalar multiplication ℬθω=bθωB_θω=b_θ\,ω . For any x∈Ux∈ U and j∈1,…,dj∈\1,…,d\ , [(ℬθω)]E(x)j=(bθω)(x)(ej(x))=bθ(x)ω(x)(ej(x))=bθ(x)[ω]E(x)j.[(B_θω)]_E(x)_j=(b_θω)(x)(e_j(x))=b_θ(x)\,ω(x)(e_j(x))=b_θ(x)\,[ω]_E(x)_j. Thus [ℬθω]E(x)=bθ(x)[ω]E(x)[B_θω]_E(x)=b_θ(x)[ω]_E(x) , and similarly [ℬθω]E′(x)=bθ(x)[ω]E′(x)[B_θω]_E (x)=b_θ(x)[ω]_E (x) . Using (3.1), [ℬθω]E′(x) [B_θω]_E (x) =bθ(x)[ω]E′(x)=bθ(x)R(x)⊤[ω]E(x) =b_θ(x)\,[ω]_E (x)=b_θ(x)\,R(x) [ω]_E(x) =R(x)⊤(bθ(x)[ω]E(x))=R(x)⊤[ℬθω]E(x). =R(x) (b_θ(x)[ω]_E(x) )=R(x) [B_θω]_E(x). (3.35) Radial nonlinearity σθ(ω)(x)=ρθ(‖ω(x)‖g)ω(x) _θ(ω)(x)= _θ(\|ω(x)\|_g)\,ω(x) . Because E and E′E are orthonormal, the coordinate vector of a 11 -form has Euclidean norm equal to the metric norm: ‖ω(x)‖g2=∑j=1d[ω]E(x)j2=‖[ω]E(x)‖22,‖ω(x)‖g2=‖[ω]E′(x)‖22.\|ω(x)\|_g^2= _j=1^d[ω]_E(x)_j^2=\|[ω]_E(x)\|_2^2, \|ω(x)\|_g^2=\|[ω]_E (x)\|_2^2. (3.36) Moreover, since R(x)∈O(d)R(x)∈ O(d) , ‖[ω]E′(x)‖2=‖R(x)⊤[ω]E(x)‖2=‖[ω]E(x)‖2.\|[ω]_E (x)\|_2=\|R(x) [ω]_E(x)\|_2=\|[ω]_E(x)\|_2. (3.37) Now compute in frame E : [σθ(ω)]E(x)=ρθ(‖ω(x)‖g)[ω]E(x).[ _θ(ω)]_E(x)= _θ(\|ω(x)\|_g)\,[ω]_E(x). (3.38) In frame E′E : [σθ(ω)]E′(x) [ _θ(ω)]_E (x) =ρθ(‖ω(x)‖g)[ω]E′(x) = _θ(\|ω(x)\|_g)\,[ω]_E (x) =ρθ(‖ω(x)‖g)R(x)⊤[ω]E(x) = _θ(\|ω(x)\|_g)\,R(x) [ω]_E(x) =R(x)⊤(ρθ(‖ω(x)‖g)[ω]E(x)) =R(x) ( _θ(\|ω(x)\|_g)\,[ω]_E(x) ) =R(x)⊤[σθ(ω)]E(x), =R(x) [ _θ(ω)]_E(x), (3.39) which proves gauge equivariance of σθ _θ . Let η:=θ,Kω+ℬθωη:=T_θ,Kω+B_θω . By linearity of coordinate representation and 3.34–3.5, [η]E′(x)=R(x)⊤[η]E(x).[η]_E (x)=R(x) [η]_E(x). (3.40) Applying 3.5 to σθ(η) _θ(η) yields [ℱ(ω)]E′(x)=[σθ(η)]E′(x)=R(x)⊤[σθ(η)]E(x)=R(x)⊤[ℱ(ω)]E(x),[F(ω)]_E (x)=[ _θ(η)]_E (x)=R(x) [ _θ(η)]_E(x)=R(x) [F(ω)]_E(x), which is exactly 3.32. Let the network be defined recursively by ω0=f _0=f and ωℓ+1=ℱℓ(ωℓ) _ +1=F_ ( _ ) for ℓ=0,…,L−1 =0,…,L-1 , where each ℱℓF_ has the form 3.31 (with its own parameters). That for each ℓ , [ωℓ+1]E′(x)=R(x)⊤[ωℓ+1]E(x)whenever[ωℓ]E′(x)=R(x)⊤[ωℓ]E(x).[ _ +1]_E (x)=R(x) [ _ +1]_E(x) [ _ ]_E (x)=R(x) [ _ ]_E(x). By induction on ℓ , this holds for all layers, and in particular for the output ωL=^θ,K(f) _L= S_θ,K(f) . Hence the full network is gauge-equivariant. 3.6 Discrete Realization and Design for Consistency While the discretization analysis is deferred to Section 6, we emphasize a design principle: the architecture should admit a discrete counterpart obtained by replacing ΔH _H with a convergent discrete Hodge Laplacian ΔH,h _H,h (e.g., from discrete exterior calculus or compatible finite elements), and replacing ⟨⋅,⋅⟩L2 ·,· _L^2 with the associated discrete inner product. Under standard spectral convergence assumptions, the truncated multiplier layer (3.4) admits a discrete implementation (θ,K(h)ωh):=∑k=0Kmθ(λk,h)⟨ωh,ψk,h⟩hψk,h, (T_θ,K^(h) _h ):= _k=0^Km_θ( _k,h) _h, _k,h _h\, _k,h, (3.41) which mirrors the continuum definition and is the key to establishing continuum–discrete commutation bounds for the learned operator. 3.7 Discussion: Relation to the Target Elliptic Solution Operator The shifted Hodge–Poisson solution operator gS_g is itself a spectral multiplier m⋆(λ)=(λ+α)−1m_ (λ)=(λ+α)^-1 . Consequently, when the network is specialized to the linear case (identity σ , ℬ≡0B≡ 0 ), the model class contains direct approximants of gS_g via learned multipliers mθm_θ on a truncated spectrum. The nonlinear extension (3.9) retains intrinsicness and gauge equivariance while enlarging expressivity to cover (i) parameterized families of elliptic operators (e.g., variable coefficients), (i) mild nonlinearity arising from modeling error or constitutive laws, and (i) practical constraints where only partial spectral information is used and residual corrections are beneficial. In the next sections we quantify approximation error in Sobolev operator norms (Section4), establish stability under metric perturbations (Section5), and prove discretization consistency for the discrete realization(3.41) (Section6). 4 Sobolev Operator Approximation This section establishes quantitative approximation bounds for the elliptic solution operator g=(ΔH+αI)−1S_g=( _H+α I)^-1 (4.1) by the gauge-equivariant intrinsic neural operator class introduced in Section 3. The analysis is carried out in Sobolev topologies on 11-forms and is stated in operator norms between Sobolev spaces. A key point is that uniform operator-norm convergence on the full space (Hs−1Ω1→Hs+1Ω1H^s-1 ^1→ H^s+1 ^1) is impossible for purely spectral truncations; convergence holds on smoother input classes (compact embeddings), which is the relevant regime for PDE solution operators and for data distributions with controlled regularity. 4.1 Sobolev Norms via the Hodge Spectrum Let ψkk≥0⊂Ω1(ℳ)\ _k\_k≥ 0⊂ ^1(M) be an L2L^2-orthonormal eigenbasis of the Hodge Laplacian ΔH _H on 11-forms, with eigenvalues 0≤λ0≤λ1≤⋯0≤ _0≤ _1≤·s, λk→∞ _k→∞: ΔHψk=λkψk. _H _k= _k _k. (4.2) For r∈ℝr , we use the standard spectral characterization of Sobolev norms on 11-forms: ‖ω‖HrΩ12:=∑k≥0(1+λk)r|⟨ω,ψk⟩L2|2,\|ω\|_H^r ^1^2:= _k≥ 0(1+ _k)^r | ω, _k _L^2 |^2, (4.3) which is equivalent to the usual chart-based definition on compact manifolds. In particular, for any bounded Borel function m:[0,∞)→ℝm:[0,∞) , the operator m(ΔH)m( _H) satisfies ⟨m(ΔH)ω,ψk⟩L2=m(λk)⟨ω,ψk⟩L2. m( _H)ω, _k _L^2=m( _k) ω, _k _L^2. (4.4) We also record the smoothing structure of the shifted resolvent m⋆(λ)=(λ+α)−1m_ (λ)=(λ+α)^-1: g=m⋆(ΔH),m⋆(λ)=1λ+α.S_g=m_ ( _H), m_ (λ)= 1λ+α. (4.5) 4.2 Truncated Functional Calculus and the Approximation Family Because implementations (and learning) typically operate on a finite spectral subspace, we introduce spectral cutoffs in a basis-independent manner. For Λ>0 >0, let Π≤Λ _≤ denote the L2L^2-orthogonal spectral projector onto the direct sum of eigenspaces with eigenvalues ≤Λ≤ : Π≤Λω:=∑λk≤Λ⟨ω,ψk⟩L2ψk. _≤ ω:= _ _k≤ ω, _k _L^2 _k. (4.6) Given a multiplier mθm_θ, define the truncated intrinsic operator ^θ,Λ:=mθ(ΔH)Π≤Λ. S_θ, :=m_θ( _H) _≤ . (4.7) This covers the linear specialization of the architecture in Section 3 (with identity nonlinearity and ℬ≡0B≡ 0), and it is the core object for Sobolev approximation of the elliptic resolvent. We decompose the approximation error into (i) spectral truncation bias and (i) multiplier approximation error: g−^θ,Λ=g(I−Π≤Λ)⏟truncation bias+(m⋆(ΔH)−mθ(ΔH))Π≤Λ⏟multiplier error.S_g- S_θ, = S_g(I- _≤ )_truncation bias+ (m_ ( _H)-m_θ( _H) ) _≤ _multiplier error. (4.8) 4.3 Mapping Properties of the Elliptic Resolvent We first state the basic boundedness of gS_g between Sobolev spaces (a standard consequence of elliptic regularity, here recorded in the spectral norm framework). Lemma 4.1 (Sobolev boundedness of gS_g). Let (ℳ,g)(M,g) be compact and without boundary, and let α>0α>0. For any s∈ℝs , the resolvent g:=(ΔH+αI)−1S_g:=( _H+α I)^-1 extends uniquely to a bounded linear operator g:Hs−1Ω1⟶Hs+1Ω1,S_g:H^s-1 ^1 H^s+1 ^1, and there exists a constant Cα<∞C_α<∞ (depending only on α) such that ‖gf‖Hs+1≤Cα‖f‖Hs−1∀f∈Hs−1Ω1.\|S_gf\|_H^s+1≤ C_α\|f\|_H^s-1 ∀ f∈ H^s-1 ^1. (4.9) Proof Since ℳM is compact and ΔH _H is an elliptic, self-adjoint, nonnegative operator on L2Ω1L^2 ^1, its spectrum is discrete and there exists an L2L^2-orthonormal basis ψkk≥0⊂Ω1(ℳ)\ _k\_k≥ 0⊂ ^1(M) with eigenvalues 0≤λ0≤λ1≤⋯0≤ _0≤ _1≤·s, λk→∞ _k→∞, such that ΔHψk=λkψk. _H _k= _k _k. Every f∈L2Ω1f∈ L^2 ^1 admits an expansion f=∑k≥0fkψkf= _k≥ 0f_k _k with fk:=⟨f,ψk⟩L2f_k:= f, _k _L^2, converging in L2Ω1L^2 ^1. We define Sobolev norms via the spectral calculus (equivalent to the standard chart-based definition on compact manifolds): ‖f‖HrΩ12:=∑k≥0(1+λk)r|fk|2,r∈ℝ.\|f\|_H^r ^1^2:= _k≥ 0(1+ _k)^r|f_k|^2, r . (4.10) This yields a Hilbert space HrΩ1H^r ^1 as the completion of Ω1(ℳ) ^1(M) under ∥⋅∥Hr\|·\|_H^r. For smooth f∈Ω1(ℳ)f∈ ^1(M) (hence in L2Ω1L^2 ^1), define gf:=∑k≥01λk+αfkψk.S_gf:= _k≥ 0 1 _k+αf_k _k. (4.11) This series converges in L2Ω1L^2 ^1 because (λk+α)−1( _k+α)^-1 is bounded by α−1α^-1 and ∑|fk|2<∞Σ|f_k|^2<∞. Let f∈Ω1(ℳ)f∈ ^1(M) and write f=∑fkψkf=Σ f_k _k. Then by 4.10 and 4.11, ‖gf‖Hs+12=∑k≥0(1+λk)s+1|1λk+α|2|fk|2.\|S_gf\|_H^s+1^2= _k≥ 0(1+ _k)^s+1 | 1 _k+α |^2|f_k|^2. (4.12) Factor the weight as (1+λk)s+11(λk+α)2=((1+λk)2(λk+α)2)(1+λk)s−1.(1+ _k)^s+1 1( _k+α)^2= ( (1+ _k)^2( _k+α)^2 )(1+ _k)^s-1. (4.13) Substituting 4.13 into 4.12 gives ‖gf‖Hs+12=∑k≥0((1+λk)2(λk+α)2)(1+λk)s−1|fk|2≤(supλ≥0(1+λ)2(λ+α)2)‖f‖Hs−12.\|S_gf\|_H^s+1^2= _k≥ 0 ( (1+ _k)^2( _k+α)^2 )(1+ _k)^s-1|f_k|^2≤ ( _λ≥ 0 (1+λ)^2(λ+α)^2 )\|f\|_H^s-1^2. (4.14) Define Cα:=(supλ≥01+λ+α).C_α:= ( _λ≥ 0 1+λ+α ). (4.15) Then 4.14 yields ‖gf‖Hs+1≤Cα‖f‖Hs−1∀f∈Ω1(ℳ).\|S_gf\|_H^s+1≤ C_α\|f\|_H^s-1 ∀ f∈ ^1(M). (4.16) It remains to check that Cα<∞C_α<∞. Consider r(λ):=1+λ+αr(λ):= 1+λ+α for λ≥0λ≥ 0. This is continuous on [0,∞)[0,∞), and limλ→∞r(λ)=1,r(0)=1α. _λ→∞r(λ)=1, r(0)= 1α. Moreover, r′(λ)=α−1(λ+α)2,r (λ)= α-1(λ+α)^2, so r is monotone increasing if α>1α>1, monotone decreasing if α<1α<1, and constant if α=1α=1. Hence supλ≥0r(λ)=max1,1α<∞, _λ≥ 0r(λ)= \1, 1α \<∞, so Cα=max1,α−1C_α= \1,α^-1 \ is a valid choice (and any larger constant also works). Therefore 4.16 holds with finite CαC_α. The space Ω1(ℳ) ^1(M) is dense in Hs−1Ω1H^s-1 ^1 by construction. Let f∈Hs−1Ω1f∈ H^s-1 ^1 and choose a sequence f(n)∈Ω1(ℳ)\f^(n)\∈ ^1(M) such that f(n)→f^(n)→ f in Hs−1H^s-1. Then by 4.16, ‖gf(n)−gf(m)‖Hs+1≤Cα‖f(n)−f(m)‖Hs−1,\|S_gf^(n)-S_gf^(m)\|_H^s+1≤ C_α\|f^(n)-f^(m)\|_H^s-1, so gf(n)\S_gf^(n)\ is Cauchy in Hs+1Ω1H^s+1 ^1 and converges to a limit we define as gfS_gf. The bound 4.9 follows by taking limits: ‖gf‖Hs+1=limn→∞‖gf(n)‖Hs+1≤limn→∞Cα‖f(n)‖Hs−1=Cα‖f‖Hs−1.\|S_gf\|_H^s+1= _n→∞\|S_gf^(n)\|_H^s+1≤ _n→∞C_α\|f^(n)\|_H^s-1=C_α\|f\|_H^s-1. Uniqueness of the extension follows because any two bounded linear extensions coincide on the dense subspace Ω1(ℳ) ^1(M). 4.4 Spectral Truncation Bias on Smoother Input Classes A key subtlety is that g(I−Π≤Λ)S_g(I- _≤ ) does not converge to zero in operator norm (Hs−1→Hs+1H^s-1→ H^s+1) as Λ→∞ →∞; one can concentrate f on arbitrarily high eigenmodes, for which truncation discards non-negligible components. Convergence holds on smoother input spaces Hs−1+γH^s-1+γ with γ>0γ>0, which is the natural regime for compactness and for many data models. Lemma 4.2 (Truncation bias in Sobolev operator norm). Let s∈ℝs and γ>0γ>0. For Λ>0 >0, let Π≤Λ _≤ be the L2L^2-orthogonal spectral projector of ΔH _H on 11-forms onto the span of eigenspaces with eigenvalues ≤Λ≤ . Then ‖g(I−Π≤Λ)‖Hs−1+γ⟶Hs+1≤Cα(1+Λ)−γ/2, \|S_g(I- _≤ ) \|_H^s-1+γ H^s+1≤ C_α(1+ )^-γ/2, (4.17) where g=(ΔH+αI)−1S_g=( _H+α I)^-1 with α>0α>0, and Cα<∞C_α<∞ depends only on α (and the choice of equivalent Sobolev norm convention). Proof Let ψkk≥0\ _k\_k≥ 0 be an L2L^2-orthonormal eigenbasis of ΔH _H on 11-forms, ΔHψk=λkψk _H _k= _k _k with 0≤λk↑∞0≤ _k \!\! ∞. For f∈Hs−1+γΩ1f∈ H^s-1+γ ^1, write f=∑k≥0fkψk,fk:=⟨f,ψk⟩L2.f= _k≥ 0f_k _k, f_k:= f, _k _L^2. By definition of Π≤Λ _≤ , (I−Π≤Λ)f=∑λk>Λfkψk.(I- _≤ )f= _ _k> f_k _k. Using the spectral representation of gS_g, g(I−Π≤Λ)f=∑λk>Λ1λk+αfkψk.S_g(I- _≤ )f= _ _k> 1 _k+αf_k _k. Using the spectral Sobolev norm convention ‖u‖Hr2=∑k≥0(1+λk)r|uk|2,\|u\|_H^r^2= _k≥ 0(1+ _k)^r|u_k|^2, we have ‖g(I−Π≤Λ)f‖Hs+12=∑λk>Λ(1+λk)s+11(λk+α)2|fk|2. \|S_g(I- _≤ )f \|_H^s+1^2= _ _k> (1+ _k)^s+1 1( _k+α)^2|f_k|^2. (4.18) Rewrite the weight in 4.18 as (1+λk)s+11(λk+α)2=((1+λk)2(λk+α)2)(1+λk)s−1.(1+ _k)^s+1 1( _k+α)^2= ( (1+ _k)^2( _k+α)^2 )(1+ _k)^s-1. Insert an additional factor (1+λk)−γ(1+λk)γ(1+ _k)^-γ(1+ _k)^γ to match the input norm: (1+λk)s−1=(1+λk)s−1+γ(1+λk)−γ.(1+ _k)^s-1=(1+ _k)^s-1+γ(1+ _k)^-γ. Therefore, (1+λk)s+11(λk+α)2=((1+λk)2(λk+α)2)(1+λk)−γ(1+λk)s−1+γ.(1+ _k)^s+1 1( _k+α)^2= ( (1+ _k)^2( _k+α)^2 )(1+ _k)^-γ(1+ _k)^s-1+γ. (4.19) Combine 4.18–4.19 to obtain ‖g(I−Π≤Λ)f‖Hs+12=∑λk>Λ((1+λk)2(λk+α)2)(1+λk)−γ(1+λk)s−1+γ|fk|2. \|S_g(I- _≤ )f \|_H^s+1^2= _ _k> ( (1+ _k)^2( _k+α)^2 )(1+ _k)^-γ(1+ _k)^s-1+γ|f_k|^2. Hence, ‖g(I−Π≤Λ)f‖Hs+12≤(supλ≥0(1+λ)2(λ+α)2)(supλ>Λ(1+λ)−γ)∑λk>Λ(1+λk)s−1+γ|fk|2. \|S_g(I- _≤ )f \|_H^s+1^2≤ ( _λ≥ 0 (1+λ)^2(λ+α)^2 ) ( _λ> (1+λ)^-γ ) _ _k> (1+ _k)^s-1+γ|f_k|^2. (4.20) Since γ>0γ>0, (1+λ)−γ(1+λ)^-γ is decreasing, so supλ>Λ(1+λ)−γ=(1+Λ)−γ. _λ> (1+λ)^-γ=(1+ )^-γ. (4.21) Moreover, the first supremum is finite because α>0α>0 and the rational function 1+λ+α 1+λ+α is bounded on [0,∞)[0,∞); explicitly, supλ≥0(1+λ)2(λ+α)2=(supλ≥01+λ+α)2≤max1,α−12. _λ≥ 0 (1+λ)^2(λ+α)^2= ( _λ≥ 0 1+λ+α )^2≤ \1,α^-1 \^2. (4.22) Finally, the sum over λk>Λ\ _k> \ is bounded by the full Sobolev norm: ∑λk>Λ(1+λk)s−1+γ|fk|2≤∑k≥0(1+λk)s−1+γ|fk|2=‖f‖Hs−1+γ2. _ _k> (1+ _k)^s-1+γ|f_k|^2≤ _k≥ 0(1+ _k)^s-1+γ|f_k|^2=\|f\|_H^s-1+γ^2. (4.23) Substituting 4.21–4.22 into 4.20 yields ‖g(I−Π≤Λ)f‖Hs+12≤Cα2(1+Λ)−γ‖f‖Hs−1+γ2, \|S_g(I- _≤ )f \|_H^s+1^2≤ C_α^2(1+ )^-γ\|f\|_H^s-1+γ^2, where Cα:=supλ≥01+λ+α≤max1,α−1C_α:= _λ≥ 0 1+λ+α≤ \1,α^-1 \. Taking square roots gives ‖g(I−Π≤Λ)f‖Hs+1≤Cα(1+Λ)−γ/2‖f‖Hs−1+γ. \|S_g(I- _≤ )f \|_H^s+1≤ C_α(1+ )^-γ/2\|f\|_H^s-1+γ. Since this holds for all f∈Hs−1+γΩ1f∈ H^s-1+γ ^1, we conclude ‖g(I−Π≤Λ)‖Hs−1+γ⟶Hs+1≤Cα(1+Λ)−γ/2. \|S_g(I- _≤ ) \|_H^s-1+γ H^s+1≤ C_α(1+ )^-γ/2. 4.5 Multiplier Approximation Error on the Truncated Subspace We now quantify the effect of approximating the true multiplier m⋆(λ)=(λ+α)−1m_ (λ)=(λ+α)^-1 on [0,Λ][0, ]. Define the uniform approximation error εΛ(θ):=sup0≤λ≤Λ|mθ(λ)−m⋆(λ)|. _ (θ):= _0≤λ≤ |m_θ(λ)-m_ (λ) |. (4.24) Lemma 4.3 (Uniform multiplier error implies Sobolev operator bound). Let s∈ℝs and Λ>0 >0. Let m⋆(λ)=(λ+α)−1m_ (λ)=(λ+α)^-1 and let mθ:[0,∞)→ℝm_θ:[0,∞) be any bounded function. Define εΛ(θ):=sup0≤λ≤Λ|mθ(λ)−m⋆(λ)|. _ (θ):= _0≤λ≤ |m_θ(λ)-m_ (λ)|. Let Π≤Λ _≤ be the L2L^2-orthogonal spectral projector of ΔH _H onto the sum of eigenspaces with eigenvalues ≤Λ≤ . Then ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λ‖Hs−1⟶Hs+1≤(1+Λ)εΛ(θ). \| (m_ ( _H)-m_θ( _H) ) _≤ \|_H^s-1 H^s+1≤(1+ ) _ (θ). (4.25) Proof Let ψkk≥0\ _k\_k≥ 0 be an L2L^2-orthonormal eigenbasis of ΔH _H on 11-forms, with eigenvalues λkk≥0\ _k\_k≥ 0. For f∈Hs−1Ω1f∈ H^s-1 ^1, write f=∑k≥0fkψk,fk:=⟨f,ψk⟩L2.f= _k≥ 0f_k _k, f_k:= f, _k _L^2. By definition of Π≤Λ _≤ , Π≤Λf=∑λk≤Λfkψk. _≤ f= _ _k≤ f_k _k. Since m(ΔH)m( _H) acts diagonally in this eigenbasis, we have (m⋆(ΔH)−mθ(ΔH))Π≤Λf=∑λk≤Λ(m⋆(λk)−mθ(λk))fkψk. (m_ ( _H)-m_θ( _H) ) _≤ f= _ _k≤ (m_ ( _k)-m_θ( _k) )f_k _k. (4.26) Using the spectral definition of Sobolev norms, ‖u‖Hs+12=∑k≥0(1+λk)s+1|uk|2,\|u\|_H^s+1^2= _k≥ 0(1+ _k)^s+1|u_k|^2, the output in 4.26 yields ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λf‖Hs+12 \| (m_ ( _H)-m_θ( _H) ) _≤ f \|_H^s+1^2 =∑λk≤Λ(1+λk)s+1|m⋆(λk)−mθ(λk)|2|fk|2. = _ _k≤ (1+ _k)^s+1 |m_ ( _k)-m_θ( _k) |^2|f_k|^2. (4.27) By definition of εΛ(θ) _ (θ), for all k with λk≤Λ _k≤ , |m⋆(λk)−mθ(λk)|≤εΛ(θ).|m_ ( _k)-m_θ( _k)|≤ _ (θ). (4.28) Substitute 4.28 into 4.27: ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λf‖Hs+12≤εΛ(θ)2∑λk≤Λ(1+λk)s+1|fk|2. \| (m_ ( _H)-m_θ( _H) ) _≤ f \|_H^s+1^2≤ _ (θ)^2 _ _k≤ (1+ _k)^s+1|f_k|^2. (4.29) Now note that for λk≤Λ _k≤ , (1+λk)s+1=(1+λk)2(1+λk)s−1≤(1+Λ)2(1+λk)s−1.(1+ _k)^s+1=(1+ _k)^2(1+ _k)^s-1≤(1+ )^2(1+ _k)^s-1. (4.30) Applying 4.30 to 4.29 gives ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λf‖Hs+12≤(1+Λ)2εΛ(θ)2∑λk≤Λ(1+λk)s−1|fk|2. \| (m_ ( _H)-m_θ( _H) ) _≤ f \|_H^s+1^2≤(1+ )^2 _ (θ)^2 _ _k≤ (1+ _k)^s-1|f_k|^2. (4.31) Since the sum over λk≤Λ\ _k≤ \ is dominated by the full sum, ∑λk≤Λ(1+λk)s−1|fk|2≤∑k≥0(1+λk)s−1|fk|2=‖f‖Hs−12. _ _k≤ (1+ _k)^s-1|f_k|^2≤ _k≥ 0(1+ _k)^s-1|f_k|^2=\|f\|_H^s-1^2. (4.32) Combining 4.31 and 4.32 yields ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λf‖Hs+12≤(1+Λ)2εΛ(θ)2‖f‖Hs−12. \| (m_ ( _H)-m_θ( _H) ) _≤ f \|_H^s+1^2≤(1+ )^2 _ (θ)^2\|f\|_H^s-1^2. Taking square roots, |(m⋆(ΔH)−mθ(ΔH))Π≤Λf|Hs+1≤(1+Λ)εΛ(θ)‖f‖Hs−1. | (m_ ( _H)-m_θ( _H) ) _≤ f |_H^s+1≤(1+ ) _ (θ)\|f\|_H^s-1. Since this holds for all f∈Hs−1Ω1f∈ H^s-1 ^1, we conclude ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λ‖Hs−1⟶Hs+1≤(1+Λ)εΛ(θ), \| (m_ ( _H)-m_θ( _H) ) _≤ \|_H^s-1 H^s+1≤(1+ ) _ (θ), which is exactly the desired bound. 4.6 Main Approximation Theorem (Bias–Approximation Decomposition) We now combine the preceding bounds using the decomposition 4.8. Theorem 4.4 (Sobolev operator approximation of the elliptic solution operator). Let s∈ℝs , γ>0γ>0, and Λ>0 >0. Let g:=(ΔH+αI)−1=m⋆(ΔH),m⋆(λ)=1λ+α.S_g:=( _H+α I)^-1=m_ ( _H), m_ (λ)= 1λ+α. Let Π≤Λ _≤ be the L2L^2-orthogonal spectral projector onto the direct sum of eigenspaces with eigenvalues ≤Λ≤ . For any parameterized multiplier mθm_θ, define the truncated intrinsic approximation ^θ,Λ:=mθ(ΔH)Π≤Λ. S_θ, :=m_θ( _H) _≤ . Define the uniform multiplier error on [0,Λ][0, ] by εΛ(θ):=sup0≤λ≤Λ|mθ(λ)−m⋆(λ)|. _ (θ):= _0≤λ≤ |m_θ(λ)-m_ (λ)|. Then there exists Cα<∞C_α<∞ depending only on α such that ‖g−^θ,Λ‖Hs−1+γ⟶Hs+1≤Cα(1+Λ)−γ/2+(1+Λ)εΛ(θ).\|S_g- S_θ, \|_H^s-1+γ H^s+1≤ C_α(1+ )^-γ/2+(1+ ) _ (θ). (4.33) Proof Using the definitions and the fact that Π≤Λ _≤ is a projection (Π≤Λ+(I−Π≤Λ)=I _≤ +(I- _≤ )=I), we write g−^θ,Λ=g(I−Π≤Λ)+(g−mθ(ΔH))Π≤Λ.S_g- S_θ, =S_g(I- _≤ )+ (S_g-m_θ( _H) ) _≤ . (4.34) Since g=m⋆(ΔH)S_g=m_ ( _H), the second term becomes (g−mθ(ΔH))Π≤Λ=(m⋆(ΔH)−mθ(ΔH))Π≤Λ. (S_g-m_θ( _H) ) _≤ = (m_ ( _H)-m_θ( _H) ) _≤ . (4.35) Thus, g−^θ,Λ=g(I−Π≤Λ)⏟AΛ+(m⋆(ΔH)−mθ(ΔH))Π≤Λ⏟Bθ,Λ.S_g- S_θ, = S_g(I- _≤ )_A_ + (m_ ( _H)-m_θ( _H) ) _≤ _B_θ, . (4.36) Let X:=Hs−1+γΩ1X:=H^s-1+γ ^1 and Y:=Hs+1Ω1Y:=H^s+1 ^1. For bounded linear operators AΛ,Bθ,Λ:X⟶YA_ ,B_θ, :X Y, ‖AΛ+Bθ,Λ‖X⟶Y≤‖AΛ‖X⟶Y+‖Bθ,Λ‖X⟶Y.\|A_ +B_θ, \|_X Y≤\|A_ \|_X Y+\|B_θ, \|_X Y. (4.37) Applying 4.37 to 4.36 yields ‖g−^θ,Λ‖Hs−1+γ⟶Hs+1≤‖g(I−Π≤Λ)‖Hs−1+γ⟶Hs+1+‖(m⋆(ΔH)−mθ(ΔH))Π≤Λ‖Hs−1+γ⟶Hs+1.\|S_g- S_θ, \|_H^s-1+γ H^s+1≤\|S_g(I- _≤ )\|_H^s-1+γ H^s+1+\| (m_ ( _H)-m_θ( _H) ) _≤ \|_H^s-1+γ H^s+1. (4.38) By Lemma 4.2, ‖g(I−Π≤Λ)‖Hs−1+γ⟶Hs+1≤Cα(1+Λ)−γ/2.\|S_g(I- _≤ )\|_H^s-1+γ H^s+1≤ C_α(1+ )^-γ/2. (4.39) Lemma 4.3 gives the bound ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λ‖Hs−1⟶Hs+1≤(1+Λ)εΛ(θ).\| (m_ ( _H)-m_θ( _H) ) _≤ \|_H^s-1 H^s+1≤(1+ ) _ (θ). (4.40) We need the operator norm on the smaller domain Hs−1+γH^s-1+γ rather than Hs−1H^s-1. Since Hs−1+γ↪Hs−1H^s-1+γ H^s-1 continuously, and with the spectral convention 4.3 one has the contractive embedding ‖f‖Hs−1≤‖f‖Hs−1+γ∀f∈Hs−1+γΩ1,\|f\|_H^s-1≤\|f\|_H^s-1+γ ∀ f∈ H^s-1+γ ^1, (4.41) because (1+λk)s−1≤(1+λk)s−1+γ(1+ _k)^s-1≤(1+ _k)^s-1+γ for all k. Therefore, for any bounded linear operator T:Hs−1⟶Hs+1T:H^s-1 H^s+1, ‖T‖Hs−1+γ⟶Hs+1=supf≠0‖Tf‖Hs+1‖f‖Hs−1+γ≤supf≠0‖Tf‖Hs+1‖f‖Hs−1=‖T‖Hs−1⟶Hs+1.\|T\|_H^s-1+γ H^s+1= _f≠ 0 \|Tf\|_H^s+1\|f\|_H^s-1+γ≤ _f≠ 0 \|Tf\|_H^s+1\|f\|_H^s-1=\|T\|_H^s-1 H^s+1. (4.42) Apply 4.42 to T=Bθ,Λ=(m⋆(ΔH)−mθ(ΔH))Π≤ΛT=B_θ, = (m_ ( _H)-m_θ( _H) ) _≤ and combine with 4.41 to obtain ‖(m⋆(ΔH)−mθ(ΔH))Π≤Λ‖Hs−1+γ⟶Hs+1≤(1+Λ)εΛ(θ).\| (m_ ( _H)-m_θ( _H) ) _≤ \|_H^s-1+γ H^s+1≤(1+ ) _ (θ). (4.43) Substitute 4.39 and 4.43 into 4.38: ‖g−^θ,Λ‖Hs−1+γ⟶Hs+1≤Cα(1+Λ)−γ/2+(1+Λ)εΛ(θ),\|S_g- S_θ, \|_H^s-1+γ H^s+1≤ C_α(1+ )^-γ/2+(1+ ) _ (θ), which is exactly the desired bound. 4.7 Discussion and Consequences (i) Interpretation of the two terms. The right-hand side of 4.33 is the sum of: • a truncation (bias) term decaying as (1+Λ)−γ/2(1+ )^-γ/2, which reflects that higher input regularity (γ>0γ>0) suppresses high-frequency components; • a multiplier approximation term growing as (1+Λ)εΛ(θ)(1+ ) _ (θ), which quantifies the difficulty of approximating the resolvent multiplier uniformly up to frequency Λ . (i) Impossibility of operator-norm convergence on Hs−1→Hs+1H^s-1→ H^s+1 via truncation. For completeness, note that for the hard truncation operator gΠ≤ΛS_g _≤ one has |g(I−Π≤Λ)|Hs−1→Hs+1↛0, |S_g(I- _≤ ) |_H^s-1→ H^s+1 → 0, since choosing f=ψkf= _k with λk>Λ _k> yields gf=(λk+α)−1ψkS_gf=( _k+α)^-1 _k while gΠ≤Λf=0S_g _≤ f=0, and the ratio |(λk+α)−1ψk|Hs+1/|ψk|Hs−1 |( _k+α)^-1 _k |_H^s+1/ | _k |_H^s-1 stays bounded away from zero as k→∞k→∞. This justifies working with smoother input classes Hs−1+γH^s-1+γ, which are compactly embedded into Hs−1H^s-1 on compact manifolds and correspond to realistic forcing distributions. (i) From frequency cutoff to mode count. If one prefers a truncation level K (number of retained modes) rather than Λ , a Weyl-type relation λK≍K2/d _K K^2/d (for fixed (ℳ,g)(M,g)) converts 4.33 into rates in K. This translation is optional for the main theory; the Λ -formulation is basis- and multiplicity-robust. (iv) Connection to the full nonlinear architecture. Theorem 4.4 is stated for the linear truncated functional calculus 4.7, which is the natural approximation class for the elliptic resolvent. The same proof strategy extends to residual/nonlinear variants provided each layer is Lipschitz on Hs+1H^s+1 and remains a function of intrinsic operators; the resulting bounds accumulate multiplicatively/additively with depth via standard stability estimates. We focus here on the linear case because it isolates the elliptic operator learning mechanism and yields the sharpest Sobolev approximation statement for (ΔH+αI)−1( _H+α I)^-1. This completes the Sobolev operator approximation result. Section 5 will quantify stability of both the true resolvent and the learned operator family under metric perturbations, and Section 6 will establish discretization consistency for convergent discrete Hodge Laplacians. 5 Stability under Metric Perturbations and Gauge Changes This section quantifies how both the true elliptic solution operator and its learned intrinsic approximation vary under perturbations of the underlying geometry. The guiding principle is that the shifted Hodge–Poisson resolvent g=(ΔH,g+αI)−1S_g=( _H,g+α I)^-1 depends continuously on the metric g in appropriate topologies, and the intrinsic neural operator family inherits a comparable stability provided its learned multipliers and nonlinearities are uniformly controlled. We separate (i) metric stability (changes in g) from (i) gauge stability (changes of local orthonormal frames). The latter is exact at the continuum level due to gauge equivariance, while the former yields quantitative bounds. 5.1 Metric Perturbation Model and Notation Let (ℳ,g)(M,g) be as in Section 2. We consider a second metric g~ g on ℳM satisfying uniform equivalence and regularity bounds. Fix an integer r≥2r≥ 2 and assume: • Uniform equivalence: there exists κ≥1κ≥ 1 such that for all x∈ℳx and v∈Txℳv∈ T_xM, κ−1gx(v,v)≤g~x(v,v)≤κgx(v,v).κ^-1g_x(v,v)≤ g_x(v,v)≤κ g_x(v,v). (5.1) • Regularity: g,g~∈Cr(ℳ)g, g∈ C^r(M) and ‖g‖Cr+‖g~‖Cr≤M\|g\|_C^r+\| g\|_C^r≤ M for some M<∞M<∞. • Small perturbation: ‖g~−g‖Cr≤ε\| g-g\|_C^r≤ , where ε is sufficiently small (depending on M,κM,κ). We denote by ΔH,g _H,g and ΔH,g~ _H, g the Hodge Laplacians on 11-forms induced by g and g~ g, and define the associated shifted operators ℒg:=ΔH,g+αI,ℒg~:=ΔH,g~+αI,α>0.L_g:= _H,g+α I, _ g:= _H, g+α I, α>0. (5.2) Their inverses g=ℒg−1S_g=L_g^-1 and g~=ℒg~−1S_ g=L_ g^-1 are bounded maps Hs−1Ω1⟶Hs+1Ω1H^s-1 ^1 H^s+1 ^1 for every s∈ℝs (Lemma 4.1). 5.2 Stability of the True Elliptic Solution Operator We first control the difference g−g~S_g-S_ g in Sobolev operator norms. The argument follows a resolvent identity coupled with bounds on the metric dependence of ΔH _H. 5.2.1 A resolvent identity For boundedly invertible operators ℒg,ℒg~L_g,L_ g on a common Banach space, one has g−g~=ℒg−1−ℒg~−1=ℒg~−1(ℒg~−ℒg)ℒg−1.S_g-S_ g=L_g^-1-L_ g^-1=L_ g^-1(L_ g-L_g)L_g^-1. (5.3) Thus, stability reduces to estimating ℒg~−ℒg=ΔH,g~−ΔH,gL_ g-L_g= _H, g- _H,g as a map between Sobolev spaces. 5.2.2 Metric dependence of the Hodge Laplacian In local coordinates, ΔH _H is a second-order elliptic differential operator on 11-forms with coefficients depending smoothly on g and finitely many derivatives of g. Under the CrC^r assumption with r≥2r≥ 2, the coefficient difference satisfies (schematically) ΔH,g~−ΔH,g=∑|β|≤2aβ(x)∂β,with‖aβ‖Cr−|β|≲‖g~−g‖Cr. _H, g- _H,g= _|β|≤ 2a_β(x)∂^β, \|a_β\|_C^r-|β| \| g-g\|_C^r. (5.4) Consequently, for each s∈ℝs with s≤r−1s≤ r-1, the operator difference extends to a bounded map ΔH,g~−ΔH,g:Hs+1Ω1⟶Hs−1Ω1,and‖ΔH,g~−ΔH,g‖Hs+1⟶Hs−1≤C‖g~−g‖Cr, _H, g- _H,g:H^s+1 ^1 H^s-1 ^1, \| _H, g- _H,g\|_H^s+1 H^s-1≤ C\| g-g\|_C^r, (5.5) where C depends on (ℳ,d)(M,d) and the a priori bounds M,κM,κ. 5.2.3 Operator-norm stability bound Combining 5.3, Lemma 4.1 applied to g and g~ g, and 5.5, yields the main stability estimate. Theorem 5.1 (Metric stability of the shifted Hodge–Poisson resolvent). Let (ℳ,g)(M,g) be a compact smooth Riemannian manifold without boundary, let α>0α>0, and let s∈ℝs . Fix an integer r≥2r≥ 2 such that r≥s+2r≥ s+2. Let g~ g be another CrC^r Riemannian metric on ℳM satisfying the uniform equivalence condition: there exists κ≥1κ≥ 1 such that κ−1gx(v,v)≤g~x(v,v)≤κgx(v,v)∀x∈ℳ,∀v∈Txℳ.κ^-1g_x(v,v)≤ g_x(v,v)≤κ g_x(v,v) ∀ x ,\ ∀ v∈ T_xM. (5.6) Assume moreover that ‖g‖Cr+‖g~‖Cr≤M\|g\|_C^r+\| g\|_C^r≤ M for some M<∞M<∞. Define ℒg:=ΔH,g+αI,g:=ℒg−1,L_g:= _H,g+α I, _g:=L_g^-1, and similarly for g~ g. Then there exists a constant C<∞C<∞, depending only on α, κ, M, ℳM, and r, such that ‖g−g~‖Hs−1⟶Hs+1≤C‖g−g~‖Cr.\|S_g-S_ g\|_H^s-1 H^s+1\ ≤\ C\|g- g\|_C^r. (5.7) Proof Since α>0α>0, both ℒgL_g and ℒg~L_ g are injective and have bounded inverses on Sobolev scales (Lemma 4.1 applied to each metric). The algebraic identity A−1−B−1=A−1(B−A)B−1A^-1-B^-1=A^-1(B-A)B^-1 with A=ℒg~A=L_ g and B=ℒgB=L_g yields g−g~=ℒg−1−ℒg~−1=g~(ℒg~−ℒg)g=g~(ΔH,g~−ΔH,g)g.S_g-S_ g=L_g^-1-L_ g^-1=S_ g(L_ g-L_g)S_g=S_ g( _H, g- _H,g)S_g. (5.8) Therefore, ‖g−g~‖Hs−1⟶Hs+1≤‖g~‖Hs−1⟶Hs+1⋅‖ΔH,g~−ΔH,g‖Hs+1⟶Hs−1⋅‖g‖Hs−1⟶Hs+1.\|S_g-S_ g\|_H^s-1 H^s+1≤\|S_ g\|_H^s-1 H^s+1·\| _H, g- _H,g\|_H^s+1 H^s-1·\|S_g\|_H^s-1 H^s+1. (5.9) Thus it suffices to (i) bound the two resolvent norms uniformly under 5.1, and (i) bound the Laplacian difference by C|g−g~|CrC |g- g |_C^r. Lemma 4.1 gives, for each fixed metric h, ‖h‖Hs−1⟶Hs+1≤Cα(h),\|S_h\|_H^s-1 H^s+1≤ C_α(h), (5.10) where the constant can be taken as Cα(h)=supλ≥01+λ+α≤max1,α−1.C_α(h)= _λ≥ 0 1+λ+α≤ \1,α^-1 \. The right-hand side is independent of h; the only subtlety is that our Sobolev norms HtΩ1H^t ^1 depend on the metric. However, under the uniform equivalence assumption 5.1 and the CrC^r bounds, the Sobolev norms defined using g and g~ g are equivalent with constants depending only on κ, M, ℳM, and r. Concretely, for each fixed t∈[−r,r]t∈[-r,r] there exists ct,Ct>0c_t,C_t>0 such that ct‖ω‖Ht(g~)≤‖ω‖Ht(g)≤Ct‖ω‖Ht(g~)∀ω.c_t\|ω\|_H^t( g)≤\|ω\|_H^t(g)≤ C_t\|ω\|_H^t( g) ∀ω. (5.11) (For compact manifolds this is standard: it follows by comparing the coordinate definitions of Sobolev norms via partitions of unity and using 5.1 plus boundedness of metric coefficients and derivatives up to order r.) Using 5.11 to move between g-Sobolev and g~ g-Sobolev norms, we may bound both resolvents by a common constant ‖g‖Hs−1⟶Hs+1+‖g~‖Hs−1⟶Hs+1≤Cα,κ,M,ℳ,r.\|S_g\|_H^s-1 H^s+1+\|S_ g\|_H^s-1 H^s+1≤ C_α,κ,M,M,r. (5.12) We will absorb these constants into a generic C. Fix a finite atlas (Uℓ,χℓ)ℓ=1L\(U_ , _ )\_ =1^L and a subordinate partition of unity ηℓ\ _ \. On each chart, the Hodge Laplacian on 11-forms is a second-order linear differential operator whose coefficients are smooth functions of the metric components and finitely many derivatives of those components. More precisely, in local coordinates one can write (ΔH,hω)i=∑|β|≤2ai,β(h)(x)∂βω(x)on Uℓ,( _H,hω)_i= _|β|≤ 2a^(h)_i,β(x)∂^βω(x) U_ , (5.13) where h∈g,g~h∈\g, g\, the index i∈1,…,di∈\1,…,d\ denotes the component of the 11-form, and each coefficient ai,β(h)a^(h)_i,β depends smoothly on hpq\h_pq\, ∂hpq\∂ h_pq\, and ∂2hpq\∂^2h_pq\ (equivalently, on h and its Christoffel symbols and curvature terms). This can be justified, for instance, by using the Weitzenböck identity on 11-forms, ΔH,h=∇h∗∇h+Rich, _H,h= _h^* _h+Ric_h, (5.14) and expanding ∇h _h in coordinates; the principal coefficients are hjkh^jk and lower-order coefficients involve Γ(h) (h) and ∂Γ(h)∂ (h), hence derivatives of h up to second order. Since g,g~∈Crg, g∈ C^r with ‖g‖Cr+‖g~‖Cr≤M\|g\|_C^r+\| g\|_C^r≤ M and r≥2r≥ 2, the coefficients satisfy ‖ai,β(g~)−ai,β(g)‖Cr−|β|(Uℓ)≤C‖g~−g‖Cr(Uℓ).\|a^( g)_i,β-a^(g)_i,β\|_C^r-|β|(U_ )\ ≤ C\| g-g\|_C^r(U_ ). (5.15) The constant C depends on M, κ, the chart, and r, but not on the particular g~,g g,g beyond those bounds. (This is a consequence of smooth dependence of the formulas for inverse metric, Christoffel symbols, and curvature on h and its derivatives, combined with the uniform ellipticity from 5.1.) Therefore, on each UℓU_ , (ΔH,g~−ΔH,g)ω=∑|β|≤2bi,β(x)∂βω,bi,β:=ai,β(g~)−ai,β(g),( _H, g- _H,g)ω= _|β|≤ 2b_i,β(x)∂^βω, b_i,β:=a^( g)_i,β-a^(g)_i,β, (5.16) with coefficient bounds 5.15. We now estimate the operator norm of the difference. Fix ℓ and localize ω by ηℓω _ ω. Since |β|≤2|β|≤ 2, ‖∂β(ηℓω)‖Hs−1(ℝd)≤C‖ηℓω‖Hs+1(ℝd).\|∂^β( _ ω)\|_H^s-1(R^d)≤ C\| _ ω\|_H^s+1(R^d). (5.17) Next we use a standard Sobolev multiplier fact on ℝdR^d: if b∈Cmb∈ C^m with integer m≥⌈s−1⌉m≥ s-1 , then multiplication by b is bounded on Hs−1H^s-1, and ‖bu‖Hs−1≤C‖b‖Cm‖u‖Hs−1.\|bu\|_H^s-1≤ C\|b\|_C^m\|u\|_H^s-1. (5.18) (One may prove 5.18 by Littlewood–Paley theory or by commutator estimates; on a compact manifold this is obtained chartwise and patched with partitions of unity. Our assumption r≥s+2r≥ s+2 ensures r−|β|≥s≥s−1r-|β|≥ s≥ s-1 in the integer case, and more generally is sufficient to apply multiplier bounds for real s.) Apply 5.18 with b=bi,βb=b_i,β and u=∂β(ηℓω)u=∂^β( _ ω), using ‖bi,β‖Cr−|β|\|b_i,β\|_C^r-|β| and r−|β|≥sr-|β|≥ s (since |β|≤2|β|≤ 2 and r≥s+2r≥ s+2). Combining with 5.17, ‖bi,β∂β(ηℓω)‖Hs−1≤C‖bi,β‖Cr−|β|‖ηℓω‖Hs+1.\|b_i,β∂^β( _ ω)\|_H^s-1≤ C\|b_i,β\|_C^r-|β|\| _ ω\|_H^s+1. (5.19) Summing over |β|≤2|β|≤ 2 and components i, and then summing over charts using the partition of unity, yields ‖(ΔH,g~−ΔH,g)ω‖Hs−1≤C(max|β|≤2‖bi,β‖Cr−|β|)‖ω‖Hs+1.\|( _H, g- _H,g)ω\|_H^s-1≤ C ( _|β|≤ 2\|b_i,β\|_C^r-|β| )\|ω\|_H^s+1. (5.20) Finally, apply 5.15 to control the coefficient differences by |g−g~|Cr |g- g |_C^r. This proves ‖ΔH,g~−ΔH,g‖Hs+1⟶Hs−1≤C‖g−g~‖Cr,\| _H, g- _H,g\|_H^s+1 H^s-1≤ C\|g- g\|_C^r, (5.21) with C depending only on κ, M, ℳM, and r. Insert 5.12 and 5.21 into 5.9: ‖g−g~‖Hs−1⟶Hs+1≤C‖g−g~‖Cr,\|S_g-S_ g\|_H^s-1 H^s+1≤ C\|g- g\|_C^r, for a constant C depending only on α, κ, M, ℳM, and r. This is exactly 5.7. Remark 5.2 (Role of the shift α). The shift α>0α>0 guarantees coercivity and removes the harmonic subspace from the kernel, ensuring the resolvent remains uniformly bounded and simplifying the perturbation analysis. Without the shift, one must work on the orthogonal complement of harmonic 1-forms and track stability of the corresponding projections. 5.3 Stability of Intrinsic Neural Approximations under Metric Perturbations We now consider the learned operator family ^θ,Λ(g) S_θ, ^(g) defined by truncated functional calculus: ^θ,Λ(g):=mθ(ΔH,g)Π≤Λ(g). S_θ, ^(g):=m_θ( _H,g) _≤ ^(g). (5.22) This is the principal approximation class for the elliptic resolvent, and it also forms the linear core of the nonlinear architecture in Section 3. Metric stability requires controlling how spectral projectors and functional calculi vary with g. We state a stability result in terms of multiplier regularity and spectral gaps. 5.3.1 Assumptions on multiplier regularity and spectral separation Assume mθ∈C1([0,Λ])m_θ∈ C^1([0, ]) and define Lθ:=sup0≤λ≤Λ|mθ′(λ)|L_θ:= _0≤λ≤ |m_θ (λ) | and Mθ:=sup0≤λ≤Λ|mθ(λ)|M_θ:= _0≤λ≤ |m_θ(λ) |. Assume further that Λ is not too close to the spectrum of ΔH,g _H,g so that the cutoff is stable: there exists δ>0δ>0 such that dist(Λ,spec(ΔH,g))≥δ,dist(Λ,spec(ΔH,g~))≥δ. dist ( , spec ( _H,g) )≥δ, dist ( , spec ( _H, g) )≥δ. (5.23) This is a standard spectral-gap condition ensuring that the projector Π≤Λ _≤ depends smoothly on the operator. 5.3.2 Stability of the learned operator Theorem 5.3 (Metric stability of truncated intrinsic neural operators). Let (ℳ,g)(M,g) be compact, without boundary. Fix s∈ℝs and an integer r≥s+2r≥ s+2. Let g~ g satisfy the perturbation model in Section 5.1, and set A:=ΔH,g,A~:=ΔH,g~.A:= _H,g, A:= _H, g. Fix Λ>0 >0 and assume the spectral separation condition: there exists δ>0δ>0 such that dist(Λ,spec(A))≥δ,dist(Λ,spec(A~))≥δ. dist ( , spec (A))≥δ, dist ( , spec ( A))≥δ. (5.24) Let mθm_θ be a real-valued multiplier on [0,Λ][0, ] with bounds Mθ:=sup0≤λ≤Λ|mθ(λ)|,Lθ:=sup0≤λ≤Λ|mθ′(λ)|.M_θ:= _0≤λ≤ |m_θ(λ) |, L_θ:= _0≤λ≤ |m_θ (λ) |. Assume in addition that mθm_θ admits a bounded complex extension mθ:→ℂm_θ:U on an open neighborhood ⊂ℂU of the contour Γ defined below, and that on Γ supz∈Γ|mθ(z)|≤CΓ(Mθ+Lθ), _z∈ |m_θ(z) |≤ C_ (M_θ+L_θ), (5.25) for a constant CΓC_ depending only on Λ,δ ,δ. Define the truncated intrinsic operator ^θ,Λ(g):=mθ(A)Π≤Λ(g),^θ,Λ(g~):=mθ(A~)Π≤Λ(g~). S^(g)_θ, :=m_θ(A) ^(g)_≤ , S^( g)_θ, :=m_θ( A) ^( g)_≤ . Then there exists C<∞C<∞ (depending on Λ,δ,α ,δ,α, ℳM, and the a priori metric bounds) such that ‖^θ,Λ(g)−^θ,Λ(g~)‖Hs−1⟶Hs+1≤C(Mθ+Lθ)‖g−g~‖Cr. \| S^(g)_θ, - S^( g)_θ, \|_H^s-1 H^s+1≤ C(M_θ+L_θ)\|g- g\|_C^r. (5.26) Proof Let Γ⊂ℂ be a positively oriented simple closed contour enclosing the real interval [0,Λ][0, ] and satisfying dist(Γ,spec(A))≥δ2,dist(Γ,spec(A~))≥δ2. dist ( , spec (A))≥ δ2, dist ( , spec ( A))≥ δ2. (5.27) Such a Γ exists by 5.24. For concreteness, one may take Γ to be the boundary of the closed δ/2δ/2-tubular neighborhood of [0,Λ][0, ] in ℂC. Denote its length by len(Γ) len ( ), which depends only on Λ and δ. For a self-adjoint operator A with discrete spectrum, the Riesz projection onto the spectral subset inside Γ is Π≤Λ(g)=12πi∫Γ(zI−A)−1z, _≤ ^(g)= 12π i _ (zI-A)^-1dz, (5.28) and similarly for A~ A. By construction, the spectrum enclosed by Γ is exactly spec(A)∩[0,Λ] spec (A)∩[0, ], hence Π≤Λ(g) _≤ ^(g) is the L2L^2-orthogonal projector onto the sum of eigenspaces with eigenvalues ≤Λ≤ . By the Dunford functional calculus (holomorphic functional calculus), for mθm_θ holomorphic on a neighborhood of Γ , mθ(A)Π≤Λ(g)=12πi∫Γmθ(z)(zI−A)−1z.m_θ(A) _≤ ^(g)= 12π i _ m_θ(z)(zI-A)^-1dz. (5.29) Indeed, the operator defined by the right-hand side equals the spectral multiplier mθm_θ applied to the part of the spectrum inside Γ and zero outside; since Γ encloses [0,Λ][0, ] and excludes the rest of the spectrum, this is precisely mθ(A)Π≤Λ(g)m_θ(A) _≤ ^(g). The same identity holds with A replaced by A~ A. Therefore, ^θ,Λ(g)−^θ,Λ(g~)=12πi∫Γmθ(z)[(zI−A)−1−(zI−A~)−1]z. S^(g)_θ, - S^( g)_θ, = 12π i _ m_θ(z) [(zI-A)^-1-(zI- A)^-1 ]dz. (5.30) For all z∈Γz∈ , both resolvents exist by 5.27. The resolvent identity gives (zI−A)−1−(zI−A~)−1=(zI−A~)−1(A~−A)(zI−A)−1.(zI-A)^-1-(zI- A)^-1=(zI- A)^-1( A-A)(zI-A)^-1. (5.31) We will estimate 5.31 as an operator from Hs−1Ω1H^s-1 ^1 to Hs+1Ω1H^s+1 ^1. Fix z∈Γz∈ . Using the eigenbasis ψk\ _k\ of A and the spectral Sobolev norms defined via A (equivalent to chart-based norms on a compact manifold), the resolvent acts diagonally: (zI−A)−1ψk=1z−λkψk.(zI-A)^-1 _k= 1z- _k _k. Hence for f=∑fkψkf=Σ f_k _k, ‖(zI−A)−1f‖Hs+12=∑k(1+λk|z−λk|)2(1+λk)s−1|fk|2.\|(zI-A)^-1f\|_H^s+1^2= _k ( 1+ _k|z- _k| )^2(1+ _k)^s-1|f_k|^2. (5.32) Therefore, ‖(zI−A)−1‖Hs−1⟶Hs+1≤supλ∈spec(A)1+λ|z−λ|.\|(zI-A)^-1\|_H^s-1 H^s+1≤ _λ∈ spec (A) 1+λ|z-λ|. (5.33) On Γ , |z−λ|≥δ/2|z-λ|≥δ/2 for all λ∈spec(A)λ∈ spec (A) by 5.27, and moreover λ↦1+λ|z−λ|λ 1+λ|z-λ| is bounded uniformly over λ≥0λ≥ 0 when z stays in a compact set away from [0,∞)[0,∞) by a fixed distance. Since Γ lies in a bounded neighborhood of [0,Λ][0, ], we may bound 1+λ1+λ on the enclosed region by 1+Λ+δ1+ +δ, while for large λ the distance |z−λ|∼λ|z-λ| λ yields a bounded ratio. Concretely, one obtains supλ≥01+λ|z−λ|≤CΛ,δ,z∈Γ, _λ≥ 0 1+λ|z-λ|≤ C_ ,δ, z∈ , (5.34) hence supz∈Γ‖(zI−A)−1‖Hs−1⟶Hs+1≤CΛ,δ. _z∈ \|(zI-A)^-1\|_H^s-1 H^s+1≤ C_ ,δ. (5.35) The same bound holds with A replaced by A~ A, with the same constant up to equivalence of Sobolev norms under the uniform metric bounds. By the coefficient-dependence estimate used in the proof of Theorem 5.1, since r≥s+2r≥ s+2, ‖A~−A‖Hs+1⟶Hs−1=‖ΔH,g~−ΔH,g‖Hs+1⟶Hs−1≤C‖g−g~‖Cr,\| A-A\|_H^s+1 H^s-1=\| _H, g- _H,g\|_H^s+1 H^s-1≤ C\|g- g\|_C^r, (5.36) where C depends only on ℳM and the a priori metric bounds. Combining 5.31, 5.35, and 5.36, for all z∈Γz∈ , ‖(zI−A)−1−(zI−A~)−1‖Hs−1⟶Hs+1≤‖(zI−A~)−1‖Hs−1⟶Hs+1⋅‖A~−A‖Hs+1⟶Hs−1⋅‖(zI−A)−1‖Hs−1⟶Hs+1≤CΛ,δ‖g−g~‖Cr.\|(zI-A)^-1-(zI- A)^-1\|_H^s-1 H^s+1≤\|(zI- A)^-1\|_H^s-1 H^s+1·\| A-A\|_H^s+1 H^s-1·\|(zI-A)^-1\|_H^s-1 H^s+1≤ C_ ,δ\|g- g\|_C^r. (5.37) (Here CΛ,δC_ ,δ absorbs the product of the two resolvent bounds and the constant in 5.36.) From 5.30, using the standard estimate for Bochner integrals of bounded operators, ‖^θ,Λ(g)−^θ,Λ(g~)‖Hs−1⟶Hs+1≤12π∫Γ|mθ(z)|‖(zI−A)−1−(zI−A~)−1‖Hs−1⟶Hs+1|dz|. \| S^(g)_θ, - S^( g)_θ, \|_H^s-1 H^s+1≤ 12π _ |m_θ(z) |\|(zI-A)^-1-(zI- A)^-1\|_H^s-1 H^s+1|dz|. (5.38) Apply 5.37 and the bound 5.25 on |mθ| |m_θ | over Γ : ‖^θ,Λ(g)−^θ,Λ(g~)‖Hs−1⟶Hs+1≤len(Γ)2πsupz∈Γ|mθ(z)|⋅CΛ,δ‖g−g~‖Cr≤C(Mθ+Lθ)‖g−g~‖Cr, \| S^(g)_θ, - S^( g)_θ, \|_H^s-1 H^s+1≤ len ( )2π _z∈ |m_θ(z) |· C_ ,δ\|g- g\|_C^r≤ C(M_θ+L_θ)\|g- g\|_C^r, (5.39) where C depends on Λ,δ ,δ, the contour choice, and the a priori metric bounds, but not on θ beyond (Mθ,Lθ)(M_θ,L_θ). Discussion. The factor Mθ+LθM_θ+L_θ expresses that stability improves when multipliers are uniformly bounded and not overly oscillatory on [0,Λ][0, ]. This motivates regularizing mθm_θ during training (e.g., spectral Lipschitz penalties). 5.3.3 Extension to the nonlinear architecture For the full network 3.9, metric stability follows by combining: • stability of each intrinsic linear layer mθℓ(ΔH,g)Π≤Λ(g)m_ _ ( _H,g) _≤ ^(g) as in 5.26, • stability of scalar pointwise multipliers bθℓ(g)b^(g)_ _ under metric pullback rules, • Lipschitz properties of the gauge-equivariant radial nonlinearities σθℓ(g) _ _ ^(g) in Sobolev norms (typically obtained by boundedness of ρθ _θ and ρθ′ _θ on the relevant range of norms). A representative bound takes the form ‖^θ,Λ(g)−^θ,Λ(g~)‖Hs−1⟶Hs+1≤Cθ‖g~−g‖Cr,\| S_θ, ^(g)- S_θ, ^( g)\|_H^s-1 H^s+1≤ C_θ\| g-g\|_C^r, (5.40) where CθC_θ grows at most polynomially/exponentially with depth depending on whether residual or contractive designs are used. In practice, enforcing uniform Lipschitz control per layer yields stable depth scaling. 5.4 Gauge Stability: Exact Equivariance at the Continuum Level We now formalize gauge stability in the sense relevant to computation: changing the local orthonormal frames used to represent 11-forms must not change the geometric output, only its coordinate representation. Let E and E′E be two local orthonormal frames on an open set U⊂ℳU , related by a smooth R:U⟶O(d)R:U O (d) as in 3.29. For any 11-form ω, its coordinate representation transforms by [ω]E′(x)=R(x)T[ω]E(x).[ω]_E (x)=R(x)^T[ω]_E(x). (5.41) Proposition 3.2 implies the following immediate corollary. Corollary 5.4 (Exact gauge equivariance implies representation stability). For the intrinsic neural operator ^θ,Λ S_θ, defined in Section 3, and for every x∈Ux∈ U, [^θ,Λ(ω)]E′(x)=R(x)T[^θ,Λ(ω)]E(x),[ S_θ, (ω)]_E (x)=R(x)^T[ S_θ, (ω)]_E(x), (5.42) i.e., the computed output transforms consistently under a change of computational gauge. Remark 5.5 (Discrete gauge error). At the discrete level (Section 6), exact gauge equivariance can be broken by discretization artifacts (e.g., approximate mass matrices, imperfect orthonormalization, or non-commuting interpolation). The continuum result 5.42 is therefore best interpreted as a target symmetry; Section 6 establishes that equivariance errors vanish under refinement under compatible discretizations. 5.5 Summary of Stability Guarantees • Metric perturbations: both gS_g and the learned intrinsic approximation vary Lipschitz-continuously with g in Sobolev operator norms, under standard regularity and spectral separation assumptions (Theorems 5.1 and 5.3). • Gauge changes: equivariance is exact at the continuum level by construction (Corollary 5.4), ensuring representation independence. These results complement Section 4 (approximation) by ensuring that the learned operator is not only accurate but also robust to geometric uncertainty and to representation choices, which is essential for cross-mesh and cross-geometry generalization. 6 Continuum–Discrete Consistency This section establishes that the intrinsic neural operators introduced in Section 3 admit discretizations on meshes that are consistent with the continuum model as the discretization is refined. In particular, we show that if the discrete Hodge Laplacian converges to its continuum counterpart in a standard sense, then the corresponding discrete truncated functional calculus converges to the continuum truncated functional calculus, and consequently the discrete neural operator converges to the continuum neural operator. We formulate the results for 1-forms using discrete exterior calculus (DEC) or compatible finite elements; the proofs require only (i) stability of discrete inner products and (i) spectral convergence of discrete Hodge Laplacians. 6.1 Discretization Model and Discrete Spaces Let (ℳ,g)(M,g) be as in Section 2. Let Thh↓0\T_h\_h 0 be a family of shape-regular triangulations of ℳM with mesh size parameter h. We assume the triangulations are quasi-uniform and geometrically consistent with ℳM (e.g., embedded triangulations with uniformly bounded aspect ratios and uniformly controlled geometry approximation). We discretize 11-forms by a finite-dimensional space Ωh1 _h^1 equipped with: • a discrete L2L^2 inner product ⟨⋅,⋅⟩h ·,· _h and norm ∥⋅∥h\|·\|_h, • a discrete Hodge Laplacian ΔH,h:Ωh1→Ωh1 _H,h: _h^1→ _h^1 that is symmetric positive semidefinite with respect to ⟨⋅,⋅⟩h ·,· _h, • a projection (interpolation) operator Πh:Ω1(ℳ)→Ωh1 _h: ^1(M)→ _h^1, and a reconstruction operator ℛh:Ωh1→Ω1(ℳ)R_h: _h^1→ ^1(M). Concrete instances include: • DEC 1-cochains with circumcentric dual and the DEC Hodge star, yielding ΔH,h=dhδh+δhdh _H,h=d_h _h+ _hd_h; • Whitney 1-forms / Nédélec-type elements with mass and stiffness matrices, yielding an equivalent discrete Hodge Laplacian. We assume the pair (Πh,ℛh)( _h,R_h) is stable in Sobolev norms and consistent on smooth forms. 6.2 Discrete Spectral Data and Discrete Truncation Let (λk,h,ψk,h)k=0Nh−1\( _k,h, _k,h)\_k=0^N_h-1 be an ⟨⋅,⋅⟩h ·,· _h-orthonormal eigenpair family for ΔH,h _H,h, ΔH,hψk,h=λk,hψk,h,⟨ψk,h,ψj,h⟩h=δkj. _H,h _k,h= _k,h _k,h, _k,h, _j,h _h= _kj. (6.1) For a cutoff Λ>0 >0, define the discrete spectral projector Π≤Λ(h)ωh:=∑λk,h≤Λ⟨ωh,ψk,h⟩hψk,h. _≤ ^(h) _h:= _ _k,h≤ _h, _k,h _h _k,h. (6.2) Given a multiplier mθm_θ, define the discrete truncated multiplier operator ^θ,Λ(h):=mθ(ΔH,h)Π≤Λ(h),wheremθ(ΔH,h)ωh:=∑kmθ(λk,h)⟨ωh,ψk,h⟩hψk,h. S_θ, ^(h):=m_θ( _H,h) _≤ ^(h), m_θ( _H,h) _h:= _km_θ( _k,h) _h, _k,h _h _k,h. (6.3) This is the exact discrete analogue of 4.7 and is the canonical discretization of the intrinsic spectral layer. 6.3 Assumptions: Stability and Spectral Convergence Our consistency theorem relies on standard approximation properties of the discrete Hodge Laplacian. We state them abstractly to cover DEC and compatible FEM. Assumption 6.1 (Uniform norm stability and quasi-optimal reconstruction). There exists C independent of h such that for all sufficiently smooth ω, ‖Πhω‖h≤C‖ω‖L2,‖ℛhωh‖L2≤C‖ωh‖h,\| _hω\|_h≤ C\|ω\|_L^2, \|R_h _h\|_L^2≤ C\| _h\|_h, (6.4) and ℛhΠhω→ωR_h _hω→ω in L2L^2 as h→0h→ 0. Assumption 6.2 (Spectral convergence below a fixed cutoff). Fix Λ>0 >0 that does not coincide with an eigenvalue of ΔH,g _H,g. Then there exists p>0p>0 and C such that for every continuum eigenpair (λk,ψk)( _k, _k) with λk<Λ _k< , there exists a discrete eigenpair (λk,h,ψk,h)( _k,h, _k,h) (with consistent indexing up to multiplicity) such that |λk,h−λk|≤Chp,‖ℛhψk,h−ψk‖L2≤Chp, | _k,h- _k |≤ Ch^p, \|R_h _k,h- _k\|_L^2≤ Ch^p, (6.5) for all sufficiently small h. Moreover, the discrete eigenspace associated with eigenvalues ≤Λ≤ converges to the corresponding continuum eigenspace in the gap sin(Θ) \!( ) metric. These assumptions are satisfied by standard compatible discretizations under regularity and shape-regularity conditions; they are widely used in numerical analysis of Hodge Laplacians on surfaces and manifolds. 6.4 Consistency of the Discrete Spectral Projector We first show that the discrete projector converges to the continuum projector under reconstruction. Lemma 6.3 (Projector consistency). Let (ℳ,g)(M,g) be compact without boundary. Fix Λ>0 >0 such that Λ∉spec(ΔH,g) ∉ spec ( _H,g). Let Π≤Λ _≤ be the L2L^2-orthogonal spectral projector of the Hodge Laplacian ΔH,g _H,g on 11-forms onto the direct sum of eigenspaces with eigenvalues ≤Λ≤ . Let ThT_h be a shape-regular family of discretizations with discrete space Ωh1 _h^1, discrete inner product ⟨⋅,⋅⟩h ·,· _h, discrete Hodge Laplacian ΔH,h _H,h, discrete projector Π≤Λ(h) ^(h)_≤ defined by 6.2, interpolation Πh _h, and reconstruction ℛhR_h. Assume: 1. (Approximation of reconstruction–interpolation) There exists p>0p>0 and s sufficiently large (depending on the method order) such that for all ω∈HsΩ1ω∈ H^s ^1, ‖ℛhΠhω−ω‖L2≤Chp‖ω‖Hs.\|R_h _hω-ω\|_L^2\ ≤\ Ch^p\|ω\|_H^s. (6.6) (For Whitney/compatible FEM/DEC on smooth manifolds, such an estimate holds with p equal to the approximation order under standard regularity assumptions.) 2. (Spectral subspace convergence below Λ ) Assumption 6.2 holds: the discrete eigenspace associated with eigenvalues ≤Λ≤ converges to the continuum eigenspace in the gap metric with order hph^p. Then there exists C independent of h such that for all ω∈HsΩ1ω∈ H^s ^1, ‖ℛhΠ≤Λ(h)Πhω−Π≤Λω‖L2≤Chp‖ω‖Hs. \|R_h ^(h)_≤ _hω- _≤ ω \|_L^2\ ≤\ Ch^p\|ω\|_H^s. (6.7) Proof Let E:=ran(Π≤Λ)⊂L2Ω1E:= ran ( _≤ )⊂ L^2 ^1 be the continuum low-frequency invariant subspace. Since Λ∉spec(ΔH,g) ∉ spec ( _H,g) and ℳM is compact, E is finite-dimensional; write m:=dim(E)m:= (E). Let Eh:=ran(Π≤Λ(h))⊂Ωh1E_h:= ran ( ^(h)_≤ )⊂ _h^1 be the discrete low-frequency eigenspace (with respect to ⟨⋅,⋅⟩h ·,· _h). Define its reconstructed counterpart in the continuum Hilbert space L2Ω1L^2 ^1 by E^h:=ℛh(Eh)⊂L2Ω1. E_h:=R_h(E_h)⊂ L^2 ^1. Let P:=Π≤ΛP:= _≤ be the L2L^2-orthogonal projector onto E, and let P^h P_h denote the L2L^2-orthogonal projector onto E^h E_h. By Assumption 6.2 (eigenspace convergence in the gap metric), for all sufficiently small h one has dim(E^h)=dim(Eh)=m ( E_h)= (E_h)=m, and there exists a constant C such that ‖P^h−P‖L2⟶L2≤Chp.\| P_h-P\|_L^2 L^2≤ Ch^p. (6.8) This is a standard equivalence: convergence of subspaces in the gap (or sin(Θ) ( )) metric is precisely operator-norm convergence of the corresponding orthogonal projectors. Define the continuum-to-continuum discrete projector (apply discrete projector to discretized input, then reconstruct): Phd:=ℛhΠ≤Λ(h)Πh:HsΩ1→L2Ω1.P_h^d:=R_h ^(h)_≤ _h:H^s ^1→ L^2 ^1. We aim to bound ‖Phdω−Pω‖L2\|P_h^dω-Pω\|_L^2. Insert P^h P_h and ℛhΠhR_h _h to obtain the decomposition Phdω−Pω=(Phdω−P^h(ℛhΠhω))⏟(I)+P^h(ℛhΠhω−ω)⏟(II)+(P^h−P)ω⏟(III).P_h^dω-Pω= (P_h^dω- P_h(R_h _hω) )_(I)+ P_h(R_h _hω-ω)_(I)+ ( P_h-P)ω_(I). (6.9) We will bound each term in L2L^2. Since P^h P_h is an orthogonal projector on L2L^2, ‖P^h‖L2⟶L2=1\| P_h\|_L^2 L^2=1. Therefore, by 6.6, ‖(II)‖L2≤‖ℛhΠhω−ω‖L2≤Chp‖ω‖Hs.\|(I)\|_L^2≤\|R_h _hω-ω\|_L^2≤ Ch^p\|ω\|_H^s. (6.10) By 6.8, ‖(III)‖L2≤‖P^h−P‖L2⟶L2‖ω‖L2≤Chp‖ω‖L2≤Chp‖ω‖Hs.\|(I)\|_L^2≤\| P_h-P\|_L^2 L^2\|ω\|_L^2≤ Ch^p\|ω\|_L^2≤ Ch^p\|ω\|_H^s. (6.11) By definition, Π≤Λ(h) _≤ ^(h) is the ⟨⋅,⋅⟩h ·,· _h-orthogonal projector onto EhE_h. Hence for any vh∈Ωh1v_h∈ _h^1, ‖vh−Π≤Λ(h)vh‖h=infwh∈Eh‖vh−wh‖h.\|v_h- _≤ ^(h)v_h\|_h= _w_h∈ E_h\|v_h-w_h\|_h. (6.12) Take vh=Πhωv_h= _hω. Then uh:=Π≤Λ(h)Πhω∈Ehu_h:= _≤ ^(h) _hω∈ E_h is the discrete best approximation to Πhω _hω in the discrete norm. To compare with P^h(ℛhΠhω) P_h(R_h _hω), we need a link between the discrete norm and the continuum L2L^2 norm after reconstruction. This link is standard for DEC/FEM mass matrices and is implicit in Assumption 6.1; we use the following consequence: there exists C such that for all sufficiently small h, ‖ℛhzh‖L2≤C‖zh‖hand‖zh‖h≤C‖ℛhzh‖L2∀zh∈Eh.\|R_hz_h\|_L^2≤ C\|z_h\|_h \|z_h\|_h≤ C\|R_hz_h\|_L^2 ∀ z_h∈ E_h. (6.13) (That is, ℛhR_h is a norm equivalence on the low-frequency space; for finite-dimensional EhE_h, this follows from stability plus uniform conditioning of the discrete mass matrix on EhE_h, ensured by shape-regularity and the spectral convergence assumption.) Let y∈E^hy∈ E_h. Then y=ℛhwhy=R_hw_h for some wh∈Ehw_h∈ E_h. Apply 6.12 and 6.13 to get ‖ℛhΠhω−ℛhuh‖L2 \|R_h _hω-R_hu_h\|_L^2 ≤C‖Πhω−uh‖h ≤ C\| _hω-u_h\|_h (6.14) =Cinfwh∈Eh‖Πhω−wh‖h =C _w_h∈ E_h\| _hω-w_h\|_h ≤Cinfwh∈Eh‖ℛhΠhω−ℛhwh‖L2 ≤ C _w_h∈ E_h\|R_h _hω-R_hw_h\|_L^2 =Cinfy∈E^h‖ℛhΠhω−y‖L2=C‖(Id−P^h)(ℛhΠhω)‖L2. =C _y∈ E_h\|R_h _hω-y\|_L^2=C\|(Id- P_h)(R_h _hω)\|_L^2. Now note that P^h(ℛhΠhω)−ℛhuh∈E^hand(Id−P^h)(ℛhΠhω)⟂E^h, P_h(R_h _hω)-R_hu_h∈ E_h (Id- P_h)(R_h _hω) E_h, hence, ‖ℛhΠhω−ℛhuh‖L22=‖(Id−P^h)(ℛhΠhω)‖L22+‖P^h(ℛhΠhω)−ℛhuh‖L22.\|R_h _hω-R_hu_h\|_L^2^2=\|(Id- P_h)(R_h _hω)\|_L^2^2+\| P_h(R_h _hω)-R_hu_h\|_L^2^2. (6.15) Combining 6.14 and 6.15 implies ‖P^h(ℛhΠhω)−ℛhuh‖L2≤C‖(Id−P^h)(ℛhΠhω)‖L2.\| P_h(R_h _hω)-R_hu_h\|_L^2≤ C\|(Id- P_h)(R_h _hω)\|_L^2. (6.16) Finally, ‖(Id−P^h)(v)‖L2≤‖v‖L2\|(Id- P_h)(v)\|_L^2≤\|v\|_L^2 for any v, and stability of ℛhΠhR_h _h yields ‖ℛhΠhω‖L2≤C‖ω‖L2≤C‖ω‖Hs\|R_h _hω\|_L^2≤ C\|ω\|_L^2≤ C\|ω\|_H^s. Therefore, ‖(I)‖L2=‖ℛhuh−P^h(ℛhΠhω)‖L2≤Chp‖ω‖Hs,\|(I)\|_L^2=\|R_hu_h- P_h(R_h _hω)\|_L^2≤ Ch^p\|ω\|_H^s, (6.17) where the final hph^p dependence is inherited from the subspace convergence (Assumption 6.2), which ensures the norm-equivalence constants and quasi-optimality transfer constants remain uniform and that E^h E_h approximates E with order hph^p. (Equivalently, one may sharpen this step by explicitly writing E^h E_h in terms of reconstructed eigenvectors and using 6.5 to bound the principal angles; this yields the same order hph^p.) Insert 6.10, 6.11, and 6.17 into the decomposition 6.9: ‖Phdω−Pω‖L2≤‖(I)‖L2+‖(II)‖L2+‖(III)‖L2≤Chp‖ω‖Hs.\|P_h^dω-Pω\|_L^2≤\|(I)\|_L^2+\|(I)\|_L^2+\|(I)\|_L^2≤ Ch^p\|ω\|_H^s. Discussion. The key point is subspace convergence of discrete eigenspaces to the continuum eigenspace below Λ , which ensures that applying the discrete projector to a discretized input selects approximately the same spectral content as the continuum projector. 6.5 Consistency of Truncated Functional Calculus We now prove convergence of the discrete multiplier operator to its continuum analogue. Theorem 6.4 (Consistency of discrete truncated multipliers). Fix Λ>0 >0 such that Λ∉spec(ΔH,g) ∉ spec ( _H,g). Let mθ∈C1([0,Λ])m_θ∈ C^1([0, ]) and set Mθ:=sup0≤λ≤Λ|mθ(λ)|,Lθ:=sup0≤λ≤Λ|mθ′(λ)|.M_θ:= _0≤λ≤ |m_θ(λ) |, L_θ:= _0≤λ≤ |m_θ (λ) |. Assume the discretization satisfies Assumptions 6.1–6.2, and moreover the interpolation/reconstruction and discrete inner product are (L2L^2)-consistent in the following standard sense: there exist p>0p>0 and s sufficiently large such that for all ω∈HsΩ1ω∈ H^s ^1 and all discrete eigenvectors ψk,h _k,h with λk,h≤Λ _k,h≤ , |⟨Πhω,ψk,h⟩h−⟨ω,ℛhψk,h⟩L2|≤Chp‖ω‖Hs. | _hω, _k,h _h- ω,R_h _k,h _L^2 |\ ≤\ Ch^p\|ω\|_H^s. (6.18) (For DEC/FEM, this is a mass-matrix/quadrature consistency bound on smooth inputs.) Then there exists C independent of h such that for all ω∈HsΩ1ω∈ H^s ^1, ‖ℛh^θ,Λ(h)Πhω−^θ,Λω‖L2≤C(Mθ+Lθ)hp‖ω‖Hs. \|R_h S^(h)_θ, _hω- S_θ, ω \|_L^2\ ≤\ C(M_θ+L_θ)h^p\|ω\|_H^s. (6.19) Proof Let E:=ran(Π≤Λ)⊂L2Ω1andEh:=ran(Π≤Λ(h))⊂Ωh1.E:= ran ( _≤ )⊂ L^2 ^1 E_h:= ran ( ^(h)_≤ )⊂ _h^1. Since Λ∉spec(ΔH,g) ∉ spec ( _H,g), E is finite-dimensional. Let m:=dim(E)m:= (E). By Assumption 6.2 (eigenspace convergence below Λ ), for sufficiently small h there exist orthonormal bases ψjj=1m⊂E(orthonormal in L2),ψj,hj=1m⊂Eh(orthonormal in ⟨⋅,⋅⟩h),\ _j\_j=1^m⊂ E (orthonormal in L^2), \ _j,h\_j=1^m⊂ E_h (orthonormal in ·,· _h), and associated eigenvalues λj≤Λ _j≤ , λj,h≤Λ _j,h≤ , such that |λj,h−λj|≤Chp,‖ℛhψj,h−ψj‖L2≤Chp. | _j,h- _j |≤ Ch^p, \|R_h _j,h- _j\|_L^2≤ Ch^p. (6.20) (If there are multiplicities, one chooses bases aligned via principal angles; 6.20 is a standard consequence of subspace convergence.) Let ω∈HsΩ1ω∈ H^s ^1 and define its continuum low-frequency component ω≤:=Π≤Λω=∑j=1mcjψj,cj:=⟨ω,ψj⟩L2. _≤:= _≤ ω= _j=1^mc_j _j, c_j:= ω, _j _L^2. (6.21) Similarly, define the discrete low-frequency coefficients of Πhω _hω in the discrete eigenbasis: cj,h:=⟨Πhω,ψj,h⟩h.c_j,h:= _hω, _j,h _h. (6.22) Then the discrete truncated multiplier applied to Πhω _hω is ^θ,Λ(h)Πhω=∑j=1mθ(λj,h)cj,hψj,h, S^(h)_θ, _hω= _j=1^mm_θ( _j,h)c_j,h _j,h, (6.23) and therefore ℛh^θ,Λ(h)Πhω=∑j=1mθ(λj,h)cj,hℛhψj,h.R_h S^(h)_θ, _hω= _j=1^mm_θ( _j,h)c_j,hR_h _j,h. (6.24) On the continuum side, since ^θ,Λ=mθ(ΔH,g)Π≤Λ S_θ, =m_θ( _H,g) _≤ , ^θ,Λω=∑j=1mθ(λj)cjψj. S_θ, ω= _j=1^mm_θ( _j)c_j _j. (6.25) Thus, the error we must bound is ‖∑j=1m(mθ(λj,h)cj,hℛhψj,h−mθ(λj)cjψj)‖L2. \| _j=1^m (m_θ( _j,h)c_j,hR_h _j,h-m_θ( _j)c_j _j ) \|_L^2. (6.26) For each j, add and subtract mθ(λj)cjℛhψj,hm_θ( _j)c_jR_h _j,h and mθ(λj,h)cjℛhψj,hm_θ( _j,h)c_jR_h _j,h to write mθ(λj,h)cj,hℛhψj,h−mθ(λj)cjψj=Tj(a)+Tj(b)+Tj(c),m_θ( _j,h)c_j,hR_h _j,h-m_θ( _j)c_j _j=T^(a)_j+T^(b)_j+T^(c)_j, (6.27) where Tj(a) T^(a)_j :=mθ(λj,h)(cj,h−cj)ℛhψj,h, =m_θ( _j,h)(c_j,h-c_j)R_h _j,h, (6.28) Tj(b) T^(b)_j :=(mθ(λj,h)−mθ(λj))cjℛhψj,h, = (m_θ( _j,h)-m_θ( _j) )c_jR_h _j,h, Tj(c) T^(c)_j :=mθ(λj)cj(ℛhψj,h−ψj). =m_θ( _j)c_j(R_h _j,h- _j). Then by the triangle inequality, ‖∑j=1m(⋅)‖L2≤∑j=1m(‖Tj(a)‖L2+‖Tj(b)‖L2+‖Tj(c)‖L2). \| _j=1^m(·) \|_L^2≤ _j=1^m (\|T^(a)_j\|_L^2+\|T^(b)_j\|_L^2+\|T^(c)_j\|_L^2 ). (6.29) We now bound each contribution. Using |mθ(λj)|≤Mθ |m_θ( _j) |≤ M_θ, Cauchy–Schwarz, and 6.20, ‖Tj(c)‖L2≤|mθ(λj)||cj|‖ℛhψj,h−ψj‖L2≤Mθ|cj|Chp.\|T^(c)_j\|_L^2≤ |m_θ( _j) | |c_j |\|R_h _j,h- _j\|_L^2≤ M_θ |c_j |Ch^p. (6.30) By the mean value theorem and the definition of LθL_θ, |mθ(λj,h)−mθ(λj)|≤Lθ|λj,h−λj|. |m_θ( _j,h)-m_θ( _j) |≤ L_θ | _j,h- _j |. (6.31) Using 6.20 and stability (|ℛhψj,h|L2≤C |R_h _j,h |_L^2≤ C) (from Assumption 6.1 and normalization), we obtain ‖Tj(b)‖L2≤|mθ(λj,h)−mθ(λj)||cj|‖ℛhψj,h‖L2≤CLθhp|cj|.\|T^(b)_j\|_L^2≤ |m_θ( _j,h)-m_θ( _j) | |c_j |\|R_h _j,h\|_L^2≤ CL_θh^p |c_j |. (6.32) We estimate cj,h−cjc_j,h-c_j. By definition and the consistency assumption 6.18, |cj,h−cj| |c_j,h-c_j | =|⟨Πhω,ψj,h⟩h−⟨ω,ψj⟩L2| = | _hω, _j,h _h- ω, _j _L^2 | (6.33) ≤|⟨Πhω,ψj,h⟩h−⟨ω,ℛhψj,h⟩L2|+|⟨ω,ℛhψj,h−ψj⟩L2| ≤ | _hω, _j,h _h- ω,R_h _j,h _L^2 |+ | ω,R_h _j,h- _j _L^2 | ≤Chp‖ω‖Hs+‖ω‖L2‖ℛhψj,h−ψj‖L2. ≤ Ch^p\|ω\|_H^s+\|ω\|_L^2\|R_h _j,h- _j\|_L^2. Using ‖ℛhψj,h−ψj‖L2≤Chp\|R_h _j,h- _j\|_L^2≤ Ch^p from 6.20 and ‖ω‖L2≤‖ω‖Hs\|ω\|_L^2≤\|ω\|_H^s (since s≥0s≥ 0 and ℳM is compact), we conclude |cj,h−cj|≤Chp‖ω‖Hs. |c_j,h-c_j |≤ Ch^p\|ω\|_H^s. (6.34) Therefore, using |mθ(λj,h)|≤Mθ |m_θ( _j,h) |≤ M_θ and ‖ℛhψj,h‖L2≤C\|R_h _j,h\|_L^2≤ C, ‖Tj(a)‖L2≤|mθ(λj,h)||cj,h−cj|‖ℛhψj,h‖L2≤CMθhp‖ω‖Hs.\|T^(a)_j\|_L^2≤ |m_θ( _j,h) | |c_j,h-c_j |\|R_h _j,h\|_L^2≤ CM_θh^p\|ω\|_H^s. (6.35) Collecting 6.30, 6.32, and 6.35 and summing over j=1,…,mj=1,…,m, we obtain ‖ℛh^θ,Λ(h)Πhω−^θ,Λω‖L2 \|R_h S^(h)_θ, _hω- S_θ, ω \|_L^2 ≤∑j=1m‖Tj(a)‖L2+∑j=1m‖Tj(b)‖L2+∑j=1m‖Tj(c)‖L2 ≤ _j=1^m\|T^(a)_j\|_L^2+ _j=1^m\|T^(b)_j\|_L^2+ _j=1^m\|T^(c)_j\|_L^2 (6.36) ≤CMθhp‖ω‖Hs+CLθhp∑j=1m|cj|+CMθhp∑j=1m|cj|. ≤ CM_θh^p\|ω\|_H^s+CL_θh^p _j=1^m |c_j |+CM_θh^p _j=1^m |c_j |. Since ψjj=1m\ _j\_j=1^m is L2L^2-orthonormal, ∑j=1m|cj|2=‖Π≤Λω‖L22≤‖ω‖L22 _j=1^m |c_j |^2=\| _≤ ω\|_L^2^2≤\|ω\|_L^2^2. Hence by Cauchy–Schwarz, ∑j=1m|cj|≤m(∑j=1m|cj|2)1/2≤m‖ω‖L2≤CΛ‖ω‖Hs, _j=1^m |c_j |≤ m ( _j=1^m |c_j |^2 )^1/2≤ m\|ω\|_L^2≤ C_ \|ω\|_H^s, (6.37) where CΛ:=mC_ := m depends only on Λ (and ℳ,gM,g, since m is the number of eigenvalues ≤Λ≤ , counted with multiplicity). Substitute 6.37 into 6.36 and absorb constants to obtain ‖ℛh^θ,Λ(h)Πhω−^θ,Λω‖L2≤C(Mθ+Lθ)hp‖ω‖Hs. \|R_h S^(h)_θ, _hω- S_θ, ω \|_L^2≤ C(M_θ+L_θ)h^p\|ω\|_H^s. Discussion. The rate is governed by (i) eigenspace convergence and (i) Lipschitz regularity of mθm_θ, which translates eigenvalue perturbations into multiplier perturbations. 6.6 Discretization Consistency of the Neural Operator We now lift the preceding consistency to the full network. We present the result in two versions: (i) linear truncated operator (the main approximation class for elliptic resolvents), and (i) nonlinear multi-layer architecture. 6.6.1 Linear model class Define the continuum operator ^θ,Λ=mθ(ΔH,g)Π≤Λ S_θ, =m_θ( _H,g) _≤ as in 4.7, and the discrete operator ^θ,Λ(h) S^(h)_θ, as in 6.3. Theorem 6.5 (Continuum–discrete commutation for the learned operator). Fix Λ>0 >0 such that Λ∉spec(ΔH,g) ∉ spec ( _H,g). Let mθ∈C1([0,Λ])m_θ∈ C^1([0, ]) and define ^θ,Λ:=mθ(ΔH,g)Π≤Λ,^θ,Λ(h):=mθ(ΔH,h)Π≤Λ(h). S_θ, :=m_θ( _H,g) _≤ , S^(h)_θ, :=m_θ( _H,h) ^(h)_≤ . Assume Assumptions 6.1–6.2, and the additional pairing consistency hypothesis 6.18 used in Theorem 6.4. Let p>0p>0 be the spectral convergence order and let s0s_0 be large enough such that Theorem 6.4 applies with s=s0s=s_0. Then there exists C independent of h such that ‖ℛh^θ,Λ(h)Πh−^θ,Λ‖Hs0⟶L2≤C(Mθ+Lθ)hp, \|R_h S^(h)_θ, _h- S_θ, \|_H^s_0 L^2\ ≤\ C(M_θ+L_θ)h^p, (6.38) where Mθ=sup[0,Λ]|mθ|M_θ= _[0, ] |m_θ | and Lθ=sup[0,Λ]|mθ′|L_θ= _[0, ] |m_θ |. Proof Let Th:=ℛh^θ,Λ(h)Πh−^θ,ΛT_h:=R_h S^(h)_θ, _h- S_θ, . We must show that the operator norm ‖Th‖Hs0⟶L2\|T_h\|_H^s_0 L^2 is O(hp)O(h^p) with the stated prefactor. By Theorem 6.4, for every ω∈Hs0Ω1ω∈ H^s_0 ^1 we have ‖Thω‖L2=‖ℛh^θ,Λ(h)Πhω−^θ,Λω‖L2≤C(Mθ+Lθ)hp‖ω‖Hs0.\|T_hω\|_L^2= \|R_h S^(h)_θ, _hω- S_θ, ω \|_L^2≤ C(M_θ+L_θ)h^p\|ω\|_H^s_0. (6.39) Here C is independent of h and ω, and depends only on the discretization family, (ℳ,g)(M,g), and the fixed cutoff Λ (through the dimension of the low-frequency space and the spectral-gap constants implicit in Assumption 6.2). By definition of the operator norm, ‖Th‖Hs0⟶L2=supω∈Hs0Ω1ω≠0‖Thω‖L2‖ω‖Hs0=sup|ω|Hs0=1‖Thω‖L2.\|T_h\|_H^s_0 L^2= _ subarraycω∈ H^s_0 ^1\\ ω≠ 0 subarray \|T_hω\|_L^2\|ω\|_H^s_0= _ |ω |_H^s_0=1\|T_hω\|_L^2. (6.40) Applying 6.39 to all ω with |ω|Hs0=1 |ω |_H^s_0=1 yields ‖Th‖Hs0⟶L2≤C(Mθ+Lθ)hp.\|T_h\|_H^s_0 L^2≤ C(M_θ+L_θ)h^p. (6.41) This is the central continuum–discrete consistency guarantee: applying the discrete learned operator to the discretized input and reconstructing converges to applying the continuum learned operator directly. 6.6.2 Nonlinear architecture Consider the depth-L network 3.9, and define its discrete counterpart by replacing each spectral layer θℓ,KT_ _ ,K with its discrete analogue 6.3, and replacing pointwise operations with their discrete fiberwise versions using the discrete metric/Hodge star at each simplex. Denote the resulting discrete network by ^θ,Λ,NL(h) S^(h)_θ, ,NL. Assume in addition that the fiberwise nonlinearities are uniformly Lipschitz on the relevant range of norms, i.e., there exists Lip(σθℓ)<∞Lip( _ _ )<∞ independent of h, and the scalar pointwise multipliers bθℓb_ _ are discretized consistently. Theorem 6.6 (Discretization consistency for the nonlinear intrinsic neural operator). Fix Λ>0 >0 with Λ∉spec(ΔH,g) ∉ spec ( _H,g). Consider the continuum depth-L network 3.9 ω0=f,ωℓ+1=ℱℓ(ωℓ):=σℓ(ℓωℓ+ℬℓωℓ),ℓ=0,…,L−1, _0=f, _ +1=F_ ( _ ):= _ (T_ _ +B_ _ ), =0,…,L-1, (6.42) where for each ℓ : • ℓ=mℓ(ΔH,g)Π≤ΛT_ =m_ ( _H,g) _≤ with mℓ∈C1([0,Λ])m_ ∈ C^1([0, ]); • ℬℓB_ is a pointwise scalar multiplication operator (ℬℓω)(x)=bℓ(x)ω(x)(B_ ω)(x)=b_ (x)ω(x) with bℓ∈L∞(ℳ)b_ ∈ L^∞(M); • σℓ _ is a fiberwise gauge-equivariant radial nonlinearity (σℓ(η))(x)=ρℓ(‖η(x)‖g)η(x)( _ (η))(x)= _ (\|η(x)\|_g)η(x) that is Lipschitz on L2Ω1L^2 ^1 on the relevant set of inputs (made precise below). Let ^θ,Λ,NL(f):=ωL S_θ, ,NL(f):= _L. Define the discrete network on Ωh1 _h^1 by ω0h=Πhf,ωℓ+1h=ℱℓh(ωℓh):=σℓh(ℓhωℓh+ℬℓhωℓh),ω^h_0= _hf, ω^h_ +1=F^h_ (ω^h_ ):=σ^h_ (T^h_ ω^h_ +B^h_ ω^h_ ), (6.43) where ℓh=mℓ(ΔH,h)Π≤Λ(h)T^h_ =m_ ( _H,h) ^(h)_≤ (discrete truncated multiplier), and ℬℓh,σℓhB^h_ ,σ^h_ are consistent discretizations of ℬℓ,σℓB_ , _ (assumptions below). Let ^θ,Λ,NL(h)(f):=ωLh S^(h)_θ, ,NL(f):=ω^h_L. Assume: • Stability of interpolation/reconstruction: There exists CRC_R independent of h such that ‖ℛhvh‖L2≤CR‖vh‖h∀vh∈Ωh1,‖Πhω‖h≤CR‖ω‖L2∀ω∈L2Ω1.\|R_hv_h\|_L^2≤ C_R\|v_h\|_h ∀ v_h∈ _h^1, \| _hω\|_h≤ C_R\|ω\|_L^2 ∀ω∈ L^2 ^1. (6.44) • Linear layer consistency: For each ℓ , there exists p>0p>0 and s0s_0 sufficiently large such that ‖ℛhℓhΠhω−ℓω‖L2≤Cℓhp‖ω‖Hs0∀ω∈Hs0Ω1, \|R_hT^h_ _hω-T_ ω \|_L^2≤ C_ h^p\|ω\|_H^s_0 ∀ω∈ H^s_0 ^1, (6.45) with Cℓ≲(Mℓ+Lℓ)C_ (M_ +L_ ), where Mℓ=sup[0,Λ]|mℓ|M_ = _[0, ] |m_ |, Lℓ=sup[0,Λ]|mℓ′|L_ = _[0, ] |m_ |. (This is Theorem 6.5 applied per layer.) • Pointwise operator stability and consistency: There exist constants Bℓ,εℓ,hB_ , _ ,h such that ‖ℬℓη‖L2≤Bℓ‖η‖L2,‖ℛhℬℓhvh−ℬℓℛhvh‖L2≤εℓ,h‖vh‖h,\|B_ η\|_L^2≤ B_ \|η\|_L^2, \|R_hB^h_ v_h-B_ R_hv_h \|_L^2≤ _ ,h\|v_h\|_h, (6.46) with εℓ,h≤Cℓ(B)hp _ ,h≤ C_ ^(B)h^p. • Nonlinearity Lipschitzness and discretization compatibility: There exists Lσ,ℓ<∞L_σ, <∞ such that for all η,ζ∈L2Ω1η,ζ∈ L^2 ^1 in the range encountered by the network, ‖σℓ(η)−σℓ(ζ)‖L2≤Lσ,ℓ‖η−ζ‖L2,\| _ (η)- _ (ζ)\|_L^2≤ L_σ, \|η-ζ\|_L^2, (6.47) and similarly the discrete nonlinearity is uniformly Lipschitz in reconstructed norm: ‖ℛh(σℓh(uh)−σℓh(vh))‖L2≤Lσ,ℓ‖ℛh(uh−vh)‖L2.\|R_h(σ^h_ (u_h)-σ^h_ (v_h))\|_L^2≤ L_σ, \|R_h(u_h-v_h)\|_L^2. (6.48) Finally, assume a consistency estimate: ‖ℛhσℓh(Πhη)−σℓ(η)‖L2≤Cℓ(σ)hp‖η‖Hs0.\|R_hσ^h_ ( _hη)- _ (η)\|_L^2≤ C_ ^(σ)h^p\|η\|_H^s_0. (6.49) Then there exists Cθ<∞C_θ<∞, depending on Lσ,ℓ,Bℓ,Mℓ,Lℓ=0L−1\L_σ, ,B_ ,M_ ,L_ \_ =0^L-1 and L but not on h, such that ‖ℛh^θ,Λ,NL(h)Πh−^θ,Λ,NL‖Hs0⟶L2≤Cθhp. \|R_h S^(h)_θ, ,NL _h- S_θ, ,NL \|_H^s_0 L^2\ ≤\ C_θh^p. (6.50) Proof For each layer ℓ , define the continuum pre-activation zℓ:=ℓωℓ+ℬℓωℓ,ωℓ+1=σℓ(zℓ),z_ :=T_ _ +B_ _ , _ +1= _ (z_ ), (6.51) and the discrete pre-activation zℓh:=ℓhωℓh+ℬℓhωℓh,ωℓ+1h=σℓh(zℓh).z^h_ :=T^h_ ω^h_ +B^h_ ω^h_ , ω^h_ +1=σ^h_ (z^h_ ). (6.52) Let the reconstructed discrete state be ω¯ℓh:=ℛhωℓh∈L2Ω1. ω^h_ :=R_hω^h_ ∈ L^2 ^1. (6.53) We will estimate the layerwise error eℓ:=‖ω¯ℓh−ωℓ‖L2.e_ :=\| ω^h_ - _ \|_L^2. (6.54) We start from eℓ+1=‖ω¯ℓ+1h−ωℓ+1‖L2=‖ℛhσℓh(zℓh)−σℓ(zℓ)‖L2.e_ +1=\| ω^h_ +1- _ +1\|_L^2=\|R_hσ^h_ (z^h_ )- _ (z_ )\|_L^2. (6.55) Add and subtract σℓ(ℛhzℓh) _ (R_hz^h_ ) and apply the triangle inequality: eℓ+1≤∥ℛhσℓh(zℓh)−σℓ(ℛhzℓh)⏟L2∥(I)ℓ+∥σℓ(ℛhzℓh)−σℓ(zℓ)⏟L2∥(II)ℓ.e_ +1≤ \|R_hσ^h_ (z^h_ )- _ (R_hz^h_ )_L^2\|_(I)_ + \| _ (R_hz^h_ )- _ (z_ )_L^2\|_(I)_ . (6.56) Bound (II)ℓ(I)_ . By the Lipschitz property 6.47, (II)ℓ≤Lσ,ℓ‖ℛhzℓh−zℓ‖L2.(I)_ ≤ L_σ, \|R_hz^h_ -z_ \|_L^2. (6.57) Bound (I)ℓ(I)_ . This term captures the discrete-vs-continuum nonlinearity mismatch. Under 6.48–6.49, applied with η=ℛhzℓhη=R_hz^h_ and using stability of Πh,ℛh _h,R_h to relate discrete arguments to continuum arguments, we obtain a bound of the form (I)ℓ≤Cℓ(σ)hp‖zℓ‖Hs0(I)_ ≤ C_ ^(σ)h^p\|z_ \|_H^s_0 (6.58) uniformly along the trajectories considered. (This is exactly 6.49 when the discrete nonlinearity is defined by applying σℓ _ pointwise after reconstruction; more general consistent constructions satisfy the same estimate.) Combining 6.56–6.58 yields eℓ+1≤Cℓ(σ)hp‖zℓ‖Hs0+Lσ,ℓ‖ℛhzℓh−zℓ‖L2.e_ +1≤ C_ ^(σ)h^p\|z_ \|_H^s_0+L_σ, \|R_hz^h_ -z_ \|_L^2. (6.59) Using 6.51–6.52, ℛhzℓh−zℓ=(ℛhℓhωℓh−ℓω¯ℓh)⏟Aℓ+(ℓω¯ℓh−ℓωℓ)⏟Bℓ+(ℛhℬℓhωℓh−ℬℓω¯ℓh)⏟Cℓ+(ℬℓω¯ℓh−ℬℓωℓ)⏟Dℓ.R_hz^h_ -z_ = (R_hT^h_ ω^h_ -T_ ω^h_ )_A_ + (T_ ω^h_ -T_ _ )_B_ + (R_hB^h_ ω^h_ -B_ ω^h_ )_C_ + (B_ ω^h_ -B_ _ )_D_ . (6.60) We bound each term in L2L^2. Term BℓB_ . Since ℓT_ is a bounded linear operator (L2⟶L2L^2 L^2) on the low-frequency range (indeed, ℓT_ has finite rank and ‖ℓ‖L2⟶L2≤Mℓ\|T_ \|_L^2 L^2≤ M_ ), we have ‖Bℓ‖L2≤‖ℓ‖L2⟶L2‖ω¯ℓh−ωℓ‖L2≤Mℓeℓ.\|B_ \|_L^2≤\|T_ \|_L^2 L^2\| ω^h_ - _ \|_L^2≤ M_ e_ . (6.61) Term DℓD_ . From 6.46, ‖Dℓ‖L2≤Bℓeℓ.\|D_ \|_L^2≤ B_ e_ . (6.62) Term AℓA_ . Insert Πhω¯ℓh _h ω^h_ and use 6.45: Aℓ=(ℛhℓhωℓh−ℛhℓhΠhω¯ℓh)+(ℛhℓhΠhω¯ℓh−ℓω¯ℓh).A_ = (R_hT^h_ ω^h_ -R_hT^h_ _h ω^h_ )+ (R_hT^h_ _h ω^h_ -T_ ω^h_ ). (6.63) The second difference is O(hp)O(h^p) by 6.45 with ω=ω¯ℓhω= ω^h_ , provided ω¯ℓh∈Hs0 ω^h_ ∈ H^s_0 along the trajectory (which is ensured if f∈Hs0f∈ H^s_0 and the layers map Hs0H^s_0 to itself; this is standard for finite-rank spectral operators plus pointwise maps with smooth scalar fields). Thus, ‖ℛhℓhΠhω¯ℓh−ℓω¯ℓh‖L2≤Cℓhp‖ω¯ℓh‖Hs0. \|R_hT^h_ _h ω^h_ -T_ ω^h_ \|_L^2≤ C_ h^p\| ω^h_ \|_H^s_0. (6.64) For the first difference in 6.63, use uniform boundedness of ℓhT^h_ and stability of ℛhR_h together with ωℓh−Πhω¯ℓhω^h_ - _h ω^h_ being a consistency residual of the reconstruction/interpolation pair; this contributes an O(hp)O(h^p) term under Assumption 6.1 (and is zero if ℛhΠhR_h _h is a projector on the discrete representation). In either case, one obtains ‖Aℓ‖L2≤Cℓ′hp‖ωℓ‖Hs0+Cℓ′eℓ,\|A_ \|_L^2≤ C_ h^p\| _ \|_H^s_0+C_ e_ , (6.65) for constants independent of h. (The eℓe_ part accounts for replacing ω¯ℓh ω^h_ by ωℓ _ inside norms.) Term CℓC_ . This is controlled by the pointwise consistency 6.46: ‖Cℓ‖L2=‖ℛhℬℓhωℓh−ℬℓω¯ℓh‖L2≤εℓ,h‖ωℓh‖h≤Chp‖ω¯ℓh‖L2≤Chp‖ωℓ‖Hs0+Chpeℓ,\|C_ \|_L^2=\|R_hB^h_ ω^h_ -B_ ω^h_ \|_L^2≤ _ ,h\|ω^h_ \|_h≤ Ch^p\| ω^h_ \|_L^2≤ Ch^p\| _ \|_H^s_0+Ch^pe_ , (6.66) where we used stability 6.44 and ‖ω¯ℓh‖L2≤‖ωℓ‖L2+eℓ\| ω^h_ \|_L^2≤\| _ \|_L^2+e_ . Putting 6.61–6.66 into 6.60 yields a bound of the form ‖ℛhzℓh−zℓ‖L2≤aℓeℓ+bℓhp‖f‖Hs0,\|R_hz^h_ -z_ \|_L^2≤ a_ e_ +b_ h^p\|f\|_H^s_0, (6.67) where aℓa_ depends on (Mℓ,Bℓ)(M_ ,B_ ) and uniform stability constants, and bℓb_ depends on the layerwise consistency constants (Cℓ,Cℓ(B))(C_ ,C_ ^(B)) and uniform bounds on the continuum trajectory in Hs0H^s_0 (which in turn depend on the network parameters and ‖f‖Hs0\|f\|_H^s_0). Insert 6.67 into 6.59: eℓ+1≤Cℓ(σ)hp‖zℓ‖Hs0+Lσ,ℓ(aℓeℓ+bℓhp‖f‖Hs0).e_ +1≤ C_ ^(σ)h^p\|z_ \|_H^s_0+L_σ, (a_ e_ +b_ h^p\|f\|_H^s_0 ). (6.68) Using uniform bounds ‖zℓ‖Hs0≤Cθ‖f‖Hs0\|z_ \|_H^s_0≤ C_θ\|f\|_H^s_0 along the continuum trajectory (a standard consequence of boundedness of the finite-rank spectral operators and pointwise multipliers on Hs0H^s_0), we obtain eℓ+1≤(Lσ,ℓaℓ)⏟=:Aℓeℓ+(Cℓ(σ)Cθ+Lσ,ℓbℓ)⏟=:Dℓhp‖f‖Hs0.e_ +1≤ (L_σ, a_ )_=:A_ e_ + (C_ ^(σ)C_θ+L_σ, b_ )_=:D_ h^p\|f\|_H^s_0. (6.69) The initial error satisfies e0=‖ℛhΠhf−f‖L2≤Chp‖f‖Hs0e_0=\|R_h _hf-f\|_L^2≤ Ch^p\|f\|_H^s_0 (6.70) by Assumption 6.1 (approximation of reconstruction–interpolation). Iterating 6.69 gives, for ℓ=0,…,L−1 =0,…,L-1, eL≤(∏j=0L−1Aj)e0+∑ℓ=0L−1(∏j=ℓ+1L−1Aj)Dℓhp‖f‖Hs0.e_L≤ ( _j=0^L-1A_j )e_0+ _ =0^L-1 ( _j= +1^L-1A_j )D_ h^p\|f\|_H^s_0. (6.71) Absorbing the finite products and sums into a constant CθC_θ depending only on the network parameters and depth (but not on h), we obtain eL≤Cθhp‖f‖Hs0.e_L≤ C_θh^p\|f\|_H^s_0. (6.72) Since eL=‖ℛh^θ,Λ,NL(h)(Πhf)−^θ,Λ,NL(f)‖L2e_L=\|R_h S^(h)_θ, ,NL( _hf)- S_θ, ,NL(f)\|_L^2, 6.72 is exactly the claimed operator bound 6.50 after taking the supremum over ‖f‖Hs0=1\|f\|_H^s_0=1. Discussion. The proof is a stability-by-composition argument: each layer has a consistency error of order hph^p, and these errors propagate through the network with amplification controlled by layerwise Lipschitz constants. Residual or contractive designs can ensure CθC_θ grows mildly with depth. 6.7 End-to-End Error Decomposition (Continuum Target vs. Discrete Learned Operator) Combining the approximation result of Section 4 with the discretization consistency above yields an end-to-end bound comparing the discrete learned operator to the true continuum PDE solution operator. For inputs in a smooth class Hs−1+γΩ1∩Hs0Ω1H^s-1+γ ^1∩ H^s_0 ^1, ‖ℛh^θ,Λ(h)Πh−g‖Hs−1+γ⟶Hs+1≤‖^θ,Λ−g‖Hs−1+γ⟶Hs+1⏟model/approximation error (Section 4)+‖ℛh^θ,Λ(h)Πh−^θ,Λ‖Hs0⟶Hs+1⏟discretization error (Section 6).\|R_h S^(h)_θ, _h-S_g\|_H^s-1+γ H^s+1≤ \| S_θ, -S_g\|_H^s-1+γ H^s+1_model/approximation error (Section 4)+ \|R_h S^(h)_θ, _h- S_θ, \|_H^s_0 H^s+1_discretization error (Section 6). (6.73) This inequality makes explicit the two independent levers for accuracy: increasing spectral capacity and improving multiplier approximation (Section 4), and refining the discretization (Section 6). 6.8 Summary The results of this section formalize a key practical requirement for operator learning on manifolds: mesh refinement should not change the learned operator in the continuum limit. Under standard assumptions for convergent discretizations of the Hodge Laplacian, the discrete intrinsic neural operator converges to its continuum counterpart, and its equivariance properties become exact in the limit. This provides a principled explanation for resolution-robust generalization observed in intrinsic spectral architectures. 7 Experiments: Algorithm and Implementation Details (Theory-Aligned) This section describes how the proposed Gauge-Equivariant Intrinsic Neural Operator (GINO) is instantiated on discrete meshes in a way that mirrors the continuum definitions in Sections 3–6. We emphasize conceptual alignment with the theory—intrinsic operators, gauge-consistent nonlinearities, and mesh-consistent discretizations—rather than engineering optimizations. 7.1 Discrete Representation of 1-Forms and Gauge Conventions 7.1.1 Discrete 1-forms On each triangulated manifold hT_h, we represent 1-forms using a compatible discretization Ωh1 _h^1 (DEC 1-cochains or Whitney 1-forms). The method provides: • a discrete inner product ⟨⋅,⋅⟩h ·,· _h approximating the L2L^2 pairing, • a discrete Hodge Laplacian ΔH,h _H,h acting on Ωh1 _h^1, • reconstruction ℛh:Ωh1→Ω1(ℳ)R_h: _h^1→ ^1(M) and interpolation Πh:Ω1(ℳ)→Ωh1 _h: ^1(M)→ _h^1. 7.1.2 Gauge for coordinate baselines When comparing to coordinate-based models that require per-vertex tangent frames, we attach an orthonormal frame EvE_v at each vertex v (e.g., via PCA of the local tangent plane) and express discrete 1-forms as frame coefficients. Gauge randomization is implemented by applying Rv∈O(d)R_v (d) to each frame, inducing the transformation av↦Rv⊤ava_v R_v a_v on coefficients. Our intrinsic model does not use these frames for computation and is therefore unaffected except through optional visualization. 7.2 Intrinsic Spectral Layers: Discrete Truncated Functional Calculus Each linear intrinsic layer is a truncated spectral multiplier of the discrete Hodge Laplacian, matching Section 4 and the discretization framework of Section 6. 7.2.1 Discrete eigensystem and truncation We compute the lowest portion of the spectrum of ΔH,h _H,h: ΔH,hψk,h=λk,hψk,h,⟨ψk,h,ψj,h⟩h=δkj. _H,h _k,h= _k,h _k,h, _k,h, _j,h _h= _kj. (7.1) Given a cutoff Λ (or equivalently a mode budget K), we form the discrete projector Π≤Λ(h)uh:=∑λk,h≤Λ⟨uh,ψk,h⟩hψk,h. ^(h)_≤ u_h:= _ _k,h≤ u_h, _k,h _h _k,h. (7.2) 7.2.2 Multiplier parameterization and application A multiplier layer is θ(h)uh:=mθ(ΔH,h)Π≤Λ(h)uh=∑λk,h≤Λmθ(λk,h)⟨uh,ψk,h⟩hψk,h.T^(h)_θu_h:=m_θ( _H,h) ^(h)_≤ u_h= _ _k,h≤ m_θ( _k,h) u_h, _k,h _h _k,h. (7.3) To align with the stability bounds in Section 5, we parameterize mθm_θ so that both Mθ=supλ∈[0,Λ]|mθ(λ)|,Lθ=supλ∈[0,Λ]|mθ′(λ)|M_θ= _λ∈[0, ] |m_θ(λ) |, L_θ= _λ∈[0, ] |m_θ (λ) | (7.4) are controlled. Concretely, we use a low-degree basis expansion on [0,Λ][0, ], mθ(λ)=∑j=0Jθjφj(λ),m_θ(λ)= _j=0^J _j _j(λ), (7.5) where φj\ _j\ are smooth functions with known derivative bounds (e.g., Chebyshev polynomials mapped to [0,Λ][0, ] or B-splines). This choice directly corresponds to Theorem 6.4, where Lipschitz regularity converts eigenvalue errors into multiplier errors. 7.3 Gauge-Equivariant Nonlinearity and Pointwise Terms Each nonlinear layer follows the intrinsic design of Section 3. 7.3.1 Pointwise scalar multiplication We include a pointwise (zeroth-order) term ℬθ(h)uh:=bθuh,B^(h)_θu_h:=b_θu_h, (7.6) where bθb_θ is a scalar field represented at vertices or elements. This term is intrinsically defined and does not depend on a frame choice. Its discretization is chosen so that ‖ℛhℬθ(h)uh−ℬθℛhuh‖L2→0as h→0,\|R_hB^(h)_θu_h-B_θR_hu_h\|_L^2→ 0 h→ 0, (7.7) consistent with the assumptions used in Theorem 6.6. 7.3.2 Fiberwise radial nonlinearity (gauge equivariant) To preserve gauge equivariance, the nonlinearity is applied radially in the fiber: σθ(u)(x)=ρθ(‖u(x)‖g)u(x). _θ(u)(x)= _θ(\|u(x)\|_g)u(x). (7.8) On the mesh, ‖u(x)‖g\|u(x)\|_g is evaluated using the discrete metric associated with the Hodge star/mass matrix. This guarantees that under any local orthonormal frame change R(x)∈O(d)R(x) (d), coefficient vectors rotate while norms remain invariant, yielding exact equivariance at the discrete level up to the discretization of the metric norm. For consistency and stability, ρθ _θ is chosen to be Lipschitz on the range of encountered norms, so that σθ _θ is Lipschitz in L2L^2, matching the hypotheses of Theorem 6.6. We parameterize ρθ _θ as a smooth bounded function (e.g., a spline or a small MLP with bounded activation), and, when needed, constrain sup|ρθ′| | _θ | to avoid exploding Lipschitz constants across depth. 7.4 Network Structure and Training Objective 7.4.1 Layer update rule A depth-L GINO is implemented as uh(0)=fh,uh(ℓ+1)=σθℓ(θℓ(h)uh(ℓ)+ℬθℓ(h)uh(ℓ)),ℓ=0,…,L−1.u^(0)_h=f_h, u^( +1)_h= _ _ (T^(h)_ _ u^( )_h+B^(h)_ _ u^( )_h ), =0,…,L-1. (7.9) This is the discrete analogue of the continuum operator in Section 3 and is precisely the model class analyzed in Sections 4–6. 7.4.2 Supervised operator regression For each geometry and forcing sample, we generate a target output uhu_h by solving the discrete PDE defining the target operator. The model prediction u^h u_h is trained by minimizing an intrinsic loss: ℒ(θ)=[‖u^h−uh‖h2],L(θ)=E [\| u_h-u_h\|_h^2 ], (7.10) optionally augmented by the energy norm associated with the operator when the target is elliptic: ℒenergy(θ)=[⟨u^h−uh,(ΔH,h+αI)(u^h−uh)⟩h].L_energy(θ)=E [ u_h-u_h,( _H,h+α I)( u_h-u_h) _h ]. (7.11) The energy norm aligns with the coercive bilinear form of the shifted elliptic operator and directly reflects the Sobolev mapping Hs−1⟶Hs+1H^s-1 H^s+1 discussed in Section 4. 7.5 Cutoffs, Mode Budgets, and Alignment with Approximation Theory The approximation bound in Theorem 4.4 decomposes error into truncation bias and multiplier approximation on [0,Λ][0, ]. Accordingly, we treat Λ (or K) as an explicit experimental axis: • Increasing Λ reduces truncation bias for sufficiently smooth inputs (Lemma 4.2). • Increasing multiplier capacity (larger J in 7.5) reduces εΛ(θ) _ (θ) on [0,Λ][0, ] (Lemma 4.3). • Cross-resolution tests directly probe the discretization-consistency rate hph^p (Theorems 6.4–6.6). We report ablations that vary (Λ,J,L)( ,J,L) to empirically validate these qualitative predictions. 7.6 Discrete–Continuum Consistency Checks To directly test Section 6, we include consistency diagnostics beyond task loss: 1. Projector consistency: evaluate ‖ℛhΠ≤Λ(h)Πhω−Π≤Λω‖L2\|R_h ^(h)_≤ _hω- _≤ ω\|_L^2 (7.12) on synthetic smooth ω, verifying convergence as h↓0h 0 (Lemma 6.3). 2. Layer consistency: evaluate the single-layer discrepancy ‖ℛhθ(h)Πhω−θω‖L2\|R_hT^(h)_θ _hω-T_θω\|_L^2 (7.13) and confirm the dependence on (Mθ,Lθ)(M_θ,L_θ) predicted by Theorem 6.4. 3. Gauge sensitivity: for frame-dependent baselines, repeat inference under multiple gauge draws and report variance of predictions; for GINO, report invariance up to numerical tolerance. 7.7 Practical Notes on Ensuring Theoretical Assumptions The stability theorems require controlling multiplier oscillations and nonlinearity Lipschitz constants. In practice we enforce this by: • Multiplier smoothness: choosing mθm_θ in a smooth basis and penalizing large derivatives on [0,Λ][0, ], approximating control of LθL_θ in 7.4. • Nonlinearity boundedness: choosing ρθ _θ bounded and Lipschitz so that each σθ _θ has a controlled Lipschitz constant in L2L^2, limiting error amplification across depth (Theorem 6.6). • Cutoff separation: avoiding cutoffs Λ near dense clusters of eigenvalues when evaluating stability to geometric perturbations (Section 5), consistent with the spectral separation logic used in Theorem 5.3. These choices are minimal and are included solely to reflect the assumptions required for the theory-to-practice correspondence. 8 Experiments We empirically validate the core theoretical claims of our Gauge-Equivariant Intrinsic Neural Operator (GINO) on controlled elliptic operator-learning problems over the 2D flat torus 2T^2 . The periodic domain admits closed-form solutions in the Fourier domain, enabling precise diagnostics of gauge equivariance, metric stability, and discretization consistency—properties directly tied to our intrinsic formulation. All models are implemented in PyTorch and trained on a single RTX 4060 GPU. 8.1 Problem Setup and Data Generation Target operator. We consider the elliptic resolvent family (Δg+αI)u=f( _g+α I)u=f , where g is a Riemannian metric on 2T^2 , Δg _g is the Laplace–Beltrami operator, and α>0α>0 ensures invertibility. Under a constant metric tensor M≻0M 0 with inverse A=M−1A=M^-1 , the solution is given spectrally by u^(k)=f^(k)k⊤Ak+α. u(k)= f(k)k Ak+α. Forcing distribution. Inputs f are sampled as 2-channel band-limited fields (interpreted as 1-forms in a local frame) via: (i) white noise in Fourier space, (i) power-law decay (1+λ)−β/2(1+λ)^-β/2 with cutoff λ≤λfλ≤ _f , and (i) per-sample RMS normalization after inverse FFT. Evaluation metrics. We report physical-space MSE, relative L2L^2 error (RelL2), and relative energy error (RelEnergy), defined as RelEnergy=‖(Λg+α)(u−u^)‖2‖(Λg+α)u‖2,RelEnergy= \|( _g+α)(u- u)\|_2\|( _g+α)u\|_2, which weights errors by the natural norm of the elliptic operator. 8.2 Models and Training Protocol GINO (ours). Our architecture applies a learnable spectral multiplier mθ(λg(k))m_θ( _g(k)) parameterized via Chebyshev polynomials on [0,Λ][0, ] , followed by a radial nonlinearity σ(u)=ρ(|u|)uσ(u)=ρ(|u|)u that enforces gauge equivariance under local frame rotations. An additional multiplier stage produces the final output. Baseline. We compare against CoordCNN—a coordinate-aware convolutional network operating directly on grid values—representing standard local, coordinate-dependent architectures lacking geometric priors. All models are optimized with AdamW and gradient clipping. Identical data generation and evaluation protocols are used across methods. 8.3 E1: Operator Learning Accuracy on Base Geometry We first evaluate GINO’s ability to learn the base resolvent g0S_g_0 under the identity metric ( M=IM=I ). As shown in Figure 1, all metrics converge rapidly within 1,000 training steps, reaching RelL2 ≈10−3≈ 10^-3 and RelEnergy ≈10−3≈ 10^-3 . A qualitative prediction (Figure 2) confirms near-perfect alignment between u u and u , with absolute error below 4×10−44× 10^-4 . This establishes that GINO’s intrinsic spectral parameterization is sufficiently expressive to approximate the target operator accurately. Figure 1: Training convergence of GINO on the base geometry (E1). Evaluation metrics (MSE, RelL2, RelEnergy) rapidly decrease and stabilize, reaching near numerical precision. Figure 2: Qualitative prediction example at step 4000 (E1). From left to right: |f||f| , |u||u| , |u^|| u| , and |u^−u|| u-u| . The model produces visually accurate predictions with small residual errors. 8.4 E2: Gauge Equivariance vs. Baseline Sensitivity To test invariance under global frame rotations R∈SO(2)R (2) , we measure whether F(R⊤f)≈R⊤F(f)F(R f)≈ R F(f) . Figure 3 reveals a stark contrast: GINO achieves an equivariance error of 1.75×10−71.75× 10^-7 (near floating-point precision), while CoordCNN exhibits large gauge dependence, with normalized standard deviation and worst-case deviation remaining at order 1 throughout training. This confirms that only GINO respects the underlying geometric symmetry. Figure 3: Gauge equivariance and sensitivity analysis (E2). Left: GINO’s equivariance error remains below 2×10−72× 10^-7 . Middle and right: CoordCNN shows strong gauge dependence with relative deviations near 1. 8.5 E3: Stability Under Metric Perturbations We evaluate generalization to anisotropic metrics M(δ)=Rdiag(1+δ,1−δ)R⊤M(δ)=Rdiag(1+δ,1-δ)R for δ∈[0,0.30]δ∈[0,0.30] , without retraining. As Figures 4 show, GINO maintains low error across all perturbations (RelL2 ∈[5×10−3,8×10−3]∈[5× 10^-3,8× 10^-3] ), whereas CoordCNN suffers from consistently high errors (RelL2 ≈0.24≈ 0.24 , RelEnergy ≈1.3≈ 1.3 ). This demonstrates that GINO learns a geometry-aware operator, while the baseline overfits to coordinate-specific representations. Figure 4: Metric perturbation stability (E3). GINO (blue) maintains low RelEnergy (left) and RelL2 (right) across ‖M−I‖F∈[0,0.4]\|M-I\|_F∈[0,0.4] , while CoordCNN (orange) shows large, invariant errors. 8.6 E4: Cross-Resolution Generalization and Discretization Consistency We assess whether learned operators commute with spectral restriction ( R ) and prolongation ( P ). In cross-resolution tests (Figures 5), GINO achieves RelL2 ≤10−2≤ 10^-2 in all train/test grid combinations, while CoordCNN degrades severely under resolution shift (e.g., coarse → fine RelL2 >0.5>0.5 ). Furthermore, GINO satisfies discretization consistency: commutation errors remain below 8×10−38× 10^-3 (Figure 6), versus >0.6>0.6 for CoordCNN. Figure 5: Cross-resolution generalization (E4). GINO generalizes across resolutions with low RelL2 (left) and RelEnergy (right); CoordCNN fails dramatically when tested on unseen grids. Figure 6: Discretization consistency diagnostic (E4). GINO achieves negligible commutation errors (<0.01<0.01), while CoordCNN errors exceed 0.6, revealing strong grid dependence. 8.7 E5: Structure-Preserving Hodge Decomposition We extend GINO to learn a regularized Helmholtz–Hodge decomposition of 1-forms into exact, coexact, and harmonic components. Final performance (Figure 7) shows RelL2 ≈3.4×10−2≈ 3.4× 10^-2 for both projections and a decomposition residual of 1.45×10−21.45× 10^-2 . Critically, gauge equivariance error remains at 3.6×10−63.6× 10^-6 (Figure 8), confirming that the operator is both structure-preserving and symmetry-respecting. Figure 7: Summary of E5 performance. Final relative errors for exact/coexact components are ≈3.4×10−2≈ 3.4× 10^-2 , with decomposition residual ≈1.45×10−2≈ 1.45× 10^-2 . Figure 8: Training dynamics for E5. Left: gauge equivariance error stays below 4×10−64× 10^-6 . Middle/right: projection errors and decomposition residual decrease rapidly during training. 8.8 E6: Ablations on Spectral Truncation and Multiplier Smoothness E6A: Spectral support. We sweep the truncation threshold Λ (Figure 9). Error metrics exhibit a U-shaped trend, with optimal performance at Λ≈100 ≈ 100 . Larger Λ increases the roughness proxy, indicating a capacity–regularization tradeoff. E6B: Smoothness–stability link. Strongly regularized models maintain flat error curves across δ (Figure 10), while unregularized (rough) models amplify RelL2 by up to 3.4× at δ=0.30δ=0.30 . The roughness proxy (Figure 10) correlates strongly with instability, providing empirical evidence that multiplier smoothness mediates geometric robustness. Figure 9: Ablation on spectral truncation (E6A). Error metrics (RelL2, RelEnergy) are minimized at intermediate Λ , while roughness increases monotonically with Λ . Figure 10: Effect of smoothness on geometric stability (E6B). Left: stronger regularization yields smoother multipliers. Right: unregularized models amplify error by up to 3.4× under metric perturbation, while smoothed models remain stable. 8.9 Summary Across six experiments, GINO consistently outperforms a standard coordinate-based baseline in accuracy, symmetry preservation, geometric stability, and discretization consistency. Crucially, ablation studies reveal that these advantages stem from principled design choices—spectral truncation and smoothness control—that align with theoretical intuitions about intrinsic operator learning. Together, these results establish GINO as a robust, geometry-aware framework for learning elliptic operators on vector-valued fields. 9 Conclusion We introduced Gauge-Equivariant Intrinsic Neural Operators (GINO) for learning elliptic PDE solution maps in a manner that is consistent with geometric structure and representation changes. By parameterizing the operator through intrinsic spectral multipliers defined on geometry-dependent spectra and combining them with gauge-equivariant nonlinearities, GINO enforces frame-consistent behavior while maintaining the expressivity required for operator learning. Experiments E1–E6 on 2T^2 provide a coherent empirical validation of the approach. GINO achieves accurate operator approximation on the base resolvent task and exhibits near numerical-precision gauge equivariance, while coordinate-based CNN baselines remain strongly gauge-sensitive. Under structured metric perturbations, GINO shows minimal degradation, and under cross-resolution evaluation it maintains low error and small commutation defects with restriction/prolongation operators, indicating reduced discretization dependence. On a regularized Helmholtz–Hodge decomposition task, the structured variant attains low projection error and small decomposition residual while preserving gauge equivariance to near floating-point tolerance. Finally, theory-aligned ablations demonstrate that multiplier smoothness strongly mediates stability: rough learned multipliers amplify perturbation errors, whereas smoothness control yields markedly improved robustness. These findings suggest that incorporating intrinsic geometry and gauge constraints is a practical route to more reliable operator surrogates for geometric PDEs. Future work includes extending the framework to spatially varying metrics and general manifolds, integrating mesh-based discretizations via discrete exterior calculus or finite element exterior calculus, and studying broader PDE families and real-world scientific benchmarks where gauge and discretization robustness are essential. Acknowledgments and Disclosure of Funding This research received no external funding. Remainder omitted in this sample. See http://w.jmlr.org/papers/ for full paper. References D. Arnold, R. Falk, and R. Winther (2010) Finite element exterior calculus: from hodge theory to numerical stability. Bulletin of the American Mathematical Society 47 (2), p. 281–354. Cited by: §1. D. N. Arnold, R. S. Falk, and R. Winther (2006) Finite element exterior calculus, homological techniques, and applications. Acta Numerica 15, p. 1–155. Cited by: §1. E. J. Bekkers (2019) B-spline cnns on lie groups. arXiv preprint arXiv:1909.12057. Cited by: §1. S. C. Brenner and L. R. Scott (2008) The mathematical theory of finite element methods. Springer New York, New York, NY. Cited by: §1. M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković (2021) Geometric deep learning: grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2104.13478. Cited by: §1. S. L. Brunton, B. R. Noack, and P. Koumoutsakos (2020) Machine learning for fluid mechanics. Annual Review of Fluid Mechanics 52 (1), p. 477–508. Cited by: §1, §1. T. S. Cohen and M. Welling (2016a) Steerable cnns. arXiv preprint arXiv:1612.08498. Cited by: §1. T. Cohen, M. Weiler, B. Kicanaoglu, and M. Geiger (2019) Gauge equivariant convolutional networks and the icosahedral cnn. In International Conference on Machine Learning (ICML), p. 1321–1330. Cited by: §1. T. Cohen and M. Welling (2016b) Group equivariant convolutional networks. In International Conference on Machine Learning (ICML), p. 2990–2999. Cited by: §1. P. de Haan, T. S. Cohen, and M. Welling (2020) Natural graph networks. Advances in Neural Information Processing Systems 33, p. 3636–3646. Cited by: §1. M. Desbrun, E. Kanso, and Y. Tong (2006) Discrete differential forms for computational modeling. In ACM SIGGRAPH 2006 Courses, p. 39–54. Cited by: §1. F. Fuchs, D. Worrall, V. Fischer, and M. Welling (2020) Se (3)-transformers: 3d roto-translation equivariant attention networks. Advances in Neural Information Processing Systems 33, p. 1970–1981. Cited by: §1. M. Geiger and T. Smidt (2022) E3n: euclidean neural networks. arXiv preprint arXiv:2207.09453. Cited by: §1. A. N. Hirani (2003) Discrete exterior calculus. California Institute of Technology. Cited by: §1. T. J. R. Hughes (2003) The finite element method: linear static and dynamic finite element analysis. Courier Corporation. Cited by: §1. J. Jost (2005) Riemannian geometry and geometric analysis. Springer Berlin Heidelberg, Berlin, Heidelberg. Cited by: §1. G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang (2021) Physics-informed machine learning. Nature Reviews Physics 3 (6), p. 422–440. Cited by: §1. J. Köhler, L. Klein, and F. Noé (2020) Equivariant flows: exact likelihood generative learning for symmetric densities. In International Conference on Machine Learning (ICML), p. 5361–5370. Cited by: §1. R. Kondor and S. Trivedi (2018) On the generalization of equivariance and convolution in neural networks to the action of compact groups. In International Conference on Machine Learning (ICML), p. 2747–2755. Cited by: §1. N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar (2023) Neural operator: learning maps between function spaces with applications to pdes. Journal of Machine Learning Research 24 (89), p. 1–97. Cited by: §1, §1. Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2020) Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895. Cited by: §1, §1. L. Lu, P. Jin, G. Pang, G. Zhang, and G. E. Karniadakis (2021) Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Intelligence 3 (3), p. 218–229. Cited by: §1, §1. M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, p. 686–707. Cited by: §1. V. G. Satorras, E. Hoogeboom, and M. Welling (2021) E (n) equivariant graph neural networks. In International Conference on Machine Learning (ICML), p. 9323–9332. Cited by: §1. K. Stachenfeld, D. B. Fielding, D. Kochkov, M. Cranmer, T. Pfaff, J. Godwin, C. Cui, S. Ho, P. Battaglia, and A. Sanchez-Gonzalez (2022) Learned simulators for turbulence. In International Conference on Learning Representations (ICLR), Cited by: §1. M. Takamoto, T. Praditia, R. Leiteritz, L. Hennigh, A. Maier, and N. Thuerey (2022) PDEBench: an extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems 35, p. 1596–1611. Cited by: §1. R. Temam (2024) Navier–stokes equations: theory and numerical analysis. AMS Chelsea Publishing, American Mathematical Society, Providence, RI. Cited by: §1. N. Thomas, T. Smidt, S. Kearnes, P. F. Riley, K. Hansen, B. McCloskey, L. Yang, and W. Li (2018) Tensor field networks: rotation-and translation-equivariant neural networks for 3d point clouds. arXiv preprint arXiv:1802.08219. Cited by: §1. F. W. Warner (1983) Foundations of differentiable manifolds and lie groups. Springer Science & Business Media. Cited by: §1. M. Weiler and G. Cesa (2019) General e (2)-equivariant steerable cnns. Advances in Neural Information Processing Systems 32. Cited by: §1. D. E. Worrall, S. J. Garbin, D. Turmukhambetov, and G. J. Brostow (2017) Harmonic networks: deep translation and rotation equivariance. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), p. 5028–5037. Cited by: §1.