\section{Introduction}
 Euler's equations on the sphere $\amsmathbb{S}^2$ describe the evolution of a homogeneous, incompressible two-dimensional fluid  constrained to move on the sphere.
The equations are given by
\begin{equation}
    \label{eq:euler}
    \begin{gathered}
        \dot v + v\cdot \nabla v = -\nabla p,\\
        \operatorname{div} v = 0
    \end{gathered}
\end{equation}
where $v$ is a divergence-free vector field for the fluid velocity and $p$ is a pressure potential.

\citet{arnold1966} showed that \Cref{eq:euler} has a remarkable geometric structure: it is a geodesic equation on the manifold of volume-preserving diffeomorphisms $\operatorname{Diff}_\mu(\amsmathbb{S}^2)$ equipped with the $\mathrm{L}^2$ metric.
Euler's equations are one among many arising this way.
By varying the manifold and the metric, a large class of equations from many different areas of physics and mathematics can be obtained, for instance the rigid body equations, the KdV equation and other shallow water equations, or equations in shape analysis (see, for example, the monograph by \citet{arnold1998} for an overview).

\Cref{eq:euler} can be rewritten as the \emph{vorticity equation},
\begin{align}
    \label{eq:vorteul}
    \begin{gathered}
        \dot \omega = \{\psi,\omega\},\\
        \psi = \Delta^{-1} \omega.
    \end{gathered}
\end{align}
Here $\omega$ is a smooth function called the \emph{vorticity function} and $\psi $ is the \emph{stream function}.
The velocity is recovered as $v = \nabla^\perp \psi $ where $\nabla^\perp$ is the skew-gradient.
The bracket $\{\cdot,\cdot\}$ denotes the Poisson bracket of smooth functions and the system \eqref{eq:vorteul} is an example of a \emph{Lie--Poisson system}.
See \citet{Marsden1999} and the references therein for more details on systems of this type.

A spatial discretization of \Cref{eq:vorteul} was introduced by \citet{zeitlin1991, zeitlin2004}.
This discretization, often called Zeitlin's model, builds on quantization theory, see \citet{hoppe1989,bordemann1991, bordemann1994,le2018}.

Zeitlin's model is a discretization preserving the underlying Lie--Poisson structure of \Cref{eq:vorteul}, resulting in a finite-dimensional isospectral flow on matrices,
\begin{equation*}
    \begin{gathered}
        \dot W = [P,W],\\
        P =  \Delta_N^{-1} W.
    \end{gathered}
\end{equation*}
Here, $W $ is a skew-Hermitian  ``vorticity matrix'' and  $P$ is a skew-Hermitian ``stream matrix''.
The operator $\Delta_N$ is the \emph{Hoppe--Yau Laplacian} \citep{hoppe1998}.
For more details on the structure of the equations and the quantization scheme, see \citet{zeitlin1991,zeitlin2004}

Zeitlin's model can be used to discretize fluid-type equations described by divergence-free vector fields.
A natural question is if this approach can be extended to equations involving general smooth vector fields.
As an example, consider the \emph{EPDiff equation}, which originates from the geodesic equation on $\operatorname{Diff}(\amsmathbb{S}^2) $ equipped with a right-invariant $\mathrm{H}^1$ metric.
It is given by
\begin{equation}
   \label{eq:eulerarnold}
    \begin{aligned}

       &\dot m =  \nabla_v^T m -\nabla_m v +\operatorname{div}(v)m,\\
       &m = (1-\Delta) v.
   \end{aligned}
\end{equation}
Here $v $ is a vector field, $\nabla$ is the covariant derivative and $\nabla^T$ its transpose with respect to $\mathrm{L}^2$.

The vector field $v$ in \Cref{eq:eulerarnold} is analogous to the velocity field $v$ in Euler's equations \eqref{eq:euler}.
An important difference, however, is that for EPDiff the vector field $v$ is not divergence-free.
The goal of this paper is to extend the quantization based discretization method of Zeitlin from divergence-free vector fields, such as in \Cref{eq:euler}, to generic vector fields, such as in \Cref{eq:eulerarnold}.

The paper is structured as follows.
In \Cref{sec:theory} we outline the theory of the discretization scheme and provide its extension to the manifold of all diffeomorphisms on the sphere.
Two computational examples are performed in \Cref{sec:computations}.
In both examples, a known vector field is quantized, and the same computations are carried out in both the non-quantized and quantized cases, and the results are qualitatively compared.
Finally, in \Cref{sec:outlook}, we discuss future directions, specifically settings for which the method presented in this paper could prove useful.

\medskip

\noindent\textbf{Acknowledgements.}
This work was supported by the Wallenberg AI,
Autonomous Systems and Software Program (WASP) funded by the Knut and Alice
Wallenberg Foundation. This work was also supported by the Swedish Research Council, grant number 2017-05040, and the Knut and Alice Wallenberg Foundation, grant number WAF2019.0201.
The authors would like to thank Paolo Cifani, Erwin Luesink and  Milo Viviani for fruitful discussions related to this work.

