Spectral Degeneracy Operators for Interpretable Turbulence Modeling

1. Introduction

The intersection of degenerate partial differential equations (PDEs), neural network symmetrization, and turbulence modeling offers a rich ground for mathematical innovation. Recent advances in inverse problems for degenerate PDEs [1] and physics-informed neural networks (PINNs) [2] have enabled the integration of adaptive diffusion mechanisms into machine learning models. Yet, a unified theoretical framework connecting these areas remains largely unexplored.

1.1. Degenerate PDEs and Inverse Problems

Degenerate PDEs naturally arise in contexts such as anisotropic diffusion [3], geometric singularities [4], and phase transitions [5]. Cannarsa et al. [1] established Lipschitz stability for reconstructing degeneracy points in parabolic equations of the form

$${\partial }_{t}w-{\partial }_x\left({|x-a|}^{\theta }{\partial }_{x}w\right)-cw=0,\theta \in [1,2),$$

(1)

from boundary measurements of ${\partial }_{x}w(1,t)$. Previous works on inverse source problems [6] and coefficient identification [7] mainly addressed scalar, one-dimensional domains, leaving multi-dimensional, vector-valued problems largely uncharted.

1.2. Neural Symmetrization

Equivariant neural networks [8] and the broader field of geometric deep learning [9] embed symmetry principles into machine learning models. Conventional group-convolution approaches [10] struggle with continuous symmetries such as $\mathrm{SO}\left(3\right)$ and anisotropic phenomena, including turbulent shear layers. Neural operators [11] and PINNs [2] offer promise for turbulence modeling, but often lack structural guarantees like rotation equivariance, energy conservation, or adaptivity to localized singularities.

1.3. Turbulence Modeling

Classical turbulence models—LES [12] and RANS [13]—rely heavily on empirical closures, which poorly capture intermittency and anisotropic dissipation. Data-driven neural closures [14,15] improve predictive performance but can violate fundamental physical constraints. Our approach enforces incompressibility via

$$\nabla ·\left({\mathbb{T}}_{\mathrm{NN}}\nabla \mathbf{u}\right)=0,$$

(2)

where ${\mathbb{T}}_{\mathrm{NN}}$ is a degeneracy-aware neural operator designed to respect the underlying PDE structure.

1.4. Contributions

We introduce spectral degeneracy operators (SDOs)—differential operators that encode both physical symmetries and adaptive singularities—and demonstrate their application to:

• Neural symmetrization, through SDO-based activation and layer design,
• Turbulence closure modeling, via data-driven calibration and spectral filtering,
• Inverse problem formulation, for reconstructing degeneracy points from sparse or boundary observations.

The key theoretical contributions of this work are:

• Generalized spectral decomposition for vector-valued SDOs (Section 2),
• Lipschitz stability results for inverse calibration in turbulence models (Section 3),
• A neural-turbulence correspondence theorem, connecting learned SDO parameters to underlying turbulent structures (Section 3).

2. Spectral Degeneracy Operators (SDOs)

2.1. Definition and Spectral Properties

Let $\Omega \subset {\mathbb{R}}^d$ be a bounded Lipschitz domain. Denote by $\mathbf{a}\in L^{\infty }(\Omega ;\Omega )$ the degeneracy points and $\mathbf{\theta }\in L^{\infty }(\Omega ;{[1,2)}^d)$ the degeneracy exponents. Define the spectral degeneracy operator (SDO) as

$${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}u\nabla ·\left({|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }}\circ \nabla u\right),$$

(3)

where

$${|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }}\prod _{i=1}^d{|x_i-a_i|}^{{\theta }_i},\circ \mathrm{denotes}\mathrm{Hadamard}\mathrm{product}.$$

(4)

The natural energy space associated with ${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}$ is

$${\mathcal{H}}_{\mathbf{\theta }}^1(\Omega )\left\{u\in L^2{(\Omega ):|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }/2}\nabla u\in L^2(\Omega ;{\mathbb{R}}^d){,u|}_{\partial \Omega }=0\right\},$$

(5)

with inner product

$${\langle u,v\rangle }_{{\mathcal{H}}_{\mathbf{\theta }}^1}{\int }_{\Omega }uvd\mathbf{x}+{\int }_{\Omega }\left({|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }/2}\nabla u\right)·\left({|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }/2}\nabla v\right)d\mathbf{x}.$$

(6)

Theorem 2.1

(Spectral Decomposition of SDOs). Let $\mathbf{a}\in \mathrm{int}(\Omega )$ and $\mathbf{\theta }\in {[1,2)}^d$. Then the operator

$${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}u:=-\nabla ·\left({\mathbf{D}}_{\mathbf{a},\mathbf{\theta }}\nabla u\right),{\mathbf{D}}_{\mathbf{a},\mathbf{\theta }}=\mathrm{diag}\left(\right|x_1-a_1{|}^{{\theta }_1},\dots ,|x_d-a_d{|}^{{\theta }_d}),$$

with domain $H_0^1(\Omega )$, isself-adjoint,positive semi-definite, and has acompact resolventon $L^2(\Omega )$.

Moreover, there exists a countable set of eigenpairs ${\left\{({\lambda }_k,{\varphi }_k)\right\}}_{k\in {\mathbb{N}}^d}$ forming an orthonormal basis of $L^2(\Omega )$ such that

