# Ensemble forecasts in reproducing kernel Hilbert space family

Benjamin Dufée, Bérenger Hug, Étienne Mémin, Gilles Tissot

► To cite this version:

Benjamin Dufée, Bérenger Hug, Étienne Mémin, Gilles Tissot. Ensemble forecasts in reproducing kernel Hilbert space family. Physica D: Nonlinear Phenomena, 2024, 459, pp.134044. 10.1016/j.physd.2023.134044. hal-04366698v2

HAL Id: hal-04366698

https://inria.hal.science/hal-04366698v2

Submitted on 31 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Ensemble forecasts in reproducing kernel Hilbert space family

Benjamin Duféeª, Bérenger Hugª, Étienne Méminª, Gilles Tissotª

ªINRIA/IRMAR Campus universitaire de Beaulieu, Rennes, 35042 Cedex, France

## Abstract

A methodological framework for ensemble-based estimation and simulation of high dimensional dynamical systems such as the oceanic or atmospheric flows is proposed. To that end, the dynamical system is embedded in a family of reproducing kernel Hilbert spaces (RKHS) with kernel functions driven by the dynamics. In the RKHS family, the Koopman and Perron-Frobenius operators are unitary and uniformly continuous. This property warrants they can be expressed in exponential series of diagonalizable bounded evolution operators defined from their infinitesimal generators. Access to Lyapunov exponents and to exact ensemble based expressions of the tangent linear dynamics are directly available as well. The RKHS family enables us the devise of strikingly simple ensemble data assimilation methods for trajectory reconstructions in terms of constant-in-time linear combinations of trajectory samples. Such an embarrassingly simple strategy is made possible through a fully justified superposition principle ensuing from several fundamental theorems.

## 1. Introduction

