Zeitlin Discretization in Fluid Dynamics

• Zeitlin discretization is a structure-preserving finite-dimensional approximation of the incompressible Euler equations that replaces the Poisson algebra with matrix Lie algebras.
• It employs sine‐bracket spectral truncation on the torus and Hoppe’s quantization on the sphere to mirror continuous Hamiltonian, Casimir, and geometric structures.
• By exactly preserving invariants, curvature, and stability properties through discrete Nambu brackets and Lie–Poisson formulations, it offers a robust framework for numerical fluid dynamics.

Zeitlin discretization is a structure-preserving finite-dimensional approximation of the incompressible Euler equations in which the Poisson algebra of functions is replaced by a sequence of matrix Lie algebras, typically $\mathfrak{su}(n)$ or $\mathfrak{u}(N)$, so that the discrete dynamics remains Hamiltonian, Lie–Poisson, and isospectral. In the flat torus setting it takes the form of a sine-bracket spectral truncation, while on the sphere it is realized through Hoppe’s quantization and a discrete Laplacian with the same spectral data as the Laplace–Beltrami operator. Its distinctive feature is that it discretizes not only the PDE but also Arnold’s geometric formulation of ideal hydrodynamics as geodesic flow on a Lie group, thereby preserving discrete energy, discrete Casimirs, coadjoint-orbit geometry, and, in several settings, curvature and Jacobi structures (Sommer et al., 2011, Modin et al., 2023).

1. Continuum origin and geometric motivation

The underlying continuum problem is the incompressible Euler equation viewed in vorticity form on a two-dimensional manifold. On the periodic torus $T^2=[0,2\pi]^2$, the barotropic vorticity equation is

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$

with

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$

The corresponding velocity is

$$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$

This dynamics can be regarded as motion on the dual of the Lie algebra of area-preserving diffeomorphisms, and it admits the Lie–Poisson bracket

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$

with Hamiltonian

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$

This structure implies infinitely many Casimirs $C[\zeta]=\int f(\zeta)\,dA$, in particular the enstrophy

$\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$

On the sphere, the same Eulerian dynamics is formulated using scalar vorticity $\mathfrak{u}(N)$0, stream function $\mathfrak{u}(N)$1, and the Poisson bracket on $\mathfrak{u}(N)$2. In Arnold’s theorem, incompressible Euler flow is geodesic flow on the group of volume-preserving diffeomorphisms with respect to a right-invariant metric. The Lie–Poisson description on vorticity space likewise yields energy and an infinite family of Casimirs. Zeitlin’s construction is motivated by the desire to preserve this full geometric package at the semi-discrete level rather than merely approximate the PDE coefficients (Modin et al., 2023).

A structure-preserving discretization in this sense keeps a Lie–Poisson bracket, conserves discrete energy and discrete analogues of Casimirs exactly, and allows algebraic constructions, such as the Nambu bracket for the two-dimensional vorticity equation, to be carried out in finite dimensions before passing to the continuum limit (Sommer et al., 2011).

2. Torus realization: sine-bracket truncation and $\mathfrak{u}(N)$3 closure

In the torus formulation, vorticity is expanded in Fourier modes,

$\mathfrak{u}(N)$4

For the full spectral system, the structure constants are

$\mathfrak{u}(N)$5

Zeitlin’s idea, following Fairlie–Zachos–Hoppe, is to replace the infinite-dimensional Lie algebra $\mathfrak{u}(N)$6 by a sequence of finite-dimensional Lie algebras $\mathfrak{u}(N)$7, approximate the Poisson bracket of Fourier modes by a sine bracket, and use an odd integer $\mathfrak{u}(N)$8 with the rectangular lattice cutoff

$\mathfrak{u}(N)$9

which yields $T^2=[0,2\pi]^2$0 modes (Sommer et al., 2011).

The discrete vorticity dynamics is

$T^2=[0,2\pi]^2$1

with trigonometric structure constants

$T^2=[0,2\pi]^2$2

Here $T^2=[0,2\pi]^2$3 enforces equality modulo $T^2=[0,2\pi]^2$4 in each component. The trigonometric factor