\section{Quantization on the sphere}
    \label{sec:theory}
    In this section we describe the quantization method.
    We begin with the geometric setting.
    The $2$-sphere $\amsmathbb{S}^2$ is a Kähler manifold: its area-form is a symplectic form $\Omega$ compatible with the standard round Riemannian metric $g$ via an integrable almost complex structure $J$, meaning that $J$ is a smooth field of automorphisms  of $T\amsmathbb{S}^2$ such that $J^2 = -I$ and $g(\cdot,\cdot) = -\Omega(J\cdot,\cdot)$.
    For more details on Kähler manifolds, see for instance \citet{silva2008}

	The gradient $\nabla f$ of a smooth function $f \in C^\infty(M)$ is defined by
	\begin{align*}
	    \,\mathrm{d} f(\cdot) =  g(\nabla f,\cdot),
	\end{align*}
	where $\,\mathrm{d}$ is the exterior derivative.
	The Hamiltonian vector field $X_f$ of $f$ is similarly defined by the relation
	\begin{align*}
	    \,\mathrm{d} f(\cdot) =\Omega(X_f,\cdot).
	\end{align*}
	Since
	\begin{align*}
	    \,\mathrm{d} f(\cdot) = g(\nabla f,\cdot) = -\Omega(J \nabla f,\cdot),
	\end{align*}
    we have that $-J\nabla f = X_f $.
    Intuitively, this means that a Hamiltonian vector field is obtained by point-wise rotation by $\pi/2$ of a gradient vector field.

    Central to the quantization procedure is the Poisson algebra of  smooth functions $(C^\infty(\amsmathbb{S}^2),\{\cdot,\cdot\})$, where $\{\cdot,\cdot\}:C^\infty(\amsmathbb{S}^2)\times C^\infty(\amsmathbb{S}^2) \to C^\infty(\amsmathbb{S}^2)$ is the \emph{Poisson bracket} defined via the symplectic form by
    \begin{align}
        \label{eq:poi}
        \{f,g\} = \Omega(X_f,X_g) = X_g(f).
    \end{align}
    This bracket is skew-symmetric, satisfies the Leibniz rule and the Jacobi identity.
    In particular, $\{\cdot,\cdot\}$ is a Lie bracket.
    Further,
    \begin{align}
        \label{eq:poi-iso}
        [X_f,X_g] = -X_{\{f,g\}},
    \end{align}
    where $[\cdot,\cdot]$ is the Lie bracket of vector fields.
    This means that $f\to X_f$ is a \emph{Poisson anti-morphism} between $(C^\infty(\amsmathbb{S}^2),\{\cdot,\cdot\})$ and $(\mathfrak{X}_{\Omega}(\amsmathbb{S}^2),[\cdot,\cdot])$, where $\mathfrak{X}_{\Omega}(\amsmathbb{S}^2)$ denotes the symplectic vector fields on $\amsmathbb{S}^2$.
    The space of all vector fields is denoted $\mathfrak{X}(\amsmathbb{S}^2)$.

    Zeitlin's model is directly applicable for $\mathfrak{X}_{\Omega}(\amsmathbb{S}^2)$.
    In order to extend it to $\mathfrak{X}(\amsmathbb{S}^2)$, we need a Poisson anti-morphism between the space of smooth vector fields and some extended Poisson algebra of functions.
    To see which space of functions this is we use standard results in Kähler geometry (cf.~\citet{le2018}): the Hodge decomposition theorem states (since $\amsmathbb{S}^2$ has trivial first coholomogy) that any $v \in \mathfrak{X}(\amsmathbb{S}^2)$ can be decomposed as
    \begin{align*}
        v = X_h + \nabla f,
    \end{align*}
    for  $f,h \in C^\infty(\amsmathbb{S}^2)$.
    As $X_h = -J\nabla h$, the two components are related by a $\pi/2$ rotation.
    Further, $Jv = \nabla h -X_f$. This means that $-f$ generates the Hamiltonian component of the rotated field while $h$ generates the gradient component.
    By letting imaginary functions generate gradient vector fields while real functions generate Hamiltonian vector fields, we can capture this behavior.
    Multiplying the vector field with $J$ corresponds to multiplying the generator with $i$.
    In other words, the \emph{real} vector field $v$ is given as the Hamiltonian vector field of a \emph{complex} function.
    To make this precise, the Hamiltonian vector field of a purely imaginary function $if \in iC^\infty(\amsmathbb{S}^2)$ is given by
    \begin{align*}
        X_{if} \coloneqq J X_f = -J^2 \nabla f = \nabla f.
    \end{align*}
    Thus, given a complex-valued function $\psi \in C^\infty(\amsmathbb{S}^2,\amsmathbb{C})$,
    \begin{align}
        \label{eq:decompos}
        X_\psi = X_{\Re \psi + i \Im \psi } =X_{\Re \psi} + \nabla \Im \psi.
    \end{align}
    Therefore,
    \begin{align*}
        JX_\psi = \nabla \Re \psi - X_{\Im \psi} = X_{-\Im \psi + i \Re \psi} =  X_{i \psi}.
    \end{align*}

    A Poisson algebra structure on $C^{\infty}(\amsmathbb{S}^2,\amsmathbb{C})$ can be obtained by complexifying $(C^\infty(\amsmathbb{S}^2),\{\})$.
    Explicitly, the bracket on $C^{\infty}(\amsmathbb{S}^2,\amsmathbb{C})$ is  given by
    \begin{align*}
         &\{{\psi_1},{\psi_2}\}_\amsmathbb{C} = \{\Re {\psi_1},\Re {\psi_2}\}-\{\Im {\psi_1},\Im {\psi_2}\}\\
            &\qquad+ i \left(\{\Re {\psi_1},\Im {\psi_2}\}+\{\Re {\psi_1},\Im {\psi_2}\}\right),
    \end{align*}
    where $\{\cdot,\cdot\}: C^\infty(M)\times C^\infty(M)\to C^\infty(\amsmathbb{S}^2)$ is the standard Poisson bracket on     $C^\infty(\amsmathbb{S}^2)$.
    By straightforward computations and the corresponding properties of $\{\cdot,\cdot\}$, one shows that $ \{\cdot,\cdot\}_\amsmathbb{C} $ is antisymmetric and satisfies the Leibniz and Jacobi identities.
    Further, it holds that $C^\infty(M,\amsmathbb{C})$ is an associative $K$-algebra where $K$ is either $\amsmathbb{C}$ or $\amsmathbb{R}$.
    Thus, $(C^\infty(M,\amsmathbb{C}),\{\}_\amsmathbb{C})$ is a Poisson algebra.

    Using the Poisson bracket we can relate the action of a Hamiltonian vector field on an imaginary function.
    As $JX_g = X_{ig}$, we have that
    \begin{align*}
        JX_g(f) = X_{ig}(f) \coloneqq i\{f,g\},
    \end{align*}
     where $f,g \in C^\infty(M)$.
    While $JX_g$ is a real vector field, applying it to a real-valued function  results in a imaginary function.
    Thus, we have for ${\psi_1},{\psi_2} \in C^\infty(\amsmathbb{S}^2,\amsmathbb{C})$,
    \begin{align*}
        &X_{\psi_2}({\psi_1}) = X_{\Re {\psi_2}} +JX_{\Im {\psi_2}} ({\psi_1}) \\
        &= X_{\Re {\psi_2}}(\Re {\psi_1})+i X_{\Re {\psi_2}}(\Im {\psi_1}) + X_{i\Im {\psi_2}}(\Re {\psi_1})+iX_{i\Im {\psi_2}}(\Im {\psi_1})  \\
        &=\{\Re {\psi_1},\Re {\psi_2}\}-\{\Im {\psi_1},\Im {\psi_2}\}\\
        &\qquad +i\left(\{\Re {\psi_1},\Im{\psi_2}\}+\{\Im {\psi_1},\Re {\psi_2} \}\right) =\{{\psi_1},{\psi_2}\}_{\amsmathbb{C}},
    \end{align*}
    and so
    \begin{equation}
        \label{eq:poidef}
        X_{{\psi_2}}({\psi_1}) = \{{\psi_1},{\psi_2}\}_\amsmathbb{C},
    \end{equation}
    a fact allowing us to  relate the Lie bracket of real vector fields with the complex Poisson bracket by the following well-known lemma.
    \begin{lemma}
        \label{th:poisiso}
        Let $\psi,{\psi_1},{\psi_2} \in C^\infty(M;\amsmathbb{C})$. Then, $\psi \mapsto X_\psi$, where $X_\psi$ is given by \Cref{eq:decompos}, is a Poisson anti-morphism, i.e,
        \begin{align*}
            -X_{\{{\psi_1},{\psi_2}\}_\amsmathbb{C}} = [X_{\psi_1},X_{\psi_2}].
        \end{align*}
    \end{lemma}
    \begin{proof}
        $X_{\psi_1}$ and $X_{\psi_2}$ are both real vector fields, so the action of $[X_{\psi_1},X_{\psi_2}]$ on a smooth complex function $h$ is also a smooth complex function. By the definition of the real Lie bracket of vector fields,
        \begin{align*}
            [X_{\psi_1},X_{\psi_2}](h) = X_{\psi_1}(X_{\psi_2}(h))-X_{\psi_2}(X_{\psi_1}(h)).
        \end{align*}
        Now, by applying \Cref{eq:poidef} twice to each term, we have that the above can be rewritten as
        \begin{align*}
            [X_{\psi_1},X_{\psi_2}](h)  =  \{{\psi_1},\{{\psi_2},h\}\}_\amsmathbb{C}+\{{\psi_2},\{h,{\psi_1}\}\}_\amsmathbb{C} = -\{h,\{{\psi_1},{\psi_2}\},h\}_\amsmathbb{C} = -X_{\{{\psi_1},{\psi_2}\}}(h).
        \end{align*}
        where the first equality is due to the Jacobi identity and the second to \Cref{eq:poidef}.
    \end{proof}

    We now derive a formulation of \Cref{eq:eulerarnold} in terms of $\{\cdot,\cdot\}_\amsmathbb{C}$.
    This corresponds to finding the vorticity formulation of the equations.
    We first let $\op{curl}m = \omega$, obtaining the equation
    \begin{align}
        \dot \omega = \op{ad}_v^T \omega.

    \end{align}

    To compute $\op{ad}_v^T \omega$, let the complex function $\omega$ act on real vector fields by the the pairing
    \begin{align*}
        \langle \omega, v\rangle = \int_{\amsmathbb{S}^2} \omega \Delta^{-1}\left( \,\mathrm{d} v^\flat + i \,\mathrm{d} (Jv)^\flat \right).
    \end{align*}
     Here $v^\flat = g(v,\cdot) = \Omega(Jv,\cdot) = \iota_{Jv} \Omega$, where $\iota$ denotes the interior product.
     By \Cref{eq:decompos} we decompose $v$ as $X_{\Re \psi_v}+\nabla \Im \psi_v$. We have that
     \begin{align*}
          \,\mathrm{d} v^\flat =  \,\mathrm{d} \iota_{-\nabla \Re \psi_v+X_{\Im \psi_v}}\Omega.
     \end{align*}
     As the symplectic form is closed, Cartan's magic formula gives that
     \begin{align*}
         \,\mathrm{d} \iota_{-\nabla \Re \psi_v+X_{\Im \psi_v}}\Omega = -\mathcal{L}_{\nabla \Re \psi_v} \Omega +\mathcal{L}_{X_{\Im \psi_v}} \Omega,
     \end{align*}
     where $\mathcal{L}$ is the Lie derivative.
     The Lie derivative of the symplectic form with respect to a symplectic vector field is zero, so the above simplifies to
     \begin{align*}
         \,\mathrm{d} v^\flat = -\mathcal{L}_{\nabla \Re \psi_v} \Omega = -\op{div} \nabla \Re \psi_v \Omega = -\Delta \Re \psi_v \Omega.
     \end{align*}
     as $\Omega$ coincides with the Riemannian volume form on $\amsmathbb{S}^2$.
     Similarly,
     \begin{align*}
          i \,\mathrm{d} (Jv)^\flat = -i \Delta \Im \psi_v \Omega
     \end{align*}
     and thus
    \begin{align*}
        \langle \omega, v\rangle = -\int_{\amsmathbb{S}^2} \omega \psi_v \Omega.
    \end{align*}
    Take now an arbitrary $u = X_{\psi_u} \in \mathfrak{X}(M)$. By  \Cref{th:poisiso},
    \begin{align*}
        &\langle \op{ad}_v^T \omega,u \rangle =  \langle  \omega,[X_{\psi_v},X_{\psi_u}] \rangle =
         -\langle  \omega,X_{\{\psi_v,\psi_u \}_{\amsmathbb{C}}} \rangle\\
        &= \int_{\amsmathbb{S}^2} \omega \{\psi_v,\psi_u\}_{\amsmathbb{C}} \Omega  = \int_{\amsmathbb{S}^2} \{\omega,\psi_v\}_{\amsmathbb{C}} \psi_u \Omega = -\langle \{\omega,\psi_v\}_{\amsmathbb{C}},u \rangle.
    \end{align*}
    The function $\psi_v = (1-\Delta)^{-1}\Delta^{-1}  \omega $ is the so-called \emph{stream function}.  Thus, dropping the subscript $v$ as the vector field was arbitrary, we have derived the vorticity formulation of
    \Cref{eq:eulerarnold},
    \begin{align}
        \label{eq:vort}
        \begin{gathered}
            \dot \omega = \{\psi,\omega\}_{\amsmathbb{C}},\\
            \psi =  (1-\Delta)^{-1}\Delta^{-1} \omega.
        \end{gathered}
    \end{align}

    We now quantize \Cref{eq:vort}.
    As a complex function consists of two real-valued functions, its real part and its imaginary part, we first turn our attention to how functions in $C^\infty(M)$ are quantized.

    The following presentation builds on \citet{hoppe1998}.
    In essence, we perform a spatial  ``truncation'' of the vorticity $\omega$ that preserves the Lie--Poisson structure of the equation.

    The vorticity equation is an equation on the dual of the vector fields, and $\mathfrak{X}(\amsmathbb{S}^2)$ is the Lie algebra of $\op{Diff}(\amsmathbb{S}^2)$.
    The quantized equation should have an equivalent structure: it should evolve on the dual of a suitable matrix Lie group.
    In order to obtain the quantization scheme, we need an
    approximation of the Laplacian.
    The quantized Laplacian allows us to build an eigenbasis of $\mathfrak{gl}(N,\amsmathbb{C})$ corresponding to the spherical harmonic basis of functions on $\amsmathbb{S}^2$.

    Explicitly, we start with the spin $\frac{N-1}{2}$ representation of  $\mathfrak{so}(3)$ in $\mathfrak{u}(N)$.
    The generators $X_1,X_2,X_3$, described by \citet{hoppe1998}, satisfy the  commutation relations of $\mathfrak{so}(3)$ up to the factor $1/N$, i.e.
    \begin{align*}
        &[X_1,X_2] = \frac{1}{N} X_3,\\
        &[X_3,X_1] = \frac{1}{N} X_2,\\
        &[X_2,X_3] = \frac{1}{N} X_1.
    \end{align*}

    As $\amsmathbb{S}^2 \subset \amsmathbb{R}^3$, there is a representation of $\mathfrak{so}(3)$ in $C^\infty(\amsmathbb{S}^2)$, with generators given by the coordinate functions $x_1,x_2,x_3$.
    They each generate rotations about their respective coordinate axis.
    This representation satisfies the  commutation relations of $\mathfrak{so}(3)$,
    \begin{align*}
        &\{x_1,x_2\} = x_3,\\
        &\{x_3,x_1\} = x_2,\\
        &\{x_2,x_3\} = x_1.
    \end{align*}

    Moreover, the Laplace-Beltrami operator on $\amsmathbb{S}^2$ is given in terms of the Poisson bracket,
    \begin{align*}
        \Delta = \sum_{k=1}^3 \{x_k,\{x_k,\cdot\}\}.
    \end{align*}
    %

    Based on the expression of $\Delta$ in terms of $\{\cdot,\cdot\}$, the quantized Laplacian on $\mathfrak{u}(N)$ should be
    \begin{align}
        \label{eq:HY}
        \Delta_N = N^2  \sum_{k=1}^3 [X_k,[X_k,\cdot]].
    \end{align}

    \citet{hoppe1998} showed that the eigenvalues of $\Delta_N$ coincide with the first $N^2$ eigenvalues of $\Delta$.
    The quantization of a function $\zeta$  is now given by a mapping from $C^\infty(\amsmathbb{S}^2)$ to $\mathfrak{u}(N)$ as follows.

    First, the expansion of  $\zeta$ in  the spherical harmonics $\{Y_{l,m}|l\geq 0, \, m = -l,\ldots,l\}$  is truncated,
     \begin{align*}
        \zeta^N  = \sum_{l=0}^N \sum_{m = -l}^l a_{l,m} Y_{l,m},
    \end{align*}
    The quantization is then performed by identifying the first $N^2$ spherical harmonics, being eigenfunctions of $\Delta$, with the eigenmatrices of $\Delta_N$, which are denoted by $T_{l,m}^N$,
    \begin{align*}
        \zeta^N \mapsto z^N = \sum_{l=0}^N \sum_{m = -l}^l a_{l,m}  T_{l,m}^N.
    \end{align*}

     It remains to describe the matrices $T_{l,m}^N$.
     The quantized Laplacian $\Delta_N$ maps $\pm m$-diagonal matrices to $\pm m$-diagonal matrices.
     A unique eigenbasis to $\Delta_N$ is given by take $T_{l,m}$ to be $\pm $-diagonal. Under this restriction, the eigenvalue problem
     \begin{align*}
         \Delta_N T_{l,m}^N = -(l+1)lT_{l,m}^N
     \end{align*}
     has a unique solution.

     This gives us a quantization of real-valued functions which we can extend to complex-valued functions by complexification.
     Just as $C^\infty(M,\amsmathbb{C}) = C^\infty(M)\otimes \amsmathbb{C}$, the quantized equivalent of complex functions is $\mathfrak{u}(N)\otimes \amsmathbb{C} \cong \mathfrak{gl}(N,\amsmathbb{C})$.
     The above quantization is thus performed twice, once for the imaginary part, and once for
     real part.

     Further, the dual of a finite-dimensional Lie algebra can be identified with the algebra itself.
     We have seen that the space of vector fields on $\amsmathbb{S}^2$ can be identified with complex valued functions, modulo a complex constant.
     From the formula for the quantized Laplacian we see that the matrices $\amsmathbb{C} I$ correspond to constants.
     Thus, the quantized equivalent of the space of vector fields is $\mathfrak{gl}(N,\amsmathbb{C})/\amsmathbb{C} I \simeq \mathfrak{sl}(N,\amsmathbb{C})$, as we identify a vector field with its generator.
     As $\mathfrak{sl}(N,\amsmathbb{C})$ is the Lie algebra of $\op{SL}(N,\amsmathbb{C})$, we conclude that quantized equivalent of $\op{Diff}(\amsmathbb{S}^2)$ is $\op{SL}(N,\amsmathbb{C})$.

     The quantized vorticity equation is given by
     \begin{align}
        \label{eq:quantvort}
         \begin{gathered}
              \dot W = [P,W]\\
              P = (1-\Delta_N)^{-1}\Delta_N^{-1} W,
         \end{gathered}
     \end{align}
     where $[\cdot,\cdot]$ denotes the matrix commutator.

    We remark that \Cref{eq:quantvort} is not fully discrete: a time-integrator is needed.
    In applications where long-term stability is required, it is important to apply a structure-preserving integrator.
    See for instance \citet{engo2001,Hairer2006,modin2019_1,modin2019} for examples of suitable integrators.

    For the convenience of the reader, we summarize the relevant groups, algebras and duals in both the quantized and classical case in Table \ref{tab:dictionary}.

    \begin{table}[ht!]
         \centering
         \begin{tabular}{l@{~~} @{}l l@{}}
             \toprule
             &\textbf{Classical} & \textbf{Quantized}          \\ \midrule
             \textbf{Group} & $\operatorname{Diff(\amsmathbb{S}^2)}$                 & $ \operatorname{SL}(N,\amsmathbb{C})$        \\
             \textbf{Algebra}&$  \mathfrak{X}(\amsmathbb{S}^2) $           & $\mathfrak{sl}(N,\amsmathbb{C})$                              \\
             \textbf{Real-valued functions} &$ C^\infty(\amsmathbb{S}^2)$       & $\mathfrak{u}(N)$ \\
             \textbf{Complex-valued functions} & $C^\infty(\amsmathbb{S}^2,\amsmathbb{C})$ & $\mathfrak{gl}(N,\amsmathbb{C})$ \\
               \textbf{Densities} &$\Omega^2(\amsmathbb{S}^2) $                & $\mathfrak{u}(N)^*$            \\ \bottomrule
             \end{tabular}
         \caption{Spaces in classical and quantized case.}
         \label{tab:dictionary}
     \end{table}