Forecasting the state of geophysical fluids has become of paramount importance for our daily life either at short time scale, for weather prediction, or at long time scale, for climate study. There is in particular a strong need to provide reliable likely scenarios for probabilistic forecasts and uncertainty quantification. To that end, methods based on linear combinations of an ensemble of numerical simulations are becoming more and more ubiquitous. Ensemble forecasting and data assimilation techniques, where the linear coefficients are further constrained through partial observations of the system and Gaussian stochastic filtering strategies, are good examples developed in meteorological centers for weather forecast operational routine use. However, although efficient in practice, this linear superposition principle of solutions remains theoretically questionable for nonlinear dynamics. Very recently, such techniques have been further coupled with machine learning to characterize surrogate dynamical models from long time-series of observations. This includes neural networks, analog forecasting, reservoir computing, kernel methods to name only a few of the works of this wide emerging research effort [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. To characterize a surrogate model of the dynamics, these techniques are built upon principles such as delay coordinates embedding, empirical basis functions, or sparse linear representations. The difficulty here is to get sufficiently long data series to represent the manifold associated to a given observable of the system as well as to exhibit an intrinsic representation of the whole system that does not depend on parameters such as the discrete time-stepping, the integration time window, or the initial conditions. Looking for such intrinsic representation often leads to focus (at least theoretically) on reversible ergodic systems.

Spectral representation of the Koopman operator [11] or of its infinitesimal generator constitutes one of the most versatile approaches to achieve this objective as they enable in theory to extract intrinsic eigenfunctions of the system dynamics. Koopman operator or its dual, the Perron-Frobenius operator, have attracted a great deal of attention in the operator-theoretic branch of ergodic theory [12] as well as in data science for the study of data driven modeling of dynamical systems. Several techniques, starting from [13] and [14], have been proposed to estimate such spectral representation. Most of them rely either on longtime harmonic averages [14, 15, 16, 17] or on finite-dimensional approximations of the Koopman operator such as the dynamic modes representation and its extensions [18, 19, 20, 21,

<footnote>*Corresponding author: Etienne.Memin@inria.fr</footnote>
22, 23, 24, 25, 26, 27, 28] or on Galerkin approximation and delay embedding [29, 30, 31]. See also the link between
singular spectrum analysis (SSA) method [32], data-adaptive harmonic decomposition [33], Hankel alternative [29]
and Koopman analysis [34, 35] for ergodic systems.

These methods enable performing spectral projection from long time series of measured observables of the system. The practical linear nature of the compositional Koopman operator, defined for any observable, $f$, of a dynamical system $\Phi_t(x)$ with initial condition $x$, as $(U_t f)(x) := f \circ \Phi_t(x)$, is nevertheless hindered by its infinite dimension. Due to this, it often exhibits a spectrum with a continuous component [36] and is not amenable to diagonalization nor to any sparsity assumption. Finite dimensional approximations and the direct use of matrix algebra numerical methods raises immediate questions on the pertinence of the associated estimation. To tackle these difficulties several techniques combining regularization, compactification of the generators, and reproducing kernel Hilbert space (RKHS) to work in spaces of smooth functions have been recently proposed [16, 17, 37]. They are possibly associated with forecast strategies [38, 39, 40] and transfer operator estimations [41, 42]. For data driven specification of dynamical systems RKHS or kernel methods have been also explored in a more direct way in terms of advanced Gaussian process regression techniques [43, 5, 44]. Kernel methods have been found particularly efficient in this context when coupled with a learning strategy for the kernel [6].

Working with RKHS brings convenient regular basis functions together with approximation convergence guar-
anties. However, the dynamical systems need still to be ergodic, and invertible to ensure that the Koopman operator
has a simple point spectrum and an associated system of orthogonal basis eigenfunctions [12]. Measure invariance and
invertibility warrants that the associated Koopman/Perron-Frobenius operators are dual unitary operators [12]. Note
nevertheless that state-of-the-art estimation techniques based on Dynamic Mode Decomposition (DMD) or Hankel
matrix [45, 29] require very long time series to converge to actual Koopman eigenvalues [34]. This is all the more
detrimental for high-dimensional state-spaces. An extended dynamical modes decomposition has also been proposed
to approximate the Koopman infinitesimal generator of deterministic or stochastic systems [46]. Like DMD, this
approach corresponds to a direct discretization of the Koopman operator.

Recent reviews on modern Koopman analysis applications and their related algorithmic developments from mea-
surements time series can be found in [47, 48, 49]. Another limitation is related to the definition of a kernel (and hence
of a RKHS) independent of the dynamics to describe the system’s observable, which de facto boils down to the strong
assumption of an invariance of the RKHS under the Koopman groups. This condition requires to define the RKHS
from the subspace spanned by the Koopman operator eigenfunctions, which excludes in general classical parametric
kernels. Note that in particular bias in the estimation of the Koopman eigenvalues have been observed for kernel-
based methods based on finite dimensional approximation and sparsity assumptions [50, 51]. A technique coupling
RKHS and the extended DMD approach of [46], enabling to provide an approximation of the Koopman generator,
was proposed in [52]. In the same way as the previous techniques, an invariance of the RKHS is also assumed.

In this study we will explore the problem from a different angle. Instead of working with long time series of data,
we will work with a set of realizations of a given dynamical system. The objective will be to learn this dynamics
locally in the phase space to synthesize in a fast way new trajectories of the system through a theoretically well
justified superposition principle. This superposition principle will in particular enable us to fully justify the linear
combination assumption employed in many ensemble data assimilation filters even though the dynamics is nonlinear.
To perform an efficient estimation inherently bound to the dynamics, we propose to embed the ensemble members in a
set of time-evolving RKHS defining a family of spaces. This setting is designed to deal with large-scale systems such
as oceanic or meteorological flows, for which it is out of the question to explore the whole attractor (if any), neither to
run very long time simulations calibrated on real-world data of a likely transient system due to climate change. The
present work relies on the fact that the kernel functions, usually called feature maps, between the native space and the
RKHS family are transported by the dynamical system. This creates, at any time, an isometry between the evolving
RKHS and the RKHS at a given initial time. Instead of being attached to a system state, the kernel is now fully
associated with the dynamics and new ensemble members embedded in the RKHS family at a given time can be very
simply forecast at a further time. For invertible dynamical systems, the Koopman and Perron-Frobenius operators in
the RKHS family are unitary, like in $L^2$. However, they are in addition uniformly continuous, with bounded generators,
and diagonalizable. As such they can be rigorously expanded in exponential form and their eigenfunctions provide
an intrinsic orthonormal basis system related to the dynamics. This set of analytical properties immediately brings
practical techniques to estimate Koopman eigenfunctions or essential features of the dynamics. Evaluations of these
Koopman eigenfunctions at the ensemble members are directly available. Finite-time Lyapunov exponents associated
with each Koopman eigenfunction are easily accessible on the RKHS family as well and can be used to get practical predictability time of the system. Expressions relating Koopman eigenfunctions to the tangent linear dynamics and its adjoint are also immediately available for data assimilation. These expressions, ensuing from the fundamental results listed previously, correspond in their simpler forms to the ensemble-based finite difference approximations intensively used in ensemble methods. Beyond that, the RKHS family finally leads to a theoretically well grounded superposition principle (for nonlinear dynamics) enabling the devise of embarrassingly simple methods for trajectory reconstructions in terms of constant-in-time linear combinations of trajectory samples.

The embedding of the trajectories in a time-evolving RKHS which follows the dynamical system can be related to concepts of local kernels [53], embedding in dynamical coordinates [54] or embedding of changing data over the parameter space as described in Coifman and Hirn [55]. In our case, we further exploit the relations with the Koopman operator together with an associated isometry, which allows us to define simple practical schemes for the setting of the time-evolving kernels and the modal information that can be extracted from them.

The main objective of this study is to propose a numerical framework for ensemble-based estimation and simulation of geophysical dynamics driven by systems of partial differential equations (PDE). In this study, we present theoretical elements for analyzing this framework. The mathematical assertions herein are based on idealized assumptions such as an infinite number of ensemble members, invertible dynamical systems and the availability of an ensemble of initial conditions. These theoretical statements will be further refined with more realistic assumptions in future works.

This paper is organized as follows. We start briefly recalling the definition and some properties of RKHS in the next section. These definitions are complemented by Appendix A which details and adapts several known results to our functional setting. The embedding of the Koopman and Perron-Frobenius operators in the RKHS family are then presented in sections 3 and 4. This is followed by the enunciation of our main results. The main theorems and properties are demonstrated in section 4 and in the appendixes. Numerical results for a quasi-geostrophic barotropic ocean model in an idealistic north Atlantic configuration will be presented in section 5 to support the pertinence of the proposed approach. This work corresponds to a presentation that has been given to the online seminar on Machine Learning for Dynamical Systems organized by the Turing institute (link).

## 2. Reproducing kernel Hilbert spaces

A RKHS, $\mathcal{H}$, is a Hilbert space of complex functions $f: E \mapsto \mathbb{C}$ defined over a non empty abstract set $E$ on which a positive definite kernel and a kernel-based inner product, $\langle \cdot, \cdot \rangle_{\mathcal{H}}$ can be defined. Throughout this work, we will consider that $E$ is a compact subset of a function vector space as we will work with an integral compact operator from which a convenient functional description can be set. In addition to enable the RKHS functions to be Gateaux differentiable, we will assume as in [56] that $E$ is defined as the closure of its nonempty interior, and we will work with RKHS defined from smooth enough kernels. RKHS possess remarkable properties, which make their use very appealing in statistical machine learning applications and interpolation problems [57, 58]. The kernels from which they are defined have a so called “reproducing property”.

**Definition 1 (Reproducing kernel)** Let $\mathcal{H}$ be a Hilbert space of $\mathbb{C}$-valued functions defined on a non-empty compact topological space $E$. A map $k: E \times E \mapsto \mathbb{C}$ is called a reproducing kernel of $\mathcal{H}$ if it satisfies the following principal properties: $\forall x \in E$,

* membership of the evaluation function $k(\cdot, x) \in \mathcal{H}$,

* reproducing property $\forall f \in \mathcal{H}, \langle f, k(\cdot, x) \rangle_{\mathcal{H}} = f(x)$.

The last property, provides an expression of the kernel as $k(x, y) = \overline{k(y, x)} = \langle k(\cdot, y), k(\cdot, x) \rangle_{\mathcal{H}}$, which is Hermitian – with $\bar{\cdot}$ denoting complex conjugate –, positive definite and associated with a continuous evaluation function $\delta_x f = \langle f, k(\cdot, x) \rangle_{\mathcal{H}} = f(x)$. The continuity of the Dirac evaluation operator is indeed sometimes taken as a definition of RKHS. By the Moore-Aronszajn theorem [59], the kernel $k$ defines uniquely the RKHS, $\mathcal{H}$, and vice versa. The set spanned by the feature maps $\text{Span}\{k(\cdot, x), x \in E\}$, is dense in $(\mathcal{H}, \|\cdot\|_{\mathcal{H}})$. We note also that useful kernel closure properties enable to define kernels through operations such as addition, Schur product, and function composition [57]. Besides, RKHS can be meaningfully characterized through integral operators, leading to an isometry with $L_{\mathbb{C}}^2(E, \nu)$ the space of square integrable functions defined on a compact metric space $E$ with finite measure $\nu$ [58].
## Integral kernel operators

Let $k: E \times E \to \mathbb{C}$ be $C^{(1,1)}(E \times E)$ (i.e. one time differentiable with respect to each argument), Hermitian, and positive definite, and let the map $\mathcal{L}_k: L_{\mathbb{C}}^2(E, \nu) \to L_{\mathbb{C}}^2(E, \nu)$ be defined as:

$$ (\mathcal{L}_k f)(x) = \int_E k(x, y)f(y) \nu(dy). \qquad (1) $$

This operator, which must be understood within the composition with the continuous inclusion $i: C^0(E, \mathbb{C}) \hookrightarrow L_{\mathbb{C}}^2(E, \nu)$, can be shown to be well defined, positive, compact and self-adjoint [58]. We will always assume that the space $L_{\mathbb{C}}^2(E, \nu)$ is infinite dimensional. The range of this operator is assumed to be dense in $L_{\mathbb{C}}^2(E, \nu)$ and hence infinite dimensional. From Mercer's theorem [60], the feature maps $k(\cdot, x)$ span a RKHS defined through the eigenpairs $(\mu_i, \varphi_i)_{i \in \mathbb{N}}$ of the kernel $k$:

$$ \mathcal{H} := \{f \in L_{\mathbb{C}}^2(E, \nu), f = \sum_{i=0}^{\infty} a_i \varphi_i : \sum_i \frac{|a_i|^2}{\mu_i} < \infty\}, \qquad (2) $$

with no null eigenvalues and an infinite sequence of eigenvalues since we have assumed that the kernel range is dense in $L_{\mathbb{C}}^2(E, \nu)$. The rank of the kernel (number of – non-zero – eigenvalues) corresponds to the dimension of $\mathcal{H}$, and will always be infinite in this work, which excludes the cases where $E$ is a finite set. The RKHS $\mathcal{H}$ is a space of smooth functions with decreasing high frequency coefficients. In fact, for all $u \in E$, there exists a constant $C > 0$ for which $\|\partial_u f\|_{L_{\mathbb{C}}^2(E, \nu)} \le C \|f\|_{\mathcal{H}}$ (Theorem 7 – Appendix A.2), where the derivative denotes the Gateau directional derivative of function $f$ in the direction $u$ defined as

$$ \forall x, u \in E, \partial_u f(x) = \lim_{\epsilon \to 0} \frac{1}{\epsilon} (f(x + \epsilon u) - f(x)). $$

In order to properly define the Gateaux derivative, a sufficient condition is to embed $E$ with a local vector space structure, which is for instance the case of differentiable manifolds. To enable such a differentiability property of the RKHS functions, we will in particular assume a weaker condition, namely that $E$ is the closure of its non empty interior [56]. We may also define uniquely a square-root symmetric isometric bijective operator $\mathcal{L}_k^{1/2}$ between $L_{\mathbb{C}}^2(E, \nu)$ and $\mathcal{H}$ (Appendix A.3). This operator enables to go from $L_{\mathbb{C}}^2(E, \nu)$ to $\mathcal{H}$ by increasing the functions regularity while its inverse lowers the function regularity by bringing them back to $L_{\mathbb{C}}^2(E, \nu)$. Both operators are bounded. The injection $j: \mathcal{H} \to L_{\mathbb{C}}^2(E, \nu)$ is continuous and compact (Theorem 8 – Appendix A.3) and $j(\mathcal{H})$ is assumed to be dense in $L_{\mathbb{C}}^2(E, \nu)$ (Prop. 11). These statements are precisely recalled in Appendix A.

**Remark 1** As mentioned above, the spaces $\mathcal{H}$ and $L_{\mathbb{C}}^2(E, \nu)$ are infinite dimensional spaces, preventing $E$ to be a finite set. Note that the denseness and the infinite dimensionality assumptions can easily be relaxed – see Remark 16 – Appendix A.3.

## 3. Dynamical systems on a RKHS family

We consider an invertible nonsingular dynamical system $X(t) = \Phi_t(X_0)$, defined from a continuous flow, $\Phi_t$ (meaning that, for any $X \in \Omega$, the mapping $t \mapsto \Phi_t(X)$ is continuous), on a compact invariant metric phase space differentiable manifold, $\Omega$, (i.e. $\Phi_t^{-1}(\Omega) = \Omega, \forall t \in \mathbb{R}^+$) of time evolving vector functions over a spatial support $\Omega_x$. The functions $X: \mathbb{R}^+ \times \Omega_x \to \mathbb{R}^d$ with $X \in C^p, p \ge 1$, are solutions of the following $d$-dimensional differential system:

$$ \left\{ \begin{array}{l} \partial_t X(t, \cdot) = \mathcal{M}(X(t, \cdot)), \text{ with } X(t, \cdot) \in \Omega, \forall t > 0, \\ X(0, \cdot) = X_0(\cdot). \end{array} \right. \qquad (3) $$

The nonlinear differential operator $\mathcal{M}: \Omega \to \Omega$ is assumed $C^1(\Omega)$, and in particular its linear tangent expression defined as the Gateaux derivative:

$$ \partial_X \mathcal{M}(X) \delta X = \lim_{\beta \to 0} \frac{1}{\beta} \left( \mathcal{M}(X(t, \cdot) + \beta \delta X(t, \cdot)) - \mathcal{M}(X(t, \cdot)) \right) $$
is such that $\sup_{X \in \Omega} \partial_X \mathcal{M}(X) < \infty$ (since $\Omega$ is compact). Let us note that such a differentiability assumption is quite common in geophysical flows, as otherwise so-called 4DVar variational assimilation strategies [61] derived from optimal control theory [62] and that are routinely used in weather and climate centers would make no sense. In the following, we consider the measure space $(\Omega_x, \mathcal{d})$ where $\mathcal{d}$ is the Lebesgue measure on $\Omega_x$. We note $L^2(\Omega_x, \mathbb{R}^d) := \{f = (f_1, \dots, f_d) : \Omega_x \to \mathbb{R}^d : \text{ for all } 1 \le i \le d \quad f_i \in L^2(\Omega_x, \mathbb{R})\}$ the vector space of square integrable functions on $\Omega_x$. We note $\|\cdot\|_{L^2}$ the norm associated with $L^2(\Omega_x, \mathbb{R}^d)$ and given by $\|f\|_{L^2}^2 := \sum_{i=1}^d \|f_i\|_{L^2(\Omega_x, \mathbb{R})}^2$ for all $f \in L^2(\Omega_x, \mathbb{R}^d)$. The set $\Omega$ is included in $L^2(\Omega_x, \mathbb{R}^d)$. The system (3) is assumed to admit a finite invariant measure $\nu$ on $\Omega$ (note that from invertibility the measure is also nonsingular with $\nu(\Phi_t^{-1})(A) = 0, \forall A \subset \Omega$ such that $\nu(A) = 0$). The system's observables are square integrable measurable complex functions with respect to measure $\nu$. They belong to the Hilbert space $L^2_\mathcal{C}(\Omega, \nu)$ with the inner product $\langle \cdot, \cdot \rangle_{L^2_\mathcal{C}(\Omega, \nu)}$ given for $f$ and $g \in L^2_\mathcal{C}(\Omega, \nu)$ by

$$ \langle f, g \rangle_{L^2_\mathcal{C}(\Omega, \nu)} := \int_{\Omega} f(y) \overline{g(y)} \nu(dy). $$

Depending on the context $X(t, \cdot)$ or $X_t$ will denote either an element of $\Omega$ or the function $X(t, \cdot) : \Omega_x \to \mathbb{R}^d$ such that $x \mapsto X(t, x)$.

In this work, the set of different initial conditions $\Omega_0$ is an infinite compact subset of $L^2(\Omega_x, \mathbb{R}^d)$. For all time $t \ge 0$, we denote by $\Omega_t \subset \Omega \subset L^2(\Omega_x, \mathbb{R}^d)$ the space defined as $\Omega_t := \Phi_t(\Omega_0)$. Furthermore, the set $\Omega_0$ of initial conditions will be assumed to be given – as it is usual in ensemble data assimilation or forecasting applications – and rich enough so that $\bigcup_{t \ge 0}(\Omega_t) = \Omega$. Hence, by this, each point of the manifold, $\Omega$, is assumed to be uniquely characterized by a given initial condition and the integration of the invertible dynamical system over a given time $t$. In other words, for any $X \in \Omega$, there exists a unique initial condition $X_0 \in \Omega_0$ and a unique time $t \in \mathbb{R}_+$ such that $X = \Phi_t(X_0)$. All the points of the phase space are tagged with a reaching time. A stationary point belongs to all sets $\Omega_t$ and a recurrent point belongs also to several sets. All the sets $\Omega_t$ will be assumed to be compact subset of a functions vector space, defined as the closure of their non-empty interior. Compactness enables the use of a classical Mercer theorem while non-empty interior allows differentiability of the associated RKHS feature maps. Compactness could be relaxed with generalized form of the Mercer theorem [63, 64]. This is however outside of the scope of this paper and we will remain in the case of compactness assumption. Let us note that an assumption of finite dimensional vector spaces for the set $\Omega_t$ could have been done (which comes immediately by the Riesz theorem, if we additionally impose a vectorial space structure together with compactness and non-empty interior). While many geophysical systems are conjectured to have finite-dimensional attractors (like the 2D Navier-Stokes equations [65, 66]), we prefer to adhere to a more general assumption here.

**Remark 2** For measure-preserving systems, the assumption $\Omega = \bigcup_{t \ge 0} \Omega_t$ is weaker than an ergodicity assumption, which is equivalent (together with Poincaré recurrence associated with measure preserving dynamical system) to the following statement [12] (Lemma 6.19): for any (non-empty) measurable set $A \subset \Omega$, then the set $\bigcup_{t \ge 0} \Phi_t^{-1}(A) = \Omega$. Ergodicity is very commonly assumed in the literature on Koopman spectral analysis [47]. We have here a weaker assumption.

Defining at each time $t$, from ensemble $\Omega_t$, a positive Hermitian kernel $k_{\Omega_t} : \Omega_t \times \Omega_t$, there exists a unique associated RKHS $\mathcal{H}_t$. In the following for the sake of concision the kernels $k_{\Omega_t}$, will be denoted by $k_t$ to refer to the dependence on the set $\Omega_t$. For all $t \ge 0$, we will use the notation $X_t = \Phi_t(X_0)$ and $(\mathcal{H}_t, \langle \cdot, \cdot \rangle_{\mathcal{H}_t}, \|\cdot\|_{\mathcal{H}_t})$ for the RKHS associated with the kernel $k_t$ defined on $\Omega_t \times \Omega_t$. The kernels are assumed to be $C^{(1,1)}(\Omega_t \times \Omega_t)$ and as a consequence as shown in appendix A, the associated feature maps have derivatives in $\mathcal{H}_t$. The RKHS $\mathcal{H}_t$ for all time $t \ge 0$ forms a family of Hilbert spaces of complex functions, each of them equipped with their own inner product $k_t(Y_t, X_t)$, for all functions $X_t, Y_t \in \Omega_t$. At any time, a measurable function of the system state, usually referred to as an observable, $f$, belonging to the RKHS $\mathcal{H}_t$ can be described as the limit of a linear combination of the feature maps $\{k_t(\cdot, X_t), X_t \in \Omega_t\}$. As it will be shown, the features maps of this RKHS family can be expressed on a time evolving orthonormal systems of basis functions, connected to each other through an exponential form and given by the eigenfunctions of the infinitesimal evolution operator of a “Koopman-like” operator defined on the RKHS family. The RKHS family is defined by $\mathcal{W} = (\mathcal{H}_t)_{t \ge 0}$. In the following, we shall present a summary of the mathematical results associated to the Koopman operator in the RKHS family.
Numerically, in practice, the setup is based on n realizations (called ensemble members) of this dynamical system, {$X_t^{(i)}, i = 1, \dots, n$}, generated from a finite set of different initial conditions {$X_0^{(i)}, i = 1, \dots, n$} $\subset \Omega_0$, are available up to time T. Still, we underline that, in the following development, the time horizon can be infinite, the sets $\Omega_t$ are infinite and the corresponding RKHS $\mathcal{H}_t$ are infinite dimensional. This setting (both practically and theoretically) is quite common for ensemble methods applied to geophysical systems. The ensemble size is small in general, while the phase space is in theory infinite dimensional (or at least of very high dimension). In the same way as it is usual in data assimilation or climate/weather forecasting, the dynamics on which we will experiment our setting will be assumed to be invertible at least on small time range. Although often derived from Euler equation, geophysical dynamics introduce in one way or the other some viscosity and are not stricto-sensu invertible. Additionally, for most of these systems, only local-in-time solutions have been demonstrated thus far. As demonstrated in the numerical section, the proposed framework yields good results for a simple geophysical PDEs system within a short-time range. As already mentioned in the introduction, the primary goal of this work is to establish a numerical framework for ensemble-based estimation and simulation. We provide below some theoretical elements for analyzing this framework. These mathematical statements rely in particular on the assumptions of an infinite number of members in the ensemble and on invertible dynamical systems, which is obviously not the case in practice. These theoretical statements will be refined with more realistic assumptions in subsequent works.

## 4. Koopman operator in the RKHS family

So far, we did not give any precise definition of the kernels associated to the RKHS family $\mathcal{H}_t$ yet. These kernels are defined from a known a priori initial kernel, $k_0 : \Omega_0 \times \Omega_0$, as:

**Definition 1** ($\mathcal{H}_t$ kernel) The kernel $k_t$ associated to the RKHS $\mathcal{H}_t$ are defined as

$$ \forall X_t, Y_t \in \Omega_t, \quad k_t(Y_t, X_t) = k_0(\Phi_t^{-1}(Y_t), \Phi_t^{-1}(X_t)), \quad (4) $$

where $k_0 : \Omega_0 \times \Omega_0$ is a given kernel.

These kernels can also be equivalently defined introducing the operators, $\mathcal{U}_t$, operating on the feature maps. An isometric property of this operator on the RKHS family, enables us to fully define the kernels along time in the same way as the previous definition. The operator $\mathcal{U}_t : \mathcal{H}_0 \to \mathcal{H}_t$ defined such that

$$ \mathcal{U}_t k_0(\cdot, X_0) = k_t(\cdot, \Phi_t(X_0)), \quad (5) $$

transports the kernel feature maps on the RKHS family by composition with the system's dynamics. This operator, and more specifically an associated infinitesimal evolution operator, will enable us to define the feature maps of $\mathcal{H}_t$ from an initial set of feature maps on $\mathcal{H}_0$. As will be fully detailed in section 4.1, we will see that the operator $\mathcal{U}_t$ is indeed directly related to the restriction on $\mathcal{H}_t$ of the adjoint of the Koopman operator $U_t$ on a bigger RKHS space, $\mathcal{H}$, associated to a fixed kernel defined on the whole phase space, $\Omega$. As $\mathcal{U}_t$ propagates forward the second argument of the feature maps, it is referred in the following to as the Koopman kernel operator in the RKHS family. Indeed, it will be pointed out that, for any $f \in \mathcal{H}$, and any $X_0 \in \Omega_0$,

$$ \begin{aligned} U_t f(X_0) &= f(X_t) = \langle R_t f, k_t(\cdot, X_t) \rangle_{\mathcal{H}_t} \\ &= \langle R_t f, \mathcal{U}_t k_0(\cdot, X_0) \rangle_{\mathcal{H}_t}, \end{aligned} $$

where $U_t$ denotes the Koopman operator operating on $L_{\mathcal{H}}^2(\Omega, \nu)$, and $R_t f$ denotes the restriction of $f$ to $\Omega_t \subset \Omega$. This expression exhibits a kernel expression of the Koopman operator definition and formally, at this point, the operator $\mathcal{U}_t$ can be thought as a kernel expression of the Koopman operator.

The global kernel $k$ (respectively the associated RKHS) is tightly bound to the time evolving kernels $k_t$ (respectively $(\mathcal{H}_t)_{t \ge 0}$). The restriction on $\mathcal{H}$ of the Koopman operator $U_t$ and its adjoint the Perron-Frobenius $P_t$ exhibit some remarkable properties. As classically, the operators $U_t$ and $P_t$ are unitary in $L_{\mathcal{H}}^2(\Omega, \nu)$ (Prop. 1), but they have the striking property to be uniformly continuous in $L_{\mathcal{H}}^2(\Omega, \nu)$ (i.e. with bounded generators – Theorem 2). As such, they can be expanded in an uniformly converging exponential series. Nevertheless, it must be outlined that the fixed kernel
$k(x, y)$ and consequently $\mathcal{H}$ are in practice only partially accessible as they are defined on the whole manifold of the dynamics and such expansion cannot be directly used. A local representation of the RKHS family $\mathcal{W}$ is on the other hand much easier to infer in practice through the time evolution of ensembles $\Omega_t$ and kernels $k_t$. As we shall see, the operator $\mathcal{U}_t$ on $\mathcal{W}$ inherits many properties from $U_t$ and in particular, a related form of exponential series expansion (Theorem 1).

The Koopman kernel operator in the RKHS family, $\mathcal{U}_t$, defines an isometry from $\mathcal{H}_0$ to $\mathcal{H}_t$ (Theorem 3)

$$ \left\langle k_t(\cdot, \Phi_t(Y_0)), k_t(\cdot, \Phi_t(X_0)) \right\rangle_{\mathcal{H}_t} = \left\langle k_0(\cdot, Y_0), k_0(\cdot, X_0) \right\rangle_{\mathcal{H}_0}, \qquad (6) \\ k_t(\Phi_t(X_0), \Phi_t(Y_0)) = k_0(X_0, Y_0). $$

This isometry ensues obviously directly from definition 1. But, it can be also guessed from definition (5) and the unitarity of $\mathcal{U}_t$ (Theorem 3, Prop.6), inherited from the unitarity in $L_{\mathbb{C}}^2(\Omega, \mathbf{v})$ of the Koopman operator and of its adjoint, the Perron-Frobenious operator. This property is of major practical interest as it allows us to define the kernels of the RKHS family from a given initial kernel fixed by the user. The kernels remain constant along the system trajectories. Alternatively, an explicit form of the feature maps can be obtained from an adjoint transport equation associated to an evolution operator in the RKHS family. Nevertheless, the isometry is far more straightforward to use to set the kernel evolution. This isometry will reveal also very useful for the data assimilation and trajectories reconstruction problems investigated in section 5. Strikingly, we have even more than this kernel isometry. An evolution operator $A_{U,t}: \mathcal{H}_t \to L_{\mathbb{C}}^2(\Omega_t, \mathbf{v})$, associated to the infinitesimal generator of Koopman operator $U_t$ can also be defined as

$$ A_{U,t} k_t(\cdot, X_t) := \partial_{\cdot, \mathcal{H}(\cdot)} k_t(\cdot, X_t), \qquad (7) $$

where $\partial_{\cdot, k_t(\cdot, X_t)}$ stands for the Gateaux directional derivative along $u \in \Omega$ of function $k(\cdot, X_t)$. This operator that will be shown to be bounded (Prop.7) and skew-symmetric (Prop.8) for the inner product of $L_{\mathbb{C}}^2(\Omega, \mathbf{v})$ plays the role of an infinitesimal generator on $\mathcal{W}$ and enables us expressing an exponential expansion of $\mathcal{U}_t$.

**Theorem 1 (The RKHS family spectral representation)** Let $\Omega$ be a compact metric differentiable manifold that is invariant by an invertible continuous flow, $\Phi_t$, defined from a dynamical system (3) admitting an invariant finite measure $\mathbf{v}$ on $\Omega$. Let $\Omega_0 \subset \Omega$ be an infinite compact subset of initial conditions and define $\Omega_t = \Phi_t(\Omega_0)$, with the assumption that $\Omega = \cup_{t \ge 0} \Omega_t$. Consider a $C^{(1,1)}$ initial kernel $k_0: \Omega_0 \times \Omega_0 \to \mathbb{C}$ that uniquely defines an infinite dimensional initial RKHS $\mathcal{H}_0$, and, the time dependent kernels $k_t$ defined in (4) together with its associated RKHS family $\mathcal{W} = (\mathcal{H}_t)_{t \ge 0}$. Then, the evolution operator $A_{U,t}: \mathcal{H}_t \to L_{\mathbb{C}}^2(\Omega_t, \mathbf{v})$ given in (7), and which is defined from the restriction of the infinitesimal generator $A_U$ of the Koopman operator in $L_{\mathbb{C}}^2(\Omega, \mathbf{v})$, can be diagonalized, at any time $t \ge 0$, by an orthonormal basis $(\psi_t^t)_{\ell}$ of $\mathcal{H}_t$ such that

$$ A_{U,t} \psi_t^t(X_t) = \lambda_t j \circ \psi_t^t(X_t), $$

where $j$ is the injection $j: \mathcal{H}_t \to L_{\mathbb{C}}^2(\Omega_t, \mathbf{v})$. We have in addition the following relation between the orthonormal basis systems along time:

$$ \forall t \ge 0 \quad \psi_t^t(X_t) = \exp(t \lambda_t) \psi_t^0(X_0), \qquad (8) $$

with $X_t = \Phi_t(X_0)$. Furthermore, the purely imaginary eigenvalues $(\lambda_t)_{\ell}$ do not depend on time.

This theorem, which constitutes our main result, provides us a time-evolving system of orthonormal bases of the RKHS family. It brings us the capability to express any observable of the system in terms of bases that are intrinsically linked to the dynamics and related to each other by an exponential relation. The eigenvalues of $A_{U,t}$ are purely imaginary since this operator is skew-symmetric in $L_{\mathbb{C}}^2(\Omega_t, \mathbf{v})$. Remarkably, the eigenvalues of each $A_{U,t}$ do not depend on time and are connected with the same covariant eigenfunctions (in the sense of (8)). These eigenfunctions correspond to restrictions of eigenfunctions of the infinitesimal generator of Koopman/Perron-Frobenius operators defined on $\mathcal{H}$. The proofs of these results are fully detailed in the next section.

The frozen in time spectrum, also referred to as iso-spectral property in the literature, connects directly the RKHS family representation with Lax pair theory of integrable system [67]. Here the Lax pair is provided by the operator $A_{U,t}$ at a given time and the evolution equation (8). The expression of the dynamical system in the RKHS
family (i.e. through the feature maps) provides a compatibility condition for the Lax pair. As the existence of a Lax pair is directly related to integrable systems, it means that the the RKHS family representation formally retains only the integrable part of a (non-necessary integrable) dynamical system. This relation between Koopman operator and integrable systems has been already put forward in several works [47, 68]. Here an ensemble characterization of such a relation is exhibited.

The kernel isometry (6) (Theorem 3) and the Koopman spectral representation in the RKHS family (Theorem 1) constitute fundamental results enabling us to build very simple ensemble-based trajectory reconstructions for new initial conditions without the requirement of resimulating the dynamical system. Amazingly, the family of RKHS together with the Koopman isometry allows to define a system's trajectory as a constant-in-time linear combination of the time varying RKHS feature maps. Several of such reconstruction techniques, based on this fully justified superposition principle, will be shown in section 5. We will in particular demonstrate the potential of this framework with a simple and very cheap data assimilation technique in the context of very scarce data in space and time. In the next section we demonstrate the different properties related to the RKHS family.

## 4.1. Proofs on the properties of the Koopman operator in the RKHS family and of the the RKHS family spectral theorem

As explained in the previous section the RKHS family, $(\mathcal{H}_t)_{t \ge 0}$, does not have a good topology to work with. We need first to start defining a “big” set with a better topology and which encompasses all the RKHS $\mathcal{H}_t$. On this big encompassing set, we shall then define a Koopman operator enabling us to study properly the Koopman operator in the RKHS family.

### 4.1.1. Construction of the “big” encompassing set $\mathcal{H}$

The phase space $\Omega$ corresponds to the set generated by the value of the dynamical system at a given time $t$. We note hence that $\Omega_t$ is a subset of $\Omega$. Each point of $\Omega$ designates a phase-space point $X = \Phi_t(X_0)$ uniquely defined from time $t$ and initial condition $X_0 \in \Omega_0$. In order to define the RKHS $\mathcal{H}$, let us define from the kernel $k_t: \Omega_t \times \Omega_t$ a symmetric positive definite map $k: \Omega \times \Omega$.

**Definition 2 (H kernel)** For all $X = \Phi_r(X_0), Y = \Phi_s(Y_0) \in \Omega$, with $X_0, Y_0 \in \Omega_0$ we define

$$k(X, Y) = k_0(X_0, Y_0) \ell(r, s) = k_t(\Phi_t(X_0), \Phi_t(Y_0)) \ell(r, s) \quad \forall t \ge 0, \quad (9)$$

where $\ell: \mathbb{R}_+ \times \mathbb{R}_+ \to \mathbb{R}$ is a symmetric kernel defined, for all $r, s \in \mathbb{R}^+$ by

$$\ell(r, s) = \varphi(r - s), \quad (10)$$

where $\varphi$ is a twice-differentiable even positive definite function such that $\varphi(0) = 1$.

The positivity and symmetry of kernel $k$ ensues from the property of kernels $k_0$ and $\ell$, which are assumed to be valid kernels. Kernel $k$ inherits the regularity conditions of $k_0$ and $\ell$ and is $C^{(1,1)}(\Omega \times \Omega)$ as well. In the trivial case where $\varphi = 1$, then comparing any pair of points on two trajectories would be equivalent to compare the initial conditions, which would result in a quite poor kernel, and degeneracy issues. In order to enrich the kernel structure, one can think of $\varphi$ as a regularized Dirac distribution, or a time Gaussian distribution, that will discriminate the points of the phase space that are reached at different times.

By the Moore-Aronszajn theorem, there exists a unique RKHS called $(\mathcal{H}, \langle \cdot, \cdot \rangle_{\mathcal{H}}, \|\cdot\|_{\mathcal{H}})$ with kernel $k$. We note that in practice the full knowledge of the phase-space is completely unreachable. Again, we therefore stress the fact that the setting of this encompassing RKHS $\mathcal{H}$ has only a theoretical purpose. The RKHS $\mathcal{H}$ can be connected to each RKHS of the family $\mathcal{W}$ through extension and restriction operators denoted $E_t$ and $R_t$ respectively, and defined as follows. For all $t \ge 0$, let

$$E_t: \begin{array}{ccc} \mathcal{H}_t & \to & \mathcal{H} \\ k_t(\cdot, X_t) & \mapsto & k(\cdot, X_t). \end{array} \quad (11)$$

As shown in Appendix B this operator is an isometry (hence continuous), and we extend this definition by linearity on $\text{Span}\{k_t(\cdot, X_t) : X_t \in \Omega_t\}$. Then, by density the function $E_t(f)$ is defined for all $f \in \mathcal{H}_t$. The restriction,

$$R_t: \begin{array}{ccc} \mathcal{H} & \to & \mathcal{H}_t \\ k(\cdot, X) = k(\cdot, \Phi_r(X_0)) & \mapsto & k(\cdot, X) \upharpoonright_{\Omega_t} = k_t(\cdot, \Phi_t(X_0)) \ell(t, r), \end{array} \quad (12)$$
which is also an isometry (Appendix B), is defined similarly for $g \in \mathcal{H}^{\text{sp}} = \text{Span}\{k(\cdot, X) : X \in \Omega\}$ and extended by density in $\mathcal{H}$ by the Moore-Aronszajn theorem [59]. The expression of $R_t(f)$ for $f \in \mathcal{H}$ is further specified in remark 19 (Appendix B). The extension map is built in such a way that each RKHS $\mathcal{H}_t$ of the family is included in the “big” encompassing RKHS $\mathcal{H}$.

In Appendix B, we list several useful properties of the restriction and extension operators. Namely, $E_t$ and $R_t$ are both isometries (Prop.12); they form an adjoint pair (Prop.13); the restriction is continuous in $L_{\mathbb{C}}^2(\Omega_t, \nu)$ (Prop.14). We define now the Koopman operator on the encompassing RKHS $\mathcal{H}$.

### 4.1.2. Koopman operators on $\mathcal{H}$

For all $t \ge 0$, we consider the Koopman operator $U_t$ defined by

$$U_t(f)(X) := f \circ \Phi_t(X) = f(\Phi_t(X)), \text{ for all } f \in \mathcal{H}. \qquad (13)$$

Since $\mathcal{H}$ is dense in $(L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$ (by proposition 11 (Appendix A.3), the operator $U_t$ can be continuously extended on $L_{\mathbb{C}}^2(\Omega, \nu)$ and to avoid notations inflation we keep noting $U_t$ this extension. We first study $U_t : L_{\mathbb{C}}^2(\Omega, \nu) \to L_{\mathbb{C}}^2(\Omega, \nu)$ with the $L_{\mathbb{C}}^2(\Omega, \nu)$ topology. The family $(U_t)_{t \ge 0}$ is a strongly continuous semi-group on $(L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$ since $t \mapsto \Phi_t(\cdot)$ is continuous on $\mathbb{R}^+$. As the feature maps are functions of $\mathcal{H}$, it can be noticed that for all $X_r \in \Omega$

$$U_t[k(\cdot, X_r)] = k(\Phi_t(\cdot), X_r), \qquad (14)$$

which justifies the stability of $\mathcal{H}$ by the operator $U_t$. This corresponds to a natural expression of the Koopman operator for any function $g \in \mathcal{H}^{\text{sp}}$, extended then by density to $\mathcal{H}$. Yet another useful equivalent expression of the Koopman operator is available for the feature maps. We have, for any points $X = \Phi_r(X_0), Y = \Phi_s(Y_0) \in \Omega$

$$U_t[k(\cdot, X)](Y) = k(\Phi_t(Y), X) = k(\Phi_{t+s}(Y_0), \Phi_r(X_0)) = k_0(Y_0, X_0)\ell(t+s, r). \qquad (15)$$

From the properties of the time kernel $\ell$, we get

$$\ell(t+s,r) = \varphi(t+s-r) = \varphi(s-(r-t)) = \ell(s,r-t), \qquad (16)$$

which leads to

$$k(\Phi_t(Y), X) = k_0(Y_0, X_0)\ell(t+s,r) = k_0(Y_0, X_0)\ell(s,r-t) = k(\Phi_s(Y_0), \Phi_{r-t}(X_0)) = k(Y, \Phi_t^{-1}(X)). \qquad (17)$$

and hence

$$U_t[k(\cdot, X)](Y) = k(Y, \Phi_t^{-1}(X)). \qquad (18)$$

As the previous equality is true for all $Y \in \Omega$, this implies that for all $X \in \Omega$

$$U_t[k(\cdot, X)] = k(\cdot, \Phi_t^{-1}(X)). \qquad (19)$$

This dual formulation of the kernel expression of the Koopman operator, which further highlights the stability of $\mathcal{H}$ by $U_t$ is intrinsically linked to the definition of the kernel $k$.

This dual expression will be of central interest in the following as it enables us to formulate the time evolution of the feature maps in terms of the Koopman operator $U_t$ and its adjoint at any time $t \ge 0$.

**Remark 3 (Transport of the kernel $k$)** For all $X = \Phi_r(X_0) \in \Omega$ and $t \ge 0$, by definition of the Koopman operator on the feature maps and (19), $U_t[k(\cdot, X)]$ has two expressions and we obtain

$$U_t[k(\cdot, X)] = k(\Phi_t(\cdot), X) = k(\cdot, \Phi_t^{-1}(X)).$$

The next remark provides a useful commutation property between $U_t$ and the kernel integral operator $\mathcal{L}_k$ defined in equation (1) or of its unique symmetric square-root $\mathcal{L}_k^{1/2}$ defined from the square-root of the kernel eigenvalues (see Appendix A.3 for a precise definition in the general case). Note that in the case of $\mathcal{H}$, the kernel integral operator is indeed a complex object which hides a time dependency.
**Remark 4 (Commutation between $\mathcal{L}_k^{1/2}$ (or $\mathcal{L}_k$) and $U_t$)** For all $X \in \Omega$, we have

$$ \mathcal{L}_k^{1/2} \circ U_t [k(\cdot, X)] = U_t \circ \mathcal{L}_k^{1/2} [k(\cdot, X)]. $$

**PROOF** By (14), we have $\mathcal{L}_k^{1/2} \circ U_t [k(\cdot, X)] = \mathcal{L}_k^{1/2} [k(\Phi_t(\cdot), X)]$. Then we obtain, for all $Y \in \Omega$

$$ (\mathcal{L}_k^{1/2} \circ U_t [k(\cdot, X)]) (Y) = (\mathcal{L}_k^{1/2} [k(\Phi_t(\cdot), X)]) (Y) = U_t [\mathcal{L}_k^{1/2} k(\cdot, X)](Y), $$

by definition of the Koopman operator since $\mathcal{L}_k^{1/2} k(\cdot, X) \in \mathcal{H}$. $\square$

<figure><img src="image_2.png" /></figure>

This commutation property that ensues directly from the compositional nature of the Koopman operator allows us to write immediately the equality

$$ \mathcal{L}_k \circ U_t [k(\cdot, X)] = U_t \circ \mathcal{L}_k [k(\cdot, X)]. \qquad (20) $$

By linearity these properties extend to all functions of $\mathcal{H}$. Let us note that the function $\mathcal{L}_k^{-1/2} k(\cdot, X)$, for all $X \in \Omega$, does not necessarily belong to $\mathcal{H}$, therefore this proof cannot be applied to $\mathcal{L}_k^{-1/2}$.

The next proposition shows the Koopman operator defined on RKHS $\mathcal{H}$ is unitary in $L_{\mathbb{C}}^2(\Omega, \nu)$, which is a classical property of the Koopman operator in $L^2$ for measure preserving invertible systems. The proof we give for this result relies on the transport kernel expression given in Remark 3.

**Proposition 1 (Unitarity of the Koopman operator in $L_{\mathbb{C}}^2(\Omega, \nu)$)** The map $U_t : (L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)}) \mapsto (L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$ is unitary for all $t \ge 0$.

**PROOF** Since the dynamical system $(\Phi_t)_{t \ge 0}$ is measure preserving, (i.e. $\forall A \in \Omega, \nu(A) = \nu(\Phi_t^{-1}(A))$, with $\Phi_t^{-1}(A)$ denoting the pre-image of $\Phi_t$) and left invariant (i.e. $\Phi_t(\Omega) = \Omega$, with $\mathrm{Ran}(\Phi_t) = \Omega$), the Koopman operator $U_t$ is an isometry: for all $f$ and $g \in L_{\mathbb{C}}^2(\Omega, \nu)$,

$$ \langle U_t(f), U_t(g) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle f, g \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. $$

Besides, it can be noticed that $U_t(\mathcal{H}^{\mathrm{sp}})$, with $\mathcal{H}^{\mathrm{sp}} = \mathrm{Span}\{k(\cdot, X) : X \in \Omega\}$, is a subset of the range of $U_t : L_{\mathbb{C}}^2(\Omega, \nu) \mapsto L_{\mathbb{C}}^2(\Omega, \nu)$, $U_t(\mathcal{H}^{\mathrm{sp}}) = \mathrm{Span}\{k(\cdot, \Phi_t^{-1}(X)) : X \in \Omega\}$ and we have the following inclusions

$$ U_t(\mathcal{H}^{\mathrm{sp}}) \subset \mathrm{Ran}(U_t) \subset L_{\mathbb{C}}^2(\Omega, \nu). $$

We show now that $U_t(\mathcal{H}^{\mathrm{sp}})$ is dense in $(L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$. Let $f \in L_{\mathbb{C}}^2(\Omega, \nu)$ be such that $\langle f, k(\cdot, \Phi_t^{-1}(X)) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = 0$ for all $X \in \Omega$. Then, from the definition of $\mathcal{L}_k$, it means that $(\mathcal{L}_k f)[\Phi_t^{-1}(X)] = 0$. The injectivity of $\mathcal{L}_k$ ensures that $f \circ \Phi_t^{-1} = 0$, so $U_{-t}f = 0$. As $U_{-t}$ is an isometry, it is in particular injective, which eventually leads to $f = 0$. By a consequence of Hahn-Banach theorem this result shows that the range of $U_t$ is dense in $L_{\mathbb{C}}^2(\Omega, \nu)$ and Proposition 1 is proved. $\square$

<figure><img src="image_3.png" /></figure>

**Proposition 1**, shows that the Koopman operator $U_t$ is invertible and that its inverse in $L_{\mathbb{C}}^2(\Omega, \nu)$ is $U_t^*$. Denoting by $P_t$ the operator defined by $P_t := U_t^*$, for all $f$ and $g \in L_{\mathbb{C}}^2(\Omega, \nu)$ we have $\langle U_t(f), g \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle f, P_t(g) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}$. The family $(P_t)_{t \ge 0}$ is a strongly continuous semi-group on $(L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$. This operator is referred to as the Perron-Frobenius operator. For all $t \ge 0$, the Perron-Frobenius operator $P_t$ verifies for all $X \in \Omega$

$$ P_t[k(\Phi_t(\cdot), X)] = k(\cdot, X). \qquad (21) $$

From Remark 1, we can write a more explicit expression of $P_t$ on the feature maps: for all $X \in \Omega$,

$$ P_t[k(\cdot, X)] = k(\cdot, \Phi_t(X)). \qquad (22) $$

Informally, if we see the function $k(\cdot, X)$ as an atom of the measure, then its expression at a future time is provided by (22), which corresponds well to the idea that the Perron-Frobenius operator advances in time the density.

As previously stated the Koopman operator $U_t$ is an isometry in $L_{\mathbb{C}}^2(\Omega, \nu)$. The next result asserts that $U_t$ is also an isometry in $\mathcal{H}$.
**Proposition 2 (Isometric relation of the kernel)** For all X and Y ∈ Ω, we have $k(\Phi_t(X), \Phi_t(Y)) = k(X, Y)$.

PROOF This results follows immediately from the kernel definition: for all $X = \Phi_r(X_0), Y = \Phi_s(Y_0) \in \Omega$, we have

$$k(\Phi_t(X), \Phi_t(Y)) = k(\Phi_{t+r}(X_0), \Phi_{t+s}(Y_0)) = k_0(X_0, Y_0)\ell(t+r, t+s) = k_0(X_0, Y_0)\ell(r, s) = k(X, Y). \quad (23)$$

Remark 3 on the application of the operator $U_t$ to $k(\cdot, X)$, with $X \in \Omega$, shows that applying the flow, $\Phi_t$, on the first variable of $k(\cdot, \cdot)$, is equivalent to applying $\Phi_t^{-1}$ to the second variable and vice-versa. Consequently, for all $f \in L_{\mathbb{C}}^2(\Omega, \nu)$ and $X \in \Omega$, we have $\langle f, k(\cdot, \Phi_t(X)) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle f, k(\Phi_t^{-1}(\cdot), X) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle f, P_t k(\cdot, X) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}$ and by definition of the adjoint, we obtain

$$\langle f, k(\cdot, \Phi_t(X)) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle U_t(f), k(\cdot, X) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. \quad (24)$$

We already knew that this equality was right for $f \in \mathcal{H}$ and for the inner product in $\mathcal{H}$, namely

$$\langle f, k(\cdot, \Phi_t(X)) \rangle_{\mathcal{H}} = f(\Phi_t(X)) = U_t f(X) = \langle U_t(f), k(\cdot, X) \rangle_{\mathcal{H}}. \quad (25)$$

Equation (24) provides a weak (in the sense that the $L_{\mathbb{C}}^2(\Omega, \nu)$-inner product against a feature map is no longer the evaluation function) formulation of the transport of any observable in $L_{\mathbb{C}}^2(\Omega, \nu)$ by the flow.

In order to further exhibit several analytical results on the Koopman operator $(U_t)_{t \ge 0}$ in $L_{\mathbb{C}}^2(\Omega, \nu)$, we introduce in the following its infinitesimal generator.

### 4.1.3. Koopman infinitesimal generator

We will note as $(A_U, \mathcal{D}(A_U))$ the infinitesimal generator of the strongly continuous semigroup $(U_t)_{t \ge 0}$ on $L_{\mathbb{C}}^2(\Omega, \nu)$ and its domain. As the Koopman and Perron-Frobenius operators are adjoint in $L_{\mathbb{C}}^2(\Omega, \nu)$ their infinitesimal generators are also adjoint of each other with possibly their own domain. The following lemma characterizes first the Perron-Frobenius infinitesimal generator $A_p$ and its domain $\mathcal{D}(A_p)$.

**Lemma 1 (Perron-Frobenius infinitesimal generator)** The Perron-Frobenius infinitesimal generator is the unbounded operator, $(A_p, \mathcal{D}(A_p))$ defined by

$$\mathcal{D}(A_p) = \{ f \in L_{\mathbb{C}}^2(\Omega, \nu) : x \mapsto \partial_{\mathcal{M}(x)} f(x) \in L_{\mathbb{C}}^2(\Omega, \nu) \} \text{ and } A_p f := -\partial_{\mathcal{M}(\cdot)} f(\cdot) \text{ for } f \in \mathcal{D}(A_p),$$

where $\mathcal{M}$ denotes the differential operator of the system dynamics (3) and $\partial_{\mathcal{M}}f$ stands for the directional derivative of $f$ along $u \in \Omega$.

PROOF See Appendix C.

For all $t \ge 0$, $P_t$ is the adjoint of $U_t$ in $L_{\mathbb{C}}^2(\Omega, \nu)$. The Koopman infinitesimal generator in $L_{\mathbb{C}}^2(\Omega, \nu)$ is consequently given by

$$\mathcal{D}(A_U) = \mathcal{D}(A_p^*) \quad \text{and} \quad A_U = A_p^*. \quad (26)$$

The Koopman and Frobenius-Perron operators are unitary in $L_{\mathbb{C}}^2(\Omega, \nu)$, as $\mathcal{D}(A_U)$ is dense in $L_{\mathbb{C}}^2(\Omega, \nu)$, by Stone's theorem, their infinitesimal generators are therefore skew symmetric for $\langle \cdot, \cdot \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}$ and we have

$$\mathcal{D}(A_p^*) = \mathcal{D}(A_p) \quad \text{and} \quad A_p^* = -A_p, \quad (27)$$

$$\mathcal{D}(A_U^*) = \mathcal{D}(A_U) \quad \text{and} \quad A_U^* = -A_U. \quad (28)$$

The adjoint has to be understood in the topology of $L_{\mathbb{C}}^2(\Omega, \nu)$. The two next propositions – well known in $L_{\mathbb{C}}^2(\Omega, \nu)$ for invertible measure preserving systems –, summarize (26), (27) and (28).

**Proposition 3 (Koopman infinitesimal generator)** The Koopman infinitesimal generator of $(U_t)_{t \in L_{\mathbb{C}}^2(\Omega, \nu)}$ is the unbounded operator, $(A_U, \mathcal{D}(A_U))$, defined by

$$\mathcal{D}(A_U) = \{ f \in L_{\mathbb{C}}^2(\Omega, \nu) : x \mapsto \partial_{\mathcal{M}(x)} f(x) \in L_{\mathbb{C}}^2(\Omega, \nu) \} \text{ and } A_U f = \partial_{\mathcal{M}(\cdot)} f(\cdot) \text{ for } f \in \mathcal{D}(A_U).$$
**Proposition 4 (Skew symmetry of the generators)** The Koopman infinitesimal generator $A_U$ is skew-symmetric in $L_{\mathbb{C}}^2(\Omega, \nu)$.

We set now a continuity property of the Koopman infinitesimal generator on a subspace of $\mathcal{D}(A_U)$. Let $f \in \mathcal{H}$, by theorem 6 (Appendix A.2), $x \mapsto \partial_{\mathcal{H}(x)}f(x) \in L_{\mathbb{C}}^2(\Omega, \nu)$, and through Cauchy-Schwarz inequality we get that $f \in \mathcal{D}(A_U)$ since the real valued function $\mathcal{M}$ belongs to $L_{\mathbb{C}}^2(\Omega, \nu)$ and consequently $\mathcal{H} \subset \mathcal{D}(A_U)$.

**Theorem 2 (Continuity of the Koopman generator on $\mathcal{H}$)** The restriction of the Koopman infinitesimal generator

$$A_U : (\mathcal{H}, \|\cdot\|_{\mathcal{H}}) \mapsto (L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$$

is continuous.

**PROOF** For all $f \in \mathcal{H}$ using the Cauchy-Schwarz inequality as above and theorem 7 (Appendix A.2) we have

$$\|A_U f\|_{L_{\mathbb{C}}^2(\Omega, \nu)} = \|\partial_{\mathcal{H}}(\cdot) f(\cdot)\|_{L_{\mathbb{C}}^2(\Omega, \nu)} \le C \|f\|_{\mathcal{H}}.$$

□

Through the dual expression of the Koopman operator (remark 3) the infinitesimal generator $A_U$ provides a dynamical system specifying the time evolution of the feature maps as, for all $X \in \Omega$, the feature map $k(\cdot, \Phi_t(X))$ verifies

$$\partial_t U_t k(\cdot, X) = \partial_t k(\Phi_t(\cdot), X) = A_U k(\Phi_t(\cdot), X). \quad (29)$$

We have the following useful differentiation formulae.

**Proposition 5 (Differentiation formulae)** For all $f \in \mathcal{H}$ (with $f = E_t g, g \in \mathcal{H}_t$) and $X_t \in \Omega_t$ we have

$$\langle \mathcal{L}_k \circ \partial_t f, k(\cdot, X_t) \rangle_{\mathcal{H}} = -\langle f, \mathcal{L}_k \circ \partial_t k(\cdot, X_t) \rangle_{\mathcal{H}}.$$

**PROOF** As $\mathcal{L}_k \circ A_U f$ belongs to $\mathcal{H}$, we have $\mathcal{L}_k \circ A_U f(X_t) = \mathcal{L}_k \circ \partial_t f(X_t)$ and

$$\langle \mathcal{L}_k \circ \partial_t f, k(\cdot, X_t) \rangle_{\mathcal{H}} = \langle \mathcal{L}_k \circ A_U f, k(\cdot, X_t) \rangle_{\mathcal{H}} = \langle j(\mathcal{L}_k^{1/2} \circ A_U f), \mathcal{L}_k^{-1/2} k(\cdot, X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}.$$

As $\mathcal{L}_k^{1/2}$ is self-adjoint in $L_{\mathbb{C}}^2(\Omega, \nu)$ and $A_U$ is skew-symmetric in $L_{\mathbb{C}}^2(\Omega, \nu)$, we obtain

$$\langle \mathcal{L}_k \circ \partial_t f, k(\cdot, X_t) \rangle_{\mathcal{H}} = \langle A_U f, k(\cdot, X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = -\langle f, A_U k(\cdot, X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = -\langle f, \partial_t k(\cdot, X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}.$$

As $\langle f, \partial_t k(\cdot, X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle f, \mathcal{L}_k \circ \partial_t k(\cdot, X_t) \rangle_{\mathcal{H}}$, this complete the proof. □

As already outlined, the whole phase space, $\Omega$, and the global embedding RKHS, $\mathcal{H}$, defined on it are both completely inaccessible for high dimensional state spaces. Instead of seeking to reconstruct this global RKHS, we will work in the following with time-varying “local” RKHS spaces built from a small (w.r.t. the phase space dimension) ensemble of initial conditions. The set of these time-varying spaces forms the RKHS family. To express the time evolution of the features maps associated to the RKHS family we now define an appropriate expression of the Koopman operator on this family.

### 4.1.4. Derivation of Koopman operator expression in $\mathcal{W}$

To fully specify the Koopman operators in the RKHS family, we rely on the family of extension, restriction mappings $(E_t)_{t \ge 0}$ and $(R_t)_{t \ge 0}$ (11, 12) relating the “big” encompassing RKHS $\mathcal{H}$ to the family of time-evolving RKHS $\mathcal{H}_t$ and use the Koopman operators $U_t : L_{\mathbb{C}}^2(\Omega, \nu) \mapsto L_{\mathbb{C}}^2(\Omega, \nu)$ for $t \ge 0$. The adjoints $P_t$ will also be very helpful as well as Remark 1 on the dual relation of the application of the flow on the global kernel $k$.

From now on, we note $X_t := \Phi_t(X_0)$ for all $X_0 \in \Omega_0$ with $X_t$ belonging to $\Omega_t$. We define the Koopman operator in the RKHS family by $\mathcal{U}_t := R_t \circ P_t \circ E_0$ for all $t \ge 0$. For all $X_0 \in \Omega_0$, we have

$$\mathcal{U}_t[k_0(\cdot, X_0)] = R_t \circ P_t[k(\cdot, X_0)] = R_t k(\cdot, X_0) = k_t(\cdot, X_0), \quad (30)$$
where the second equality is due to (22), and the third equality holds true from the definition of $R_t$ and the fact that $\ell(t,t) = 1$. The Koopman operator in the RKHS family, $\mathcal{U}_t : \mathcal{H}_0 \mapsto \mathcal{H}_t$, transports the kernel feature maps on the RKHS family by composition with the system's dynamics. It inherits some of the nice properties of the Perron-Frobenius operator defined on the encompassing global RKHS. As shown by the following theorem, proposition 2 remains valid for the family of kernels $(k_t)_{t \ge 0}$ and $\mathcal{U}_t$ is still unitary in the sense of the following theorem.

**Theorem 3 (Koopman RKHS isometry)** The Koopman operator on the RKHS family defines an isometry from $\mathcal{H}_0$ to $\mathcal{H}_t$: for all $X_0$ and $Y_0 \in \Omega_0$

$$ \langle \mathcal{U}_t k_0(\cdot, Y_0), \mathcal{U}_t k_0(\cdot, X_0) \rangle_{\mathcal{H}_t} = \langle k_0(\cdot, Y_0), k_0(\cdot, X_0) \rangle_{\mathcal{H}_0}. $$

The range of $\mathcal{U}_t : \mathcal{H}_0 \mapsto \mathcal{H}_t$ is dense in $\mathcal{H}_t$.

**PROOF** Upon applying equations (B.2 of Prop. 12, Appendix B), which states that $R_t : \mathcal{H} \to \mathcal{H}_t$ is an isometry, we have

$$ \langle \mathcal{U}_t k_0(\cdot, Y_0), \mathcal{U}_t k_0(\cdot, X_0) \rangle_{\mathcal{H}_t} = \langle P_t k(\cdot, Y_0), P_t k(\cdot, X_0) \rangle_{\mathcal{H}} = \langle k(\cdot, \Phi_t(Y_0)), k(\cdot, \Phi_t(X_0)) \rangle_{\mathcal{H}} = k(\Phi_t(X_0), \Phi_t(Y_0)). $$

With proposition 2 we obtain

$$ \langle \mathcal{U}_t k_0(\cdot, Y_0), \mathcal{U}_t k_0(\cdot, X_0) \rangle_{\mathcal{H}_t} = k(X_0, Y_0) = k_0(X_0, Y_0) = \langle k_0(\cdot, Y_0), k_0(\cdot, X_0) \rangle_{\mathcal{H}_0}. $$

Besides, the image of $\text{Span}\{k_0(\cdot, Y_0) : Y_0 \in \Omega_0\}$ under $\mathcal{U}_t$ is the set $\text{Span}\{k_t(\cdot, Y_t) : Y_t \in \Omega_t\}$ which is dense in $(\mathcal{H}_t, \|\cdot\|_{\mathcal{H}_t})$. The range of $\mathcal{U}_t : \mathcal{H}_0 \mapsto \mathcal{H}_t$ is thus dense in $\mathcal{H}_t$ and the operator $\mathcal{U}_t$ is, in that sense, unitary. $\square$

Let us now determine the adjoint of the Koopman operator in the RKHS family. Let $\mathcal{P}_t := R_0 \circ U_t \circ E_t$ for all $t \ge 0$ and let $X_t = \Phi_t(X_0) \in \Omega_t$. With the same arguments as for $\mathcal{U}_t$, we have

$$ \mathcal{P}_t[k_t(\cdot, X_t)] = k_0(\cdot, X_0). \qquad (31) $$

The mapping $\mathcal{P}_t : \mathcal{H}_t \mapsto \mathcal{H}_0$ with $t \ge 0$ constitutes the Perron-Frobenius family of operators in the RKHS family. The mapping $\mathcal{P}_t : \mathcal{H}_t \mapsto \mathcal{H}_0$ is unitary for the RKHS family topology (isometry from $\mathcal{H}_t$ to $\mathcal{H}_0$ and the range of $\mathcal{P}_t$ is dense in $\mathcal{H}_0$). The next proposition justifies that $\mathcal{U}_t$ and $\mathcal{P}_t$ have inverse action on the feature maps.

**Proposition 6 (Koopman Perron-Frobenius duality)** For all $X_0 \in \Omega_0$ and $Y_t = \Phi_t(Y_0) \in \Omega_t$, we have

$$ \langle \mathcal{U}_t k_0(\cdot, X_0), k_t(\cdot, Y_t) \rangle_{\mathcal{H}_t} = \langle k_0(\cdot, X_0), \mathcal{P}_t[k_t(\cdot, Y_t)] \rangle_{\mathcal{H}_0}. $$

**PROOF** By the equations (30) and (31), we have

$$ \langle \mathcal{U}_t k_0(\cdot, X_0), k_t(\cdot, Y_t) \rangle_{\mathcal{H}_t} = k_t(Y_t, \Phi_t(X_0)) \quad \text{and} \quad \langle k_0(\cdot, X_0), \mathcal{P}_t[k_t(\cdot, Y_t)] \rangle_{\mathcal{H}_0} = k_0(Y_0, X_0), $$

and both terms are equal by construction because of Definition 1. $\square$

In order to derive the Koopman and Perron-Frobenius operators' spectral representation in the RKHS family, we exhibit now two propagation operators that will allow us to express the evolution of the feature maps in $\mathcal{W}$.

### 4.1.5. The RKHS family spectral representation

We specify hereafter a family of operator $A_{U,t}$, related to the Koopman infinitesimal generator $A_U$, and that give rise to an evolution equation on $\Omega_t$ akin to (29). They will play, in that sense, the role of infinitesimal generators on the RKHS family.

For all $t \ge 0$, let $A_{U,t}$ be defined by $A_{U,t} := R_t \circ A_U \circ E_t$ with $E_t : \mathcal{H}_t \mapsto \mathcal{H}$ and $R_t : L_{\mathbb{C}}^2(\Omega_t, \nu) \mapsto L_{\mathbb{C}}^2(\Omega_t, \nu)$ (B.4 Appendix B).

**Proposition 7 (Continuity of $A_{U,t}$)** The mapping $A_{U,t} : (\mathcal{H}_t, \|\cdot\|_{\mathcal{H}_t}) \mapsto (L_{\mathbb{C}}^2(\Omega_t, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega_t, \nu)})$ is well-defined and continuous.
PROOF Upon applying (B.1 of Prop. 12, Appendix B), the continuity of $A_U$ (Theorem 2) and of $R_t$ (Proposition 14, Appendix B) we obtain the proof of the proposition. □

In a very similar way as the infinitesimal generator on $\mathcal{H}$, the operator $A_{U,t}$ can be understood as an evolution equation of the feature maps defined on $\Omega_t$ and associated to $\mathcal{H}_t$. For $X_t \in \Omega_t$, by proposition 3 it can be noticed that

$$A_{U,t} k_t(\cdot, X_t) = R_t \circ A_U [k(\cdot, X_t)] = R_t [\partial_{\mathcal{H}(\cdot)} k(\cdot, X_t)]$$

and in particular on $\Omega_t$

$$A_{U,t} k_t(\cdot, X_t) = \partial_{\mathcal{H}(\cdot)} k_t(\cdot, X_t). \qquad (32)$$

**Remark 5** Operator $A_{U,t}$ is alike the Koopman infinitesimal generator $A_U$. As a matter of fact by proposition 3 and theorem 2 we have for all $X \in \Omega$

$$A_U k(\cdot, X) = \partial_{\mathcal{H}(\cdot)} k(\cdot, X),$$

and for all $X_t \in \Omega_t$

$$A_{U,t} k_t(\cdot, X_t) = \partial_{\mathcal{H}(\cdot)} k_t(\cdot, X_t).$$

Notice that these two operators are different and act on different domains.

**Remark 6** Operators $A_{U,t}$ can be understood as an evolution operator for the feature maps. As a matter of fact through (29) and the definition of $A_{U,t}$, we have, for all $X \in \Omega_0$

$$\partial_t k(\cdot, \Phi_t(X))|_{\Omega_t} = A_{U,t} k_t(\cdot, \Phi_t(X)). \qquad (33)$$

Through the above remark the operator $A_{U,t}$ inherits the properties of operator $A_U$ defined on the global encompassing RKHS $\mathcal{H}$. As shown in the next proposition it remains in particular skew-symmetric.

**Proposition 8 (Skew-symmetry of $A_{U,t}$)** The operator $A_{U,t}$ is skew-symmetric in $L_{\mathcal{H}}^2(\Omega_t, \nu)$: for all $f$ and $g \in \mathcal{H}_t$

$$\langle A_{U,t} f, j \circ g \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)} = - \langle f, A_{U,t} g \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)}.$$

where $j$ is the injection $j: \mathcal{H}_t \mapsto L_{\mathcal{H}}^2(\Omega_t, \nu)$.

For the need of this proof only, we introduce the following notation. For a function $h$ defined on $\Omega_t$, we note $\tilde{h}$ the extension on $\Omega$ defined by

$$\tilde{h} := \begin{cases} h & \text{on } \Omega_t \\ 0 & \text{else.} \end{cases}$$

PROOF As $g = E_t(g)$ holds on $\Omega_t$, we have

$$\langle A_{U,t} f, j \circ g \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)} = \langle A_U \circ E_t(f) |_{\Omega_t}, E_t(g) \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)} = \langle A_U \circ \widetilde{E_t}(f), E_t(g) \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)}.$$