$T^2=[0,2\pi]^2$5

is what turns the truncated mode set into a closed Lie algebra isomorphic to $T^2=[0,2\pi]^2$6. In the limit $T^2=[0,2\pi]^2$7, the sine recovers the linear cross product, and the discrete structure constants converge to the continuum ones,

$T^2=[0,2\pi]^2$8

The corresponding finite-dimensional Lie–Poisson bracket is

$T^2=[0,2\pi]^2$9

which is bilinear, antisymmetric, and satisfies the Jacobi identity. The truncated vorticity equations can therefore be written as

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$0

The discrete Hamiltonian is the truncated kinetic energy,

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$1

and the quadratic Casimir coincides, up to normalization, with the truncated enstrophy,

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$2

Thus

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$3

The discrete system has $$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$4 independent Casimirs, but the quadratic one is central in the Nambu construction and in the correspondence with continuum enstrophy (Sommer et al., 2011).

This formulation distinguishes Zeitlin discretization from naive spectral truncation. In the latter, the truncated Poisson bracket may fail Jacobi, or Casimirs may not be preserved exactly; by contrast, the sine-bracket construction yields a genuine finite-dimensional Lie–Poisson system with exact invariant conservation (Sommer et al., 2011).

3. Sphere formulation: quantization, matrix hydrodynamics, and Euler–Arnold structure

On $$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$5, Zeitlin’s model is formulated through Hoppe’s explicit quantization of the sphere. Let $$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$6 be spherical harmonics with

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$7

Hoppe constructs matrices $$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$8 such that

$$\partial_t \zeta + J(\psi,\zeta)=0,\qquad \zeta=\Delta\psi,$$9

form an orthonormal basis of $$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$0 with respect to

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$1

and the truncated Fourier–Toeplitz quantization map is

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$2

The Lie algebra structure is scaled by

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$3

Bordemann–Meinrenken–Schlichenmaier show that this family is an $$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$4-approximation of the Poisson algebra, and Charles–Polterovich give the estimate

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$5

Real-valued functions correspond to skew-Hermitian matrices, and zero-mean functions correspond to traceless skew-Hermitian matrices $$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$6 (Modin et al., 2023).

The sphere model introduces a quantized Laplacian

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$7

so that, up to truncation, $$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$8 shares eigenvalues with the Laplace–Beltrami operator. There is also the Hoppe–Yau representation-theoretic identity

$$J(a,b):=\frac{\partial a}{\partial x_1}\frac{\partial b}{\partial x_2} -\frac{\partial a}{\partial x_2}\frac{\partial b}{\partial x_1}.$$9

The Euler–Zeitlin equations on $$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$0 are

$$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$1

where $$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$2 is the vorticity matrix and $$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$3 is the stream matrix. This is a finite-dimensional Lie–Poisson system on $$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$4 with Hamiltonian

$$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$5

The corresponding right-invariant Riemannian metric on $$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$6 is induced by

$$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$7

and the reconstruction equation is

$$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$8

Accordingly, the sphere model is not only a discretization of the vorticity equation but also a finite-dimensional realization of Arnold’s geodesic picture (Modin et al., 2023).

A parallel matrix formulation is used in later work on the sphere. There, for each integer $$\mathbf{v}=\left(-\partial_{x_2}\psi,\ \partial_{x_1}\psi\right).$$9,

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$0

with scaled Frobenius inner product

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$1

The Hamiltonian is

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$2

and the system is again an Euler–Arnold equation: geodesic flow on $$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$3 with a right-invariant metric (Melzi et al., 11 Mar 2026).

4. Invariants, Nambu structure, and exact discrete geometry

A defining property of Zeitlin discretization is exact preservation of invariants associated with Lie–Poisson geometry. On the torus, conservation of the Hamiltonian and the quadratic Casimir follows directly from the discrete bracket. On the sphere, the Euler–Zeitlin system is isospectral: the eigenvalues of $$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$4 are constant in time. This yields discrete Casimirs

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$5

which converge as $$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$6 to the continuous Casimirs $$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$7. In the sphere formulation used for Arnold stability, the conserved quantities are

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$8