\section{Numerical example}
\label{sec:computations}
In this section we provide numerical examples. Consider the vector field on the sphere given by
\begin{align*}
    v = \alpha(n\times\nabla +\nabla) Y_{1,0}
\end{align*}
where $n$ is the normal vector of the sphere and $\alpha$ is a scalar constant ensuring simple expressions.

Further, $v$ is a linear combination of \emph{vector spherical harmonics}.
It holds that  $\operatorname{div} (n\times\nabla Y_{1,0} )= 0$ and that $\operatorname{curl}( \nabla Y_{1,0}) =0$ \citep{barrera1985}.
In the language of \Cref{sec:theory},
\begin{align}
    v = X_{-\alpha Y_{1,0}- i \alpha Y_{1,0}}.
\end{align}
Explicitly, the vector field has the Cartesian coordinate expression
\begin{align*}
    v(x,y,z) = ( y-xz,-x-yz,1-z^2).
\end{align*}
The flow of $v$ starting at $y_0$ is
\begin{align}
    \label{eq:igs}
    \Phi_{y_0}(t) = \begin{bmatrix}
        A e^t \frac{\cos(t)}{C + e^{2t}} + Be^t \frac{\sin(t)}{C + e^{2t}} \\
        -A e^t \frac{\sin(t)}{C + e^{2t}} + Be^t \frac{\cos(t)}{C + e^{2t}}\\
        \frac{e^{2t}-C}{C+e^{2t}}\\
    \end{bmatrix},