By proposition 4 and since $E_t(f) |_{\Omega_t} = f$ on $\Omega_t$ we infer that

$$\langle A_{U,t} f, j \circ g \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)} = - \langle E_t(f), A_U \circ \widetilde{E_t}(g) \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)} = - \langle f, A_{U,t}(g) \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)} = - \langle f, A_{U,t}(g) \rangle_{L_{\mathcal{H}}^2(\Omega_t, \nu)}.$$

□

We are now ready to prove the RKHS family spectral representation theorem, which states that the bounded operator $A_{U,t}: \mathcal{H}_t \mapsto L_{\mathcal{H}}^2(\Omega_t, \nu)$ is diagonalizable for all $t \ge 0$.
### 4.1.6. Proof of Theorem 1 [The RKHS family spectral representation]

The full proof, organized in two main steps, is thoroughly detailed in Appendix D. For this proof, we consider the restriction of the Koopman infinitesimal generator on $\mathcal{H}$, $A_U : \mathcal{H} \to L_{\mathbb{C}}^2(\Omega, \nu)$, which is connected to each $A_{U,t}$ for $t \ge 0$ through the restriction operator $R_t$. In the first step, the diagonalization of operator $A_U : \mathcal{H} \to L_{\mathbb{C}}^2(\Omega, \nu)$ in $L_{\mathbb{C}}^2(\Omega, \nu)$ is first performed. To that end we introduce an intermediate (approximating) operator denoted $\tilde{A}_U$, directly related to $A_U$ and whose inverse is shown to be compact and self-adjoint. The second step of the proof consists in deducing the diagonalization of each $A_{U,t}$ from the diagonalization of $A_U$ obtained at step 1. These two steps are thoroughly presented in Appendix D.

## 4.2. Tangent linear dynamics

A result of practical interest concerns the establishment of a rigorous ensemble expression of the tangent linear dynamics operator. Recall that we note $X_t = \Phi_t(X_0)$ for all $t \ge 0$. We define $\delta X(t, \cdot) := \Phi_t[X_0(\cdot) + \delta X_0(\cdot)] - \Phi_t[X_0(\cdot)]$ for all $t \ge 0$, and where $\delta X_0(x)$ is a perturbation of the initial condition at point $x$. The function $\delta X(t, \cdot) \in L^2(\Omega_x, \mathbb{R}^d)$ is the perturbation of the flow at time $t$ with respect to the initial condition $X_0$. We have

$$ \Phi_t(X_0 + \delta X_0) - \Phi_t(X_0) = \int_0^t \mathcal{M}[\Phi_s(X_0 + \delta X_0)] - \mathcal{M}[\Phi_s(X_0)] ds + o(\|\delta X_0\|) = \int_0^t d(\mathcal{M} \circ \Phi_s(X_0)) \delta X_0 ds + o(\|\delta X_0\|) $$

and we obtain also that

$$ \Phi_t(X_0 + \delta X_0) - \Phi_t(X_0) = \int_0^t d\mathcal{M}[\Phi_s(X_0)] [d\Phi_s(X_0)\delta X_0] ds + o(\|\delta X_0\|) = \int_0^t \partial_{\mathcal{X}}\mathcal{M}[\Phi_s(X_0)] \delta X(s, \cdot) ds + o(\|\delta X_0\|). $$

The variation of the flow verifies $\delta X(t, \cdot) = \int_0^t \partial_{\mathcal{X}}\mathcal{M}[X_s]\delta X(s, \cdot) ds$ almost everywhere on $\Omega_x$. Recall that $\mathcal{M}$ is assumed $C^1$, in particular $\sup_{X \in \Omega} \partial_{\mathcal{X}}\mathcal{M}(X) < \infty$ (since $\Omega$ is compact). The function $\delta X(t, \cdot)$ belongs to $L^2(\Omega_x, \mathbb{R}^d)$ and verifies for all $t \ge 0$

$$ \partial_t \delta X(t, \cdot) = \partial_{\mathcal{X}}\mathcal{M}[X_t]\delta X(t, \cdot). \quad (34) $$

Each component $[\delta X(t, \cdot)]_i$ of $\delta X(t, \cdot)$ belongs to $L^2(\Omega_x, \mathbb{R})$. Let $[g_x^\epsilon(\delta X_t)]_i : \Omega \to \mathbb{C}$ be the representation of $[\delta X(t, \cdot)]_i$ in $L_{\mathbb{C}}^2(\Omega, \nu)$ defined for all $x \in \Omega_x$ by

$$ [g_x^\epsilon(\delta X_t)]_i := [\delta X_t(x)]_i. $$

For any $x \in \Omega_x$, the family of functions $g_x^\epsilon(X_t)$ correspond to vector observables of the dynamical system and we may immediately write for all $t \ge 0$

$$ \partial_t g_x^\epsilon(\delta X_t) = \partial_{\mathcal{X}}\mathcal{M}(X_t)g_x^\epsilon(\delta X_t). \quad (35) $$

We are now ready to exhibit a kernel expression of the tangent linear operator $\partial_{\mathcal{X}}\mathcal{M}(X_t)$. For all $1 \le i \le d$, the function $\mathcal{L}_k[g_x^\epsilon]_i \in \mathcal{H}$ verifies for all $\delta X_t \in \Omega$

$$ \langle \mathcal{L}_k \circ \partial_t \mathcal{L}_k [g_x^\epsilon]_i, k(\cdot, \delta X_t) \rangle_{\mathcal{H}} = \mathcal{L}_k (\partial_t \mathcal{L}_k [g_x^\epsilon]_i)(\delta X_t). \quad (36) $$

Upon applying the differentiation formulae of proposition 5 and the evolution equation (29), we have

$$ \begin{aligned} \langle \mathcal{L}_k \circ \partial_t \mathcal{L}_k [g_x^\epsilon]_i, k(\cdot, \delta X_t) \rangle_{\mathcal{H}} &= -\langle \mathcal{L}_k [g_x^\epsilon]_i, \mathcal{L}_k \circ \partial_t k(\cdot, \delta X_t) \rangle_{\mathcal{H}} \\ &= -\langle \mathcal{L}_k [g_x^\epsilon]_i, \mathcal{L}_k [A_U k(\cdot, \delta X_t)] \rangle_{\mathcal{H}} \\ &= -\langle \mathcal{L}_k^{1/2} [g_x^\epsilon]_i, \mathcal{L}_k^{1/2} [A_U k(\cdot, \delta X_t)] \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. \end{aligned} $$

The operator $\mathcal{L}_k^{1/2}$ being self-adjoint in $L_{\mathbb{C}}^2(\Omega, \nu)$ and $A_U$ skew-symmetric for the inner product of $L_{\mathbb{C}}^2(\Omega, \nu)$, we have

$$ \begin{aligned} \langle \mathcal{L}_k \circ \partial_t \mathcal{L}_k [g_x^\epsilon]_i, k(\cdot, \delta X_t) \rangle_{\mathcal{H}} &= \langle A_U \mathcal{L}_k [g_x^\epsilon]_i, k(\cdot, \delta X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} \\ &= \langle \mathcal{L}_k^{1/2} A_U \mathcal{L}_k [g_x^\epsilon]_i, \mathcal{L}_k^{-1/2} k(\cdot, \delta X_t) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} \\ &= \langle \mathcal{L}_k A_U \mathcal{L}_k [g_x^\epsilon]_i, k(\cdot, \delta X_t) \rangle_{\mathcal{H}} \\ &= \mathcal{L}_k A_U \mathcal{L}_k [g_x^\epsilon]_i(\delta X_t). \end{aligned} $$
Combining the right-hand side of the above expression with (36), we have

$$ \mathcal{L}_k(\partial_t \mathcal{L}_k[g_x^\epsilon]_i)(\delta X_t) = \mathcal{L}_k(A_U \mathcal{L}_k[g_x^\epsilon]_i)(\delta X_t). \quad (37) $$

As $\mathcal{L}_k$ is injective, the kernel of the tangent linear operator in the RKHS $\mathcal{H}$ reads:

$$ (\partial_t \mathcal{L}_k[g_x^\epsilon]_i)(\delta X_t) = A_U \mathcal{L}_k[g_x^\epsilon]_i(\delta X_t). $$

As $\Omega$ is bounded, the function $t \mapsto \delta X(t, \cdot)$ belongs to $L^\infty([0, T])$ and therefore $\delta X_t \mapsto g_x^\epsilon(\delta X_t)$ belongs to $L^\infty([0, T], L^2(\Omega, \nu))$. By (35), we have also that $\delta X_t \mapsto \partial_t g_x^\epsilon(\delta X_t)$ belongs to $L^\infty([0, T], L^2(\Omega, \nu))$ since $\sup_{Y \in \Omega} \partial_X \mathcal{M}(Y) < \infty$. We have therefore on $\Omega$

$$ \partial_t \mathcal{L}_k[g_x^\epsilon]_i(\delta X_t) = \int_{\Omega} k(\delta X_t, z) \partial_t [g_x^\epsilon]_i(z) \nu(dz) = \mathcal{L}_k(\partial_t [g_x^\epsilon]_i)(\delta X_t), $$

and thus,

$$ \partial_t [g_x^\epsilon]_i(\delta X_t) = \mathcal{L}_k^{-1} A_U \mathcal{L}_k[g_x^\epsilon]_i(\delta X_t). $$

By the commutation property of remark 4, together with (35) and as $\mathcal{L}_k^{1/2}[g_x^\epsilon]_i \in \mathcal{H}$, the following equalities hold for all $\delta X_t \in \Omega$

$$ \begin{gathered} \partial_t [g_x^\epsilon]_i(\delta X_t) = \mathcal{L}_k^{-1/2} A_U \mathcal{L}_k^{1/2} [g_x^\epsilon]_i(\delta X_t), \\ \partial_X \mathcal{M}(X_t) g_x^\epsilon(\delta X_t) = \mathcal{L}_k^{-1/2} A_U \mathcal{L}_k^{1/2} g_x^\epsilon(\delta X_t). \end{gathered} \quad (38) $$

Let's specify now the kernel expression of the tangent linear dynamics in $\Omega_s$ for all $s \ge 0$. The function $\mathcal{L}_k^{1/2}[g_x^\epsilon]_{\Omega_s}$ belongs to $\mathcal{H}_s$. We have as well $\mathcal{L}_k^{1/2}(\partial_t g_x^\epsilon) = \mathcal{L}_{k_s}^{1/2}(\partial_t g_x^\epsilon)$ on $\Omega_s$. For all $Y_s$ and $X'_s \in \Omega_s$ we have that $A_U k(\cdot, X'_s)(Y_s) = A_{U,s} k_s(\cdot, X'_s)(Y_s)$. In particular, we get $A_U \mathcal{L}_k^{1/2}[g_x^\epsilon](Y_s) = A_{U,s} \mathcal{L}_{k_s}^{1/2}[g_x^\epsilon](Y_s)$. We obtain hence the kernel expression for all $Y_s \in \Omega_s$

$$ \partial_X \mathcal{M}(X_s) g_x^\epsilon(Y_s) = \mathcal{L}_{k_s}^{-1/2} A_{U,s} \mathcal{L}_{k_s}^{1/2} g_x^\epsilon(Y_s). \quad (39) $$

Note that the domains of the kernel expressions of the tangent linear (38) and (39) are different. The right-hand side of (39) provides a convenient kernel expression of the tangent linear operator, enabling us to evaluate the tangent linear dynamics from an ensemble of feature maps. The adjoint of the tangent linear dynamics is straightforwardly given by:

$$ \partial_X \mathcal{M}^*(X_s) g_x^\epsilon(Y_s) = - \mathcal{L}_{k_s}^{-1/2} A_{U,s} \mathcal{L}_{k_s}^{1/2} g_x^\epsilon(Y_s). \quad (40) $$