$${\varphi }_k(\mathbf{x};\mathbf{a},\mathbf{\theta })=\prod _{i=1}^d{\varphi }_{k_i}^{\left(i\right)}(x_i;a_i,{\theta }_i),$$

(7)

where each ${\varphi }_{k_i}^{\left(i\right)}$ satisfies the 1D Bessel-type Sturm–Liouville problem

$$-\left(\right|x_i-a_i{{|}^{{\theta }_i}{\varphi }_{k_i}^{\left(i\right)\prime })}^{\prime }={\lambda }_{k_i}{\varphi }_{k_i}^{\left(i\right)},x_i\in (0,1),{\varphi }_{k_i}^{\left(i\right)}\left(0\right)={\varphi }_{k_i}^{\left(i\right)}\left(1\right)=0.$$

(8)

The eigenvalues admit the asymptotic behavior

$${\lambda }_k(\mathbf{a},\mathbf{\theta })\sim \sum _{i=1}^d{\left({\displaystyle \frac{j_{{\nu }_i,k_i}}{|1-a_i|}}\right)}^2,{\nu }_i={\displaystyle \frac{{\theta }_i-1}{2-{\theta }_i}},$$

(9)

where $j_{{\nu }_i,k_i}$ denotes the $k_i$-th positive zero of the Bessel function $J_{{\nu }_i}$ of the first kind.

Proof. 

Assume $u\left(\mathbf{x}\right)={\prod }_{i=1}^d{u}_i\left(x_i\right)$. Then the operator acts as

$${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}u\left(\mathbf{x}\right)=\sum _{i=1}^d\left[\prod _{j\ne i}{|x_j-a_j|}^{{\theta }_j}\right]\left(\right|x_i-a_i{{|}^{{\theta }_i}{u}_i^{\prime })}^{\prime }.$$

Thus, the d-dimensional spectral problem reduces to d decoupled 1D weighted Sturm–Liouville problems (8).

Each 1D operator

$${\mathcal{L}}_i{u}_i:=-\left(\right|x_i-a_i{{|}^{{\theta }_i}{u}_i^{\prime })}^{\prime },u_i\in H_0^1(0,1),$$

is self-adjoint with respect to the weighted inner product

$${\langle u_i,v_i\rangle }_{L^2((0,1),|x_i-a_i{|}^{{\theta }_i}d{x}_i)}:={\int }_0^1{u}_i\left(x_i\right)v_i\left(x_i\right){|x_i-a_i|}^{{\theta }_i}d{x}_i.$$

The operator is positive semi-definite because

$$\langle {\mathcal{L}}_i{u}_i,u_i\rangle ={\int }_0^1|x_i-a_i{|}^{{\theta }_i}{\left|u_i^{\prime }\right|}^{2}d{x}_i\ge 0.$$

Compactness of the resolvent follows from standard embeddings $H_0^1\hookrightarrow L^2$ for singular-weighted Sobolev spaces [1].

By tensorizing the 1D eigenfunctions ${\varphi }_{k_i}^{\left(i\right)}$, we obtain

$${\varphi }_k\left(\mathbf{x}\right)=\prod _{i=1}^d{\varphi }_{k_i}^{\left(i\right)}\left(x_i\right),$$

with eigenvalues

$${\lambda }_k=\sum _{i=1}^d{\lambda }_{k_i}.$$

Orthonormality in $L^2(\Omega )$ follows from the product structure and orthonormality of each 1D set.

For each 1D Sturm–Liouville problem with weight $|x_i-a_i{|}^{{\theta }_i}$, standard singular Bessel function theory [16] implies

$${\lambda }_{k_i}\sim {\left({\displaystyle \frac{j_{{\nu }_i,k_i}}{|1-a_i|}}\right)}^2,{\nu }_i={\displaystyle \frac{{\theta }_i-1}{2-{\theta }_i}}.$$

Summing over i yields the asymptotic formula (9).

Compact resolvent and orthonormality guarantee that ${\left\{{\varphi }_k\right\}}_{k\in {\mathbb{N}}^d}$ forms a complete orthonormal basis of $L^2(\Omega )$, establishing the spectral decomposition. □

2.2. Neural Symmetrization via SDOs

Define SDO layers by

$$\left[{\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}\mathbf{u}\right]\left(\mathbf{x}\right)=\nabla ·\left({|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }}\circ \nabla \mathbf{u}\left(\mathbf{x}\right)\right),\mathbf{x}\in \Omega ,$$

(10)

with trainable $\mathbf{a}\in \Omega$ and $\mathbf{\theta }\in {[1,2)}^d$. Let $G(\mathbf{x},\mathbf{y};\mathbf{a},\mathbf{\theta })$ be the Green’s function:

$${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}{G(·,\mathbf{y})=\delta (·-\mathbf{y}),G|}_{\partial \Omega }=0.$$

(11)

Definition 2.2

(SDO-Net). AnSDO-Netcontains layers of the form