\end{align}
where $A,B$ and $C$ are constants determined by $y_0$.

\begin{figure}[ht!]
      \centering
     \begin{subfigure}[t]{0.49\textwidth}
         \centering
         \includegraphics[trim=100 50 100 50,clip,width=\textwidth]{gradpart.pdf}
         \caption{The gradient part $\alpha \nabla Y_{1,0}$. It flows from the south pole towards the north pole. It is irrational, and has nonzero divergence.}
         \label{fig:gradpart}
     \end{subfigure}
     \hfill
     \begin{subfigure}[t]{0.49\textwidth}
         \centering
         \includegraphics[trim=100 50 100 50,clip,width=\textwidth]{hampart.pdf}
         \caption{The Hamiltonian part $\alpha (n \times \nabla) Y_{1,0}$. Note that it is a pure rotation, and has zero divergence.}
         \label{fig:hampart}
     \end{subfigure}
     \hfill

        \begin{subfigure}[b]{0.52\textwidth}
         \centering
         \includegraphics[trim=100 50 100 50,clip, width=\textwidth]{vec.pdf}
         \caption{The vector field $v$. Note how it spirals upwards from the south pole towards the north pole. It is neither divergence-free nor irrotational.}
         \label{fig:vec}
     \end{subfigure}

        \caption{The vector field $v$ together with its Hamiltonian and gradient components.  }

        \label{fig:vecfields}