**Remark 7 (Projection observables)** The point observable functions $g_x^\epsilon$ used above can be extended to other functions defined from basis $(\Psi_j)_{j \ge 0}$ of $L^2(\Omega_x, \mathbb{R}^d)$, with

$$ \begin{aligned} g_x^\psi X_t &:= \sum_j^\infty \langle X_t, \Psi_j \rangle_{L^2(\Omega_x)} \Psi_j(x), \\ g_x^\psi : \Omega &\mapsto \sum_j^\infty \langle \cdot, \Psi_j \rangle_{L^2(\Omega_x)} \Psi_j(x) \in \mathbb{C}. \end{aligned} $$

**Remark 8 (Dependence on $X_t$)** It should be noted that in (39) and (40) the right-hand side does not depend on $X_s$ whereas the tangent linear operator or its adjoint on the left-hand side does. Function $X_s$ indicates around which function of $\Omega_s$ the nonlinear system is linearized. The equivalent of the infinitesimal generator of the Koopman operator, representing the dynamics' linear tangent operator on the RKHS family, depends necessarily also on this function. This dependence is here implicit and induced by the considered sampled functions used to define the RKHS $\mathcal{H}_t$. If the set of members $(X_t^{(i)})_{i \ge 1}$ are centered around a particular function $X_s$, the operators $A_{U,t}$ can be interpreted as a representation of the tangent linear operator around function $X_s$. In ensemble methods, $X_s$, is in general taken as the ensemble mean, and $(X_s^{(i)})_{i \ge 1}$ is an ensemble of time-dependent perturbations around this mean.

**Remark 9 (More regularity on $\delta X_t$)** For all $1 \le i \le d$, if we suppose that the function $[g_x^\epsilon]_i$ belongs to $\mathcal{H} \subset L^2_{\mathbb{C}}(\Omega, \nu)$, the proof can be simplified and (39) is replaced by

$$ \partial_X \mathcal{M}(X_s) g_x^\epsilon(\delta X_s) = A_{U,s} g_x^\epsilon(\delta X_s). \quad (41) $$
It can be pointed out that the expression above corresponds to the approximation of the tangent linear dynamics used in ensemble method if we work in a finite dimensional space such that $\Omega \subset \mathbb{R}^n$ and assume that $k_t(\cdot, X_t)$ is defined as $(p-1)^{-1/2}\langle(X_t - \bar{X}_t), \cdot - \bar{X}_t\rangle_{\mathbb{R}^n}$, with $\bar{X}_t$ the empirical ensemble mean. With that definition, we have $A_{U,t}k_t(\cdot, X^{(i)}) = \sum_j \partial_{X_j} k_t(\cdot, X_t^{(i)}).\mathcal{M}(\cdot)_j = \langle(\delta_{X_t^{(i)}} - \frac{1}{p}\sum_l \delta_{X_t^{(l)}}).\mathcal{M}(\cdot), \cdot - \bar{X}_t\rangle_{\mathbb{R}^n}$ that reads $(p-1)^{-1/2}\langle\mathcal{M}(X_t^{(i)}) - \frac{1}{p}\sum_j \mathcal{M}(X_t^{(j)}), \cdot - \bar{X}_t\rangle_{\mathbb{R}^n}$, for $i=1, \dots, p$ and for which, when associated to the $\mathbb{R}^n$ Euclidean inner product on a resolution grid of size $n$, the left-hand side of this latter expression corresponds to the so-called $(p \times n)$ anomaly matrix built from $p$ ensemble members of the dynamical system. The tangent linear approximation provided by ensemble methods can be thus immediately interpreted as a particular instance of feature maps together with a given choice of specific inner product to define the reproducing kernel. Keeping a finite dimensional approximation but working without assuming that the functions $[g_x^\epsilon]$ belong to $\mathcal{H}$, and thus with now expression (38) for the ensemble tangent linear expression, corresponds to the case in which a localization procedure identified to the square root operator $\mathcal{L}_{k_t}^{1/2}$ has been considered. These two choices embed the problem within a particular RKHS family of functions. The relation between the tangent linear dynamics and the anomaly matrix is in our case exact and does not correspond to a finite difference approximation as classically presented in ensemble methods. The RKHS family can be seen as a way of linearizing locally a nonlinear system in a convenient sequence of spaces of smooth functions.

### 4.3. Finite time Lyapunov exponents

The kernel of the Koopman operator provides also a direct access to the finite time Lyapunov exponents. Recalling from (35) that for any punctual observable $g_x$, as defined previously, we have

$$g_x(\delta X_t) = \int_0^t \partial_{X_t} \mathcal{M}(X_s) g_x(\delta X_s) ds + o(\|\delta X_0\|). \quad (42)$$

With the expression of the tangent linear operator in terms of the evolution operator (38) on $\Omega_s$, we have then

$$g_x(\delta X_t) = \int_0^t \mathcal{L}_{k_s}^{-1/2} A_{U,s} \mathcal{L}_{k_s}^{1/2} g_x(\delta X_s) \, ds. \quad (43)$$

For all $s \in [0,t]$, we consider a perturbation $g_x(\delta X_s) = \mathcal{L}_{k_s}^{-1/2} \psi_\ell^s(\delta X_s)$ along a Koopman generator eigenfunction associated to the eigenvalue of maximal modulus $|\lambda_\ell|$. By Theorem 1, we have

$$\partial_t g_x(\delta X_s) = \lambda_\ell \mathcal{L}_{k_s}^{-1/2} \psi_\ell^s = \lambda_\ell g_x(\delta X_s).$$

Therefore we get

$$|g_x(\delta X_t)| = e^{\lambda_\ell |t|} |g_x(\delta X_0)|, \quad (44)$$

and the finite time Lyapunov exponent is consequently defined as

$$\sigma = |\lambda_{\ell}|. \quad (45)$$

For regular perturbations $g_x \in \mathcal{H}_s$ with unitary perturbation $g_x(\delta X_s) = \psi_\ell^s(\delta X_s)$, the derivation is even simpler as we obtain from remark 9

$$\partial_t g_x(\delta X_s) = A_{U,t} g_x(\delta X_s),$$

which yields directly to expression (44) and to the same expression for the Lyapunov exponent.

The modulus of the larger Koopman eigenvalue in the RKHS family provides thus an estimate of the Lyapunov exponent. It can be outlined that the computation of Lyapunov exponents for large scale systems is computationally very demanding as it requires the construction of the linear tangent dynamics operator and the solution of an eigenvalue problem of very big dimension. The construction of the exact numerical tangent linear operator is in general a tedious task when expressed in $L^2(\Omega_x)$, as in equation (34). The ensemble based method provided by our formalism is on the contrary very simple by expressing the tangent linear operator in $L^2_{\mathbb{C}}(\Omega_t, \nu)$, as in equation (35). It can be noticed that, by this change of norm in the definition, the computed values are not the same.

Three distinct values can then be defined for practical computations. First, the Lyapunov spectrum expressed in $L^2_{\mathbb{C}}(\Omega_t, \nu)$ can be determined by computing the singular values of $\mathcal{L}_{k_t}^{-1/2} A_{U,t} \mathcal{L}_{k_t}^{1/2}$. The time integral is dropped since
the evaluations are constant along trajectories. It can be viewed as an advantage of working in $L^2_{\mathbb{C}}(\Omega_t, \nu)$ instead of $L^2(\Omega_x)$ as performed classically. The time independence is due to the fact that the Koopman operator is intrinsic to the dynamical system. However, from a numerical point of view, as the computation is performed in practice through an ensemble with a limited number of members, the learned spectrum is representative only of the local dynamics at the time ($t = t_0$) at which the kernel has been evaluated. As an alternative, modal Lyapunov exponents can be defined by the square root of the first singular values of $\mathcal{L}_{k_t}^{-1/2} \psi_\ell |\lambda_\ell|^2 \psi_\ell^* K_t^{-1} \mathcal{L}_{k_t}^{1/2} = \mathcal{L}_{k_t}^{-1/2} \psi_\ell |\lambda_\ell|^2 \psi_\ell^* \mathcal{L}_{k_t}^{-1/2}$ with $K_t(i, j) = k_t(X_t^{(i)}, X_t^{(j)})$. We call these singular values the *Koopman modal Lyapunov exponents* (KMLE). Finally, as a third option, equation (45) can simply be considered to evaluate modal exponents. We can notice that the two modal Lyapunov exponents definitions are very similar; the former being expressed in $L^2_{\mathbb{C}}(\Omega_t, \nu)$ and the latter in $\mathcal{H}_t$. An example of the estimation of the two types of modal exponents will be presented for a quasi-geostrophic dynamics in the numerical section.

### 4.4. Practical considerations

Let us stress again that in practice, we only have access to the mappings $k_t$ and $A_{U,t}$ with $t \ge 0$. The mapping $k$ and $A_U$ are completely inaccessible for high-dimensional systems as they require the complete knowledge of the phase space or at least of a long enough orbit with a density assumption in the whole phase space. This last assumption is associated to strong requirements of the dynamical system and is not necessarily valid for a given time series of a particular observable. Instead of working with an infinite (dense) trajectory, Theorem 1 enables us to estimate the eigenvalues and eigenfunctions of the Koopman operator locally in the RKHS family, which can locally conveniently be accessed from an ensemble of finite time trajectories. As it will be described in the following, operators $A_{U,t}$, can be discretized as an ensemble matrix – itself related, as we saw it previously, to the tangent linear dynamics operator. This matrix is then diagonalized to get access to Koopman eigenvalues and their associated eigenfunctions. In theory, the diagonalization of $A_{U,t}$ needs to be performed only once, at a given time, to access the Koopman eigenpairs $(\psi_\ell^t)_\ell$ and $(\lambda_\ell)_\ell$. However, the exponential relation between distinct instants allows us also to consider averaging strategies to eventually robustify the estimation in practice. This capability will be exploited in the numerical section for the data assimilation of time series.

#### Diagonalization in practice

For all $t \ge 0$, let $m_t$ be the kernel expression of the operator $A_{U,t}$ given by $m_t(X_t, Y_t) := A_{U,t}^*[k_t(\cdot, Y_t)](X_t)$ for all $X_t$ and $Y_t \in \Omega_t$. By Proposition 8, we have $m_t(X_t, Y_t) = -A_{U,t}[k_t(\cdot, Y_t)](X_t)$.

Let us denote by $\{X_t^{(i)} : 1 \le i \le n\}$ an ensemble of members generated by the dynamical system and by $\{k_t(\cdot, X_t^{(j)}) : 1 \le i \le n\}$ the $n$ associated feature maps. For all time $t \ge 0$, these $n$ feature maps enable us to build a kernel expression of operator $A_{U,t}$ as the $n \times n$ matrix $M_t = (m_t(X_t^{(i)}, X_t^{(j)}))_{1 \le i,j \le n}$ with:

$$ (M_t)_{ij} := -A_{U,t}[k_t(\cdot, X_t^{(j)})](X_t^{(i)}) = [-\partial_{\mathcal{H}}(\cdot) k_t(\cdot, X_t^{(j)})](X_t^{(i)}). $$

As shown in the following, this matrix enables us to access to the Koopman generator eigenvalues and to the evaluation of the eigenfunctions at the ensemble members.

By definition of $A_{U,t}$, we have $A_{U,t}[k_t(\cdot, X_t^{(j)})] = A_U[k(\cdot, X_t^{(j)})]$. Otherwise, as $-\mathcal{L}_k^{1/2} \circ A_U[k(\cdot, X_t^{(j)})]$ belongs to $\mathcal{H}$, by the proof (Step 1) (Appendix D) of Theorem 1 and denoting $(\psi_\ell)_\ell$ an orthonormal basis of $\mathcal{H}$ sets from the eigenfunctions of $A_U$, we get

$$ \begin{aligned} -\mathcal{L}_k^{1/2} \circ A_U[k(\cdot, X_t^{(j)})] &= -\sum_{\ell=0}^{\infty} \langle \mathcal{L}_k^{1/2} \circ A_U[k(\cdot, X_t^{(j)})], \psi_\ell \rangle_{\mathcal{H}} \psi_\ell, \\ &= -\sum_{\ell=0}^{\infty} \langle A_U k(\cdot, X_t^{(j)}), \mathcal{L}_k^{-1/2} \psi_\ell \rangle_{L^2_{\mathbb{C}}(\Omega, \nu)} \psi_\ell. \end{aligned} $$

Noting $\beta_\ell$ the eigenvalues of $\mathcal{L}_k$, we have also (step 1 Appendix D theorem 1) $\mathcal{L}_k j(\psi_\ell) = \beta_\ell j(\psi_\ell)$ and

$$ -\mathcal{L}_k^{1/2} \circ A_U[k(\cdot, X_t^{(j)})] = -\sum_{\ell=0}^{\infty} \beta_\ell^{-1/2} \langle A_U k(\cdot, X_t^{(j)}), j(\psi_\ell) \rangle_{L^2_{\mathbb{C}}(\Omega, \nu)} \psi_\ell. $$
By the skew symmetry of the generators (proposition 4) we have that

$$ -\mathcal{L}_k^{1/2} \circ A_U [k(\cdot, X_t^{(j)})] = \sum_{\ell=0}^{\infty} \lambda_\ell \beta_\ell^{-1/2} \langle k(\cdot, X_t^{(j)}), j(\psi_\ell) \rangle_{L^2(\Omega, \nu)} \psi_\ell = \sum_{\ell=0}^{\infty} \lambda_\ell \beta_\ell^{-1/2} \overline{\mathcal{L}_k j(\psi_\ell)(X_t^{(j)})} \psi_\ell, $$

which leads to

$$ -\mathcal{L}_k^{1/2} \circ A_U [k(\cdot, X_t^{(j)})] = \sum_{\ell=0}^{\infty} \lambda_\ell \beta_\ell^{1/2} \overline{\psi_\ell(X_t^{(j)})} \psi_\ell, $$

and, applying $\mathcal{L}_k^{-1/2}$ on both sides, we get

$$ -A_U [k(\cdot, X_t^{(j)})] = \sum_{\ell=0}^{\infty} \lambda_\ell \beta_\ell^{1/2} \overline{\psi_\ell(X_t^{(j)})} \mathcal{L}_k^{-1/2} \psi_\ell = \sum_{\ell=0}^{\infty} \lambda_\ell \overline{\psi_\ell(X_t^{(j)})} j(\psi_\ell). $$

By the restriction expression (remark 19, Appendix B), and keeping only $n$ eigen functions to represent the values at the $n$ members, we finally obtain the following equality for all $t \ge 0$

$$ m_t(X_t^{(i)}, X_t^{(j)}) = \sum_{\ell=0}^{n} \lambda_\ell \psi_\ell^t(X_t^{(i)}) \overline{\psi_\ell^t(X_t^{(j)})}, $$

which shows that the diagonalization of $M_t$ provides a set of the Koopman generator eigenpairs in the RKHS family.

In practice, the skew-symmetric matrix $M_t := (m_t(X_t^{(i)}, X_t^{(j)}))_{1 \le i,j \le n}$ is assembled from the definition of $A_{U,t}$ and a given choice of the kernel. As explained in the previous section this matrix corresponds to a kernel expression of $A_{U,t}$ with

$$ M_t(X_t^{(i)}, X_t^{(j)}) := -A_{U,t}[k_t(\cdot, X_t^{(j)})](X_t^{(i)}) = [-\partial_{\mathcal{H}(\cdot)} k_t(\cdot, X_t^{(j)})](X_t^{(i)}). $$

This matrix can be interpreted in the RKHS setting as resulting from the matrix multiplication:

$$ M_t(X_t^{(i)}, X_t^{(j)}) := \sum_{\ell} \mathbb{F}_t(j, \ell) K_t(\ell, i), $$

with $\mathbb{F}_t(i, \ell) = -A_{U,t}[k_t(X_t^{(i)}, X_t^{(j)})]$. This indeed corresponds to a discretization of the kernel expression of operator $A_{U,t}$ through the empirical Dirac measure. Numerically, instead of working with matrix $M_t$ (i.e. the evaluation of $A_{U,t}k_t(\cdot, X)$ at several discrete points), we will work directly with matrix $\mathbb{F}_t = M_t K_t^{-1}$. This has the advantage of directly working with an implicit discretization of operator $A_{U,t}$, and to relax somewhat its dependency on the kernel choice. The skew-symmetric matrix $M_t K_t^{-1}$ is then diagonalized through a direct numerical procedure (using LAPACK library and working numerically on the anti-symmetric part of $M_t K_t^{-1}$) and can be written $M_t K_t^{-1} = V_t \Lambda V_t^T$, with $V_t$ a unitary matrix and $\Lambda$ a diagonal matrix. The matrix $V_t$ gathers eigenvectors $V_t^j$ of $\mathbb{F}$, which is a discretization of $A_{U,t}k_t(\cdot, X)$ giving access to the values of the Koopman eigenfunctions $\psi_t^j = \mathcal{L}_k^{1/2} V_t^j$ at the $n$ ensemble members points $X_t^{(i)}$. The matrix $\Lambda$ is composed of Koopman eigenvalues with conjugate pairs of pure imaginary eigenvalues.

As previously mentioned this diagonalization can be performed at a single time or at several instants accompanied with an averaging procedure. Theorem 1 and (8), give access to the eigenvectors evaluation along a trajectory for all time instants.

We provide below as examples, expressions of the evaluation of $M_t(X_t^{(i)}, X_t^{(j)}) = (-\partial_{\mathcal{H}(\cdot)} k_t(\cdot, X_t^{(j)}))(X_t^{(i)})$ for the empirical covariance kernel ($k_E$) and the Gaussian kernel ($k_G$). The empirical covariance kernel is defined through the kernel isometry property as

$$ k_t(X_t^{(i)}, X_t^{(j)}) = k_E(X_0^{(i)}, X_0^{(j)}) = \langle X_0^{(i)}, X_0^{(j)} \rangle_{\Omega_x}, $$

with $(\cdot, \cdot)_{\Omega_x}$ the inner product of $L^2(\Omega_x, \mathbb{R}^d)$ and where $\Omega_x$ denotes the physical domain of the considered dynamics. We obtain

$$ M_t(X_t^{(i)}, X_t^{(j)}) = \left\langle \frac{\partial}{\partial t} X_0^{(j)}, X_0^{(i)} \right\rangle_{\Omega_x}. \qquad (46) $$
In this expression, we see that the time derivative of the ensemble members at the initial time is required. Similarly, the Gaussian kernel is defined as

$$k_t(X_t^{(i)}, X_t^{(j)}) = k_G(X_0^{(i)}, X_0^{(j)}) = \exp\left(-\frac{1}{\ell_G^2} \|X_0^{(i)} - X_0^{(j)}\|_{\Omega_x}^2\right),$$

This leads to

$$M_t(X_t^{(i)}, X_t^{(j)}) = -\frac{2}{\ell_G^2} \left\langle \frac{\partial}{\partial t} X_0^{(j)}, (X_0^{(i)} - X_0^{(j)}) \right\rangle_{\Omega_x} \exp\left(-\frac{1}{\ell_G^2} \|X_0^{(i)} - X_0^{(j)}\|_{\Omega_x}^2\right). \quad (47)$$

Thanks to the isometry property, this matrix needs to be evaluated only at a single time.

# 5. Application to ensemble methods

## 5.1. Quasi-geostrophic model

Based on a set of p ensemble members, the ability of the proposed method to predict a trajectory associated with a new initial condition is demonstrated on a quasi-geostrophic barotropic model of a double-gyre [69] in a “rectangular-box ocean”. This idealistic model of North-Atlantic large-scale circulation considers the transport of the ageostrophic state variable components by a velocity in geostrophic balance. The adimensioned barotropic vorticity equation written for the potential vorticity q and the stream function ψ in Ωx = ([0 : L] × [-L : L]), with L = 1, is

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial q}{\partial t} + J(\psi, q) = f - (\delta/L)^5 \nabla^4 \omega \\[1.5ex] \Delta\psi = \omega \quad ; \quad q = R_o\omega + y_c, \end{array} \right. \qquad (48)$$

where ω is the vorticity and $x_c = (x_c, y_c)$ are the latitude and the longitude increasing eastward and northward, respectively, while $J(f,g) = \partial_x f \partial_y g - \partial_x g \partial_y f$ stands for the Jacobian (with x and y denoting horizontal and vertical coordinates of $\mathbb{R}^2$). In the model (48), the state of the system can be uniquely determined from the potential vorticity q. The set of all possible states of the system is then $\Omega \subset L^2(\Omega_x, \mathbb{R})$ (i.e. $d=1$). The dimensionless Rossby number $R_o = 0.0036$ measures the ratio between inertial forces and the Coriolis force, and a steady wind-forcing is defined as $f = \sin(\pi y_c)$. The parameter $\delta/L = 0.032$ controls the hyper-viscosity model, with $\delta$ the Munk boundary layer width. Scales smaller than $\delta$ will not be resolved. Numerical details and links with dimensional variables associated with a mid-latitude ocean basins, such as North Atlantic, as well as the way the initial conditions of the ensemble have been generated are given in Appendix E. Let us note that this system is ergodic when associated to an uncorrelated random forcing [70, 71, 72]. Here the forcing is stationary, and the system is not ergodic. As a consequence, such a system cannot be used for all the techniques relying on an ergodic assumption to extract a spectral representation of the Koopman operator.

## 5.2. Reconstructions using RKHS

In this study we focus only on two specific kernels: the empirical covariance kernel ($k_E$), and the Gaussian Kernel ($k_G$). They are defined respectively by

$$k_E(X_t^{(i)}, X_t^{(j)}) = \left\langle q^{(i)}(\cdot, t_0), q^{(j)}(\cdot, t_0) \right\rangle_{\Omega_x},$$

and

$$k_G(X_t^{(i)}, X_t^{(j)}) = \exp \frac{1}{\ell_G^2} \|q^{(i)}(\cdot, t_0) - q^{(j)}(\cdot, t_0)\|_{\Omega_x}^2,$$

where $\langle \cdot, \cdot \rangle_{\Omega_x}$ stands for the standard $L^2$ inner product in $\Omega_x$ and $q^{(i)}(x,t)$ is the solution of equation (48) of the $i^{th}$ ensemble member at time $t$ at the space location $x$. The former has the advantage to be simple, and to be directly related to standard ensemble methods employed in data assimilation. However, it is strongly rank deficient and introduces spurious correlations in the physical domain $\Omega_x$ between two far apart locations. Localization techniques [73] have been introduced to tackle this problem and may be interpreted in our setting as embeddings in specific RKHS. This
interpretation may open new computationally efficient generalisations related for instance to localization techniques in the model space as explained after remark 9. The parallel done between operators $\mathcal{L}_{k_t}^{1/2}$ and localization procedure should in particular enable to define time-varying and flow dependant localization techniques that are notably difficult to built in practice. Gaussian kernels are widely used in machine learning for its high regularity and non-compact support in $L^2(\Omega, \nu)$ [74]. It can be also interpreted as a way to perform localization. In this study, compared to other tested – but not presented – kernels (polynomial and exponential for instance), it led to better performances.