$${\mathbf{u}}_{l+1}=\sigma \left({\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}{\mathbf{u}}_l+{\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l\right),$$

(12)

with Lipschitz activation σ, weight ${\mathbb{W}}_l$, bias ${\mathbf{b}}_l$, and trainable ${\mathbf{a}}_l,{\mathbf{\theta }}_l$.

Lemma 2.3

(Weighted SDO Solve). Let $\mathbf{a}\in \mathrm{int}(\Omega )$ and $\mathbf{\theta }\in {[1,2)}^d$. For any $f\in L^2(\Omega )$, the boundary-value problem

$${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}{u=f,u|}_{\partial \Omega }=0,$$

(13)

admits a unique solution $u\in {\mathcal{H}}_{\mathbf{\theta }}^1(\Omega )$, where the weighted Sobolev space is defined as

$${\mathcal{H}}_{\mathbf{\theta }}^1(\Omega ):=\left\{v\in L^2(\Omega ):\sum _{i=1}^d{\int }_{\Omega }|x_i-a_i{|}^{{\theta }_i}|{\partial }_{x_i}v{{|}^{2}d\mathbf{x}<\infty ,v|}_{\partial \Omega }=0\right\}.$$

Moreover, there exists a constant $C_{\mathbf{\theta }}>0$ such that

$${\parallel u\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1(\Omega )}\le C_{\mathbf{\theta }}{\parallel f\parallel }_{L^2(\Omega )}.$$

(14)

Proof. 

Consider the bilinear form

$$a(u,v):={\int }_{\Omega }\sum _{i=1}^d{|x_i-a_i|}^{{\theta }_i}{\partial }_{x_i}u{\partial }_{x_i}vd\mathbf{x},u,v\in {\mathcal{H}}_{\mathbf{\theta }}^1(\Omega ).$$

The weak form of (13) is: find $u\in {\mathcal{H}}_{\mathbf{\theta }}^1(\Omega )$ such that

$$a(u,v)={\int }_{\Omega }fvd\mathbf{x},\forall v\in {\mathcal{H}}_{\mathbf{\theta }}^1(\Omega ).$$

(15)

By definition,

$$\left|a(u,v)\right|\le {\parallel u\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1}{\parallel v\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1}.$$

Also, $a(u,u)={\parallel u\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1}^2$ shows coercivity.

The Lax–Milgram theorem then guarantees a unique $u\in {\mathcal{H}}_{\mathbf{\theta }}^1(\Omega )$ solving (15).

From (15) with $v=u$ and Cauchy–Schwarz inequality:

$${\parallel u\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1}^2=a(u,u)={\int }_{\Omega }fud\mathbf{x}\le {\parallel f\parallel }_{L^2}{\parallel u\parallel }_{L^2}.$$

Using a weighted Poincaré inequality [1]:

$${\parallel u\parallel }_{L^2}\le C_P{\parallel u\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1},$$

we obtain

$${\parallel u\parallel }_{{\mathcal{H}}_{\mathbf{\theta }}^1}\le C_{\mathbf{\theta }}{\parallel f\parallel }_{L^2},$$

which proves (14).

Alternatively, using the spectral decomposition $\left\{({\lambda }_k,{\varphi }_k)\right\}$ from Theorem 2.1, one can write

$$u=\sum _{k\in {\mathbb{N}}^d}{\displaystyle \frac{{\langle f,{\varphi }_k\rangle }_{L^2}}{{\lambda }_k}}{\varphi }_k,$$

which also satisfies the same norm estimate due to orthonormality and positivity of ${\lambda }_k$. □

Theorem 2.4

(Well-Posedness of SDO Layers). Let ${\mathbf{u}}_l\in {\mathcal{H}}_{{\mathbf{\theta }}_l}^1(\Omega )$, and let the SDO layer be defined by

$${\mathbf{u}}_{l+1}=\sigma \left({\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}^{-1}({\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l)\right),$$

(16)

where σ is Lipschitz continuous with constant $L_{\sigma }$, ${\mathbb{W}}_l$ is a linear operator, and ${\mathbf{b}}_l\in L^2(\Omega )$. Then:

• Existence and uniqueness:  There exists a unique ${\mathbf{u}}_{l+1}\in {\mathcal{H}}_{{\mathbf{\theta }}_l}^1(\Omega )$ satisfying (16).
• Lipschitz bound:

  $$\parallel {\mathbf{u}}_{l+1}{\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\le L_{\sigma }{C}_{{\mathbf{\theta }}_l}\left(\parallel {\mathbf{u}}_l{\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}+{\parallel {\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l\parallel }_{L^2}\right),$$

  (17)
  where $C_{{\mathbf{\theta }}_l}$ is the constant from Lemma 2.3.
• Continuous dependence on parameters:

  $$\parallel {\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}^{-1}-{\mathcal{L}}_{{\tilde{\mathbf{a}}}_l,{\tilde{\mathbf{\theta }}}_l}^{-1}{\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^{-1}\to {\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\to 0\mathrm{as}{\mathbf{a}}_l\to {\tilde{\mathbf{a}}}_l,{\mathbf{\theta }}_l\to {\tilde{\mathbf{\theta }}}_l.$$

  (18)

Proof. 

From Lemma 2.3, for any $f\in L^2(\Omega )$, the problem

$${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}{\mathbf{u}=f,\mathbf{u}|}_{\partial \Omega }=0,$$

admits a unique solution $\mathbf{u}\in {\mathcal{H}}_{{\mathbf{\theta }}_l}^1(\Omega )$ with

$${\parallel \mathbf{u}\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\le C_{{\mathbf{\theta }}_l}{\parallel f\parallel }_{L^2}.$$

Setting $f={\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l$, we obtain existence and uniqueness of ${\mathbf{u}}_{l+1}$ before activation.

Applying the Lipschitz continuous activation $\sigma$ with constant $L_{\sigma }$ gives

$$\parallel {\mathbf{u}}_{l+1}{\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}=\parallel \sigma \left({\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}^{-1}({\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l)\right){\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\le L_{\sigma }{\parallel {\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}^{-1}({\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l)\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}.$$

Using Lemma 2.3:

$$\parallel {\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}^{-1}({\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l){\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\le C_{{\mathbf{\theta }}_l}{\parallel {\mathbb{W}}_l{\mathbf{u}}_l+{\mathbf{b}}_l\parallel }_{L^2},$$

which yields (17).

Let $\mathbf{u}={\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}^{-1}f$ and $\tilde{\mathbf{u}}={\mathcal{L}}_{{\tilde{\mathbf{a}}}_l,{\tilde{\mathbf{\theta }}}_l}^{-1}f$. Then

$${\mathcal{L}}_{{\tilde{\mathbf{a}}}_l,{\tilde{\mathbf{\theta }}}_l}(\mathbf{u}-\tilde{\mathbf{u}})=({\mathcal{L}}_{{\tilde{\mathbf{a}}}_l,{\tilde{\mathbf{\theta }}}_l}-{\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l})\mathbf{u}.$$

Taking the ${\mathcal{H}}_{{\mathbf{\theta }}_l}^1$ norm and applying the stability estimate from Lemma 2.3:

$$\parallel \mathbf{u}-\tilde{\mathbf{u}}{\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\le C_{{\mathbf{\theta }}_l}{\parallel ({\mathcal{L}}_{{\tilde{\mathbf{a}}}_l,{\tilde{\mathbf{\theta }}}_l}-{\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l})\mathbf{u}\parallel }_{L^2}.$$

Since ${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}$ depends continuously on ${\mathbf{a}}_l$ and ${\mathbf{\theta }}_l$, we have

$$\parallel \mathbf{u}-\tilde{\mathbf{u}}{\parallel }_{{\mathcal{H}}_{{\mathbf{\theta }}_l}^1}\to 0\mathrm{as}{\mathbf{a}}_l\to {\tilde{\mathbf{a}}}_l,{\mathbf{\theta }}_l\to {\tilde{\mathbf{\theta }}}_l,$$

which establishes (18). □

Lemma 2.5

(Spectral Decomposition of SDO Layer). Let ${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}$ be the SDO at layer l as in Theorem 2.1, and let ${\mathbf{u}}_l\in L^2(\Omega )$. Then ${\mathbf{u}}_l$ admits the spectral decomposition

$${\mathbf{u}}_l\left(\mathbf{x}\right)=\sum _{k\in {\mathbb{N}}^d}{\langle {\mathbf{u}}_l,{\varphi }_k\rangle }_{L^2(\Omega )}{\varphi }_k(\mathbf{x};{\mathbf{a}}_l,{\mathbf{\theta }}_l),$$

(19)

and the action of the operator satisfies

$${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}{\mathbf{u}}_l\left(\mathbf{x}\right)=\sum _{k\in {\mathbb{N}}^d}{\lambda }_k({\mathbf{a}}_l,{\mathbf{\theta }}_l){\langle {\mathbf{u}}_l,{\varphi }_k\rangle }_{L^2(\Omega )}{\varphi }_k(\mathbf{x};{\mathbf{a}}_l,{\mathbf{\theta }}_l),$$

(20)

where $({\lambda }_k,{\varphi }_k)$ are the eigenpairs of ${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}$ as in Theorem 2.1.

Proof. 

From Theorem 2.1, $-{\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}$ is self-adjoint, positive semi-definite, and has compact resolvent in $L^2(\Omega )$. Therefore, the spectral theorem for compact self-adjoint operators applies, ensuring the existence of a countable orthonormal basis of eigenfunctions ${\left\{{\varphi }_k\right\}}_{k\in {\mathbb{N}}^d}$ with corresponding eigenvalues ${\left\{{\lambda }_k\right\}}_{k\in {\mathbb{N}}^d}$.

Since $\left\{{\varphi }_k\right\}$ forms an orthonormal basis of $L^2(\Omega )$, any ${\mathbf{u}}_l\in L^2(\Omega )$ can be expanded as

$${\mathbf{u}}_l\left(\mathbf{x}\right)=\sum _{k\in {\mathbb{N}}^d}{\langle {\mathbf{u}}_l,{\varphi }_k\rangle }_{L^2(\Omega )}{\varphi }_k(\mathbf{x};{\mathbf{a}}_l,{\mathbf{\theta }}_l),$$

with convergence in $L^2(\Omega )$. This proves (19).

Applying ${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}$ to the expansion termwise gives

$$\begin{array}{cc}\hfill {\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}{\mathbf{u}}_l\left(\mathbf{x}\right)& =\sum _{k\in {\mathbb{N}}^d}{\langle {\mathbf{u}}_l,{\varphi }_k\rangle }_{L^2(\Omega )}{\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}{\varphi }_k(\mathbf{x};{\mathbf{a}}_l,{\mathbf{\theta }}_l)\hfill \\ \hfill & =\sum _{k\in {\mathbb{N}}^d}{\lambda }_k({\mathbf{a}}_l,{\mathbf{\theta }}_l){\langle {\mathbf{u}}_l,{\varphi }_k\rangle }_{L^2(\Omega )}{\varphi }_k(\mathbf{x};{\mathbf{a}}_l,{\mathbf{\theta }}_l).\hfill \end{array}$$

(21)

where we used the eigenvalue equation ${\mathcal{L}}_{{\mathbf{a}}_l,{\mathbf{\theta }}_l}{\varphi }_k={\lambda }_k{\varphi }_k$. Convergence is in $L^2(\Omega )$.

The decomposition allows efficient representation of SDO layers and is crucial for spectral analysis of deep networks. The dependence of ${\varphi }_k$ and ${\lambda }_k$ on ${\mathbf{a}}_l$ and ${\mathbf{\theta }}_l$ is continuous due to the parameter-dependence properties of SDOs (see Theorem 2.1 and Theorem 2.4). This completes the proof. □

Corollary 2.6

(Spectral Stability under Training). Let $\Omega \subset {\mathbb{R}}^d$ be a bounded Lipschitz domain. For each layer l, denote by

$${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}u:=-\nabla ·\left(D_{\mathbf{a},\mathbf{\theta }}\nabla u\right)\mathrm{with}{D}_{\mathbf{a},\mathbf{\theta }}\left(x\right)=diag\left(|x_1-a_1{|}^{{\theta }_1},\dots ,{|x_d-a_d|}^{{\theta }_d}\right),$$

acting on $H_0^1(\Omega )$, and let ${\left\{{\lambda }_k(\mathbf{a},\mathbf{\theta })\right\}}_{k\ge 1}$ be the eigenvalues enumerated in nondecreasing order with associated orthonormal eigenfunctions ${\left\{{\varphi }_k(·;\mathbf{a},\mathbf{\theta })\right\}}_{k\ge 1}$. Assume that between two training steps t and $t+1$ the parameters change by

$$\delta :={∥{\mathbf{a}}^{(t+1)}-{\mathbf{a}}^{\left(t\right)}∥}_{{\mathbb{R}}^d}+{∥{\mathbf{\theta }}^{(t+1)}-{\mathbf{\theta }}^{\left(t\right)}∥}_{{\mathbb{R}}^d}.$$

Then, for each fixed $k\in \mathbb{N}$,

$${\displaystyle \frac{{\lambda }_k\left({\mathbf{a}}^{(t+1)},{\mathbf{\theta }}^{(t+1)}\right)}{{\lambda }_k\left({\mathbf{a}}^{\left(t\right)},{\mathbf{\theta }}^{\left(t\right)}\right)}}=1+{\mathcal{O}}_k\left(\delta \right)\mathrm{as}\delta \to 0.$$

(22)

In particular, the spectrum of each mode remains stable under sufficiently small parameter updates during training.

Remark 2.7

(Physics-Informed Kernels and Modewise Filtering). Using the Green’s function representation, each layer admits the expansion

$${\mathbf{u}}_{l+1}\left(\mathbf{x}\right)=\sum _{k\in {\mathbb{N}}^d}{\displaystyle \frac{{\langle {\mathbf{u}}_l,{\varphi }_k(·;{\mathbf{a}}_l,{\mathbf{\theta }}_l)\rangle }_{L^2(\Omega )}}{{\lambda }_k({\mathbf{a}}_l,{\mathbf{\theta }}_l)}}{\varphi }_k(\mathbf{x};{\mathbf{a}}_l,{\mathbf{\theta }}_l).$$

(23)

Thus each SDO layer acts as a modewise filter whose anisotropy adapts to the learned parameters. When the eigenvalues vary smoothly with $({\mathbf{a}}_l,{\mathbf{\theta }}_l)$, the gains $1/{\lambda }_k$ of the dominant modes change only by $\mathcal{O}\left(\delta \right)$ for small updates, ensuring stable physics-informed filtering across layers.

3. Main Results

3.1. Universality for Turbulence Closure

Theorem 3.1

(Universality of SDO-Nets). Let $\Omega \subset {\mathbb{R}}^d$ be a bounded Lipschitz domain, and let $\mathcal{T}:L^2(\Omega ;{\mathbb{R}}^d)\to L^2(\Omega ;{\mathbb{R}}^{d\times d})$ be a turbulence closure operator mapping velocity fields to stress tensors. Then, for any compact set $\mathcal{U}\subset L^2(\Omega ;{\mathbb{R}}^d)$ and any $\epsilon >0$, there exists an SDO-Net ${\mathcal{T}}_{\mathrm{NN}}$ such that

$$\underset{\mathbf{u}\in \mathcal{U}}{sup}{∥\mathcal{T}\left(\mathbf{u}\right)-{\mathcal{T}}_{\mathrm{NN}}\left(\mathbf{u}\right)∥}_{L^2(\Omega )}<\epsilon ,$$

(24)

with ${\mathcal{T}}_{\mathrm{NN}}\left(\mathbf{u}\right)$ weakly divergence-free:

$$\nabla ·{\mathcal{T}}_{\mathrm{NN}}\left(\mathbf{u}\right)=0\mathrm{in}{\mathcal{D}}^{\prime }(\Omega ).$$

(25)

Proof. 

By Theorem 2.1, the set of eigenfunctions ${\left\{{\varphi }_k\right\}}_{k\in {\mathbb{N}}^d}$ of the spectral degeneracy operator ${\mathcal{L}}_{\mathbf{a},\mathbf{\theta }}$ forms a complete orthonormal basis of $L^2(\Omega )$. Thus, for each component ${\mathcal{T}}_{ij}\left(\mathbf{u}\right)$ of the closure tensor, and for any $\mathbf{u}\in \mathcal{U}$, we can approximate

$${\mathcal{T}}_{ij}\left(\mathbf{u}\right)\approx \sum _{\left|k\right|\le N}{c}_k^{\left(ij\right)}\left(\mathbf{u}\right){\varphi }_k(\mathbf{x};\mathbf{a},\mathbf{\theta }),$$

(26)

with coefficients $c_k^{\left(ij\right)}\left(\mathbf{u}\right)\in \mathbb{R}$ and error arbitrarily small in $L^2$ norm by choosing N sufficiently large.

Each coefficient map $\mathbf{u}\mapsto c_k^{\left(ij\right)}\left(\mathbf{u}\right)$ is a continuous functional on the compact set $\mathcal{U}\subset L^2(\Omega ;{\mathbb{R}}^d)$. By the **universal approximation theorem** for feedforward networks [17], for any $\delta >0$, there exists a neural network ${\mathcal{N}}_k^{\left(ij\right)}$ such that

$$\underset{\mathbf{u}\in \mathcal{U}}{sup}\left|c_k^{\left(ij\right)}\left(\mathbf{u}\right)-{\mathcal{N}}_k^{\left(ij\right)}\left(\mathbf{u}\right)\right|<\delta .$$

(27)

Let $\mathbf{F}\left(\mathbf{x}\right)\in L^2(\Omega ;{\mathbb{R}}^{d\times d})$ be any tensor field. By the Hodge decomposition (or Helmholtz projection) in bounded Lipschitz domains, there exists a unique decomposition

$$\mathbf{F}={\mathbf{F}}^{\mathrm{div}}+\nabla \blacksquare ,$$

(28)

with $\nabla ·{\mathbf{F}}^{\mathrm{div}}=0$ and potential $\blacksquare$ vanishing on $\partial \Omega$. Applying this projection to the neural approximation of the spectral expansion (26), we obtain a divergence-free SDO-Net:

$${\mathcal{T}}_{\mathrm{NN}}\left(\mathbf{u}\right)={\mathbb{P}}_{\mathrm{div}}\sum _{\left|k\right|\le N}{\mathcal{N}}_k\left(\mathbf{u}\right){\varphi }_k\left(\mathbf{x}\right),$$

(29)

where ${\mathbb{P}}_{\mathrm{div}}$ denotes the $L^2$-orthogonal projection onto divergence-free tensor fields.

By triangle inequality,

$$\begin{array}{cc}\hfill {∥\mathcal{T}\left(\mathbf{u}\right)-{\mathcal{T}}_{\mathrm{NN}}\left(\mathbf{u}\right)∥}_{L^2}& \le {∥\mathcal{T}\left(\mathbf{u}\right)-\sum _{\left|k\right|\le N}{c}_k\left(\mathbf{u}\right){\varphi }_k∥}_{L^2}\hfill \\ \hfill & +{∥\sum _{\left|k\right|\le N}(c_k\left(\mathbf{u}\right)-{\mathcal{N}}_k\left(\mathbf{u}\right)){\varphi }_k∥}_{L^2}\hfill \\ \hfill & +{∥\sum _{\left|k\right|\le N}{\mathcal{N}}_k\left(\mathbf{u}\right){\varphi }_k-{\mathbb{P}}_{\mathrm{div}}\sum _{\left|k\right|\le N}{\mathcal{N}}_k\left(\mathbf{u}\right){\varphi }_k∥}_{L^2}.\hfill \end{array}$$

(30)

Each term can be made smaller than $\epsilon /3$ by choosing sufficiently large N and accurate neural networks ${\mathcal{N}}_k$. Hence ${\mathcal{T}}_{\mathrm{NN}}$ satisfies (24) and is weakly divergence-free by construction (29).

□

Remark 3.2

(Interpretation). This theorem shows that **SDO-Nets provide a universal approximator** for turbulence closure operators while exactly enforcing incompressibility. The spectral degeneracy basis allows adaptive representation of anisotropic and localized structures, critical for turbulent flows.

3.2. Inverse Calibration of Degeneracy Points

Theorem 3.3

(Lipschitz Stability of Degeneracy Points). Let $\Omega \subset {\mathbb{R}}^d$ be a bounded Lipschitz domain, and let ${\mathbf{u}}_1,{\mathbf{u}}_2$ be solutions to the degenerate Navier–Stokes system

$${\partial }_t\mathbf{u}+(\mathbf{u}·\nabla )\mathbf{u}+\nabla p-\nabla ·\left({|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }}\nabla \mathbf{u}\right)=\mathbf{f},\nabla ·\mathbf{u}=0,$$

(31)

with degeneracy points ${\mathbf{a}}_1,{\mathbf{a}}_2\in \Omega$ and identical initial/boundary conditions.

Assume that boundary measurements satisfy

$${∥\left(|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}\nabla {\mathbf{u}}_1-{|\mathbf{x}-{\mathbf{a}}_2|}^{\mathbf{\theta }}\nabla {\mathbf{u}}_2\right)·\mathbf{n}∥}_{L^2(\partial \Omega \times (0,T))}\le \delta .$$

(32)

Then there exist constants $C>0$ and $\gamma \in (0,1]$, depending on Ω, $\mathbf{\theta }$, and the spectral gap of the SDO, such that

$${∥{\mathbf{a}}_1-{\mathbf{a}}_2∥}_{{\mathbb{R}}^d}\le C{\delta }^{\gamma }.$$

(33)

Proof. 

We outline a rigorous four-step argument:

Set $\mathbf{w}={\mathbf{u}}_1-{\mathbf{u}}_2$. Subtracting the two PDEs (31) gives

$${\partial }_t\mathbf{w}+({\mathbf{u}}_1·\nabla )\mathbf{w}+(\mathbf{w}·\nabla ){\mathbf{u}}_2+\nabla (p_1-p_2)-\nabla ·\left(|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}\nabla \mathbf{w}\right)=\nabla ·\left(\left(\right|\mathbf{x}-{\mathbf{a}}_2{|}^{\mathbf{\theta }}-|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }})\nabla {\mathbf{u}}_2\right).$$

(34)

Multiply (34) by $\mathbf{w}$ and integrate over $\Omega$:

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{∥\mathbf{w}∥}_{L^2}^2+{\int }_{\Omega }|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}{|\nabla \mathbf{w}|}^{2}dx\le {\int }_{\Omega }\left|\right(|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}-|\mathbf{x}-{\mathbf{a}}_2{|}^{\mathbf{\theta }}\left)\right||\nabla {\mathbf{u}}_2\left|\right|\nabla \mathbf{w}|dx+\mathcal{N},$$

(35)

where $\mathcal{N}={\int }_{\Omega }|(\mathbf{w}·\nabla ){\mathbf{u}}_2·\mathbf{w}|dx$ is the nonlinear advection term.

Apply Cauchy–Schwarz and Young inequalities:

$${\int }_{\Omega }\left|\right(|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}-|\mathbf{x}-{\mathbf{a}}_2{|}^{\mathbf{\theta }}\left)\right||\nabla {\mathbf{u}}_2\left|\right|\nabla \mathbf{w}|dx\le ϵ{∥\nabla \mathbf{w}∥}_{L^2\left(\right|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}dx)}^2+C_{ϵ}{∥{\mathbf{a}}_1-{\mathbf{a}}_2∥}_{{\mathbb{R}}^d}^{2\theta }{∥\nabla {\mathbf{u}}_2∥}_{L^2}^2.$$

(36)

By Theorem 2.1, the SDO has a positive first eigenvalue ${\lambda }_1>0$:

$${\int }_{\Omega }|\mathbf{x}-{\mathbf{a}}_1{|}^{\mathbf{\theta }}{|\nabla \mathbf{w}|}^{2}dx\ge {\lambda }_1{∥\mathbf{w}∥}_{L^2}^2.$$

(37)

Combining (35)–(37) and bounding the nonlinear term via ${∥\nabla {\mathbf{u}}_2∥}_{L^{\infty }}{∥\mathbf{w}∥}_{L^2}^2$, we obtain a Grönwall inequality:

$${\displaystyle \frac{d}{dt}}{∥\mathbf{w}∥}_{L^2}^2+{\lambda }_1{∥\mathbf{w}∥}_{L^2}^2\le C{∥{\mathbf{a}}_1-{\mathbf{a}}_2∥}^{2\theta }+C{∥\mathbf{w}∥}_{L^2}^2.$$

(38)

Using the boundary measurement condition (32) and standard elliptic estimates (or Carleman inequalities for SDOs [1]), we can control the interior norm ${∥\mathbf{w}∥}_{L^2(\Omega )}$ in terms of $\delta$. Thus, solving (38) gives

$${∥{\mathbf{a}}_1-{\mathbf{a}}_2∥}_{{\mathbb{R}}^d}\le C{\delta }^{\gamma },$$

(39)

for some $\gamma \in (0,1]$, completing the proof. □

Remark 3.4

(Interpretation). This result guarantees Lipschitz-type stability for inverse calibration of degeneracy points: small changes in boundary fluxes lead to controlled shifts in the inferred centers ${\mathbf{a}}_l$. It forms the theoretical foundation for learning adaptive attention centers in SDO-Nets from partial or boundary measurements.

3.3. Neural-Turbulence Correspondence

Theorem 3.5

(Neural-Turbulence Correspondence – Refined). Let ${\mathcal{T}}_{\mathrm{NN}}:L^2(\Omega ;{\mathbb{R}}^d)\to L^2(\Omega ;{\mathbb{R}}^{d\times d})$ be an SDO-Net trained to minimize the residual energy functional

$${\mathcal{E}}_N\left({\mathcal{T}}_{\mathrm{NN}}\right)={∥{\partial }_t\overline{\mathbf{u}}+(\overline{\mathbf{u}}·\nabla )\overline{\mathbf{u}}+\nabla \overline{p}{-\nabla ·\left(\right|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }}\nabla \overline{\mathbf{u}})-\nabla ·{\mathcal{T}}_{\mathrm{NN}}\left(\overline{\mathbf{u}}\right)∥}_{L^2(\Omega )}.$$

(40)

Assume:

• The dataset ${\left\{{\overline{\mathbf{u}}}_i\right\}}_{i=1}^N$ is dense in the function space of resolved velocities as $N\to \infty$.
• The SDO-Net satisfies the Lipschitz stability property from Theorem 3.3.
• The loss functional ${\mathcal{E}}_N$ is equi-coercive and lower semicontinuous with respect to the degeneracy points ${\mathbf{a}}_N$.

Then, as $N\to \infty$, the learned degeneracy points ${\mathbf{a}}_N$ converge to the true turbulence structures ${\mathbf{a}}^{*}$ in $L^1$:

$$\underset{N\to \infty }{lim}{∥{\mathbf{a}}_N-{\mathbf{a}}^{*}∥}_{L^1(\Omega )}=0.$$

(41)

Proof. 

The proof is based on a three-step argument combining consistency, stability, and variational convergence:

By the Germano identity [18], the true subgrid stress satisfies

$$\nabla ·{\mathcal{T}}^{*}\left(\overline{\mathbf{u}}\right)=\underset{N\to \infty }{lim}\nabla ·{\mathcal{T}}_{\mathrm{NN}}\left({\overline{\mathbf{u}}}_N\right),$$

(42)

ensuring that the SDO-Net can approximate the exact residual as $N\to \infty$.

Applying Theorem 3.3 to the mapping

$$\mathbf{a}\mapsto {\nabla ·\left(\right|\mathbf{x}-\mathbf{a}|}^{\mathbf{\theta }}\nabla \overline{\mathbf{u}})+\nabla ·{\mathcal{T}}_{\mathrm{NN}}\left(\overline{\mathbf{u}}\right)$$

guarantees that small residual errors in ${\mathcal{E}}_N$ induce controlled deviations in the inferred ${\mathbf{a}}_N$:

$${∥{\mathbf{a}}_N-{\mathbf{a}}^{*}∥}_{{\mathbb{R}}^d}\le C{\mathcal{E}}_N{\left({\mathcal{T}}_{\mathrm{NN}}\right)}^{\gamma }.$$

(43)

Let ${\mathcal{E}}_{\infty }\left(\mathbf{a}\right)={lim}_{N\to \infty }{\mathcal{E}}_N\left(\mathbf{a}\right)$. By equi-coercivity and lower semicontinuity,

$${\mathbf{a}}_N\to \underset{\mathbf{a}}{arg\; min}{\mathcal{E}}_{\infty }\left(\mathbf{a}\right)={\mathbf{a}}^{*}.$$

(44)

Combining (43) and (44) yields

$${∥{\mathbf{a}}_N-{\mathbf{a}}^{*}∥}_{L^1(\Omega )}\le C{\mathcal{E}}_N{\left({\mathcal{T}}_{\mathrm{NN}}\right)}^{\gamma }\to 0\mathrm{as}N\to \infty .$$

(45)

Hence, the SDO-Net learns degeneracy points that asymptotically align with true turbulent structures. □

Remark 3.6

(Practical Implication). The theorem rigorously justifies the use of trainable degeneracy points in SDO-Nets: during large-data training, these points act as learned attention centers, automatically localizing key turbulence structures in a physically consistent way.

4. Results

We summarize here the main theoretical and computational contributions developed in the preceding sections:

• Spectral theory of SDOs. We proved self-adjointness, compact resolvent, and a complete tensor-product basis of eigenfunctions with Bessel-type asymptotics, enabling efficient spectral representations of SDO layers.
• Stability and inverse calibration. Using min–max principles and Carleman-type estimates we derived Lipschitz-type stability bounds for recovering degeneracy points from boundary measurements in degenerate Navier–Stokes systems.
• Universal approximation for turbulence closure. We established that SDO-Nets approximate arbitrary divergence-free closure operators on compact subsets of $L^2$, combining spectral expansions with neural coefficient maps.
• Neural–turbulence correspondence. We proved convergence of trainable degeneracy points to true turbulent structures as the amount of training data increases, rigorously justifying their interpretation as learned attention centers.

These results jointly demonstrate that SDO-based architectures retain the expressive power of modern neural operators while inheriting structural guarantees from the underlying PDEs.

5. Conclusions

We have introduced and analyzed Spectral Degeneracy Operators as a bridge between degenerate PDE theory and physics-informed neural networks. Our analysis shows that SDO layers act as stable, adaptive, and interpretable spectral filters capable of representing anisotropic and intermittent features of turbulence without violating conservation laws. The inverse stability results provide a theoretical foundation for data-driven calibration of degeneracy points, and the universality theorem establishes the approximation power of SDO-Nets for turbulence closures. Future work will extend these ideas to fully three-dimensional, time-dependent turbulent flows, incorporate stochastic parameterizations, and explore SDO-based architectures in other domains such as geophysical fluid dynamics and nuclear-reactor thermo-hydraulics.