\end{figure}
The vector field and its two components are illustrated in \Cref{fig:vecfields}.
With this vector field we  perform two numerical examples. To achieve this, we remark that the quantized equivalent of $v$, or rather, of its generator, is the stream matrix given by
\begin{align*}
    P = -P_{\,\alpha Y_{10}} - iP_{\,\alpha Y10} \in \mathfrak{gl}(N,\amsmathbb{C}),
\end{align*}
where $P_{\,\alpha Y_{1,0}}$ is the diagonal stream matrix corresponding to the function $\alpha Y_{1,0}$.

\subsection*{Point transportation}
We consider the problem of  moving points with $v$.
The vector fields transports a  point $y_0 \in \amsmathbb{S}^2$ by the ODE
\begin{align*}
    \dot y(t) = v(y(t)),\, t \in [0,T],\, y(0) = y_0.
\end{align*}
\Cref{eq:igs} is the solution to this ODE, so we avoid numerical integration of the flow.

As this is a simple illustrative first example, consider the point $y_0 = (-1,0,0)$.
A natural first question is how to represent  $y_0$ in the quantzied world.
To this end, consider the Dirac delta function centered at $y_0$, denoted by $\delta_{y_0}$.
As we cannot hope to exactly find a quantized delta function, we instead approximate it by first creating a matrix $B_{\operatorname{north}} \in \operatorname{u}(N)$ given by
\begin{align*}
    [W]_{mn} = \begin{cases}
        &i \text{ if $m = n = N$},\\
        &0 \text{ else}.
    \end{cases}