Once the kernel matrix $K_{ij} = k_t(X_t^{(i)}, X_t^{(j)})$ for $(i, j) \in [1 \dots p] \times [1 \dots p]$ (equal for all $t$ – kernels isometry, Theorem 3) is evaluated, a very simple reconstruction of the whole trajectory can be performed. Denoting $X_0'$ a new ensemble member initial condition, the linear combination vector $\beta$ associated to the reconstructed trajectory $\hat{X}_t' = \sum_{i=1}^p \beta_i X_t^{(i)}$ can be estimated by kernel regression from the evaluation $k(X_0', X_0^j)$, for $j = [1 \dots p]$, and the reproducing property as

$$ \beta = K^{-1}k(X_0', X_0^{[1 \dots p]}). \quad (49) $$

Strikingly, Koopman isometry implies that these coefficients are constant along time, and brings hence a direct reconstruction of the whole trajectory.

### 5.3. Evaluation of the Koopman eigenfunctions

Performing the $\mathcal{H}_t$ inner product of equation (7) with the ensemble members' feature maps leads to the matrix

$$ \mathbb{M}(X_t^{(i)}, X_t^{(j)}) = (-\partial_{\mathcal{H}(\cdot)} k_t(\cdot, X_t^{(j)})) (X_t^{(i)}). \quad (50) $$

As explained previously eigen-elements of $MK^{-1}$, noted $(\lambda_\ell, V_\ell(X_t^{(i)}))$, give access to the exact evaluation of Koopman eigenfunctions on the RKHS family at the ensemble members $\psi_\ell'(X_t^{(i)}) = \mathcal{L}_k^{1/2}V_\ell(X_t^{(i)})$, and to the associated eigenvalues $\lambda_\ell$. We note that $\partial_{\mathcal{H}(\cdot)}k_t(\cdot, X_t^{(j)})$ is obtained by differentiating the kernel expression (see equations 46, 47 for the two selected kernels). Skew-symmetry of $MK^{-1}$, inherited from the skew-symmetry of $A_{U,t}$, is only approached numerically. Evaluations are thus enhanced by considering only the skew-symmetric part of the matrix $MK^{-1}$, which is dominant in the numerical evaluation of (50).

### 5.4. Lyapunov times

The Lyapunov times evaluated using the two kernels $k_E$ and $k_G$ are compared. Three methods to compute the Lyapunov exponents $\sigma$ are detailed in section 4.3. The associated Lyapunov times are defined by $T = \sigma^{-1} \ln C$, with, $C = 10^3$.

The Lyapunov times associated with exponents expressed in $L_\mathcal{C}^2(\Omega_t, \nu)$ are shown in figure 1(a). The spectrum associated with all singular values, truncated such that $\sigma > 10^{-5}$, is displayed. Due to this truncation, very slow modes are not well captured. We can see that the Gaussian kernel $k_G$ estimates in average larger Lyapunov times.

In figure 1(b), the Koopman modal Lyapunov times are shown as a function of the absolute value of the Koopman eigenvalue. The Lyapunov times expressed with the RKHS norm equation (45) is directly $T_\ell = 1/|\lambda_\ell| \ln C$, which explains the perfect decaying slope in logarithmic scale. We observe that the computation of the modal exponents in $L_\mathcal{C}^2(\Omega_t, \nu)$ leads simply to a vertical shift of the curves due to the scaling factor induced by the change of norm. We can notice that with the Gaussian kernel $k_G$, slower Koopman eigenfunctions are evaluated compared to $k_E$, which leads to longer Lyapunov times. This explains longer predictions, which will be presented in section 5.5. As it will be seen, the modal exponents and their associated Lyapunov times will reveal very useful, as they will allow us to filter out sequentially the contributions of the Koopman eigenfunctions beyond their predictability time.

### 5.5. Reconstructions

The ability of trajectories reconstruction is illustrated in figure 2, where a full trajectory is estimated from the embedding of a new initial condition in the RKHS family. We compare the anomaly (figure 2(a)) of a given new ensemble member at an advanced time (3.7 days), with reconstructions using the Gaussian kernel $k_G$ (figure 2(c)). Figure 2(b) shows the best possible reconstruction as a linear combination of ensemble members, with regards to the $L^2(\Omega_x)$-norm of the projection error.
<figure><img src="image_4.png" /><figcaption>(a) Lyapunov times spectrum expressed in $L_C^2(\Omega_t, \nu)$ for the Gaussian kernel $k_G$ and the empirical covariance $k_E$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Modal Lyapunov times associated with individual Koopman modes. Comparison with expression in $L_C^2(\Omega_t, \nu)$ and in $H_t$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2: Comparison between the reference potential vorticity of a test ensemble member at $t = 0.07$ (corresponding to 3.7 days of prediction) (a) and the predicted reconstruction using the Gaussian kernel ($l_G = 1$) (c). The best possible reconstruction using a linear combination of ensemble members with respect to the $L^2$ norm is shown in (b).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3: Mean estimation error using an empirical covariance kernel (short red dashes), a Gaussian kernel (green dash-dots), the reconstructed Gaussian kernel using all Koopman eigenfunctions (orange squares) and the Lyapunov filtered kernel (blue dot-dashes). “Ensemble average” (long black dashes) is the estimation error of the ensemble average. The projection error (purple solid line) corresponds to the optimal constant-in-time linear combination of ensemble members w.r.t. the L²(Ωₓ)-norm. The vertical width of the stripes regions over the mean error is equal to the standard deviation of the error. One dimensionless time t corresponds to 92 days. Errors are averaged over the 100 members of the test ensemble.</figcaption></figure>

Reconstruction errors are quantified by averaging the square of the L²(Ωₓ)-norm over the whole test ensemble composed of 100 new members (i.e. not used in the training phase). The standard deviations of the error are also displayed on this figure. The time evolution of this error is presented in figure 3. It is compared to the projection error defined as the minimal linear combinations of ensemble members with a constant-in-time coefficient w.r.t. L²(Ωₓ)-norm. The reconstruction error performed by the ensemble average is also shown. This error is obviously high for small times, but estimations are difficult to be more precise than the ensemble average estimation once the predictability time of the system has passed. The reconstructions performed with $k_E$, grow and overtake the ensemble average reconstruction around $t = 0.2$ (i.e. 18.4 days). On this basis, the use of $k_G$ ($l_G = 1.0$) leads to better results. Indeed, despite a slight reconstruction penalty at the initial time due to a stronger smoothing (see figure 2), the reconstruction errors become quickly significantly lower, and the error never reaches the ensemble anomaly error. This small loss of performance in L²(Ωₓ)-norm at the initial time allows to have a gain in robustness in terms of long-time prediction. Let us note that the lengthscale parameter $l_G$ is sensitive. In this work, we fixed it a priori. However, techniques that aim at learning such a parameter have been proposed in [5, 6].

The reconstructions based on (49) are valid for all times. Koopman isometry warrants the use of the same linear combination for the whole trajectory. Nevertheless, in practice the ensemble are of limited (small) size, and the eigenpairs of the infinitesimal generator close to zero (representing almost stationary functions) are difficult to estimate accurately. Modal predictability time associated to Lyapunov exponents (44) can be used to robustify the forecast of new initial conditions. The kernel can be expanded in the Koopman eigenfunctions basis as

$$k(X_t^{(i)}, X_t^{(j)}) = \sum_{\ell=1}^{\infty} \psi_{\ell}^{t}(X_t^{(i)}) \overline{\psi_{\ell}^{t}}(X_t^{(j)}) = \sum_{\ell=1}^{\infty} \psi_{\ell}^{0}(X_0^{(i)}) \overline{\psi_{\ell}^{0}}(X_0^{(j)}). \quad (51)$$

With the numerical procedure proposed here, the $p$ evaluated Koopman eigenfunctions allow to perform an exact reconstruction of the kernel. An appealing potential of equation (51) for faithful predictions consists in filtering the kernel by switching-off components that have passed their predictability time. We set this time from the finite-time Lyapunov exponents (44) as: $T_\ell = |\lambda_\ell|^{-1} \ln C$, with, $C = 10^3$. The filtered kernel $\tilde{k}(X_t^{(i)}, X_t^{(j)})$ and its associated $p \times p$ kernel matrix $\tilde{K}$ allows us to define reconstructions, with the modified linear coefficients $\tilde{\beta} = \tilde{K}^+ \tilde{K} \tilde{\beta}$. The superscript $\bullet^+$ stands for the Moore-Penrose pseudo-inverse. It can be remarked that the operator $\tilde{K}^+ \tilde{K}$ is a projector performing filtering directly on the vector of coefficients $\tilde{\beta}$.

In figure 3, it is assessed that the kernel is exactly reconstructed by the Koopman eigenfunctions (superimposition of the green dash-dots curve and the orange squares). When the kernel is Lyapunov filtered, the reconstruction almost fits at short time the previous predictions (superimposition of the blue dot-dashes curve, the green dash-dots curve and the orange squares). For longer time it significantly reduces the error and clearly separates the spread of the ensemble from the ensemble average estimation. For very long times (not shown in the figure), filtering acts on all modes,
which reduces to the ensemble average estimation. This strategy takes advantage of the modal predictability time. Slow eigenfunctions are used for long time predictions, and fast ones only for short time forecasts. As a final remark, we note that very fast modes (with large Koopman eigenvalues modulus) does not seem to contribute significantly to the reconstruction, since at short times the filtered and full reconstructions match.

## 5.6. Data assimilation

A possible application of the proposed procedure is ensemble-based data assimilation. In the context of satellite data, ensemble optimal interpolation (EnOI) is widely used to reconstruct spatial fields [75, 76]. It is for instance employed in operational gridded altimetry reconstructions [77]. EnOI consists in expressing at a given time the solution as a linear combination of ensemble members, which minimises an objective cost function. As shown in [78, 79], incorporating information from other times together with an adequate dynamical model significantly improves the estimations. This is obviously at the price of forecasting the whole ensemble by the dynamics.

Koopman isometry provided by the RKHS family allows us to incorporate time-series of observations without any new simulation of the dynamics. The spatial reconstruction is embarrassingly simply given as the constant-in-time best linear combination of ensemble members trajectories.

Based on the same ensemble as in the previous sections, we propose two demonstration tests considering i) a time-series of extremely sparse (two satellite swaths like) observations ii) an even sparser case with the same swaths but at a single time only. To that end, we consider a synthetic observation mimicking two along tracks measurements of an orbiting satellite, idealized by two vertical lines. In the first case these lines are assumed to be known every 0.9 days $\Delta t = 0.01$, while in the second case only day 5.5 is available ($t = 0.06$). The observations are perturbed by a centered Gaussian white noise with the value of the selected standard deviation equal to 10% of the signal root-mean-square. These 2 benchmarks are carried out for each of the 100 members of the test ensemble. As previously, the members of the test ensemble have not been used in the learning phase.

In this numerical experiment, we search for the constant in time coefficients $\beta$ (thanks to isometry), which minimises

$$ \mathcal{J} = \frac{1}{2T} \int_{0}^{T} \| \mathbb{H}(\hat{X}_{t}) - \mathcal{Y} \|_{R^{-1}}^{2} dt + \frac{\alpha^{2}}{2} \| \beta - \bar{\beta} \|^{2} \text{ s.t. } \mathbb{H}(\hat{X}_{t}) = \sum_{i=1}^{p} \beta_{i} \mathbb{H}(X_{t}^{(i)}), \quad (52) $$

where $\mathbb{H}$ is the observation operator associated with the satellite measurements $\mathcal{Y}$, $\hat{X}_t$ is the estimated state, and $\bar{\beta} = \frac{1}{p}(1, \dots, 1)^T$. The observation operator is assumed to have its range in the RKHS family and hence provides observables belonging to the RKHS family, which can be interpolated using the RKHS family kernels. It can be noticed that the expression of the constraint is due to kernel interpolation, and not to some linear property of $\mathbb{H}$. This treatment is similar to the one performed for instance in discrete empirical interpolation (DEIM) [80]. In the case of a single observation, the time integral is dropped, and the cost function corresponds to a standard EnOI cost function. The matrix $R$ is the observation covariance matrix, defined by the time average of the empirical covariance, localized by performing a Schur product with matrix $\exp(-r_{ij}^2/\ell_{loc}^2)$, where $r_{ij}$ is the distance between measurements $i$ and $j$, and $\ell_{loc} = 0.05$. The regularization parameter has been set to $\alpha = 10$. For this problem we selected the simple empirical kernel, defined from a non localized potential vorticity matrix as feature maps of the ensemble. Working with the isometry property justifies the estimation of constant in time coefficients for the linear combination. This simplistic method could be easily extended by introducing more elaborated kernels, Koopman eigenfunctions representation to project any system's observable and forecast them, as well as the Lyapunov filtering scheme presented in section (5.5).

Figure 4 shows the potential vorticity trajectory (panel a), the single-time observation (panel b), the trajectory reconstruction associated to this single observation (panel c) and the one obtained from the swaths time-series (panel d). Three consecutive times spaced of 0.02 (1.8 days) are displayed. We can see that despite extremely sparse measurements, a single observation already provides a fairly good estimation, together with forward and backward in time forecasts based only on the estimated linear combination of the ensemble members (i.e. without any simulation). As expected, a time-series of swaths improves further the performances.

In figure 5, the average of the errors obtained over all 100 members of the test ensemble are plotted. As a reference, the best possible reconstruction with respect to the $L^2$ norm with constant in time linear coefficients $\beta$ is shown, as well as the estimation by the ensemble average. The error performed with the time-series of satellite observations is almost optimal in terms of average and standard deviation. We note that for the second test-case, near the observation time,
<figure><img src="image_8.png" /><figcaption>Figure 4: Estimation of potential vorticity at times $t = 0.04$ (top), $t = 0.06$ (middle) and $t = 0.08$ (bottom). Panel (a): True anomaly. Panel (b): Observation associated with the single-time swath observation. Panel (c): Estimation knowing the two satellite swaths only at day 5.5 ($t = 0.06$). Panel (d): Estimation knowing the two swaths every 0.9 days ($\Delta t = 0.01$).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 5: Estimation errors averaged over the whole 100 members of the test ensemble. The best reconstruction error with constant coefficients (purple), estimation by the ensemble average (black), estimation with the time-series of satellite swaths observations (yellow) and only at day 5.5 ($t = 0.06$) (blue).</figcaption></figure>

the error is lower than the one obtained for the whole time-series. Compared to the time series of observations it may be interpreted as an over-fitting at the observation time since only a single observation is used to infer the ensemble members' coefficients. As expected, this error grows, but reasonably, as the estimation is propagated forward and backward in time. As a final remark, let us emphasize that considering time series of the sparse observables (in space and time) such as the satellite swaths considered in this study, which constitute typical altimetric observations of the ocean surface, would require an extremely long series to reconstruct the Koopman eigenpairs (see [34]) at the scale of a large basin such as the North Atlantic basin. For geophysical flows this can be highly problematic in the nowadays state of the fast climate change we are observing. The proposed technique, based on an ensemble of realizations instead of time series, constitutes a promising alternative to face these difficulties. We assessed here the technique with an idealistic dynamics model and more realistic models need to be tested. This will be the subject of future works.

## 6. Concluding remarks

We proposed a theoretical setting to estimate the Koopman eigen-pairs of compact dynamical systems in em-bedding the system in a RKHS family. Beyond the Koopman operator spectral representation, we showed that the RKHS family provides access to finite time Lyapunov exponents, as well as to the tangent linear dynamics. These features enabled us to propose very simple strategies for data assimilation of very sparse observations in time and space. Obviously more efficient data assimilation techniques could be used instead of the simple regression strategy explored here. As the RKHS family boils down to a linear setting it should be perfectly adapted to ensemble Kalman filters. Extension of classical ensemble data assimilation filters on the RKHS family are currently in development and will be the subject of future publications. The technique lies on the sequence of time-evolving RKHS to repre-sent the ensemble dynamics and on the superposition principle attached to the RKHS family. This study provides a theorem-based setting to estimate the Koopman eigenvalues and eigenfunctions, as well as exact ensemble-based representation of the tangent linear dynamics and of the flow Lyapunov exponents. This work relies on the assumption of an invertible deterministic dynamical system defined on a compact manifold. In future studies, we will focus on the generalisation to stochastic fluid flow dynamics with transport noise models [81, 82, 83, 84] and the introduction of more sophisticated ensemble data assimilation techniques [85, 86].

## Acknowledgments