and in particular

$$\{F_1,F_2\}_{\text{PDE}} := \int_{T^2} dA\, \zeta\, J\!\left(\frac{\delta F_1}{\delta\zeta}, \frac{\delta F_2}{\delta\zeta}\right),$$9

is the squared Frobenius norm, the discrete analogue of enstrophy. On $$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$0 there is also an $$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$1 symmetry, and the discrete angular momentum

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$2

is conserved; equivalently, the projection of vorticity onto the first eigenspace is frozen (Modin et al., 2023, Melzi et al., 11 Mar 2026).

In the torus case, the semi-simple Lie algebra structure allows an explicit algebraic construction of a Nambu bracket. For a semi-simple Lie algebra with structure constants $$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$3 and Killing form $$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$4, define

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$5

Applying this to the Zeitlin algebra yields the scaled Nambu tensor

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$6

and the discrete Nambu bracket

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$7

It is trilinear, totally antisymmetric, and satisfies the Leibniz rule. Most importantly,

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$8

so the discrete vorticity dynamics may be written in Nambu form,

$$H=\frac12\int_{T^2}dA\,|\nabla\psi|^2 =\frac12\int_{T^2}dA\,\zeta\,\psi.$$9

Under convergence of the discrete functionals, this bracket converges to the continuum Nambu bracket

$C[\zeta]=\int f(\zeta)\,dA$0

and inserting the enstrophy as the third argument reproduces the continuous Lie–Poisson bracket (Sommer et al., 2011).

A recurrent misconception is that the existence of a Nambu representation here implies the full generalized Jacobi identity. It does not: the resulting Nambu bracket, both discrete and continuous, does not satisfy Takhtajan’s identity, and explicit counterexamples are given. Its significance is instead algebraic and structural: the bracket arises algorithmically from a structure-preserving finite-dimensional approximation and recovers the Lie–Poisson formulation through the quadratic Casimir (Sommer et al., 2011).

Exact structure preservation also underlies recent model-reduction work. A time-dependent low-rank factorization of the vorticity matrix on the sphere keeps the approximate flow isospectral and Lie–Poisson, and the error in the solution, in the approximation of the Hamiltonian, and of the Casimir functions only depends on the approximation of the vorticity matrix at the initial time (Pagliantini, 2024).

5. Curvature, Jacobi equations, stability, and rigidity

One of the main reasons Zeitlin discretization is used in geometric hydrodynamics is that it carries over Riemannian structures that are absent in standard discretizations. On $C[\zeta]=\int f(\zeta)\,dA$1, sectional curvature and Jacobi fields encode Lagrangian and Eulerian stability. In the discrete model on $C[\zeta]=\int f(\zeta)\,dA$2, the same right-invariant-metric framework applies, and Arnold’s curvature formula produces an exact matrix analogue of the continuum curvature expression: $C[\zeta]=\int f(\zeta)\,dA$3 For $C[\zeta]=\int f(\zeta)\,dA$4, the corresponding curvature converges linearly in $C[\zeta]=\int f(\zeta)\,dA$5: $C[\zeta]=\int f(\zeta)\,dA$6 Thus, for sufficiently large $C[\zeta]=\int f(\zeta)\,dA$7, discrete and continuous curvatures have the same sign (Modin et al., 2023).

The Jacobi equation also admits a discrete analogue. For stationary base flows, the split Jacobi equations in the matrix model are

$C[\zeta]=\int f(\zeta)\,dA$8

For matching initial data, the embedded discrete Jacobi fields converge in $C[\zeta]=\int f(\zeta)\,dA$9 to the continuous Jacobi fields, uniformly on bounded time intervals. This provides a direct link between Eulerian and Lagrangian stability in the Euler equations and in the finite-dimensional matrix model (Modin et al., 2023).

Arnold’s nonlinear stability method has also been carried over to the matrix setting on $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$0. For a steady state $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$1 satisfying

$\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$2

the second variation is expressed by

$\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$3

A key sub-diagonal formula gives

$\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$4