\end{align*}
$B_{\operatorname{north}}$ corresponds to a ``blob'' at the north pole.
$B_{\operatorname{north}}$ is then rotated with an appropriate rotation matrix to center the blob at the desired point $y_0$.
See \Cref{fig:point} and \Cref{fig:blob} for an illustration of the $y_0$ and its quantized equivalent.
A blob matrix is analogous to an approximation of a delta function with bump funcions.

Another issue is how the stream matrix $P$ acts on blobs.
As a first approach, note that the blobs correspond to densities.

To understand what the action on densities must be, note that $P \in \mathfrak{gl}(N,\amsmathbb{C})$ generates a path of ``quantized diffeomorphisms'' $F$ by the equation
\begin{align}
    \begin{gathered}
        \dot F = P F,\\
         F(0) = I.
    \end{gathered}
    \label{eq:floweq}
\end{align}

The solution to \Cref{eq:floweq} is $F(t) = \exp(P t)$.
The action of a quantized diffeomorphism on a quantized density, such as a blob matrix $B$ is
\begin{align*}
    F.B = F B F^*.
\end{align*}
Indeed, this action on positive definite Hermitian matrices $B$ captures the properties of the action of diffeomorphisms on densities: (i) the trace in conserved (the total volume is conserved); (ii) if $B=I$ and $F\in \mathrm{SU}(N)$ then the action degenerates (the action of a volume preserving diffeomorphisms on the volume form does nothing); (iii) the orbit consists of all positive definite Hermitian matrices with the same trace (the orbit of the action of diffeomorphisms on densities is transitive, a result by~\citet{Mo1965}). For more details on this analogy, we refer to \citet{Mo2017}.

In \Cref{fig:blobpoint} we see the results of moving a blob starting at the point $(-1,0,0)$, compared to how the is moved by the non-quantized vector field.
Note that the quantization behaves as expected: the center of the blob is moved to the point $\Phi_{y_0}(t)$.