The authors acknowledge the support of the ERC project 856408-STUOD and the French National program LEFE (Les Enveloppes Fluides et l'Environnement).
## Appendix A. Reproducing kernel Hilbert spaces

In the following appendices we provide complementary elements for the detailed proofs of the results listed in the paper. We first recall some classical theorems on RKHS and present their adaptations to our functional setting. In particular we present useful results on differentiability of RKHS functions and reproducing kernels. In the following $L^2(\Omega_x)$ stands for the Hilbert space of Lebesgue measurable $\mathbb{R}^d$-valued of square-integrable functions on the set $\Omega_x \subset \mathbb{R}^p$ $p = 2, 3$ endowed at a fixed time, $t$, with the inner product

$$ \langle f, g \rangle_{L^2(\Omega_x)} = \int_{\Omega_x} f(y) g(y) \, dy. $$

Let $E \subset L^2(\Omega_x)$ denote a compact subset equipped with the norm of $L^2(\Omega_x)$, by the Aronszajn Theorem, [59], the reproducing kernel Hilbert spaces (RKHS) of $\mathbb{C}$-valued functions defined on $E$ are uniquely characterized through a reproducing kernel $k: E \times E \to \mathbb{C}$. This kernel is a Hermitian, positive definite function. In the following we will assume furthermore it is $C^{(1,1)}(E \times E)$, which means here it is twice differentiable with respect to its two variables. The (Gateaux) derivative with respect to the first and second variable in the direction $v \in E$ are denoted $\partial_v^{(1)}k$, and $\partial_v^{(2)}k$, respectively. The reproducing kernel defines an evaluation function $k(\cdot, x) \in \mathcal{H}$ with a reproducing property (Def. 1):

$$ \forall x \in E, \forall f \in \mathcal{H}, \langle f, k(\cdot, x) \rangle_{\mathcal{H}} = f(x). $$

The RKHS $\mathcal{H}$ corresponds to a space of functions with a strong uniform convergence property. The elements of $\mathcal{H}$ are in addition continuous functions and we have for all $f \in \mathcal{H}$ that

$$ \forall x \in E \quad |f(x)| \le \sup_{y \in E} k(y, y)^{1/2} \|f\|_{\mathcal{H}}. \qquad (A.1) $$

### Appendix A.1. Functional description of $\mathcal{H}$

We recall in the following the Mercer's theorem which gives a functional description of the RKHS $\mathcal{H}$ and of its associated reproducing kernel $k$.

Let $L_{\mathbb{C}}^2(E, v)$ be the space of measurable complex valued functions defined on the compact set $E$ for which the square modulus is integrable with respect to measure $v$, $L_{\mathbb{C}}^2(E, v)$ is a Hilbert space equipped with the inner product $\langle \cdot, \cdot \rangle_{L_{\mathbb{C}}^2(E, v)}$ given for $f$ and $g \in L_{\mathbb{C}}^2(E, v)$ by

$$ \langle f, g \rangle_{L_{\mathbb{C}}^2(E, v)} := \int_E f(y) \overline{g(y)} v(dy), $$

so that it is linear in its first argument and antilinear in its second. Let us first define $\mathcal{L}_k: L_{\mathbb{C}}^2(E, v) \mapsto L_{\mathbb{C}}^2(E, v)$ the integral operator with kernel $k \in C^{(1,1)}(E \times E)$ as:

$$ (\mathcal{L}_k f)(x) = \int_E k(x, y) f(y) v(dy). \qquad (A.2) $$

We give below a useful remark relating the integral operator $\mathcal{L}_k$ and the inner product $\langle \cdot, \cdot \rangle_{L_{\mathbb{C}}^2(E, v)}$.

**Remark 10 (Links between $\mathcal{L}_k$ and $\langle \cdot, \cdot \rangle_{L_{\mathbb{C}}^2(E, v)}$)** As $k$ is Hermitian, it can be noticed that for all $x \in E$, $(\mathcal{L}_k f)(x) = \int_E f(y) \overline{k(y, x)} v(dy)$, and we have thus $(\mathcal{L}_k f)(x) = \langle f, k(\cdot, x) \rangle_{L_{\mathbb{C}}^2(E, v)}$.

This kernel integral operator can be shown to be well-defined and compact (by Arzelà Ascoli theorem). It is immediate to see that as $k$ is positive definite and Hermitian, $\mathcal{L}_k$ is positive and self-adjoint on $L_{\mathbb{C}}^2(E, v)$. By the spectral theorem of compact self-adjoint operators, there exists an orthonormal basis $(\phi_i)_i$ of $L_{\mathbb{C}}^2(E, v)$ consisting of eigenvectors of $\mathcal{L}_k$ with corresponding decreasing real eigenvalues $(\mu_i)_i: \mu_i \ge \mu_{i+1} \ge 0$ for all integers $i$. These results are proved in [58]. We consider in the following that $\mu_i \ne 0$ for all $i$ and $\mathcal{L}_k$ is injective. This denseness assumption could be relaxed (see remark 16) but simplifies significantly the presentation. The eigenpairs $(\mu_i, \phi_i)_{i \in \mathbb{N}}$ enables to formulate the following theorems (with proofs available for instance in [58] and [60]). Elaborating on remark 10 and the kernel eigenfunctions (often referred to as Mercer's eigenfunctions due to theorem 4) we get immediately relations between the eigenfunction evaluation and the $L^2(E, v)$ inner product.
**Remark 11 (Links between $(\varphi_i)_{i \ge 0}$ and $\langle \cdot, \cdot \rangle_{L^2_{\mathbb{C}}(E, v)}$)** For all $x \in E$, by Remark 10, we note that

$$ \langle \varphi_i, k(\cdot, x) \rangle_{L^2_{\mathbb{C}}(E, v)} = \mathcal{L}_k(\varphi_i)(x) = \mu_i \varphi_i(x) \quad \text{and} $$

$$ \langle k(\cdot, x), \varphi_i \rangle_{L^2_{\mathbb{C}}(E, v)} = \overline{\langle \varphi_i, k(\cdot, x) \rangle}_{L^2_{\mathbb{C}}(E, v)} = \mu_i \overline{\varphi_i(x)}. $$

**Remark 12** As $k \in C^{(1,1)}(E \times E)$, the range of $\mathcal{L}_k$ is included in $C^1(E)$ and the eigenfunctions of $\mathcal{L}_k$ are $C^1(E)$.

**Theorem 4 (Mercer decomposition)** For all $x$ and $y \in E$, $k(x,y) = \sum_{i=0}^{\infty} \mu_i \varphi_i(x) \overline{\varphi_i(y)}$ where the convergence is absolute (for each $x,y \in E \times E$) and uniform on $E \times E$.

In particular, the series $\sum \mu_i$ is convergent as

$$ \int_E k(y, y) v(dy) = \int_E \sum_{i=0}^{\infty} \mu_i \varphi_i(y) \overline{\varphi_i(y)} v(dy) = \sum_{i=0}^{\infty} \mu_i < \infty. \quad (A.3) $$

**Theorem 5 (Mercer representation of RKHS)** The RKHS $(\mathcal{H}, \langle \cdot, \cdot \rangle_{\mathcal{H}}, \|\cdot\|_{\mathcal{H}})$ is the Hilbert space characterized by

$$ \mathcal{H} = \{ f \in L^2_{\mathbb{C}}(E, v), f = \sum_{i=0}^{\infty} a_i \varphi_i : \sum_i \frac{|a_i|^2}{\mu_i} < \infty \}, $$

$$ \langle f, g \rangle_{\mathcal{H}} = \sum_{i=0}^{\infty} \frac{1}{\mu_i} \langle f, \varphi_i \rangle_{L^2_{\mathbb{C}}(E, v)} \overline{\langle g, \varphi_i \rangle}_{L^2_{\mathbb{C}}(E, v)} $$

$$ \text{and} \quad \|f\|_{\mathcal{H}}^2 = \sum_{i=0}^{\infty} \frac{|\langle f, \varphi_i \rangle_{L^2_{\mathbb{C}}(E, v)}|^2}{\mu_i}. $$

The inner product $\langle \cdot, \cdot \rangle_{\mathcal{H}}$ is sesqui-linear (linear in its first argument and antilinear in its second). The proof of Mercer's theorem warrants that the linear space spanned by $\{k(\cdot, x) : x \in E\}$ is dense in $L^2_{\mathbb{C}}(E, v)$. If $f = \sum_{i=0}^{\infty} a_i \varphi_i \in \mathcal{H}$ then the series converges absolutely and uniformly to $f$.

The following remark shows that the feature map's conjugate as well as its real and imaginary parts belong to $\mathcal{H}$.

**Remark 13 (Consequence of the Mercer representation)** For all $y \in E$, the functions $\overline{k(\cdot, y)}$, $\text{Re}[k(\cdot, y)]$ and $\text{Im}[k(\cdot, y)]$ belong to $\mathcal{H}$. By Theorems 5 and 4, we have indeed that

$$ \sum_{i=0}^{\infty} \frac{1}{\mu_i} |\langle \overline{k(\cdot, y)}, \varphi_i \rangle_{L^2_{\mathbb{C}}(E, v)}|^2 = \sum_{i=0}^{\infty} \frac{1}{\mu_i} |\mathcal{L}_k \varphi_i(y)|^2 = \sum_{i=0}^{\infty} \mu_i |\varphi_i(y)|^2 < \infty, \quad (A.4) $$

and $\overline{k(\cdot, y)} \in \mathcal{H}$. By operation, $\text{Re}[k(\cdot, y)]$ and $\text{Im}[k(\cdot, y)]$ belong thus to $\mathcal{H}$.

**Remark 14 (Comparison of the norms $\mathcal{H}$ and $L^2_{\mathbb{C}}(E, v)$)** For all $f \in \mathcal{H}$, by Theorem 5 we have for all $x \in E$, the inequality

$$ f(x) \le |f(x)| = \left| \sum_{i=0}^{\infty} a_i \varphi_i(x) \right| \le \left( \sum_{i=0}^{\infty} \frac{|a_i|^2}{\mu_i} \right)^{1/2} \left( \sum_{i=0}^{\infty} \mu_i |\varphi_i(x)|^2 \right)^{1/2}, $$

with theorem 4 we obtain

$$ \|f\|_{L^2_{\mathbb{C}}(E, v)} \le \left( \int_E |k(x, x)|^2 v(dx) \right)^{1/2} \|f\|_{\mathcal{H}}. \quad (A.5) $$
**Remark 15** For all $f \in L_{\mathbb{C}}^2(E, \nu)$, $\mathcal{L}_k(f)$ belongs to $\mathcal{H}$, in fact we have

$$ \sum_{i=0}^{\infty} \frac{1}{\mu_i} \left| \langle \mathcal{L}_k(f), \varphi_i \rangle_{L_{\mathbb{C}}^2(E, \nu)} \right|^2 = \sum_{i=0}^{\infty} \mu_i \left| \langle f, \varphi_i \rangle_{L_{\mathbb{C}}^2(E, \nu)} \right|^2 \le \left( \sum_{i=0}^{\infty} \mu_i \right) \|f\|_{L_{\mathbb{C}}^2(E, \nu)}^2 < \infty. $$

Moreover, by Remark 14 the mapping $\mathcal{L}_k : \mathcal{H} \to \mathcal{H}$ is continuous : for all $f \in \mathcal{H}$ we have

$$ \| \mathcal{L}_k(f) \|_{\mathcal{H}}^2 \le \left( \sum_{i=0}^{\infty} \mu_i \right) \| f \|_{L_{\mathbb{C}}^2(E, \nu)}^2 \le C \| f \|_{\mathcal{H}}^2. $$

With a $C^{(1,1)}(E \times E)$ regularity condition on the kernel, we will prove that the RKHS $\mathcal{H}$ is a subspace of smoother functions than $C^0(E, \mathbb{C})$ (theorem 6). In the following appendix we discuss differentiability property of the feature maps.

## Appendix A.2. Differentiability of feature maps and RKHS functions $\mathcal{H}$

In this appendix we extend results given in [87, 88, 56] for $E \subset \mathbb{R}^d$ to the case where $E$ is a compact subset of a function vector space. In the following, aligned with the assumptions enunciated in the paper's preamble, $E$ designates a compact set included in $L^2(\Omega_x, \mathbb{R}^d)$, the space of measurable functions that are square integrable on $\Omega_x$, and like [56], to enable differentiability of the RKHS function, we will assume that $E$ is the closure of its nonempty interior $\mathring{E}$. We assume also that $\nu$ is a finite measure on $E$, (i.e. $\nu(E) < \infty$). Denoting by $\partial_u^{(2)} k(\cdot, y)$ the kernel functional derivative in the direction $u$ with respect to the second variable $y$ – while the upper-script $\partial_u^{(1)}$ indicates a directional derivative in the first variable –, we have first the following lemma on the feature maps derivative.

**Lemma 2 (Feature maps differentiability)** The kernel derivatives belong to $\mathcal{H}$. More precisely, let $k \in C^{(1,1)}(E \times E)$, then for all $y \in E$ and $u \in E$, we have $\partial_u^{(2)} k(\cdot, y) \in \mathcal{H}$.

**PROOF** Let $y \in \mathring{E}$ and $u \in E$. Let $(h_n)_{n \ge 0}$ be a sequence of $\mathbb{C}\setminus\{0\}$ which converges to 0 and such that for all $n$, $y+h_nu \in \mathring{E}$ (since $\mathring{E}$ is an open set). For all integers $n$, let

$$ d_n(\cdot, y) := k(\cdot, y + h_n u) - k(\cdot, y) \text{ and } t_n(\cdot, y) := \frac{d_n(\cdot, y)}{h_n}. $$

We first notice that $(t_n(\cdot, y))_{n \ge 0}$ converges pointwise to $\partial_u^{(2)} k(\cdot, y)$ since $k$ is $C^{(1,1)}(E \times E)$. Then, we now show that $(t_n(\cdot, y))_n$ is a Cauchy sequence of $(\mathcal{H}, \|\cdot\|_{\mathcal{H}})$. For all integers $p$ and $q$, we have

$$ \|t_p(\cdot, y) - t_q(\cdot, y)\|_{\mathcal{H}}^2 = \|t_p(\cdot, y)\|_{\mathcal{H}}^2 + \|t_q(\cdot, y)\|_{\mathcal{H}}^2 - 2\mathrm{Re}[\langle t_p(\cdot, y), t_q(\cdot, y) \rangle_{\mathcal{H}}]. \quad (A.6) $$

By the reproducing property of $\mathcal{H}$ and as $k$ is $C^{(1,1)}(E \times E)$, we have for all $x \in E$ and for some $p$ and $q$ large enough,

$$ \begin{align*} \langle t_p(\cdot, y), t_q(\cdot, x) \rangle_{\mathcal{H}} &= \frac{1}{h_p h_q} \langle k(\cdot, y + h_p u) - k(\cdot, y), k(\cdot, x + h_q u) - k(\cdot, x) \rangle_{\mathcal{H}} \\ &= \frac{1}{h_p h_q} (k(x + h_q u, y + h_p u) - k(x, y + h_p u) - k(x + h_q u, y) + k(x, y)) \\ &= \frac{1}{h_p h_q} \left( \partial_u^{(1)} k(x, y + h_p u) \times h_q - \partial_u^{(1)} k(x, y) \times h_q + o(h_q) \right) \\ &= \partial_u^{(2)} \partial_u^{(1)} k(x, y) + o_{p,q}(1) \end{align*} $$

We have also for some $p$ and $q$ large enough (with $y = x$) and $\varepsilon > 0$,

$$ \left| \langle t_p(\cdot, y), t_q(\cdot, y) \rangle_{\mathcal{H}} - \partial_u^{(2)} \partial_u^{(1)} k(y, y) \right| \le \varepsilon, $$
and by (A.6) we obtain hence $\|t_p(\cdot, y) - t_q(\cdot, y)\|_{\mathcal{H}}^2 \le 4\epsilon$, which means $(t_n(\cdot, y))_n$ is a Cauchy sequence in $\mathcal{H}$. The sequence $(t_n(\cdot, y))_n$ converges to a function $t(\cdot, y)$ in $\mathcal{H}$ and by (A.1) the convergence is pointwise. Consequently, we obtain also $\partial_u^{(2)} k(\cdot, y) = t(\cdot, y) \in \mathcal{H}$ for $y \in \mathring{E}$.

Let $y \in E \setminus \mathring{E}$. There exists a sequence $(y^i)$ of $\mathring{E}$ converging to $y$. For all integers $i$ and $j$, we have,

$$ \| \partial_u^{(2)} k(\cdot, y^i) - \partial_u^{(2)} k(\cdot, y^j) \|_{\mathcal{H}} = \partial_u^{(1)} \partial_u^{(2)} k(y^i, y^j) + \partial_u^{(1)} \partial_u^{(2)} k(y^j, y^i) - 2\mathrm{Re}(\partial_u^{(1)} \partial_u^{(2)} k(y^i, y^j)). $$

As $(y^i)$ is a Cauchy sequence of $E$ and $\partial_u^{(1)} \partial_u^{(2)} k$ is continuous on $E \times E$, we infer that $(\partial_u^{(2)} k(\cdot, y^i))_i$ is a Cauchy sequence of $\mathcal{H}$. In particular, $(\partial_u^{(2)} k(\cdot, y^i))_i$ converge to a function $h$. For all $x \in E$, we have

$$ h(x) = \langle h, k(\cdot, x) \rangle_{\mathcal{H}} = \lim_{i \to \infty} \langle \partial_u^{(2)} k(\cdot, y^i), k(\cdot, x) \rangle_{\mathcal{H}} = \partial_u^{(2)} k(x, y^i). $$

We obtain thus that $h = \partial_u^{(2)} k(\cdot, y) \in \mathcal{H}$ for all $y \in E \setminus \mathring{E}$, and consequently, $\partial_u^{(2)} k(\cdot, y) \in \mathcal{H}$ for all $y \in E$. $\square$

As the sequence $(t_n(\cdot, x))_n$ converges to $\partial_u^{(2)} k(\cdot, x)$ in $\mathcal{H}$, we obtain also

$$ \langle f, \partial_u^{(2)} k(\cdot, x) \rangle_{\mathcal{H}} = \lim_{p \to \infty} \left\langle f, \frac{k(\cdot, x+h_p) - k(\cdot, x)}{h_p} \right\rangle_{\mathcal{H}} = \lim_{p \to \infty} \frac{f(x+h_p u) - f(x)}{h_p} = \partial_u f(x). \quad (A.7) $$

We provide in theorem 6 an expression of the directional derivative along direction $u \in E$ of an RKHS function $f$. To that end we first introduce useful remarks together with an intermediate lemma.

As $k$ is $C^{(1,1)}(E \times E)$, the Mercer eigenfunctions $\varphi_i$ are $C^1(E, \mathbb{C})$. Similarly to remark 11, we get a computational rule relating the kernel operator derivatives and the inner product $\langle \cdot, \cdot \rangle_{L^2_{\mathcal{H}}(E,v)}$. For all $y \in E$ and $u \in E$, we note that

$$ \langle \partial_u^{(2)} k(\cdot, y), \varphi_i \rangle_{L^2_{\mathcal{H}}(E,v)} = \int_E \partial_u^{(2)} k(x, y) \overline{\varphi_i(x)} v(dx) = \partial_u \int_E \overline{k(y, x)} \overline{\varphi_i(x)} v(dx) = \partial_u \overline{\mathcal{L}_k(\varphi_i)(y)} = \mu_i \partial_u \overline{\varphi_i(y)}. $$

We obtain for all $y \in E$ and $u \in E$ that

$$ \langle \partial_u^{(2)} k(\cdot, y), \varphi_i \rangle_{L^2_{\mathcal{H}}(E,v)} = \mu_i \partial_u \overline{\varphi_i(y)} \quad \text{and} \quad \langle \varphi_i, \partial_u^{(2)} k(\cdot, y) \rangle_{L^2_{\mathcal{H}}(E,v)} = \mu_i \partial_u \varphi_i(y). \quad (A.8) $$

These remarks allow us to prove the next lemma.

**Lemma 3** The series $\sum \mu_i |\partial_u \varphi_i(y)|^2$ converges uniformly on $E$.

**PROOF** The functions $\mu_i \varphi_i$ are $C^1(E, \mathbb{C})$. By definition of $\|\cdot\|_{\mathcal{H}}$ and the previous remark, we have

$$ \| \partial_u^{(2)} k(\cdot, y) \|_{\mathcal{H}}^2 = \sum_{i=0}^{\infty} \frac{|\langle \partial_u^{(2)} k(\cdot, y), \varphi_i \rangle_{L^2_{\mathcal{H}}(E,v)}|^2}{\mu_i} = \sum_{i=0}^{\infty} \mu_i |\partial_u \varphi_i(y)|^2. $$

The series $\sum \mu_i |\partial_u \varphi_i(y)|^2$ converges pointwise to the continuous function $\|\partial_u^{(2)} k(\cdot, y)\|_{\mathcal{H}}^2$. By Dini's theorem (the sequence of partial sums is monotonically increasing as $\mu_i$ is positive), the convergence is also uniform. $\square$

**Theorem 6 (Mercer's decomposition of derivatives)** Let $u \in E$. For all $f = \sum_{i=0}^{\infty} a_i \varphi_i \in \mathcal{H}$ (with $a_i = \langle f, \varphi_i \rangle_{L^2_{\mathcal{H}}(E,v)}$), the directional derivative of the function $f$ in the direction $u$ exists and we have

$$ \partial_u f = \sum_{i=0}^{\infty} a_i \partial_u \varphi_i, $$

where the convergence is uniform on $E$.
PROOF (THEOREM 6) Let $f = \sum_{i=0}^{\infty} a_i \varphi_i \in \mathcal{H}$. For all $i$, the function $a_i \varphi_i$ belongs to $C^1(E, \mathbb{C})$ and the series $\sum_i a_i \varphi_i$ converges pointwise to $f$ on $E$. As for the convergence of the derivative series, we have, by Cauchy-Schwarz inequality, that for all $n \in \mathbb{N}$ and $y \in E$,

$$ \left| \sum_{i=n}^{\infty} a_i \partial_u(\varphi_i)(y) \right| \le \left( \sum_{i=n}^{\infty} \frac{|a_i|^2}{\mu_i} \right)^{1/2} \left( \sum_{i=n}^{\infty} \mu_i |\partial_u(\varphi_i)(y)|^2 \right)^{1/2}. $$

By Lemma 3, the term $\sum_{i=n}^{\infty} \mu_i |\partial_u(\varphi_i)(y)|^2$ converges uniformly on $E$ to $0$. Since $f \in \mathcal{H}$, we obtain also that the series $\sum a_i \partial_u(\varphi_i)$ converges uniformly on $E$. Then, through term by term differentiation the result is proved. $\square$

We remark immediately that as functions of the RKHS the kernel feature maps $k(\cdot, x)$ can be differentiated to obtain a term by term Mercer's type decomposition of the kernel derivatives.

**Proposition 9 (Mercer Kernel derivatives decomposition)** For all $x, y \in E$ and for all $u, w \in E$, we have

$$ \partial_u^{(1)} \partial_w^{(2)} k(x, y) = \sum_{i=0}^{\infty} \mu_i \partial_u \varphi_i(x) \partial_w \overline{\varphi_i(y)} \quad (A.9) $$

where the convergence is uniform on $E \times E$.

PROOF The functions $(x, y) \mapsto \mu_i \varphi_i(x) \overline{\varphi_i(y)}$ are $C^{(1,1)}(E \times E)$. For all $u$ and $w \in E$, we apply theorem 6: the series $\sum \mu_i \partial_u \varphi_i(x) \overline{\varphi_i(y)}$ converges uniformly to $\partial_u^{(1)} k(x, y)$ and $\sum \mu_i \varphi_i(x) \partial_w \overline{\varphi_i(y)}$ converges uniformly to $\partial_w^{(2)} k(x, y)$ on $E \times E$. For all $w \in E$, theorem 6 can be applied to the function $x \mapsto \partial_w^{(2)} k(x, y) \in \mathcal{H}$ and for all $u \in E$ and all $x, y \in E$,

$$ \partial_u^{(1)} \partial_w^{(2)} k(x, y) = \sum_{i=0}^{\infty} a_{i,w}(y) \partial_u \varphi_i(x), $$

with $a_{i,w}(y) := \langle \partial_w^{(2)} k(\cdot, y), \varphi_i \rangle_{L^2(E,v)}$. By (A.8), we obtain that $a_{i,w}(y) = \mu_i \partial_w \overline{\varphi_i(y)}$ and

$$ \partial_u^{(1)} \partial_w^{(2)} k(x, y) = \sum_{i=0}^{\infty} \mu_i \partial_u \varphi_i(x) \partial_w \overline{\varphi_i(y)}. $$

Finally, by Cauchy-Schwarz, we have that for all $N \in \mathbb{N}$

$$ \left| \sum_{i=0}^{N} \mu_i \partial_u \varphi_i(x) \partial_w \overline{\varphi_i(y)} \right|^2 \le \left( \sum_{i=0}^{N} \mu_i |\partial_u \varphi_i(x)|^2 \right) \left( \sum_{i=0}^{N} \mu_i |\partial_w \varphi_i(y)|^2 \right), $$

and by Lemma 3 the convergence of $\sum \mu_i \partial_u \varphi_i(x) \partial_w \overline{\varphi_i(y)}$ to $\partial_u^{(1)} \partial_w^{(2)} k(x, y)$ is uniform on $E \times E$. $\square$

These results can be easily extended to higher order derivatives provided a kernel with enough regularity. They correspond, in slightly different forms, to results given in [87, 88] for $\mathbb{R}^d$. Note that the existence of such higher order derivatives pertains to differentiability conditions of the Mercer eigenfunctions, $(\varphi_i)_{i \in \mathbb{N}}$, associated with kernel $k$.

By Theorem 6, the functions of $\mathcal{H}$ admit a directional derivative along any vector of $E$. The $L^2_{\mathbb{C}}(E, v)$ norm of the directional derivatives of functions of $\mathcal{H}$ can then be controlled by the RKHS norm.

**Theorem 7 (Continuity in $L^2_{\mathbb{C}}(E, v)$ of RKHS derivatives)** There exists a constant $C > 0$ for which:

$$ \| \partial_u f \|_{L^2_{\mathbb{C}}(E, v)} \le C \| f \|_{\mathcal{H}} \quad \text{for all } f \in \mathcal{H} \text{ and } u \in E. $$
PROOF Let $u \in E$. We prove the existence of $C > 0$ such that $\|\partial_u f\|_{L^2_{\mathbb{C}}(E, v)} \le C \|f\|_{\mathcal{H}}$ for all $f = \sum_{i=0}^{\infty} a_i \varphi_i \in \mathcal{H}$.

By theorem 6 and Fubini's theorem, we have

$$ \partial_u f = \sum_{i=0}^{\infty} a_i \partial_u \varphi_i = \sum_{i=0}^{\infty} a_i \sum_{\ell=0}^{\infty} \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \varphi_\ell = \sum_{\ell=0}^{\infty} \left( \sum_{i=0}^{\infty} a_i \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \right) \varphi_\ell. $$

With Cauchy-Schwarz inequality we infer that

$$ \sum_{\ell=0}^{\infty} \left| \sum_{i=0}^{\infty} a_i \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \right|^2 \le \sum_{\ell=0}^{\infty} \left( \sum_{i=0}^{\infty} \frac{|a_i|^2}{\mu_i} \right) \left( \sum_{i=0}^{\infty} \mu_i \left| \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \right|^2 \right) \le \|f\|^2_{\mathcal{H}} \left( \sum_{\ell=0}^{\infty} \sum_{i=0}^{\infty} \mu_i \left| \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \right|^2 \right). $$

The last sum is finite as we have, using Fubini's theorem,

$$ \sum_{\ell=0}^{\infty} \sum_{i=0}^{\infty} \mu_i \left| \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \right|^2 = \sum_{i=0}^{\infty} \mu_i \sum_{\ell=0}^{\infty} \left| \langle \partial_u(\varphi_i), \varphi_\ell \rangle_{L^2_{\mathbb{C}}(E, v)} \right|^2 \le \sum_{i=0}^{\infty} \mu_i \int_E |\partial_u(\varphi_i)(x)|^2 v(dx). $$

Upon applying the monotone convergence theorem (since $\mu_i > 0$) and proposition 9 we obtain

$$ \sum_{i=0}^{\infty} \mu_i \int_E |\partial_u(\varphi_i)(x)|^2 v(dx) = \int_E \sum_{i=0}^{\infty} \mu_i |\partial_u(\varphi_i)(x)|^2 v(dx) = \int_E \partial_u^{(2)} \partial_u^{(1)} k(x,x) v(dx), $$

and $(\int_E \partial_u^{(2)} \partial_u^{(1)} k(x,x) v(dx))^{1/2} \le (v(E) \sup_{x \in E} |k(x,x)| \sup_{u \in E} |u|^2)^{1/2} < \infty$, where continuity of the kernel has been used on both arguments and noting that $v$ is a finite measure on $E$. We denote by $C$ the last constant. $\square$

Appendix A.3. Isometry between $L^2_{\mathbb{C}}(E, v)$ and $\mathcal{H}$

Let us now focus on the links between $L^2_{\mathbb{C}}(E, v)$ and $\mathcal{H}$. Let us recall that integral kernel operator $\mathcal{L}_k$ defined eq.(1) is positive definite, compact and self-adjoint, with a range that is assumed to be dense in $L^2_{\mathbb{C}}(E, v)$. As a consequence, it admits a summable non-increasing sequence of strictly positive eigenvalues $(\mu_i)_{i \in \mathbb{N}}$ with corresponding eigenfunctions, $(\varphi_i)_{i \in \mathbb{N}}$ forming an orthonormal basis of $L^2_{\mathbb{C}}(E, v)$. Since $\mu_i > 0$ for all integers $i$, the integral operator $\mathcal{L}_k$ admits a unique symmetric square root operator denoted by $\mathcal{L}_k^{1/2}$ and defined for all integers $i$ as $\mathcal{L}_k^{1/2}(\varphi_i) = \mu_i^{1/2} \varphi_i$. The operator $\mathcal{L}_k^{1/2} : L^2_{\mathbb{C}}(E, v) \to \mathcal{H}$ is a bijective isometry: for all $f, g \in L^2_{\mathbb{C}}(E, v)$

$$ \langle \mathcal{L}_k^{1/2} f, \mathcal{L}_k^{1/2} g \rangle_{\mathcal{H}} = \langle f, g \rangle_{L^2_{\mathbb{C}}(\Omega, v)}. \qquad (A.10) $$

This operator must be understood as $\mathcal{L}_k^{1/2} \circ \mathcal{L}_k^{1/2} = \mathcal{L}_k$. The inverse of $\mathcal{L}_k^{1/2}$, noted $\mathcal{L}_k^{-1/2} : \mathcal{H} \mapsto L^2_{\mathbb{C}}(E, v)$, and defined for all integers $i \ge 0$ by $\mathcal{L}_k^{-1/2}(\varphi_i) = \mu_i^{-1/2} \varphi_i$ is an isometry as well: for all $f, g \in \mathcal{H}$

$$ \langle \mathcal{L}_k^{-1/2} f, \mathcal{L}_k^{-1/2} g \rangle_{L^2_{\mathbb{C}}(E, v)} = \langle f, g \rangle_{\mathcal{H}}. \qquad (A.11) $$

It can be checked easily that $\mathcal{L}_k^{1/2}$ and $\mathcal{L}_k^{-1/2}$ are Hermitian for the inner product of $L^2_{\mathbb{C}}(E, v)$: for all $f, g \in L^2_{\mathbb{C}}(E, v)$

$$ \langle j(\mathcal{L}_k^{1/2}(f)), g \rangle_{L^2_{\mathbb{C}}(E, v)} = \langle f, j(\mathcal{L}_k^{1/2}(g)) \rangle_{L^2_{\mathbb{C}}(E, v)}, \qquad (A.12) $$

with $j : \mathcal{H} \to L^2_{\mathbb{C}}(E, v)$ the continuous injection (by comparison of the norms) and for all $f, g \in \mathcal{H}$

$$ \langle \mathcal{L}_k^{-1/2}(f), g \rangle_{L^2_{\mathbb{C}}(E, v)} = \langle f, \mathcal{L}_k^{1/2}(g) \rangle_{L^2_{\mathbb{C}}(E, v)}. \qquad (A.13) $$

By construction, the family of Mercer eigenfunctions $(\varphi_i)_{i \in \mathbb{N}}$ is an orthonormal basis of $L^2_{\mathbb{C}}(E, v)$ and $\mathcal{L}_k^{1/2}$ enables to get an orthonormal basis of $\mathcal{H}$. We prove that in the next proposition.
**Proposition 10 (Orthonormal basis of $\mathcal{H}$)** The sequence of rescaled kernel eigenfunctions $\left(\mathcal{L}_k^{1/2}(\varphi_i)\right)_{i \ge 0}$ is an orthonormal basis of $\mathcal{H}$.

**PROOF** By (A.10), the family $\left(\mathcal{L}_k^{1/2}(\varphi_i)\right)_i$ is orthonormal in $\mathcal{H}$. Furthermore, we have for all $f \in \mathcal{H}$

$$ \sum_{i=0}^{\infty} \left| \left\langle f, \mathcal{L}_k^{1/2}(\varphi_i) \right\rangle_{\mathcal{H}} \right|^2 = \sum_{i=0}^{\infty} \left| \left\langle \mathcal{L}_k^{-1/2} f, \varphi_i \right\rangle_{L_{\mathcal{H}}^2(E, \nu)} \right|^2 = \left\| \mathcal{L}_k^{-1/2}(f) \right\|_{L_{\mathcal{H}}^2(E, \nu)}^2 = \|f\|_{\mathcal{H}}^2, $$

by (A.13) and (A.11). We obtain the result through Bessel characterization.

<figure><img src="image_10.png" /></figure>

The next theorem and proposition provide us useful results for the following

**Theorem 8 (Compact injection)** The injection $j: \mathcal{H} \to L_{\mathcal{H}}^2(E, \nu)$ is continuous and compact.

**PROOF** The comparison of the norms of $\mathcal{H}$ and $L_{\mathcal{H}}^2(E, \nu)$ and (A.3) justifies the continuity of $j$ (see remark 14). By proposition 10 and (A.3), we have

$$ \sum_{i=0}^{\infty} \left\| j(\mathcal{L}_k^{1/2} \varphi_i) \right\|_{L_{\mathcal{H}}^2(E, \nu)}^2 = \sum_{i=0}^{\infty} \mu_i < \infty; $$

$j$ is also an Hilbert-Schmidt operator and therefore a compact operator.

<figure><img src="image_11.png" /></figure>

**Proposition 11 (Compact injection range denseness)** The range of $j$, $j(\mathcal{H})$, is dense in $L_{\mathcal{H}}^2(E, \nu)$, more precisely the linear set spanned by $\{k(\cdot, X) : X \in E\}$ is dense in $L_{\mathcal{H}}^2(E, \nu)$.

**PROOF** Let $f \in L_{\mathcal{H}}^2(E, \nu)$ such that $\langle f, g \rangle_{L_{\mathcal{H}}^2(E, \nu)} = 0$ for all $g \in \mathcal{H}$. In particular, for all $X \in E$ we have $\langle f, k(\cdot, X) \rangle_{L_{\mathcal{H}}^2(E, \nu)} = 0$, which means that $(\mathcal{L}_k f)(X) = 0$. Since $\mathcal{L}_k$ is injective, we obtain that $f = 0$. □

From now on, we note $\mathcal{H}^{\text{sp}} := \text{Span}\{k(\cdot, X) : X \in E\}$. Proposition (11) justifies that $\mathcal{H}^{\text{sp}}$ is dense in $(L_{\mathcal{H}}^2(E, \nu), \|\cdot\|_{L_{\mathcal{H}}^2(E, \nu)})$. Furthermore, by the Moore-Aronszajn theorem [59], $\mathcal{H}^{\text{sp}}$ is also dense in $(\mathcal{H}, \|\cdot\|_{\mathcal{H}})$.

**Remark 16 (Relaxation of $\mathcal{L}_k$ injectivity)** The injectivity property of $\mathcal{L}_k$ (equivalent to its positive definite character) is crucial to show the denseness of $j(\mathcal{H})$. If some eigenvalues of $\mathcal{L}_k$ are zero, the embedding space $L_{\mathcal{H}}^2(E, \nu)$ we are working with should be replaced by the linear subspace, $\mathcal{V} = \text{Ker}(\mathcal{L}_k)^\perp \subset L_{\mathcal{H}}^2(E, \nu)$, which is the actual closure of $j(\mathcal{H})$, and satisfies $\mathcal{V} = \text{Span}\{\varphi_i : i \ge 0 \text{ with } \mu_i \ne 0\}$.

**Remark 17 (Higher-order RKHS differentiability)** Let $p \ge 1$ be an integer. Assuming $k \in C^{(p,p)}(E \times E)$, the result of Lemma 2 and theorem 6 can be extended to higher order differentiation. As a matter fact, for all $u \in E$, all $x \in E$ and all integer $0 \le m \le p$, we can prove that the $m$th order derivative $\partial_u^{(2),m} k(\cdot, x) \in \mathcal{H}$ and $\mathcal{H} \subset C^p(E, \mathbb{C})$. This corresponds, in slightly different forms, to results provided in [87, 88] for $\mathbb{R}^d$.

**Remark 18 (Links between basis of $L_{\mathcal{H}}^2(E, \nu)$ and $\mathcal{H}$)** Proposition 10 proves that if $(e_i)_i$ is an orthonormal basis of $L_{\mathcal{H}}^2(E, \nu)$ then $(\mathcal{L}_k^{1/2} e_i)_i$ is an orthonormal basis of $\mathcal{H}$. The converse statement is also true and can be shown with the same arguments: if $(f_i)_i$ is an orthonormal basis of $\mathcal{H}$ then $(\mathcal{L}_k^{-1/2} f_i)_i$ is a Hilbert basis of $L_{\mathcal{H}}^2(E, \nu)$.

## Appendix B. Properties of the Extension and Restriction operators $E_t$ and $R_t$

In this appendix we elaborate on the properties of the restriction and extension operators (namely that they are isometries; form an adjoint pairs in $\mathcal{H}_t$; provide a representation of $R_t$ in $\mathcal{H}$ and state the continuity of $R_t$ in $L^2(\Omega_t, \nu)$).

**Proposition 12 (Extension/Restriction isometry)** For all $t \ge 0$, the mapping $E_t: \mathcal{H}_t \to \mathcal{H}$ is an isometry:

$$ \forall f, g \in \mathcal{H}_t \qquad \langle E_t(f), E_t(g) \rangle_{\mathcal{H}_t} = \langle f, g \rangle_{\mathcal{H}_t} \qquad (B.1) $$

and $R_t: \mathcal{H} \to \mathcal{H}_t$ is an isometry as well:

$$ \forall f, g \in \mathcal{H} \qquad \langle R_t(f), R_t(g) \rangle_{\mathcal{H}_t} = \langle f, g \rangle_{\mathcal{H}}. \qquad (B.2) $$
**PROOF (FOR THE EXTENSION)** Let $f$ and $g \in \operatorname{Span}\{k_t(\cdot, X_t) : X_t \in \Omega_t\} : f = \sum_{i=1}^{N_1} \alpha_i k_t(\cdot, X_t^{(i)})$ and $g = \sum_{j=1}^{N_2} \beta_j k_t(\cdot, Y_t^{(j)})$

with $X_t^{(i)}$ and $Y_t^{(j)} \in \Omega_t$. By definition, we have, for $X_t = (X_t, t) \in \Omega$ and $Y_t = (Y_t, t) \in \Omega$

$$ \langle E_t(f), E_t(g) \rangle_{\mathcal{H}} = \sum_{i=1}^{N_1} \sum_{j=1}^{N_2} \alpha_i \bar{\beta}_j \langle k(\cdot, X_t^{(i)}), k(\cdot, Y_t^{(j)}) \rangle_{\mathcal{H}} = \sum_{i=1}^{N_1} \sum_{j=1}^{N_2} \alpha_i \bar{\beta}_j k(Y_t^{(j)}, X_t^{(i)}). $$

By construction of $k$, we have $k(Y_t^{(j)}, X_t^{(i)}) = k_t(Y_t^{(j)}, X_t^{(i)})$ for all $i$ and $j$, hence

$$ \langle E_t(f), E_t(g) \rangle_{\mathcal{H}} = \sum_{i=1}^{N_1} \sum_{j=1}^{N_2} \alpha_i \bar{\beta}_j k_t(Y_t^{(j)}, X_t^{(i)}) = \sum_{i=1}^{N_1} \sum_{j=1}^{N_2} \alpha_i \bar{\beta}_j \langle k_t(\cdot, X_t^{(i)}), k_t(\cdot, Y_t^{(j)}) \rangle_{\mathcal{H}} = \langle f, g \rangle_{\mathcal{H}_t}. $$

By density of $\operatorname{Span}\{k_t(\cdot, X_t) : X_t \in \Omega_t\}$ in $(\mathcal{H}_t, \|\cdot\|_{\mathcal{H}_t})$, this result extends to all $f$ and $g \in \mathcal{H}_t$ and (B.1) is proved. $\square$

With similar arguments for $R_t$, (B.2) is also proved. By (11) and (12), the mappings $E_t$ and $R_t$ have inverse action on the mapping $k$ and $k_t$. The next proposition justifies that the adjoint of $E_t$ is $R_t$.

**Proposition 13** For all $t \ge 0$, the mappings $E_t$ and $R_t$ are such that for all $f \in \mathcal{H}_t$ and $g \in \mathcal{H}$

$$ \langle E_t(f), g \rangle_{\mathcal{H}} = \langle f, R_t(g) \rangle_{\mathcal{H}_t}. \qquad (B.3) $$

**PROOF** Let us first show the property for feature maps: for $X_t = \Phi_t(X_0) \in \Omega_t$ and $Y = \Phi_s(Y_0) \in \Omega$, we have

$$ \langle E_t[k_t(\cdot, X_t)], k(\cdot, Y) \rangle_{\mathcal{H}} = k(Y, X_t) = k_0(Y_0, X_0)\ell(s, t) $$

and

$$ \langle k_t(\cdot, X_t), R_t[k(\cdot, Y)] \rangle_{\mathcal{H}_t} = \langle k_t(\cdot, X_t), k_t(\cdot, \Phi_t(Y_0))\ell(t, s) \rangle_{\mathcal{H}_t} = k_t(\Phi_t(Y_0), X_t)\ell(t, s) = k_0(Y_0, X_0)\ell(s, t), $$

where the first equality comes from the definition of the restriction in equation (12), and the last equality follows from Definition 1 and the properties of the time kernel in Definition 2. By linearity we obtain that $\langle E_t(f), g \rangle_{\mathcal{H}} = \langle f, R_t(g) \rangle_{\mathcal{H}_t}$ for all $f \in \operatorname{Span}\{k_t(\cdot, X_t) : X_t \in \Omega_t\}$ and $g \in \operatorname{Span}\{k(\cdot, X) : X \in \Omega\}$.

Let $f \in \mathcal{H}_t$ and $g \in \mathcal{H}$. There exists sequences $(f_n)_n$ of $\operatorname{Span}\{k_t(\cdot, X_t) : X_t \in \Omega_t\}$ and $(g_n)_n$ of $\operatorname{Span}\{k(\cdot, X) : X \in \Omega\}$ such that $\|f - f_n\|_{\mathcal{H}_t} \to 0$ and $\|g - g_n\|_{\mathcal{H}} \to 0$. We have for all $n$ the equality

$$ \langle E_t(f_n), g_n \rangle_{\mathcal{H}} = \langle f_n, R_t(g_n) \rangle_{\mathcal{H}_t}. $$

Furthermore, by Cauchy-Schwarz and (B.1) we have that

$$ \begin{aligned} |\langle E_t(f_n), g_n \rangle_{\mathcal{H}} - \langle E_t(f), g \rangle_{\mathcal{H}}| & \le |\langle E_t(f_n) - E_t(f), g \rangle_{\mathcal{H}}| + |\langle E_t(f), g_n - g \rangle_{\mathcal{H}}| \\ & \le \|f_n - f\|_{\mathcal{H}_t} \|g\|_{\mathcal{H}} + \|f\|_{\mathcal{H}} \|g_n - g\|_{\mathcal{H}} \end{aligned} $$

and $\langle E_t(f_n), g_n \rangle_{\mathcal{H}} \to \langle E_t(f), g \rangle_{\mathcal{H}}$. With similar arguments, it can be readily proved that $\langle f_n, R_t(g_n) \rangle_{\mathcal{H}} \to \langle f, R_t(g) \rangle_{\mathcal{H}_t}$ by using (B.2). We obtain hence (B.3). $\square$

**Remark 19 (Expression of $R_t(f)$ for $f \in \mathcal{H}$)** The function $R_t(f)$ belongs to $\mathcal{H}_t$. By proposition 13, we have for all $X_t \in \Omega_t$

$$ R_t(f)(X_t) = \langle R_t(f), k_t(\cdot, X_t) \rangle_{\mathcal{H}} = \langle f, k(\cdot, X_t) \rangle_{\mathcal{H}} = f(X_t). $$

More generally, for all $t \ge 0$ the mapping $R_t$ can be extended on $L_{\mathcal{H}}^2(\Omega, \nu)$ (keeping the same notation) as:

$$ R_t: \quad L_{\mathcal{H}}^2(\Omega, \nu) \quad \begin{array}{rcl} \mapsto & L_{\mathcal{H}}^2(\Omega_t, \nu) & \\ f \mapsto & f \upharpoonright_{\Omega_t} & . \end{array} \qquad (B.4) $$

In particular, on $\mathcal{H} \subset L_{\mathcal{H}}^2(\Omega, \nu)$ this boils down to the definition given in (12). As $\Omega_t \subset \Omega$, we have directly the next proposition.

**Proposition 14 (Continuity of the restriction)** The mapping $R_t : (L_{\mathcal{H}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathcal{H}}^2(\Omega, \nu)}) \to (L_{\mathcal{H}}^2(\Omega_t, \nu), \|\cdot\|_{L_{\mathcal{H}}^2(\Omega_t, \nu)})$ is continuous.
## Appendix C. Proof of Lemma 1 (Perron-Frobenius infinitesimal generator)

In this appendix we provide the proof on the Perron-Frobenius infinitesimal generator expression and of its domain, $\mathcal{D}(A_p)$. We recall that $\mathcal{H}^{\text{sp}}$ is the linear span of the set $\{k(\cdot, x) : x \in \Omega\}$.

**PROOF** Let $f \in \mathcal{D}(A_p)$ and $\psi \in \mathcal{H}^{\text{sp}}$, we have for all $t \ge 0$,

$$ \left\langle \frac{P_t \circ f - f}{t}, \psi \right\rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} = \left\langle f, \frac{U_t \circ \psi - \psi}{t} \right\rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} \quad (C.1) $$