where $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$5 and $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$6 are simultaneously diagonalized by $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$7. Using spectral-slope bounds and the fact that the operator norm of $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$8 on the subspace orthogonal to the first spherical harmonic eigenspace is $\mathcal{E}=\frac12\int_{T^2}dA\,\zeta^2.$9, one obtains the principal Lyapunov criterion: if

$\mathfrak{u}(N)$00

then the steady state is Lyapunov stable in the Frobenius norm. In the functional subclass

$\mathfrak{u}(N)$01

this holds whenever $\mathfrak{u}(N)$02 everywhere. The result mirrors the Constantin–Germain threshold for Euler on $\mathfrak{u}(N)$03 (Melzi et al., 11 Mar 2026).

The same paper proves rigidity statements. Under the same $\mathfrak{u}(N)$04 spectral-slope condition, there exists a rotation $\mathfrak{u}(N)$05 such that $\mathfrak{u}(N)$06 is diagonal. Under the stronger bound

$\mathfrak{u}(N)$07

necessarily $\mathfrak{u}(N)$08. These results are the matrix analogues of the continuum zonal-flow rigidity picture: Arnold-stable steady states are highly constrained, and the stronger lower bound eliminates all nontrivial steady states (Melzi et al., 11 Mar 2026).

6. Extensions, reduction strategies, and scope

A major limitation of the classical formulation is domain dependence. On the torus, the construction uses a planar Fourier basis and a sine bracket. On the sphere, it depends on spherical harmonics, Toeplitz quantization, and the Hoppe–Yau Laplacian. Earlier work emphasized that similar structure-preserving truncations were not readily available for three-dimensional Euler or shallow-water equations. A recent extension addresses part of that gap: for axisymmetric solutions of the Euler equations on $\mathfrak{u}(N)$09, Zeitlin’s approach has been extended to a finite-dimensional Euler–Arnold system on an Abelian extension of matrix Lie algebras, providing the first discretization of the 3-D Euler equations that fully preserves the geometric structure, albeit restricted to axisymmetric solutions (Modin et al., 2024).

In that axisymmetric $\mathfrak{u}(N)$10 setting, the reduced continuum system on $\mathfrak{u}(N)$11 is

$\mathfrak{u}(N)$12

Its discrete analogue is formulated on $\mathfrak{u}(N)$13 with bracket

$\mathfrak{u}(N)$14

metric

$\mathfrak{u}(N)$15

and Euler–Arnold equations

$\mathfrak{u}(N)$16

The model preserves the Hamiltonian

$\mathfrak{u}(N)$17

and discrete Casimirs

$\mathfrak{u}(N)$18

for any real analytic $\mathfrak{u}(N)$19. Because the semi-discrete system is finite-dimensional, curvature and Jacobi equations are again well-defined, and explicit Ricci and conjugate-point calculations are available in low-dimensional cases (Modin et al., 2024).

At the computational level, geometric low-rank approximation has been developed for the sphere model. If

$\mathfrak{u}(N)$20

then evolving

$\mathfrak{u}(N)$21

produces a rank-$\mathfrak{u}(N)$22 approximation that remains skew-Hermitian, isospectral, and Lie–Poisson. The exact discrete Casimir error is

$\mathfrak{u}(N)$23

and the Hamiltonian error is likewise frozen in time. The computational complexity of solving the approximate model scales quadratically with the order of the vorticity matrix and linearly if a further approximation of the stream function is introduced; in the full model the cost per time step is $\mathfrak{u}(N)$24, while the low-rank and truncated variants reduce this to $\mathfrak{u}(N)$25-type and $\mathfrak{u}(N)$26-type regimes under the assumptions stated in the paper (Pagliantini, 2024).

These developments clarify the scope of Zeitlin discretization. It is not a generic-purpose discretization for arbitrary fluid PDEs, and its implementation is algebraically more involved than simple spectral truncation. Its strength lies elsewhere: it is the only known discretization that preserves the rich geometric structure in the two-dimensional setting considered on the sphere, and it has proved sufficiently robust to support convergence of curvature, convergence of Jacobi equations, Lyapunov stability theory, rigidity results, Nambu constructions, and structure-preserving low-rank reduction within a single finite-dimensional framework (Modin et al., 2023).