\begin{figure}[ht!]
      \centering
     \begin{subfigure}[t]{0.49\textwidth}
         \centering
         \includegraphics[trim=100 50 100 50,clip,width=\textwidth]{blob.pdf}
         \caption{Initial point. }
         \label{fig:point}
     \end{subfigure}
     \hfill
       \begin{subfigure}[t]{0.49\textwidth}
         \centering
         \includegraphics[trim=100 50 100 50,clip,width=\textwidth]{blob_def.pdf}
         \caption{Transported point. }
         \label{fig:point_warp}
     \end{subfigure}
     \hfill
     \begin{subfigure}[t]{0.49\textwidth}
         \centering
         \includegraphics[trim=50 50 40 50,clip,width=\textwidth]{quantblob.pdf}
         \caption{ Initial blob}
         \label{fig:blob}
     \end{subfigure}
      \begin{subfigure}[t]{0.49\textwidth}
         \centering
         \includegraphics[trim=50 50 40 50,clip,width=\textwidth]{quantblob_def.pdf}
         \caption{Transported blob}
         \label{fig:blob_warp}
     \end{subfigure}
        \caption{Upper row: the sphere with the initial point $y_0 = (-1,0,0)$ (A) and the transported point $\Phi_{y_0}(0.5)$ (B). Lower row: gridlines in spherical coordinates projected with the Hammer projection, with the initial blob at $y_0$ (C) and the transported blob (D). Note how the blob moves as the point.}
        \label{fig:blobpoint}
\end{figure}

One issue to be addressed when moving blob density matrices in this manner is that the magnitude of the nonzero values in the matrix may grow very large.
This is due to the irrotational part of the vector field contracting the blob as it moves towards the north pole.
This to be expected.

To address this issue, we now suggest another approach to transporting blobs, in which we move the center point rather  than transporting the blob as a quantized density.

Let us begin with the well-known fact that for any smooth functions $f,g$ on $\amsmathbb{S}^2$, it holds that
\begin{align*}
    \int_{\amsmathbb{S}^2} \Delta g f \,\mathrm{d} x = -\int_{\amsmathbb{S}^2} \nabla g \cdot \nabla f \,\mathrm{d} x.
\end{align*}
Analogously, for any $F,G \in \mathfrak{gl}(N,\amsmathbb{C})$ it holds that
\begin{align*}
    \operatorname{Tr}(\Delta_N F G ) = N^2\sum_{i=1}^3 \operatorname{Tr}([X_k,[X_k,F]] G ) = -N^2 \sum_{i=1}^3 \operatorname{Tr}([X_k,F][X_k,G]).
\end{align*}
Thus, the $k$th component of the quantized equivalent of the gradient of the stream function is  matrix
\begin{align}
    \label{eq:grad}
    (\nabla_N P)_k = N[X_k,P], \, k = 1,2,3.
\end{align}
where $P$ is the stream matrix and \Cref{eq:HY} was used in the first step.

Having found the quantized gradient, we move on the quantized Hamiltonian vector field.
The Hamiltonian vector field from a gradient by a rotation, and the same idea holds for the quantized equivalent, which  has exactly components of the form
\begin{align}
    \label{eq:hamilt}
    (X_P^N)_k = [N[X_k,P],B].
\end{align}
It can be verified by the commutation relations of $X_1,X_2$ and $X_3$ that this behaves analogously to the cross product of the gradient with a fixed vector.
As $B$ is centered around a point $x$, taking the Lie bracket with $B$ corresponds to the cross product of the gradient with the normal vector in $x$.
The ``quantized vector" field $V$ of $P$ thus has components given by \Cref{eq:grad} and \Cref{eq:hamilt}.

Given a function $f \in C^\infty(M)$, we can extract its value at the point $x$ by
\begin{align*}
    f(x) = \int_{\amsmathbb{S}^2} f(y) \delta_{x}(y) \,\mathrm{d} x.
\end{align*}
As blobs are analogs of delta functions, the value of the $k$th component of the quantized vector field $V$ is taken to be
\begin{align*}
    a_k(B) = \operatorname{Tr}(V_k B),
\end{align*}
and as the $k$th component of $B$ is $[X_k,B]$, the ODE governing the evolution of $B$, now understood as a smooth curve of blob matrices starting at the initial blob $B_0$, is
\begin{align}
    \label{eq:blobLP}
    \begin{gathered}
        \dot B = \sum_{k=1}^3 a_k(B)[X_k,B], \\
        B(0) = B_0
    \end{gathered}
\end{align}

\Cref{eq:blobLP} is a Lie--Poisson system and can therefore be integrated using for instance Lie--Poisson Euler, see \citet{engo2001}.

In short, given $B_n$, the following value of the blob is given by
\begin{align*}
    B_{n+1} = e^{h\sum_{k=1}^3 a_k(B) X_k } B_n  e^{-h\sum_{k=1}^3 a_k(B) X_k },
\end{align*}
where $h$ is a predetermined step size.
The result of running this recursion for $n= 200$ steps with $h = 1$ is plotted in \Cref{fig:blob_warp_center}.

\begin{figure}[ht!]
    \centering
    \includegraphics[width=0.75\textwidth]{blob_warp_nb.pdf}
    \caption{Transported blob}
    \label{fig:blob_warp_center}
\end{figure}

Note that while the blob is not centered in the same point as in \Cref{fig:blob_warp}, the blob moves in the same direction.
Further, transporting blobs in this manner does not result in the same type of blow-up that can be observed when letting the quantized diffeomorphism generated by \Cref{eq:floweq} act on the blob as if it were a density.