For all $x \in \Omega$, we have

$$ \frac{U_t(\psi)(x) - \psi(x)}{t} \underset{t \to 0^+}{\longrightarrow} \frac{d}{dt} [U_t(\psi)(x)] \Big|_{t=0} = \partial_{\mathcal{H}(x)} \psi(x). $$

We pass to the limit $[t \to 0^+]$ in (C.1), and recalling that $\mathcal{M}$ is real valued, divergence free with $\mathcal{M} \upharpoonright_{\partial\Omega} = 0$, and measure invariant we obtain

$$ \int_{\Omega} A_p f(x) \overline{\psi(x)} \nu(dx) = \int_{\Omega} -\partial_{\mathcal{H}(x)} f(x) \overline{\psi(x)} \nu(dx) $$

thus, $x \mapsto \partial_{\mathcal{H}(x)} f(x)$ belongs to $L_{\mathcal{C}}^2(\Omega, \nu)$ and

$$ \mathcal{D}(A_p) \subset \{f \in L_{\mathcal{C}}^2(\Omega, \nu) : x \mapsto \partial_{\mathcal{H}(x)} f(x) \in L_{\mathcal{C}}^2(\Omega, \nu)\} $$

$$ \text{and } A_p f(x) = -\partial_{\mathcal{H}(x)} f(x) \text{ for } f \in \mathcal{D}(A_p) \text{ and } x \in \Omega. $$

Conversely, considering now $f \in L_{\mathcal{C}}^2(\Omega, \nu)$ such that $x \mapsto \partial_{\mathcal{H}(x)} f(x) \in L_{\mathcal{C}}^2(\Omega, \nu)$, we show that $f \in \mathcal{D}(A_p)$. For all $\psi \in \mathcal{H}^{\text{sp}}$, we have

$$ \langle P_t(f) - f, \psi \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} = \langle f, U_t(\psi) - \psi \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} = \\ \langle f, \int_0^t \frac{dU_s(\psi)}{ds} ds \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} = \langle f, \int_0^t \partial_{\mathcal{H}(\phi_s(\cdot))} \psi[\Phi_s(\cdot)] ds \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)}. $$

For all $\psi \in \mathcal{H}^{\text{sp}}$, we have $\langle -\partial_{\mathcal{H}(\cdot)} f(\cdot), \psi \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} = \langle f, -\partial_{\mathcal{H}(\cdot)} \psi(\cdot) \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)}$ (since $\mathcal{H}$ is divergence free) and we obtain also the following equality

$$ \left\langle \frac{P_t(f) - f}{t} - \partial_{\mathcal{H}(\cdot)} f(\cdot), \psi \right\rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} = \left\langle f, \frac{1}{t} \int_0^t \partial_{\mathcal{H}(\phi_s(\cdot))} \psi[\Phi_s(\cdot)] ds - \partial_{\mathcal{H}(\cdot)} \psi(\cdot) \right\rangle_{L_{\mathcal{C}}^2(\Omega, \nu)}. $$

Since $t \mapsto \int_0^t \partial_{\mathcal{H}(\phi_s(\cdot))} \psi[\Phi_s(\cdot)] ds$ is differentiable in $0^+$, we obtain in particular

$$ \left\langle \frac{P_t(f) - f}{t} - \partial_{\mathcal{H}(\cdot)} f(\cdot), \psi \right\rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} \underset{t \to 0^+}{\longrightarrow} 0, $$

for all function $\psi \in \mathcal{H}^{\text{sp}}$. By density of $\mathcal{H}^{\text{sp}}$ in $L_{\mathcal{C}}^2(\Omega, \nu)$ we get

$$ \left\| \frac{P_t f - f}{t} - \partial_{\mathcal{H}(\cdot)} f(\cdot) \right\|_{L_{\mathcal{C}}^2(\Omega, \nu)} \underset{t \to 0^+}{\longrightarrow} 0 $$

and hence $f \in \mathcal{D}(A_p)$.

□
## Appendix D. Proof of step1 and step2 of the RKHS family spectral theorem

This appendix gathers step 1 and step 2 of the proof of the RKHS family representation theorem. Step 1 Uses Lions-Lax-Milgram theorem to show that an approximation $\tilde{A}_U$ of the restriction on $\mathcal{H}$ of the infinitesimal generator $A_U$ is bijective. The inverse operator $\tilde{A}_U^{-1}$ seen as an operator of $L_{\mathcal{C}}^2(\Omega, \nu) \mapsto L_{\mathcal{C}}^2(\Omega, \nu)$ through the composition with the compact injection $j: \mathcal{H} \mapsto L_{\mathcal{C}}^2(\Omega, \nu)$ (Theorem 8, Appendix A.3) is shown to be compact and self-adjoint (Lemmas 5 and 6). The diagonalization of $A_U$ is then obtained from the diagonalization of $\tilde{A}_U^{-1}$ and $\mathcal{L}_k$ (Prop.15). Step 2 finally constructs from the eigenfunctions of $A_U$ an orthonormal basis of $\mathcal{H}_t$ that diagonalizes the operator $A_{U,t}$ (Prop.16).

### Step 1

Let $\tilde{A}_U$ be the linear mapping defined by

$$ \tilde{A}_U : (\mathcal{H}, \|\cdot\|_{\mathcal{H}}) \mapsto (L_{\mathcal{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathcal{C}}^2(\Omega, \nu)}) \\ f \mapsto [iA_U + \lambda \mathcal{L}_k^{-1/2}] (f) \qquad (D.1) $$

with $\lambda$ a fixed real such that

$$ |\lambda| > \|A_U\| \left( \int_{\Omega} |k(x,x)|^2 \nu(dx) \right)^{1/2}. $$

By Theorem 7 (Appendix A.2), the norm of $A_U$ verifies

$$ \|A_U\| \le \left( \nu(\Omega) \sup_{x \in \Omega} \|k(x,x)\| \sup_{x \in \Omega} |\mathcal{M}(x)|^2 \right)^{1/2}, $$

and to lighten the notations in the following, we note $M := \|A_U\|$. We start analyzing the properties of $\tilde{A}_U$. First of all, the mapping $\tilde{A}_U$ is continuous. We have for all $f \in \mathcal{H}$ the inequality

$$ \|\tilde{A}_U f\|_{L_{\mathcal{C}}^2(\Omega, \nu)} \le \|iA_U f\|_{L_{\mathcal{C}}^2(\Omega, \nu)} + |\lambda| \|\mathcal{L}_k^{-1/2} f\|_{L_{\mathcal{C}}^2(\Omega, \nu)} = \|A_U f\|_{L_{\mathcal{C}}^2(\Omega, \nu)} + |\lambda| \|\mathcal{L}_k^{-1/2} f\|_{L_{\mathcal{C}}^2(\Omega, \nu)}. $$

The identity of (A.11) and Theorem 2 allow us to write

$$ \|\tilde{A}_U f\|_{L_{\mathcal{C}}^2(\Omega, \nu)} \le (M + |\lambda|) \|f\|_{\mathcal{H}}, $$

which shows $\tilde{A}_U$ is continuous on $\mathcal{H}$. As shown below, the constant $\lambda$ is chosen to ensure the inverse of $\tilde{A}_U$ is well defined.

**Lemma 4 (Bijectivity of $\tilde{A}_U$)** The operator $\tilde{A}_U : (\mathcal{H}, \|\cdot\|_{\mathcal{H}}) \mapsto (L_{\mathcal{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathcal{C}}^2(\Omega, \nu)})$ is bijective.

**PROOF** Let $a$ be the bilinear form defined by

$$ a: \mathcal{H} \times \mathcal{H} \mapsto \mathbb{C} \\ (f, g) \mapsto \langle \tilde{A}_U(f), \mathcal{L}_k^{-1/2}(g) \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)}. $$

The bilinear form $a$ is continuous on $\mathcal{H}$. By (A.11), for all $f \in \mathcal{H}$ we have

$$ a(f, f) = \langle iA_U(f), \mathcal{L}_k^{-1/2}f \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)} + \lambda \|f\|_{\mathcal{H}}^2. $$

In particular, we get

$$ |a(f, f)| \ge |\lambda| \|f\|_{\mathcal{H}}^2 - |\langle A_U(f), \mathcal{L}_k^{-1/2} f \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)}| \ge |\lambda| \|f\|_{\mathcal{H}}^2 - |\langle A_U(f), \mathcal{L}_k^{-1/2} f \rangle_{L_{\mathcal{C}}^2(\Omega, \nu)}|. $$
Using Cauchy-Schwarz inequality we have

$$-|\langle A_U(f), \mathcal{L}_k^{-1/2} f \rangle_{L^2(\Omega, v)}| \geq -\|A_U(f)\|_{L^2(\Omega, v)} \|\mathcal{L}_k^{-1/2} f\|_{L^2(\Omega, v)}.$$

Through the continuity of $A_U : \mathcal{H} \mapsto L^2_{\mathcal{H}}(\Omega, v)$, remark 14 and (A.11) we infer that

$$-|\langle A_U(f), \mathcal{L}_k^{-1/2} f \rangle_{L^2(\Omega, v)}| \geq -M \left( \int_{\Omega} |k(x,x)|^2 v(dx) \right)^{1/2} \|f\|_{\mathcal{H}}^2.$$

We obtain hence that $|a(f, f)| \geq C \|f\|_{\mathcal{H}}^2$ for all $f \in \mathcal{H}$ with the constant

$$C := \left( |\lambda| - M \left( \int_{\Omega} |k(x,x)|^2 v(dx) \right)^{1/2} \right) > 0,$$

which is positive by definition of $\lambda$. The bilinear form $a$ verifies thus a coercivity condition in $\mathcal{H}$.
To show that $\tilde{A}_U : \mathcal{H} \mapsto L^2_{\mathcal{H}}(\Omega, v)$ is bijective, let us consider $f \in L^2_{\mathcal{H}}(\Omega, v)$.

The linear mapping $g \mapsto \langle f, \mathcal{L}_k^{-1/2}(g) \rangle_{L^2_{\mathcal{H}}(\Omega, v)}$ is continuous on $\mathcal{H}$. By the Lions-Lax-Milgram theorem, there exists thus a unique $x_f \in \mathcal{H}$ such that

$$a(x_f, g) = \langle f, \mathcal{L}_k^{-1/2}(g) \rangle_{L^2_{\mathcal{H}}(\Omega, v)} \quad \text{for all } g \in \mathcal{H}.$$

As the operator $\mathcal{L}_k^{-1/2} : \mathcal{H} \mapsto L^2_{\mathcal{H}}(\Omega, v)$ is bijective, we obtain also

$$\langle \tilde{A}_U(x_f), h \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = \langle f, h \rangle_{L^2_{\mathcal{H}}(\Omega, v)} \quad \text{for all } h \in L^2_{\mathcal{H}}(\Omega, v). \qquad (D.2)$$

In particular, we obtain hence $\tilde{A}_U(x_f) = f$.

□

This lemma justifies the existence of the operator $\tilde{A}_U^{-1} : (L^2_{\mathcal{H}}(\Omega, v), \|\cdot\|_{L^2_{\mathcal{H}}(\Omega, v)}) \mapsto (\mathcal{H}, \|\cdot\|_{\mathcal{H}})$. By the bounded inverse theorem $\tilde{A}_U^{-1}$ is also bounded. In order to diagonalize $A_U$, we will diagonalize the inverse of $\tilde{A}_U$ in $L^2_{\mathcal{H}}(\Omega, v)$. To that end, we consider $j \circ \tilde{A}_U^{-1} : L^2_{\mathcal{H}}(\Omega, v) \mapsto L^2_{\mathcal{H}}(\Omega, v)$ with $j$ the injection of $\mathcal{H}$ in $L^2_{\mathcal{H}}(\Omega, v)$. As the injection is compact (Theorem 8), Appendix A.3), we obtain directly that the mapping

$$j \circ \tilde{A}_U^{-1} : (L^2_{\mathcal{H}}(\Omega, v), \|\cdot\|_{L^2_{\mathcal{H}}(\Omega, v)}) \mapsto (L^2_{\mathcal{H}}(\Omega, v), \|\cdot\|_{L^2_{\mathcal{H}}(\Omega, v)})$$

is continuous and compact. By abuse of notation denoting the composition $j \circ \tilde{A}_U^{-1}$ as $\tilde{A}_U^{-1} : L^2_{\mathcal{H}}(\Omega, v) \mapsto L^2_{\mathcal{H}}(\Omega, v)$ we have the next lemma.

**Lemma 5 (Compactness of $\tilde{A}_U^{-1}$)** The operator $\tilde{A}_U^{-1} : L^2_{\mathcal{H}}(\Omega, v) \mapsto L^2_{\mathcal{H}}(\Omega, v)$ is compact.

**Lemma 6 (Self-adjointness $\tilde{A}_U^{-1}$)** The operator $\tilde{A}_U^{-1}$ is self-adjoint in $L^2_{\mathcal{H}}(\Omega, v)$.

**PROOF** For all $f$ and $g \in L^2_{\mathcal{H}}(\Omega, v)$, we have

$$\langle \tilde{A}_U^{-1}(f), g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = \langle j \circ x_f, g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = \langle j \circ x_f, \tilde{A}_U x_g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} \qquad (D.3)$$

where we used the notation of the proof of Lemma 4 (and the elements $x_f$ and $x_g$ belong to $\mathcal{H}$ with $j$ the injection $\mathcal{H} \mapsto L^2_{\mathcal{H}}(\Omega, v)$). By definition, we have

$$\langle j \circ x_f, \tilde{A}_U x_g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = -i \langle j \circ x_f, A_U(x_g) \rangle_{L^2_{\mathcal{H}}(\Omega, v)} + \lambda \langle j \circ x_f, \mathcal{L}_k^{-1/2}(x_g) \rangle_{L^2_{\mathcal{H}}(\Omega, v)}.$$

Skew symmetry of the generator $A_U$ (Proposition 4), symmetry of $\mathcal{L}_k^{-1/2}$ (A.13) and as $\lambda$ is real allow us to write

$$\langle j \circ x_f, \tilde{A}_U x_g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = i \langle A_U(x_f), j \circ x_g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} + \lambda \langle \mathcal{L}_k^{-1/2}(x_f), j \circ x_g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = \langle \tilde{A}_U x_f, j \circ x_g \rangle_{L^2_{\mathcal{H}}(\Omega, v)}.$$

With similar arguments as (D.3), we conclude that $\langle \tilde{A}_U^{-1}(f), g \rangle_{L^2_{\mathcal{H}}(\Omega, v)} = \langle f, \tilde{A}_U^{-1}(g) \rangle_{L^2_{\mathcal{H}}(\Omega, v)}$ for all $f$ and $g \in L^2_{\mathcal{H}}(\Omega, v)$. □
With lemmas 5 and 6, we may apply the spectral theorem of compact self-adjoint operator and consequently there exists a orthonormal basis $(e_n)_{n \ge 0}$ of $(L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$ and a sequence of real $(\alpha_n)_{n \in (\sigma_p(\tilde{A}_U^{-1}) - \{0\})^\mathbb{N}}$ such that

$$ \lim_{n \to \infty} \alpha_n = 0 \quad \text{and} \quad \tilde{A}_U^{-1}(e_n) = \alpha_n e_n \quad \text{for all } n \in \mathbb{N}. \qquad (\text{D.4}) $$

Moreover, for all integers $n$, the real eigenvalue $\alpha_n$ has finite multiplicity. It can be noticed that $\alpha_n \ne 0$ for all $n$ since $\tilde{A}_U^{-1}$ is invertible and the eigenfunctions $e_n$ belong to the RKHS $\mathcal{H}$. Finally, we have for all $n \in \mathbb{N}$

$$ \tilde{A}_U(e_n) = \left[ iA_U + \lambda \mathcal{L}_k^{-1/2} \right] (e_n) = \alpha_n^{-1} e_n. $$

We now are in position to fully specify the eigenvalues of $A_U$. If all the eigenvalues of $\tilde{A}_U^{-1}$ are simple, it is easier to relate them to the eigenvalues of $A_U$. This is explained in Remark 20. In the other case, we need first to establish a commutativity property between $\tilde{A}_U^{-1}$ and $\mathcal{L}_k$.

**Lemma 7 (Commutation between $A_U$ and $\mathcal{L}_k$)** The operators $A_U$ and $\mathcal{L}_k$ are commutative on $\mathcal{H}$. For all $f \in \mathcal{H}$ we have the identity: $A_U \circ \mathcal{L}_k(f) = \mathcal{L}_k \circ A_U(f)$.

Before demonstrating this result, we first detail the commutation between $\mathcal{L}_k$ and $U_t$ on $\mathcal{H}^{\text{sp}}$. Considering $f \in \mathcal{H}^{\text{sp}}$ the element $\mathcal{L}_k(f)$ belongs to $\mathcal{H}$. We apply Proposition 2 and we have for all $X \in \Omega$ and all $t \ge 0$

$$ \mathcal{L}_k[U_t(f)](X) = \int_{\Omega} k(X,y) U_t(f)(y) \nu(dy) = \int_{\Omega} k(\Phi_t(X), \Phi_t(y)) f \circ \Phi_t(y) \nu(dy). $$

Since the dynamical system $(\Phi_t)_{t \ge 0}$ is measure preserving and $\mathcal{L}_k(f)$ belongs to $\mathcal{H}$, we obtain

$$ \mathcal{L}_k[U_t(f)](X) = \int_{\Omega} k(\Phi_t(X), y) f(y) \nu(dy) = [\mathcal{L}_k f](\Phi_t(X)) = U_t[\mathcal{L}_k f](X). $$

**PROOF (LEMMA 7)** First, we prove the result for $f \in \mathcal{H}^{\text{sp}}$. By the previous equality we have for all $g \in L_{\mathbb{C}}^2(\Omega, \nu)$ and $t \ge 0$

$$ \left\langle \frac{U_t[\mathcal{L}_k f] - \mathcal{L}_k f}{t}, g \right\rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \left\langle \frac{\mathcal{L}_k[U_t f] - \mathcal{L}_k f}{t}, g \right\rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. $$

The operator $\mathcal{L}_k$ being self-adjoint on $L_{\mathbb{C}}^2(\Omega, \nu)$ we infer that

$$ \left\langle \frac{U_t[\mathcal{L}_k f] - \mathcal{L}_k f}{t}, g \right\rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \left\langle \frac{U_t f - f}{t}, \mathcal{L}_k(g) \right\rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. $$

Taking the limit as $[t \to 0^+]$ in the left-hand side, as $\mathcal{L}_k(f)$ belongs to $\mathcal{H} \subset \mathcal{D}(A_U)$, we have

$$ \lim_{t \to 0^+} \left\langle \frac{U_t[\mathcal{L}_k f] - \mathcal{L}_k f}{t}, g \right\rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle A_U \mathcal{L}_k(f), g \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. $$

Doing the same in the right-hand side, as $f$ belongs to $\mathcal{H} \subset \mathcal{D}(A_U) = \mathcal{D}(A_P)$, we also have

$$ \lim_{t \to 0^+} \left\langle \frac{U_t f - f}{t}, \mathcal{L}_k(g) \right\rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle A_U f, \mathcal{L}_k(g) \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle \mathcal{L}_k \circ A_U f, g \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}. $$

We deduce that $\langle A_U \circ \mathcal{L}_k(f), g \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} = \langle \mathcal{L}_k \circ A_U(f), g \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)}$ for all $g \in L_{\mathbb{C}}^2(\Omega, \nu)$. The lemma is proved on $\mathcal{H}^{\text{sp}}$. Then, we suppose that $f \in \mathcal{H}$. There exists a sequence $(f_n)_n$ of $\mathcal{H}^{\text{sp}}$ which converges to $f$ in $\mathcal{H}$.

The continuity of $A_U : \mathcal{H} \to L_{\mathbb{C}}^2(E, \nu)$ and Remark 15 justifies for all integers $n$

$$ \|A_U \circ \mathcal{L}_k(f_n) - A_U \circ \mathcal{L}_k(f)\|_{L_{\mathbb{C}}^2(E, \nu)} \le C \|\mathcal{L}_k(f_n - f)\|_{\mathcal{H}} \le C \|f_n - f\|_{\mathcal{H}}. $$
As $L_k : L_{\mathbb{C}}^2(E, \nu) \to L_{\mathbb{C}}^2(E, \nu)$ is continuous, we infer

$$ \| \mathcal{L}_k \circ A_U(f_n) - \mathcal{L}_k \circ A_U(f) \|_{L_{\mathbb{C}}^2(E, \nu)} \le C \| A_U(f_n - f) \|_{L_{\mathbb{C}}^2(E, \nu)} \le C \|f_n - f\|_{\mathcal{H}}. $$

For all integers $n$, we have $A_U \circ \mathcal{L}_k(f_n) = \mathcal{L}_k \circ A_U(f_n)$. We pass to the limit to obtain $A_U \circ \mathcal{L}_k(f) = \mathcal{L}_k \circ A_U(f)$. $\square$

**Lemma 8 (Commutation between $\tilde{A}_U^{-1}$ and $\mathcal{L}_k$.)** The operators $\tilde{A}_U^{-1}$ and $\mathcal{L}_k$ are commutative on $L_{\mathbb{C}}^2(\Omega, \nu)$.

**PROOF** We apply the lemma 7 and we infer for all $f \in \mathcal{H}$

$$ \mathcal{L}_k \circ \tilde{A}_U(f) = i \mathcal{L}_k \circ A_U(f) + \lambda \mathcal{L}_k \circ \mathcal{L}_k^{-1/2}(f) = i A_U \circ \mathcal{L}_k(f) + \lambda \mathcal{L}_k^{-1/2} \circ \mathcal{L}_k(f) = \tilde{A}_U \circ \mathcal{L}_k(f). $$

As $\tilde{A}_U : \mathcal{H} \mapsto L_{\mathbb{C}}^2(\Omega, \nu)$ is bijective, we obtain that $\mathcal{L}_k \circ \tilde{A}_U^{-1}(g) = \tilde{A}_U^{-1} \circ \mathcal{L}_k(g)$ for all $g \in L_{\mathbb{C}}^2(\Omega, \nu)$.

**Proposition 15 (Diagonalization of $\tilde{A}_U^{-1}$ and $\mathcal{L}_k$)** The operators $\tilde{A}_U^{-1}$ and $\mathcal{L}_k$ are diagonalizable in the same orthonormal basis of $L_{\mathbb{C}}^2(\Omega, \nu)$.

**PROOF** For $n \in \mathbb{N}$, we note $E_{\alpha_n} \subset L_{\mathbb{C}}^2(\Omega, \nu)$ the eigenspace of $\tilde{A}_U^{-1}$ associated with eigenvalue $\alpha_n$. By Lemma 8, we have for all $f \in E_{\alpha_n}$ that

$$ [\tilde{A}_U^{-1} - \alpha_n I] \circ \mathcal{L}_k(f) = \mathcal{L}_k \circ [\tilde{A}_U^{-1} - \alpha_n I](f) = 0, $$

and $E_{\alpha_n}$ is thus an invariant subspace of $\mathcal{L}_k$. The operator $\mathcal{L}_k : E_{\alpha_n} \mapsto E_{\alpha_n}$ can be diagonalized. There exists an orthonormal basis $(f_{n,p})_p$ of $E_{\alpha_n}$ and reals $(\beta_{n,p})_p$ (with $1 \le p \le \dim E_{\alpha_n}$) such that

$$ \mathcal{L}_k(f_{n,p}) = \beta_{n,p} f_{n,p} \quad \text{and} \quad \tilde{A}_U^{-1}(f_{n,p}) = \alpha_n f_{n,p}. $$

The family $(f_{n,p})_p$ diagonalizes both the operators $\mathcal{L}_k$ and $\tilde{A}_U^{-1}$ on $E_{\alpha_n}$. It can be noticed that $\beta_{n,p} > 0$ since $\mathcal{L}_k$ is positive definite. As $\alpha_n \ne 0$ for all $n$, the functions $f_{n,p}$ belong to $\mathcal{H}$. The concatenation of the vectors $(f_{n,p})_{n,p}$ for all $n \in \mathbb{N}$ and $1 \le p \le \dim E_{\alpha_n}$ forms a orthonormal basis of $L_{\mathbb{C}}^2(\Omega, \nu)$. As a matter of fact, we have by construction that, $\text{Span}\{f_{n,p} : n \in \mathbb{N}$ and $1 \le p \le \dim E_{\alpha_n}\} = \text{Span}\{e_n : n \in \mathbb{N}\}$ and the family $(f_{n,p})_{n,p}$ is also complete in $L_{\mathbb{C}}^2(\Omega, \nu)$. For all $n$ the family $(f_{n,p})_p$ is orthonormal in $L_{\mathbb{C}}^2(\Omega, \nu)$ (by construction). The eigenspaces of $\tilde{A}_U^{-1}$ are in addition orthogonal: the vectors $f_{n,p}$ and $f_{m,q}$ are orthogonal in $L_{\mathbb{C}}^2(\Omega, \nu)$ for all $n \ne m$, $1 \le p \le \dim E_{\alpha_n}$ and $1 \le q \le \dim E_{\alpha_n}$. The family $(f_{n,p})_{n \in \mathbb{N}, 1 \le p \le \dim E_{\alpha_n}}$ is thus orthonormal in $L_{\mathbb{C}}^2(\Omega, \nu)$ and constitutes a orthonormal basis of $L_{\mathbb{C}}^2(\Omega, \nu)$. $\square$

The result of step 1 can now be obtained. Keeping the same notations for the eigenfunctions of $\tilde{A}_U$, for all $n \in \mathbb{N}$ and $1 \le p \le \dim E_{\alpha_n}$, we have $\tilde{A}_U(f_{n,p}) = \alpha_n^{-1} f_{n,p}$ and $\mathcal{L}_k^{-1/2}(f_{n,p}) = \beta_{n,p}^{-1/2} f_{n,p}$. With the expression of $A_U$ in (D.1) we obtain

$$ A_U(f_{n,p}) = i(\lambda \beta_{n,p}^{-1/2} - \alpha_n^{-1})f_{n,p}. \qquad (D.5) $$

Denoting for the sake of clarity, $(V_\ell)_\ell$ the orthonormal basis $(f_{n,p})_{n,p}$ of $L_{\mathbb{C}}^2(\Omega, \nu)$ and $(\lambda_\ell)_\ell$ the associated eigenvalues of $A_U$ (given (D.5) with $\lambda_\ell \in i\mathbb{R}$) and upon applying $\mathcal{L}_k^{1/2}$ on the orthonormal basis $(V_\ell)_{\ell \in \mathbb{N}}$ of $L_{\mathbb{C}}^2(\Omega, \nu)$, we get also an orthonormal basis of $\mathcal{H}$ (Remark 18) which diagonalizes $A_U$:

$$ A_U(\mathcal{L}_k^{1/2} V_\ell) = \lambda_\ell \mathcal{L}_k^{1/2} V_\ell. \qquad (D.6) $$

This completes the proof of step 1.

It can be noticed that for all $g \in \mathcal{H}$ we have

$$ A_U(g) = \sum_{\ell=0}^{\infty} \lambda_{\ell} \langle g, V_{\ell} \rangle_{L_{\mathbb{C}}^2(\Omega, \nu)} V_{\ell} \quad \text{and} \quad A_U(g) = \sum_{\ell=0}^{\infty} \lambda_{\ell} \langle g, \mathcal{L}_k^{1/2}(V_{\ell}) \rangle_{\mathcal{H}} \mathcal{L}_k^{1/2}(V_{\ell}). $$
**Remark 20 (Particular case: $\tilde{A}_U$ has simple eigenvalues)** A simpler case arises when the eigenvalues of $\tilde{A}_U$ are simple. This corresponds to the case of ergodic dynamical systems [12]. For all $n$, by (D.4) and Lemma 7, the eigenspace $E_{\alpha_n}$ of $\tilde{A}_U$ is an invariant space of $\mathcal{L}_k^{-1/2}$. As $\dim E_{\alpha_n} = 1$ and $E_{\alpha_n}$ is an invariant subspace of $\mathcal{L}_k^{-1/2}$, there exists reals $\beta_n$ such that $\mathcal{L}_k^{-1/2}e_n = \beta_n e_n$ and we obtain

$$A_U(e_n) = i(\lambda \beta_n - \alpha_n^{-1})e_n,$$

with the orthonormal basis $(e_n)_n$ of $L_{\mathbb{C}}^2(\Omega, \nu)$ obtained from the diagonalization of $\tilde{A}_U^{-1}$. We have as well for the associated orthonormal basis of $\mathcal{H}$

$$A_U(\mathcal{L}_k^{1/2}e_n) = i(\lambda \beta_n - \alpha_n^{-1})(\mathcal{L}_k^{1/2}e_n).$$

**Remark 21 (Exponential form of the Koopman semi-group)** The range of $\tilde{A}_U^{-1}$ is included in $\mathcal{H}$. It can be also noticed that $V_\ell$ belongs to $\mathcal{H}$ and the following inclusions hold

$$\mathrm{Span}\{V_\ell : \ell \in \mathbb{N}\} \subset \mathcal{H} \subset L_{\mathbb{C}}^2(\Omega, \nu).$$

The set $\mathrm{Span}\{V_\ell : \ell \in \mathbb{N}\}$ is dense in $(L_{\mathbb{C}}^2(\Omega, \nu), \|\cdot\|_{L_{\mathbb{C}}^2(\Omega, \nu)})$ and is an invariant space of the Koopman infinitesimal generator $A_U$. For all $f$ belonging to this set and all $t \ge 0$, as $A_U$ is bounded, the Koopman operator is uniformly continuous and can be written in terms of an exponential series as

$$U_t(f) = \exp(tA_U)(f),$$

with convergence in the operator norm. Its adjoint, the Perron-Frobenius operator, can be written in the same way as

$$P_t(f) = \exp(-tA_U)(f).$$

Denoting for the sake of simplicity, $\psi_\ell = \mathcal{L}_k^{1/2}(V_\ell)$ for all integers $\ell$. The family $(\psi_\ell)_\ell$ is an orthonormal basis of $\mathcal{H}$. We have $U_t(\psi_\ell) = \exp(t\lambda_\ell)\psi_\ell$ and we obtain hence

$$\psi_\ell(X_t) = \exp(t\lambda_\ell) \psi_\ell(X_0). \qquad (D.7)$$

### Step 2

From this orthonormal basis $(\psi_\ell)_\ell$ of $\mathcal{H}$, we now build a orthonormal basis of $\mathcal{H}_t$ that diagonalizes the operator $A_{U,t}$ for all time $t \ge 0$.

**Proposition 16 (Hilbert basis of $\mathcal{H}_t$)** For all $t \ge 0$, the family $(R_t(\psi_\ell))_\ell$ is an orthonormal basis of $(\mathcal{H}_t, \|\cdot\|_{\mathcal{H}_t})$.

**PROOF** For all $\ell$ and $\ell'$, upon applying Proposition 12 we have $\langle R_t(\psi_\ell), R_t(\psi_{\ell'}) \rangle_{\mathcal{H}_t} = \langle \psi_\ell, \psi_{\ell'} \rangle_{\mathcal{H}} = \delta_{\ell, \ell'}$, consequently the family $(R_t(\psi_\ell))_\ell$ is orthonormal in $\mathcal{H}_t$. Besides, the family $\mathrm{Span}\{R_t(\psi_\ell) : \ell \in \mathbb{N}\}$ is dense in $\mathcal{H}_t$. Let $h \in \mathcal{H}_t$ be such that $\langle h, R_t(\psi_\ell) \rangle_{\mathcal{H}_t} = 0$ for all integers $\ell$. By Proposition 13, we have $\langle h, R_t(\psi_\ell) \rangle_{\mathcal{H}_t} = \langle E_t(h), \psi_\ell \rangle_{\mathcal{H}} = 0$. As $(\psi_\ell)_\ell$ is an orthonormal basis of $\mathcal{H}$, we conclude that $E_t(h) = 0$ and thus $h = 0$ on $\mathcal{H}_t$. $\square$

The orthonormal basis $(R_t(\psi_\ell))_\ell$ of $\mathcal{H}_t$ diagonalizes the operator $A_{U,t}$. As a matter of fact, for all $\ell$,

$$A_{U,t}[R_t \psi_\ell] = R_t \circ A_U \circ E_t \circ R_t(\psi_\ell) = \lambda_\ell [R_t \psi_\ell].$$

For conciseness purpose, we note in the following $\psi_\ell^t := R_t(\psi_\ell)$. Through Remark 19 and (D.7), it can be noticed that for all $X_0 \in \Omega_0$

$$\psi_\ell^t(X_t) = \exp(t\lambda_\ell) \psi_\ell^0(X_0),$$

which ends the proof of step 2, and completes the proof of the spectral representation of Koopman operator in the RKHS family.

**Remark 22** The set $\mathrm{Span}\{\psi_\ell^t : \ell \in \mathbb{N}\}$ is an invariant space of $A_{U,t}$. As a matter of fact, considering remark 21 and proposition 12, for all $f \in \mathrm{Span}\{\psi_\ell^t : \ell \in \mathbb{N}\}$ we have $\exp(tA_{U,t})(f) = R_t \circ \exp(tA_U) \circ E_t(f)$.
## Appendix E. Numerical details of the QG implementation

The QG model is discretized on a staggered grid of size $N_x \times N_y$, where $N_x = 64$ and $N_y = 128$, with order-2 finite differences. Advection term is computed using the Arakawa discretization [89] and centered differences are considered for the other terms. The Poisson equation is solved using spectral method by projection on a sine basis. Time integration is performed with a fourth-order Runge-Kutta scheme. An ensemble of $p = 100$ members is considered for the training, and another test ensemble of $p_2 = 100$ members is generated for predictions assessment. Starting from an established state of a particular trajectory, a set of initial conditions is generated by perturbing this state by the potential vorticity of $N_{\text{POD}} = 50$ proper orthogonal decomposition (POD) modes [90] weighted by a Gaussian random perturbation $N(0, 10\lambda_i/\lambda_1)$, where $\lambda_i$ is the $i^{\text{th}}$ POD eigenvalue. Then, a spin-up time integration of $t = 0.03$ is performed to evacuate eventual nonphysical flow features. The training and testing sets are generated independently using the same rule.

Following the dimensioning given in [91], $R_o = \frac{U_r}{\beta L_r^2} = 0.0036$ with $U_r$ and $L_r$ the reference velocity and length scale, respectively. We give as well the gradient of the Coriolis parameter $\beta = 1.75 \cdot 10^{-11}$, typical of mid-latitudes at the center of the basin. Setting $L_r = 2000$ km, leads to $U_r = 0.252 \text{ m.s}^{-1}$, and then to a reference time $T_r = \frac{L_r}{U_r} = 92$ days.

## References

[1] M. Bocquet, J. Brajard, A. Carrassi, L. Bertino, Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization, F. of Data Science 2 (2020) 55–80.

[2] S. Brunton, J. Proctor, N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (2016) 3932–3937.

[3] R. Fablet, B. Chapron, L. Drumetz, E. Mémin, O. Pannekoucke, F. Rousseau, Learning variational data assimilation models and solvers, Journal of Advances in Modeling Earth Systems 13 (2021) e2021MS002572. E2021MS002572 2021MS002572.

[4] G. A. Gottwald, S. Reich, Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation, Physica D: Nonlinear Phenomena 423 (2021) 132911.

[5] B. Hamzi, H. Owhadi, Learning dynamical systems from data: A simple cross-validation perspective, part i: Parametric kernel flows, Physica D: Nonlinear Phenomena 421 (2021) 132817.

[6] H. Owhadi, G. R. Yoo, Kernel flows: From learning kernels from data into the abyss, Journal of Computational Physics 389 (2019) 22–47.

[7] J. Pathak, Z. Lu, B. Hunt, M. Girvan, E. Ott, Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (2017) 121102.

[8] J. Pathak, B. Hunt, M. Girvan, Z. Lu, E. Ott, Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach, Phys. Rev. Lett. 120 (2018) 024102.

[9] M. Raissi, P. Perdikaris, G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.

[10] Z. Zhao, D. Giannakis, Analog forecasting with dynamics-adapted kernels, Nonlinearity 29 (2016) 2888–2939.

[11] B. O. Koopman, Hamiltonian systems and transformation in Hilbert space, Proceedings of the National Academy of Sciences 17 (1931) 315–318.

[12] T. Eisner, B. Farkas, M. Haase, R. Nagel, Operator Theoretic Aspects of Ergodic Theory, Graduate Texts in Mathematics, Springer, 2015.

[13] M. Dellnitz, O. Junge, On the approximation of complicated dynamical behavior, SIAM Journal on Numerical Analysis 36 (1999) 491–515.

[14] I. Mezic, Spectral properties of dynamical systems, model reduction and decomposition, Nonlinear Dynamics 41 (2005) 309–325.

[15] M. Budišić, R. Mohr, I. Mezić, Applied Koopmanism, Chaos: An Interdisciplinary Journal of Nonlinear Science 22 (2012) 047510.

[16] S. Das, D. Giannakis, Delay-coordinate maps and the spectra of Koopman operators, J. Stat. Phys. 175 (2019) 1107–1145.

[17] S. Das, D. Giannakis, Koopman spectra in reproducing kernel Hilbert spaces, Appl. Comput. Harmon. Anal. 49 (2020) 573–607.

[18] C. Rowley, I. Mezic, S. Bagheri, P. Schlatter, D. Henningson, Spectral analysis of nonlinear flows, J. Fluid Mech. 641 (2009) 115–127.

[19] P. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech. 656 (2010) 5–28.

[20] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, J. N. Kutz, On dynamic mode decomposition: Theory and applications, Journal of Computational Dynamics 1 (2014) 391–421.

[21] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25 (2015) 1307–1346.

[22] M. O. Williams, C. W. Rowley, I. G. Kevrekidis, A kernel-based method for data-driven Koopman spectral analysis, Journal of Computational Dynamics 2 (2015) 247–265.

[23] N. Kutz, X. Fu, S. Brunton, Multiresolution dynamic mode decomposition, SIAM Journal on Applied Dynamical Systems 15 (2016) 713–735.

[24] A. M. DeGennaro, N. M. Urban, Scalable extended dynamic mode decomposition using random kernel approximation, SIAM Journal on Scientific Computing 41 (2019) A1482–A1499.

[25] G. Buza, S. Jain, G. Haller, Using spectral submanifolds for optimal mode selection in nonlinear model reduction, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477 (2021) 20200725.

[26] S. Pan, N. Arnold-Medabalimi, K. Duraisamy, Sparsity-promoting algorithms for the discovery of informative Koopman-invariant subspaces, Journal of Fluid Mechanics 917 (2021) A18.
[27] M. J. Colbrook, A. Townsend, Rigorous data-driven computation of spectral properties of Koopman operators for dynamical systems, arXiv preprint https://arxiv.org/abs/2111.14889 (2021).

[28] P. J. Baddoo, B. Herrmann, B. J. McKeon, S. L. Brunton, Kernel learning for robust dynamic mode decomposition: linear and nonlinear disambiguation optimization, Proceedings of the Royal Society A 478 (2022) 20210830.

[29] S. Brunton, B. Brunton, J. Proctor, E. Kaiser, J. Kutz, Chaos as an intermittently forced linear system, Nat. Commun. 8 (2017).

[30] D. Giannakis, J. Slawinska, Z. Zhao, Spatiotemporal feature extraction with data-driven Koopman operators, in: Journ. of Mach. Learn. Res.: Workshop and Conference Proceedings, volume 44, 2015, pp. 103–115.

[31] D. Giannakis, Data-driven spectral decomposition and forecasting of ergodic dynamical systems, Appl. Comput. Harmon. Anal. https://doi.org/10.1016/j.acha.2017.09.001 (2017).

[32] R. Vautard, M. Ghil, Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series, Physica D: Nonlinear Phenomena 35 (1989) 395–424.

[33] D. Kondrashov, E. A. Ryzhov, P. Berloff, Data-adaptive harmonic analysis of oceanic waves and turbulent flows, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2020), 061105.

[34] Y. Zhen, B. Chapron, E. Mémin, L. Peng, Eigenvalues of autocovariance matrix: A practical method to identify the Koopman eigenfrequencies, Phys. Rev. E 105 (2022) 034205.

[35] Y. Zhen, B. Chapron, E. Mémin, Bridging Koopman operator and time-series auto-correlation based Hilbert–Schmidt operator, in: B. Chapron, D. Crisan, D. Holm, E. Mémin, A. Radomska (Eds.), Stochastic Transport in Upper Ocean Dynamics, Springer International Publishing, 2023, pp. 301–316.

[36] I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics 41 (2005) 309–325.

[37] S. Das, D. Giannakis, J. Slawinska, Reproducing kernel Hilbert space compactification of unitary evolution groups, Applied and Computational Harmonic Analysis 54 (2021) 75–136.

[38] F. Gilani, D. Giannakis, J. Harlim, Kernel-based prediction of non-Markovian time series, Physica D: Nonlinear Phenomena 418 (2021) 132829.

[39] R. Alexander, D. Giannakis, Operator-theoretic framework for forecasting nonlinear time series with kernel analog techniques, Physica D: Nonlinear Phenomena 409 (2020) 132520.

[40] D. Burov, D. Giannakis, K. Manohar, A. Stuart, Kernel analog forecasting: Multiscale test problems, Multiscale Model. Simul. 19 (2021) 1011–1040.

[41] S. Klus, I. Schuster, K. Muandet, Eigendecompositions of transfer operators in reproducing kernel Hilbert spaces, Journal of Nonlinear Science 30 (2020) 283–315.

[42] N. Santitissadeekorn, E. M. Bollt, Ensemble-based method for the inverse Frobenius–Perron operator problem: Data-driven global analysis from spatiotemporal “movie” data, Physica D: Nonlinear Phenomena 411 (2020) 132603.

[43] J. Lee, E. De Brouwer, B. Hamzi, H. Owhadi, Learning dynamical systems from data: A simple cross-validation perspective, part iii: Irregularly-sampled time series, Physica D: Nonlinear Phenomena 443 (2023) 133546.

[44] B. Hamzi, R. Maulik, H. Owhadi, Simple, low-cost and accurate data-driven geophysical forecasting with learned kernels, Proceedings of the Royal Society A 477 (2021).

[45] H. Arbabi, I. Mezic, Ergodic theory, dynamic mode decomposition and computation of spectral properties of the Koopman operator, SIAM Journal on Applied Dynamical Systems 16 (2017) 2096–2126.

[46] S. Klus, F. Nüske, S. Peitz, J.-H. Niemann, C. Clementi, C. Schütte, Data-driven approximation of the Koopman generator: Model reduction, system identification, and control, Physica D: Nonlinear Phenomena 406 (2020) 132416.

[47] S. L. Brunton, M. Budišić, E. Kaiser, J. N. Kutz, Modern Koopman theory for dynamical systems, SIAM Review 64 (2022) 229–340.

[48] P. Bevanda, S. Sosnowski, S. Hirche, Koopman operator dynamical models: Learning, analysis and control, Annual Reviews in Control 52 (2021) 197–212.

[49] C. R. S. Otto, Koopman operators for estimation and control of dynamical systems, Annual Reviews of Control, Robotics, and Autonomous systems 4 (2021) 59–87.

[50] V. Kostic, P. Novelli, A. Maurer, C. Ciliberto, L. Rosasco, M. Pontil, Learning dynamical systems via Koopman operator regression in reproducing kernel Hilbert spaces., in: Advances in Neural Information Processing Systems, 2022.

[51] V. Kostic, K. Lounici, P. Novelli, M. Pontil, Koopman operator learning: Sharp spectral rates and spurious eigenvalues, arXiv, https://arxiv.org/abs/2302.02004, 2023.

[52] S. Klus, F. Nüske, B. Hamzi, Kernel-based approximation of the Koopman generator and schrödinger operator, Entropy 22 (2020).

[53] T. Berry, T. Sauer, Local kernels and the geometric structure of data, Applied and Computational Harmonic Analysis 40 (2016) 439–469.

[54] R. Banisch, P. Koltai, Understanding the geometry of transport: Diffusion maps for Lagrangian trajectory data unravel coherent sets, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (2017) 035804.

[55] R. R. Coifman, M. J. Hirn, Diffusion maps for changing data, Applied and Computational Harmonic Analysis 36 (2014) 79–107.

[56] D.-X. Zhou, Derivative reproducing properties for kernel methods in learning theory, Journal of Computational and Applied Mathematics 220 (2008) 456–463.

[57] T. Berlinet, C. Thomas-Agnan, Reproducing kernel Hilbert spaces in Probability and Statistics, Kluwer Academic Publishers, 2001.

[58] F. Cucker, S. Smale, On the mathematical foundation of learning, Bul. of the Amer. Math. Soc. 39 (2001) 1–49.

[59] N. Aronszajn, Theory of reproducing kernels, Trans. of the American Math. Society 68 (1950).

[60] H. König, Eigenvalue distribution of compact operators with application to integral operators, Linear Algebra and its application 84 (1986) 111–122.

[61] F.-X. Le Dimet, O. Talagrand, Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects, Tellus 38A (1986) 97–110.

[62] J. L. Lions, Optimal control of systems governed by PDEs, Springer-Verlag, New-York, 1971.

[63] C. Carmeli, E. De Vito, A. Toigo, Vector valued reproducing kernel Hilbert spaces of integrable functions and Mercer theorem, Analysis and Applications 04 (2006) 377–408.
[64] I. Steinwart, C. Scovel, Mercer's theorem on general domains: On the interaction between measures, kernels, and rkhss, Constructive Approximation 35 (2012) 363-417.

[65] O. A. Ladyzhenskaya, Finite-dimensionality of bounded invariant sets for navier-stokes systems and other dissipative systems, Zapiski Nauchnykh Seminarov POMI 115 (1982) 137-155.

[66] P. Constantin, C. Foiş, R. Temam, Attractors representing turbulent flows, volume 314, American Mathematical Soc., 1985.

[67] P. Lax, Integrals of nonlinear equations of evolution and solitary waves, Communications on pure and applied mathematics 21 (1968) 467-490.

[68] Y. C. Li, Ergodic isospectral theory of the Lax pairs of Euler equations with harmonic analysis flavor, Proc. of the Amer. Math. Soc. 133 (2005) 2681-2687.

[69] R. J. Greatbatch, B. T. Nadiga, Four-gyre circulation in a barotropic model with double-gyre wind forcing, Journal of Physical Oceanography 30 (2000) 1461-1471.

[70] J. R. Brannan, J. Duan, T. Wanner, Dissipative quasi-geostrophic dynamics under random forcing, Journal of Mathematical Analysis and Applications 228 (1998) 221-233.

[71] J. Duan, B. Goldys, Ergodicity of stochastically forced large scale geophysical flows, International Journal of Mathematics and Mathematical Sciences 28 (2001) 313-320.

[72] L. Yang, X. Pu, Ergodicity of large scale stochastic geophysical flows with degenerate gaussian noise, Applied Mathematics Letters 64 (2017) 27-33.

[73] T. M. Hamill, J. S. Whitaker, C. Snyder, Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter, Monthly Weather Review 129 (2001) 2776-2790.

[74] H. Q. Minh, Some properties of Gaussian reproducing kernel Hilbert spaces and their implications for function approximation and learning theory, Constructive Approximation 32 (2010) 307-338.

[75] P. Y. Le Traon, F. Nadal, N. Ducet, An improved mapping method of multisatellite altimeter data, Journal of Atmospheric and Oceanic Technology 15 (1998) 522-534.

[76] N. Ducet, P. Y. Le Traon, G. Reverdin, Global high-resolution mapping of ocean circulation from topex/poseidon and ers-1 and -2, Journal of Geophysical Research: Oceans 105 (2000) 19477-19498.

[77] AVISO, SSALTO/DUACS User Handbook: MSLA and (M)ADT Near-Real Time and Delayed Time Products, Technical Report CLS-DOS-NT-06-034, SALP-MU-P-EA-21065-CLS, 5rev0, CNES, 2016. Available online at https://www.aviso.altimetry.fr/fileadmin/documents/data/tools/hdbk_duacs.pdf.

[78] C. Ubelmann, B. Cornuelle, L.-L. Fu, Dynamic mapping of along-track ocean altimetry: Method and performance from observing system simulation experiments, Journal of Atmospheric and Oceanic Technology 33 (2016) 1691-1699.

[79] F. Guillou, S. Metref, E. Cosme, C. Ubelmann, M. Ballarotta, J. Le Sommer, J. Verron, Mapping altimetry in the forthcoming swot era by back-and-forth nudging a one-layer quasigeostrophic model, Journal of Atmospheric and Oceanic Technology 38 (2021) 697-710.

[80] S. Chaturantabut, D. C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal of Scientific Computing 32 (2010) 2737-2764.

[81] E. Mémin, Fluid flow dynamics under location uncertainty, Geophys. & Astro. Fluid Dyn. 108 (2014) 119-146.

[82] W. Bauer, P. Chandramouli, B. Chapron, L. Li, E. Mémin, Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: A stochastic approach, Journal of Physical Oceanography 50 (01 Apr. 2020) 983-1003.

[83] B. Chapron, P. Dérrian, E. Mémin, V. Resseguier, Large-scale flows under location uncertainty: a consistent stochastic framework, QJRMS 144 (2018) 251-260.

[84] A. Debussche, B. Hug, E. Mémin, A consistent stochastic large-scale representation of the Navier-Stokes equations, Journal of Mathematical Fluid Mechanics 25 (2023) 19.

[85] B. Dufée, E. Mémin, D. Crisan, Observation-Based Noise Calibration: An Efficient Dynamics for the Ensemble Kalman Filter, in: Stochastic Transport in Upper Ocean Dynamics, volume 10 of Mathematics of Planet Earth, Springer International Publishing, 2023, pp. 43-56. URL: https://hal.science/hal-03910764. doi:10.1007/978-3-031-18988-3_4.

[86] B. Dufée, E. Mémin, D. Crisan, Stochastic parametrization: An alternative to inflation in ensemble Kalman filters, Quarterly Journal of the Royal Meteorological Society 148 (2022) 1075-1091.

[87] C.-J. Simon-Gabriel, B. Schölkopf, Kernel distribution embeddings: Universal kernels, characteristic kernels and kernel metrics on distributions, Journal of Machine Learning Research 19 (2018) 1-29.

[88] I. Steinwart, A. Christmann, Support Vector Machines. Information Science and Statistics, Springer, 2008.

[89] A. Arakawa, Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I, Journal of Computational Physics 1 (1966) 119-143.

[90] J. L. Lumley, The structure of inhomogeneous turbulent flows, Atmospheric turbulence and radio wave propagation 1 (1967) 166-178.

[91] O. San, A. E. Staples, Z. Wang, T. Iliescu, Approximate deconvolution large eddy simulation of a barotropic ocean circulation model, Ocean Modelling 40 (2011) 120-132.