\subsection*{Grid deformation}
    As a second experiment, let us consider the ``icosasphere'' triangulation of the sphere into a triangular mesh, see \Cref{fig:triangulation}.

    \begin{figure}[ht]
        \centering
        \includegraphics[trim=100 50 100 50,clip]{icosphere.pdf}
        \caption{Icosahedral triangulation of the sphere. Created using \url{https://github.com/meshpro/meshzoo}. The triangulation is built up from successive refinements of coarser grids. In this image, $10$ refinements were used.}
        \label{fig:triangulation}
    \end{figure}
    The vector field generates a curve of diffeomorphisms by the differential equation
    \begin{gather*}
        \dot \gamma(t)  = v\circ \gamma(t), \\
        \gamma(0) = \operatorname{Id}.
    \end{gather*}

    The vector field $v$ acts on the vertices of the mesh as in the previous example.
    The faces of the triangulation are thereby deformed.
    They remain triangles, but their area changes.
    By coloring the faces of the deformed triangulation by the quotient of the deformed area with the undeformed area, we obtain an indication of how the diffeomorphism $\gamma(T)$ deforms the uniform density on $\amsmathbb{S}^2$.
    See \Cref{fig:def_triangulation} for the coloring when $T = 1$.
    Comparing with \Cref{fig:vec}, we see that this is expected.
    As $v$ spins points towards the north pole, the triangles already close to the north pole become smaller, and those closer to the south pole more drawn out.

    \begin{figure}[ht!]
        \centering
        \includegraphics{icosphere_def_col.pdf}
        \caption{The colored and deformed icosahedral sphere. Note that the triangles closer to the north pole are contracted, but the area of those closer to the south pole increases. }
        \label{fig:def_triangulation}
    \end{figure}

    The quantized diffeomorphism $F$ acts on the uniform density matrix  $ I$
    by
    \begin{align*}
        F.I =  F IF^* = FF^*.
    \end{align*}
    We plot the resulting density in \Cref{fig:quantdef}.
    We see that the behavior is qualitatively similar: larger values are seen close to the south pole, whereas the density is smaller around the north pole.
     \begin{figure}[ht]
        \centering
        \includegraphics[width=0.75\textwidth]{quant_col.pdf}
        \caption{Plot of the action of $F$ on $I$. Note that while the values are generally small, the larger values are seen around the south pole.}
        \label{fig:quantdef}
    \end{figure}

\section{Outlook}
\label{sec:outlook}
In this paper we have described how Zeitlin's model can be extended to problems on the entire diffeomorphism group.
We have performed two numerical experiments.
In order to provide the reader with some outlook we consider the following possible applications, to be pursued in future work.
\subsection*{Compressible flows}
    There is an extension of Arnold's description from incompressible to compressible fluids (see survey by~\cite{KhMiMo2020}).
    This extension corresponds to replacing the group of volume-preserving diffeomorphisms with all diffeomorphisms.
    Thus, geometric discretizations of compressible fluid is a direct application of the technique developed in this paper.
    For compressible fluids one also need a potential function (in addition to the metric defining the kinetic energy).
\subsection*{Optimal mass transport}
    Optimal mass transport is closely related to compressible fluids via \emph{Otto calculus} (cf.~\citet{BeBr2000}, \citet{Ot2001}). It is natural to use the quantization framework presented here to develop a quantized version of optimal mass transport, such that the transport maps are replaced by group elements in $\mathrm{SL}(N,\amsmathbb{C})$.
\subsection*{Diffeomorphic density matching}
    A possible use case from the realm of shape matching is the problem finding an an optimal diffeomorphism that transforms one initial density to a target density.
    This in turn could have potential uses in sampling \citep{bauer2017}, computer graphics \citep{dominitz2010} and medical image analysis \citep{gorbunova2012,haker2004,rottman2015}.

    An approach to the inexact density matching problem suggested in \citet{bauer2015} is as follows:  given two elements
    \begin{align*}
        I_0,I_1 \in \operatorname{Dens}(\amsmathbb{S}^2) \coloneqq \left\{I \in C^\infty(\amsmathbb{S}^2) | I\geq 0, \int_{\amsmathbb{S}^2} I \,\mathrm{d} x = 1\right\},
    \end{align*}
    the problem consists of finding a diffeomorphism $\varphi \colon \amsmathbb{S}^2 \to \amsmathbb{S}^2$ warping the initial density $I_0$ by  $\varphi \cdot I_0 = |D \varphi| I_0 \circ \varphi$.
    Here $|D\phi|$ denotes the Jacobian determinant of $\varphi$.

    The diffeomorphism $\varphi$ is given as the endpoint of a smooth curve $\gamma(t) \in \operatorname{Diff}(\amsmathbb{S}^2)$.
    It should be optimal in the sense that it minimizes the \emph{inexact matching energy functional}
    \begin{align}
        \label{eq:densmatch}
        E(\varphi) = \int_{\amsmathbb{S}^2} |\varphi \cdot I_0 -I_1|^2 \,\mathrm{d} x + \int_{\amsmathbb{S}^2} v \cdot Lv \mathrm{d} x
    \end{align}
    where $v=\dot{\gamma} \circ \gamma^{-1}$ and $L= (1-\alpha \Delta)^k$ for some $k \in \amsmathbb{N} \cup \{0\} $ and $\alpha>0$.
    Here $\Delta$ denotes the Laplace--Beltrami operator acting on vector fields.

    The first term of $E(\phi)$ is a similarity measure, penalizing dissimilarity between $\varphi \cdot I_0$ and $I_1$.
    The second term is a regularizer penalizing irregular transformations.
    It should be noted that this setting is similar to the large diffeomorphic metric mapping \citep{Bruveris2013}.

    The matching problem admits a dynamical formulation in the sense that the \emph{Euler--Arnold} framework applies and the matching problem is reduced to an equation similar to \Cref{eq:eulerarnold}.
    This equation can in turn be quantized with the method described in \Cref{sec:theory} to obtain a spatial discretization of the inexact density matching problem.

\bibliographystyle{amsplainnat.bst}