Providing efficient and accurate parametrizations for model reduction is a key goal in many areas of science and technology. Here we present a strong link between data-driven and theoretical approaches to achieving this
goal. Formal perturbation expansions of the Koopman operator allow us to derive general stochastic parametrizations of weakly coupled dynamical systems. Such parametrizations yield a set of stochastic integro-differential equations with explicit noise and memory kernel formulas to describe the effects of unresolved variables. We show that the perturbation expansions involved need not be truncated when the coupling is additive.
The unwieldy integro-differential equations can be recast as a simpler multilevel Markovian model, and we establish an intuitive connection with a
generalized Langevin equation. This connection helps setting up a parallelism between the top-down, equations-based methodology herein and the well-established empirical model reduction (EMR) methodology that has been shown to provide efficient dynamical closures to partially observed systems. Hence, our findings support, on the one hand, the physical basis and robustness of the EMR methodology and, on the other hand, illustrate the practical relevance of the perturbative expansion used for deriving the parametrizations.

Parametrizations aim to reduce the complexity of high-dimensional dynamical systems. Here, a theory-based and a data-driven approach
for the parametrization of coupled systems are compared, showing that both yield the same stochastic multilevel structure. The results provide a very strong support to the use of empirical methods in model reduction and
clarify the practical relevance of the proposed theoretical framework.

§ INTRODUCTION AND MOTIVATION

Multiscale systems are typically characterized by the presence of significant variability over a large range of spatial and temporal scales. This multiscale character is due to a combination of the following factors: the nature of the external forcings, the inhomogeneity of the properties of
the system's various components, the complexity of the coupling mechanisms between the components, and the variety of instabilities, dissipative processes, and feedbacks acting at different scales. In many cases, both the theoretical understanding of such systems and the formulation of numerical models for simulating their properties are based on focusing upon a reduced range of large spatial and long temporal scales of interest, and upon devising an efficient way to capture effectively the impact of the faster dynamical processes acting predominantly in the neglected smaller spatial scales <cit.>.

Given a high-dimensional dynamical system, we are thus interested in reformulating it in such a way that only the variables of interest are resolved.  A first-guess would be to ignore altogether the unresolved variables
and consider uniquely the filtered evolution laws valid for the targeted range of scales. However, this is well-known to be inadequate, because nonlinearity guarantees that the unresolved variables have an impact on the
resolved ones, both in terms of the detailed dynamics and of its statistical properties. Therefore, the problem of constructing accurate and efficient reduced-order models — or, equivalently, of defining the coarse-grained dynamics — is an essential and fundamental aspect of studying multiscale systems, both theoretically and through numerical simulations.

For the sake of concreteness, let us consider, the case of climate science. Is It is well-nigh impossible, therefore, given our current scientific
knowledge and our available or even foreseeable technological capabilities, to create a numerical model able to directly simulate the climate system in all details for all the relevant time scales, which span a range of
over 15 orders of magnitude <cit.>. Hence, one has to focus on a specific range of scales through suitably
developed, approximate evolution equations that provide the basis for the
numerical  modeling. Such equations are derived from the fundamental laws
of climate dynamics through systematic asymptotic expansions that are based on imposing an approximate balance between the forces acting on geophysical flows. These balance relations lead to removing small-scale, fast processes that are assumed to play a minor role at the scales of interest by filtering out the corresponding waves  <cit.>.

In  climate science, parametrization schemes have been traditionally formulated in such a way that one expresses the net impact on the scales of interest of processes occurring within the unresolved scales via deterministic functions of the resolved variables, as in the pioneering work on the parametrization of convective processes by Arakawa and Schubert <cit.>. More recently, it has been recognized, mostly on empirical grounds, that parametrizations should involve  stochastic and non-Markovian components <cit.>.  Machine learning methods have been proposed as the next frontier of data-driven parametrizations <cit.>, able to deliver a new generation of Earth system models <cit.>; see, though, the caveats discussed by <cit.>.

§.§ The Projection Operator Formalism of Mori and Zwanzig

Let us reformulate the problem of constructing parameterizations as a projection of the dynamics onto the set of resolved variables. By working at
the level of observables, Mori <cit.> and Zwanzig <cit.> showed that the evolution laws for the projected dynamics incorporate a deterministic term that would be obtained by neglecting altogether the impact of the unresolved variables, to which a stochastic and a non-Markovian correction had to be added. A. J. Chorin and coauthors <cit.> played an important role in
developing further these ideas and applying them to several important problems. We briefly recapitulate below the Mori-Zwanzig projection operator

Formally, let $\Phi$ denote a generic observable defined on a state space
viewed as the product of two finite-dimensional spaces $\mathcal{X}\times
\mathcal{Y}$, with variables $\xx \in \cX$ and $\yy \in \cY$ being the resolved and unresolved variables, respectively. Next, let us define $\mathbb{P}$ to be a projector onto functions depending only on the target variables in $\mathcal{X}$, with the complementary projector on the unresolved variables being defined by $\mathbb{Q}=\mbox{Id} - \mathbb{P}$.

Given a smooth flow $\psi_t$ on $\mathcal{X}\times \mathcal{Y}$, we consider its action on smooth observables $\Phi = \Phi(\xx, \yy)$ defined by

U_t Φ(,)=Φ(ψ_t(,)),   (, ) ∈×,   t≥0.
The operator $U_t$ is the Koopman operator and the family $\{U_t\}_{t\geq
0}$ of Koopman operators indexed by time forms a semigroup; see Sec. <ref> below for further details. The Koopman operator describes how functions on the phase space change under the action of the flow; its time evolution obeys the following equation:
∂_t (U_tΦ)=ℒ (U_t Φ).
The linear operator $\mathcal{L}$ gives the instantaneous rate of change of $\Phi$ under the action of the flow $\psi_t$. Representing $U_t$ formally as the exponential of $\LLL$, $U_t=e^{t\LLL}$, is both useful and fully justified by operator semigroup theory <cit.>. With this representation, $\LLL$ satisfies the following identities

\begin{align}
\partial _t \big(U_t\Phi\big) & = \mathcal{L} e^{t\LLL}\Phi = e^{t\LLL}\mathcal{L}\Phi = e^{t\LLL}\left(\mathbb{P}+\mathbb{Q}\right)\mathcal{L}\Phi   \\
&= e^{t\LLL}\mathbb{P}\LLL\Phi + e^{t\mathbb{Q}\LLL}\mathbb{Q}\LLL\Phi
+ \int_{0}^{t}e^{(t-s)\LLL}\mathbb{P}\LLL e^{s\mathbb{Q}\LLL}\mathbb{Q}\LLL \Phi  \dd s, \label{mori zwanzig 2}
\end{align}

where we have employed the Dyson identity <cit.>
to obtain Eq. (<ref>), as will be discussed in Sect. <ref>
The first term in Eq. (<ref>) is the contribution of the resolved variables $\xx$ alone to the instantaneous rate of change of $\Phi$. The second term models the fluctuating effects of the unresolved $\yy$-variable by itself, while the third and last term represents, via an integral, the time-delayed influence upon $\xx$ of its interactions with $\yy$.

This formal calculation suggests that any closed model for the $\xx$-variables should incorporate a fluctuating term to account for the $\yy$ contributions and a memory or integral term for the $\yy$–$\xx$ interactions. Unfortunately, the Mori-Zwanzig equation (<ref>)— also known as the Generalized Langevin Equation (GLE) <cit.> — does not provide explicit analytic formulas to
determine each of the three summands in the right-hand side (RHS) of Eq. (<ref>). Hence, we need efficient ways to approximate such
an equation.

In the limit of perfect time scale separation between the $\xx$- and  $\yy$-variables, the non-Markovian term drops out and the fluctuating term can be represented as a — possibly multiplicative — white-noise term, thus recovering the basic results obtained via homogenization theory <cit.>. When no such separation exists, however, one has to resort to finding an integral kernel beyond the
abstract formulation of Eq. (<ref>); see, for instance, the theoretical ansatz based on perturbation expansion presented in refs. <cit.> and discussed later in the paper, and  refs. <cit.> for concrete applications.

In parallel with the theoretical approaches to approximate the Mori-Zwanzig equation (<ref>), data-driven methods have been proposed to model fluctuations and memory effects arising as a result of projecting a large state space onto (much) smaller subspaces. To this end, the work in <cit.> provides a rigorous connection between the Mori-Zwanzig equation (<ref>) and multilevel regression
models that were initially introduced for climatological purposes in <cit.>. More recently, ref. <cit.> proposed Nonlinear Auto-Regressive Moving Average with eXogenous input (NARMAX) models as a data-driven methodology that is comparable with the Mori-Zwanzig formalism and applied such models to the deterministic Kuramoto-Sivashinsky equation and to a stochastic Burgers equation. The complementarity of theory-based and data-driven model reduction methods in the absence of time scale separation is very well documented in <cit.> as well. A highly complementary recent contribution aimed at finding a common thread between data-driven methods and the Mori-Zwanzig theoretical framework is <cit.>.

Efforts at using ingeniously selected basis functions as a stepping stone in data-driven methods for model reduction, effective simulation with partial data, and even prediction are multiple and flourishing. Thus,

the eigenvalues and eigenfunctions of the Koopman and transfer operators have been used to capture the modes of variability of the underlying flow, regardless of the latter being deterministic or stochastic; see, respectively, <cit.> or <cit.>.
Dynamic mode decomposition <cit.> allows one to reconstruct from observations the eigenvectors and eigenvalues of the Koopman operator for observables of interest even in high-dimensional dynamical systems  <cit.>. The latter approach is complementary to the one presented herein because we shall use the  eigenvectors of the Koopman operator to build the projected dynamics for the observables of interest, which can then be re-written as a multilevel Markovian
stochastic model.

Examples of other types of selection of dynamically interesting and effective bases are multichannel singular-spectrum analysis analysis <cit.> and data-adaptive harmonic decomposition <cit.>.
Both methodologies have been used extensively and successfully in the simulation, as well as the prediction, of complex phenomena <cit.>.

Our main theoretical result, Theorem <ref> below, establishes how such spectral elements determine in turn the constitutive elements of data-driven non-Markovian closure of partially observed complex systems, when rewritten as a multilevel Markovian stochastic model. This theorem highlights, in particular, new bridges with Koopman modes and the DMD
decomposition <cit.>, as well with other kinds of projections onto spectral bases <cit.>. A more thorough discussion of the complementary approaches involved is beyond
the scope of this paper.

§.§ Goals of This Paper

Many approaches to constructing theoretically rigorous parametrizations have been devised. These can be broadly divided into top-down and data-driven approaches: Top-down approaches aim at deriving the parametrizations by applying suitable approximations to the equations describing the dynamics of the whole systems <cit.>, while data-driven parametrizations are built by constructing a statistical-dynamical model of the
impact of unresolved scales on the scales of interest. In fact, partial observations of a time-evolving system can be used to deduce the fluctuating and delayed effects of the unobserved processes, as shown in <cit.>.

In this paper, we will discuss and compare the properties of the Wouters-Lucarini (WL) top-down parametrization <cit.> and of the Empirical Model Reduction (EMR) data-driven parametrization <cit.>. We will also see when and how the integro-differential equation occurring in the WL parametrization can be recast into a set of Markovian stochastic differential equations (SDEs). In other words, we investigate the quasi-Markovianity of the latter parametrization <cit.>.

The two aforementioned methodologies are conceptually and practically distinct, even though the ultimate goal is to provide a  computationally practical approximation for the Mori-Zwanzig or GLE integro-differential equation. In other words, both approaches — top-down and data-driven — provide fluctuations in the form of stochastic noise and memory effects determined by an integral kernel. On the one hand, the WL approach assumes prior knowledge about the decoupled hidden dynamics but no information about the statistical properties of the coupled system. The empirical approach, on the other hand, samples the observed variables evolving according
to the latter. The structure of the multilevel stochastic models (MSMs) that generalize EMRs <cit.> allows one, moreover,
to derive explicit formulas for the fluctuating and memory correction terms that parametrize the influence of hidden processes.

The overall goal of this paper is to provide a conceptual and analytical link between these two approaches, aiming, on the one hand, to buttress the practical relevance of the WL perturbative approach and, on the other,
to provide further insight into the well-documented robustness of the EMR
method. Moreover, we will clarify how multilevel systems arise from both the top-down (WL) and the bottom-up (EMR) approaches. The paper explores the complete set of boxes and explains all the arrows in Figure <ref>. The diagram in the figure shows that, starting from
the top box, one can arrive at a memory equation, like (<ref>), via top-down or bottom-up methods, as indicated by the left and right sequence of arrows, respectively.

The paper is structured as follows. Section <ref> revisits the derivation of the WL parametrization method by applying the  Dyson expansion to the Koopman operator associated with two
weakly coupled dynamical systems. We show that such an expansion need not
be truncated for additively coupled models and consider more general coupling laws than those in <cit.>.

Furthermore, we study the problem of finding Markovian representations of the memory equation in the WL approach based on the spectral decomposition of the Koopman semigroup. Specifically, Theorem <ref> shows, in the case of a scalar observable, how to recast the stochastic integro-differential equation arising in the WL parametrization as a multilevel Markovian stochastic system involving explicitly the spectral elements of the (uncoupled) Koopman operator, and we point out in Remark <ref> how such a Markovianization extends to the multidimensional setting. Section <ref> provides new insights into the Markovian representation adopted in the MSM framework; these insights help one to determine, in particular, the number of levels required for EMR to converge.

Section <ref> presents a comparison of the data-driven and top-down parametrization approaches using a simple conceptual stochastic climate model. Finally, we discuss the conclusions obtained from this investigation in Section <ref>.

Appendices are included in order to avoid diverting the attention of the reader from the main message of the paper. Appendix <ref> provides the proof of the key Theorem <ref>. Appendix <ref>  briefly revisits the spectral decomposition of correlation functions and provides a criterion to quantify the loss of the semigroup property that is useful in identifying non-Markovian effects. Appendix <ref> discusses the stochastic Itô integration of the elementary form of
an MSM. Appendix <ref> shows how EMR approaches can capture the dynamical properties of partially observed systems; in it, we consider a simple climate model, obtained by coupling the Lorenz atmospheric model <cit.> and the Lorenz convection model <cit.>.

Schematic view of the two complementary approaches
studied in this paper.
The arrows on the left-hand side indicate top-down, perturbative parametrizations;
on the right, they refer to bottom-up, empirical parametrizations.

§ REVISITING THE WEAK-COUPLING–LIMIT PARAMETRIZATION

To study dynamical systems in which one can separate the variables into two groups, with weak coupling between the two, one often resorts to so-called parametrizations of the effects of one group on the other. In the weak-coupling limit, the coupling itself can be treated as a perturbation of the main dynamics <cit.>. Granted such an assumption and some degree of structural stability of the system, one can apply response theory to derive explicit stochastic and memory terms to describe the impact of the variables we want to neglect on the variables of interest, in the Mori-Zwanzig spirit. Note that, to do so, no assumption on time-scale separation between the two groups of variables is necessary. This point is particularly relevant in fields like the climate sciences, where no clear time-scale separation is observed, so that asymptotic expansions of the kind used in homogenization theory are of
limited utility.

§.§ Deriving the WL Approximation for the Mori-Zwanzig Formalism

Here, using a perturbative approach, we review the derivation of the parametrization presented in refs. <cit.>. Formally, we want to couple two dynamical systems generated independently by two vector fields $\FFF : \mathcal{X}\subseteq \mathbb{R}^{d_1} \longrightarrow \mathcal{X}$ and $\GGG : \mathcal{Y}\subseteq \mathbb{R}^{d_2}\longrightarrow \mathcal{Y}$ with possibly $d_1 \neq d_2$ and typically $d_1\ll d_2$.
We study a broad class of systems of the form:

\begin{align}
\dot{\mathbf{x}}(t) &= \FFF(\mathbf{x}(t)) + \epsilon \CCC_{\xx}^{\xx}\left(\xx(t) \right) : \CCC_{\yy}^{\xx}\left(\yy(t) \right),
\label{mle1}\\
\dot{\mathbf{y}}(t) &= \GGG(\mathbf{y}(t)) + \epsilon \CCC_{\xx}^{\yy}\left(\xx(t) \right):\CCC_{\yy}^{\yy}\left(\yy(t) \right).
\label{mle2}
\end{align}

The operation indicated by the colon $\xx:\yy$ denotes the Hadamard product that multiplies vectors or matrices component-wise. Here, four new vector fields have been introduced to model the coupling law, namely $\CCC_{\xx}^{\xx}:\mathcal{X} \longrightarrow \mathcal{X}$, $\CCC_{\xx}^{\yy}: \mathcal{X} \longrightarrow \mathcal{Y}$, $\CCC_{\yy}^{\xx}:\mathcal{Y} \longrightarrow \mathcal{X}$ and $\CCC_{\yy}^{\yy}:\mathcal{Y} \longrightarrow \mathcal{Y}$.

The real parameter $\epsilon$  controls the strength of the coupling between the two groups of variables, $\xx(t)$ and $\yy(t)$, so that the $\xx$-  and $\yy$-variables are uncoupled for $\epsilon = 0$. We assume that
the vector fields $ \FFF$ and $\GGG$, as well as the coupling laws
in Eqs. (<ref>,<ref>) are such that system (<ref>) possesses a global attractor. Furthermore, we assume throughout this article
that this global attractor supports an physical invariant probability measure $\mu$ that describes the distribution of trajectories onto the global attractor.

The WL parametrization views the coupling as an $\epsilon$-perturbation of the otherwise independent $\xx$- and $\yy$-processes, with $\xx$ the observed and $\yy$ the hidden variables. One next assumes that the impacts of perturbations applied to these processes can be addressed using response theory <cit.>, so that response formulas can be used to derive an effective equation for the $\xx$-variables.

Taking the Mori-Zwanzig <cit.> point of view, we wish to calculate the evolution of observables that depend on the observed variables $\xx$ alone, $\Phi = \Phi (\xx(t))$. The idea, following <cit.>, is to perform a perturbative expansion
of the differential operator $\LLL$ governing the evolution of $\Phi(\xx(t))$ under the action of the flow associated with Eq. (<ref>). Denoting by $u(\xx,\yy,t)$ the time evolution of a smooth observable $\Phi$ in $\mathcal{C}^\infty(\cX\times \cY)$,  the first step of this
Dyson-like operator
expansion reads as follows:

∂_t u = u =(_0 +ϵ_1) u.
Here $ \LLL_0$ and $\LLL_1$ account for the advective effects of the uncoupled and coupling terms, respectively, that compose the RHS of Eq. (<ref>), namely

\begin{align}
\LLL_0 &= \begin{bmatrix}
\FFF (\xx) \\ \GGG (\yy)
\end{bmatrix}\cdot \begin{bmatrix}
\nabla_{\xx} \\ \nabla_{\yy}
\end{bmatrix}, \label{WL_LLL_0} \\
\LLL_1 &=  \begin{bmatrix}
\CCC_{\xx}^{\xx}\left(\xx \right) : \CCC_{\yy}^{\xx}\left(\yy \right) \\ \CCC_{\xx}^{\xx}\left(\xx \right) : \CCC_{\yy}^{\xx}\left(\yy \right)
\end{bmatrix}\cdot \begin{bmatrix}
\nabla_{\xx} \\ \nabla_{\yy}
\end{bmatrix};\label{WL_LLL_1}
\end{align}

in Eqs. (<ref>), $\nabla_{\xx}$ and $\nabla_{\yy}$ denote the vector differential operators with respect to the variables $\xx$ and $\yy$.

Recalling Eq. (<ref>), the solution operator of Eq. (<ref>) is the Koopman operator. Formally, its dual acts on densities and it is the so-called transfer operator <cit.>. Equation (<ref>) is thus a transport equation, where the physical quantity or observable is advected by the vector field on the RHS of Eq. (<ref>).

Note that the operator formalism presented here in the deterministic dynamical systems setting — and the associated semigroup theory — extends
to Markov diffusion processes driven by a stochastic forcing <cit.>.  In the latter case, the transport equation (<ref>) becomes the so-called backward Kolmogorov equation that describes
the evolution of the expected value of observables. Loosely speaking, the
corresponding extension amounts to adding a Laplacian-like operator to the advection operator $\LLL$ <cit.>. See <cit.> for the appropriate context in testing the applicability of response theory to the independent  $\xx$  and $\yy$ processes, when

More precisely, one associates to the solution $u(\xx,\yy,t)$ of Eq. (<ref>) unfolding from an observable $\Phi = \Phi(\xx,\yy)$ at time $t=0$, a family of linear Koopman operators indexed by time
$\{U_t\}_{t\geq 0}$ such that $u(\xx,\yy,t)=U_t \Phi (\xx, \yy)$, for any $t\geq 0$ and $(\xx, \yy)$ in $\cX \times \cY$. These operators are defined — as mentioned already in connection with introducing the GLE (<ref>) — as exponentials of the operator $\LLL$, i.e., $U_t=e^{t\LLL}$. This notation is formal as the operator $\LLL$ is unbounded; it is, however, useable as $\{U_t\}_{t\geq 0}$ satisfies  the semigroup property, i.e. $U_{t+s}=U_t U_s$, $t,s\geq 0$, as for a standard exponential. Over the appropriate function space[ Such a space can be chosen, for instance, as $D_p=\{ \Phi \in L^{p}_{\mu} (\cX \times \cY) \, | \, A\Phi=\lim_{t\rightarrow 0} \ t^{-1}(U_t \Phi -\Phi) \textrm{ exists}\}$ for some $ p \in [1,\infty]$, with $\mu$ denoting a relevant
invariant measure of the system (<ref>) supported by the global attractor, while the limit is taken in the sense of strong convergence <cit.>.] of observables $\Phi$, this family actually forms
a strongly continuous contracting semigroup <cit.>.

The action of the flow on an observable $\Phi$ becomes thus more transparent thanks to the operator $U_t$, according to the equation

\begin{equation}\label{Def_Koopman}
U_t \Phi(\xx_0, \yy_0) = e^{t\LLL}\Phi (\xx_0, \yy_0) = \Phi (\xx(t;\xx_0),\yy(t;\yy_0)),
\end{equation}

where $(\xx(t;\xx_0),\yy(t;\yy_0))$ denotes the system's solution at time
$t$ emanating from the initial state $(\xx_0, \yy_0)$ at time $t=0$.  In what follows, we omit the subscript $0$ in $(\xx_0, \yy_0)$ but still take it as an initial state.

The semigroup $\{U_t\}_{t\geq 0}$ is known as the Koopman semigroup and for each $t$, $U_t$ is the Koopman operator mentioned above; see also <cit.>. When the coupling parameter $\epsilon$ in system (<ref>) is small, one can use formal perturbation expansions of the Koopman semigroup to better isolate and assess the coupling effects at the level of observables. To do so, we follow here the perturbation expansion first introduced by Freeman J. Dyson in the context of quantum electrodynamics <cit.> and later formulated rigorously in mathematical terms in <cit.>. Formally, this expansion reads as follows:

\begin{align}
%\Phi(\xx(t),\yy(t)) &=
%%%%%MC--->VC & MSG: useless to recall the relation between  \Phi(\xx(t),\yy(t))  and e^{t\LLL}\Phi(\xx,\yy) here
U_t\Phi (\xx,\yy) = e^{t\LLL}\Phi(\xx,\yy) & = e^{t\LLL_0+t\epsilon \LLL_1}  \Phi(\xx,\yy) \\ \label{pertexp1} &= e^{t\LLL_0}\Phi(\xx,\yy)+\epsilon\int_{0}^te^{s\LLL}\LLL_1e^{(t-s)\LLL_0}\Phi(\xx,\yy)\dd s\\ \label{pertexp2}
&= e^{t\LLL_0}\Phi(\xx,\yy)+\epsilon\int_{0}^te^{(t-s)\LLL_0}\LLL_1e^{s\LLL}\Phi(\xx,\yy)\dd s,
\end{align}

and it yields the following expansion of the Koopman operator in $\epsilon$:

U_tΦ(,) = e^t_0Φ(,)+ϵ∫_0^te^(t-s)_0_1e^s_0Φ(,)s + 𝒪(ϵ^2).
This identity shows that the evolution of a generic observable can be described as an  $\epsilon$-perturbation of its decoupled evolution according to ${\LLL_0}$.
We note that these expansions are purely formal and in particular it is not  clear in which sense this expansion might converge. For a bounded perturbation operator $\LLL_1$, it would be straightforward to prove boundedness of the resulting perturbed semigroup. However, $\LLL_1$ here is a differential linear operator, for which direct estimates are more laborious. Leaving alone the functional analysis framework that would make such an expansion rigorously convergent, we shall use nevertheless the expansion (<ref>) throughout this article.

The objective now is, using this operator expansion, to derive an effective reduced-order model for the evolution of the $\xx$-variable without having to resolve the $\yy$-process. We start observing the system at $t=0$, but assume that it has already attained a steady state. Since we are only concerned with observables depending solely on the $\xx$-variables, we formulate now an evolution equation for such observables. To do so, we
consider first the Liouville equation for a generic $\yy$-independent observable $\Phi$, this is, $\Phi(\xx,\yy)=\Phi(\xx)$, for every $\xx$ and

For such an observable, at the time we start observing the coupled system, Eq. (<ref>) reduces to

∂_t (U_t Φ) |_t=0=[ () + ϵ_^():_^()  ]·∇_Φ,
where $\cdot$ denotes an inner product in $\mathcal{X}$. Equation (<ref>) illustrates the trivial fact that the time evolution in Eq. (<ref>) of an $\xx$-dependent physical quantity is also affected by the $\yy$-variables.

Following <cit.>, the decoupled equations are assumed to have been evolving for some time prior to the coupling. Hence, we
have to formally parametrize the evolution of the $\CCC_{\yy}^{\xx}\left(\yy \right)$-contribution to the vector field  which is, ultimately, a vector-valued observable.

We do so by introducing an extended version of the Koopman operators that
act on vectors component-wise, rather than just on real-valued observables. Consider $\mathbf{v}: \mathcal{X}\times \mathcal{Y}\longrightarrow \mathbb{R}^d$, for some positive integer $d$, and define the action of the Koopman operator $e^{t\LLL}$ on $\mathbf{v}$ as:

[e^t𝐯(,)]_i = e^t[𝐯(,)]_i
for every $i=1,\ldots,d$. The definition (<ref>) will allow us to use the semigroup notation for observables of possibly different dimensions, all of which take their inputs in the phase space $\mathcal{X} \times \mathcal{Y}$ . Ultimately, this is a component-wise evaluation of our extended Koopman operator family, and its generator can be obtained analogously. As mentioned above, we have to model the effects of the coupling vector field $\CCC _{\yy}^{\xx}(\yy)$, whose state at time $t=0$ is the
product of the evolution from time $-t$ to $0$. We then have, with the dynamics starting at time $-t$ and initial state $(\xx_0,\yy_0)$,

\begin{equation}
\CCC_{\yy}^{\xx}\left(\yy \right) = e^{t\LLL}\CCC_{\yy}^{\xx}\left(\yy_0,-t \right)=e^{t\LLL}\CCC_{\yy}^{\xx}\left(\yy_0 \right).
\end{equation}

Now, by using the perturbative expansions in Eqs. (<ref>)-(<ref>), we obtain:

\begin{align}
\CCC_{\yy}^{\xx}\left(\yy \right) &= e^{t\LLL}\CCC_{\yy}^{\xx}\left(\yy_0 \right)=e^{t\LLL_0+t\epsilon \LLL_1}	\CCC_{\yy}^{\xx}\left(\yy_0 \right) \\ \label{exact expansion1}
&= e^{t\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)+\epsilon\int_{0}^te^{s\LLL}\LLL_1e^{(t-s)\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\dd s\\ \label{exact expansion2}
&= e^{t\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)+\epsilon\int_{0}^te^{(t-s)\LLL_0}\LLL_1e^{s\LLL}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\dd s \\\label{lw approx}
&= e^{t\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)+\epsilon\int_{0}^te^{(t-s)\LLL_0}\LLL_1e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\dd s + \mathcal{O}\left(\epsilon^2\right).
\end{align}

Plugging the indentity in Eq. (<ref>) into Eq. (<ref>, we find the following expression:

\begin{equation}\label{time zero liouville}
\partial _t \big(U_t \Phi\big) |_{t=0} = \left[ \FFF(\xx) + \epsilon \CCC_{\xx}^{\xx}\left(\xx \right):\Big\{ e^{t\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)+\epsilon\int_{0}^te^{s\LLL}\LLL_1e^{(t-s)\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\dd s \Big\}  \right]\cdot \nabla_{\xx} \Phi.
\end{equation}

This equation is an exact reformulation of the problem induced by Eq. (<ref>). This reformulation demonstrates that memory effects enter at second order in powers of the coupling parameter. Notice, though, that even if Eq. (<ref>) reduces the dimensionality
of the problem from $d_1 + d_2$ to $d_1$, it does not constitute an approximation for the evolution of $\Phi$ as an observable of $\xx$ alone, since it depends on the evolution of the $\yy$-variables in the coupled regime by means of the action of $e^{s\LLL}$ onto $\LLL_1$.

Therefore, we need to perform a further approximation by considering Eq. (<ref>) instead, which leads to:

\begin{equation}\label{dysonapprox}
\partial _t \big(U_t \Phi\big) |_{t=0} \simeq \left[ \FFF(\xx) + \epsilon \CCC_{\xx}^{\xx}\left(\xx \right):\Big\{ e^{t\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)+\epsilon\int_{0}^te^{(t-s)\LLL_0}\LLL_1e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\dd s \Big\}  \right]\cdot \nabla_{\xx} \Phi ,
\end{equation}

where the terms of order $\epsilon ^3$ have been dropped. Equation (<ref>)  is our equivalent of the Dyson approximation in quantum electrodynamics; it approximates the evolution of the $\xx$-variables with no need for the evolution of the $\yy$-variables in the coupled regime.

This result amounts to saying that — by observing only the statistical properties of the decoupled dynamics of the $\yy$-process, obtained with $\epsilon=0$ — one can
construct a Markovian contribution
_^():{e^t_0_^(_0 )},
and a non-Markovian contribution
_^():{ ∫_0^te^(t-s)_0_1e^s_0_^(_0 )s },
to the dynamics of the $\xx$-variables that is able to describe the impact the the coupling.

Expanding the kernel $\tilde{\mathcal{K}}$ of the memory contribution, we

\begin{align}
\tilde{\mathcal{K}}(t,s,\xx_0,\yy_0):&=e^{(t-s)\LLL_0}\LLL_1e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right) \\ \label{kernel2}
&= e^{(t-s)\LLL_0}\left(\left[ \CCC_{\xx}^{\xx}(\xx_0):\CCC_{\yy}^{\xx}(\yy_0) \right]\cdot \nabla_{\xx} + \left[\CCC_{\xx}^{\yy}(\xx_0):\CCC_{\yy}^{\yy}(\yy_0)  \right]\cdot \nabla_{\yy}\right)e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\\
&= \left[e^{(t-s)\LLL_0} \left(\CCC_{\xx}^{\yy}(\xx_0):\CCC_{\yy}^{\yy}(\yy_0)  \right)\right]\cdot \nabla_{\yy}e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right).
\end{align}

Note that the leading-order Koopman operator $e^{s\LLL_0}$ models the evolution of the observables in the uncoupled regime. Since there is no prior knowledge on initializing the coupled system at time $-t$, the initial state $\yy_0$ in the hidden variables should be drawn from an ensemble, according to a probability density function. At this stage, there is freedom in the choice of such a prior. However, since we are assuming that the
coupled system was initialized at time $-t$, it is natural to draw $\yy_0$ according to the invariant measure $\nu$ associated with the dynamical system generated by the vector field $\GGG$ from Eq. (<ref>).

We wish to sample initial conditions from the coupled steady state, but do not assume any prior knowledge of the coupled statistics. As discussed in <cit.>, we can take advantage of response theory to address this situation. Indeed, for any sufficiently smooth observable $\Psi$, we have:

$$ 	\langle \Psi \rangle_{\epsilon} = \langle  \Psi \rangle_{\epsilon=0} + \sum_{k=1}^{\infty}\epsilon ^{k}\delta_k [\Psi],$$

where $\langle \Psi \rangle_{\epsilon}$ is the expectation value of $\Psi$ in the coupled system (<ref>),
$\langle \Psi \rangle_{\epsilon=0}$ is the expectation value of $\Psi$ according to the statistics generated by the uncoupled $\yy$ process obtained by setting $\epsilon=0$ in Eq. (<ref>), and $\epsilon ^{k}\delta_k [\Psi]$ is the $k^{\mathrm{th}}$-order response. In what follows, we remove the subscripts for the averages when $\epsilon = 0$. Therefore, we have that the expected value of the coupling function reads as:

\begin{equation}
\bigg\langle  \CCC_{\yy}^{\xx} \bigg \rangle_{\epsilon} = \bigg\langle
\CCC_{\yy}^{\xx}\bigg \rangle + \sum_{k=1}^{\infty}\epsilon ^{k}\delta_k [\CCC^{\xx}_{\yy}].
\end{equation}

Likewise, we can calculate the average of such function at time $t$:

\begin{equation}
\bigg\langle  e^{t\LLL_0}\CCC_{\yy}^{\xx} \bigg \rangle_{\epsilon} = \bigg\langle  e^{t\LLL_0}\CCC_{\yy}^{\xx} \bigg \rangle + \sum_{k=1}^{\infty}\epsilon ^{k}\delta_k [e^{t\LLL_0}\CCC^{\xx}_{\yy}].
\end{equation}

Now, by letting $\tilde{\eta}(t,\yy_0)= e^{t\LLL_0}\CCC_{\yy}^{\xx}(\yy_0)$, we find that in order for the approximate statistics to agree up to
second order in $\epsilon$ with the exact one, we only need to impose the
following conditions upon the first two moments of the parametrized noisy
fluctuations (see also <cit.>):

\begin{align}
\bigg \langle \tilde{\eta}(t,\yy_0) \bigg \rangle &= \int \nu\left(\dd
\yy_0\right)\CCC_{\yy}^{\xx}\left(\yy_0 \right), \label{moments1} \\
\bigg \langle \tilde{\eta}(t,\yy_0)\tilde{\eta}^{\top}(0,\yy_0) \bigg \rangle &= \int \nu(\dd \yy_0) e^{t\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)\left(\CCC_{\yy}^{\xx}\left(\yy_0 \right)\right)^{\top}; \label{moments2}
\end{align}

here $(\cdot)^{\top}$ stands for the transpose of a vector or a matrix. It follows that any stochastic noise $\eta(t)$ that satisfies the two conditions above will be suitable for parametrizing the fluctuations in the $\yy$-dynamics tied to the lack of knowledge in the initial state. Each of
the entries in the correlation matrix given by Eq. (<ref>) is the correlation function between the components of the vector field $\CCC^{\xx}_{\yy}$ and these will become explicit provided a suitable spectral decomposition is at hand. Such a decomposition will be provided later in Sect. <ref>, although the reader is referred at this point to Appedix <ref> for clarity.

In the memory term, though, we neglect $\epsilon$-corrections to its statistics since memory effects are of order $\epsilon ^2$ already. Thus, we have, by averaging the kernel $\tilde{\mathcal{K}}(t,s,\xx_0,\yy_0)$ in Eq. (<ref>) with respect to the $\yy$-variables,

\begin{align}
\mathcal{K}(t,s,\xx):&=\int \nu(\dd \yy_0) \left[e^{(t-s)\LLL_0} \left(\CCC_{\xx}^{\yy}(\xx_0):\CCC_{\yy}^{\yy}(\yy_0)  \right)\right]\cdot \nabla_{\yy}e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right) \\
& = \int \nu(\dd \yy_0)  \left[e^{(t-s)\LLL_0} \CCC_{\xx}^{\yy}(\xx_0): e^{(t-s)\LLL_0}\CCC_{\yy}^{\yy}(\yy_0)\right]  \cdot \nabla_{\yy}e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right) \\
& = \int \nu(\dd \yy_0)\left[ \CCC_{\xx}^{\yy}(\xx(t-s)): e^{(t-s)\LLL_0}\CCC_{\yy}^{\yy}(\yy_0)\right]  \cdot \nabla_{\yy}e^{s\LLL_0}\CCC_{\yy}^{\xx}\left(\yy_0 \right)
\end{align}

This way, the memory kernel only depends on the $\xx$-variables. Hence, we find a self-consistent evolution of the $\xx$-variables, subject to the
influence of unobserved variables $\yy$, in the form of a stochastic integro-differential equation (SIDE) resembling the
GLE (<ref>):

\begin{equation}\label{lw parametrization}
\dot{\xx}(t)=\FFF(\xx) + \epsilon \CCC_{\xx}^{\xx}\left(\xx \right): \eta(t)+\epsilon^2\CCC_{\xx}^{\xx}(\xx)\cdot \int_{0}^t\mathcal{K}(t,s,\xx)\dd s,
\end{equation}

where $\eta(t)$ is a stochastic forcing that agrees with the mean and correlation properties stated in Eq. (<ref>).

We emphasize that the solution $\xx(t)$ of the original system of ordinary differential equations (<ref>) does not satisfy Eq. (<ref>): it is just the proposed reduced-order model for the $\xx$-variables. The closure provided by expressing the corrections in the second and third term on the RHS of Eq. (<ref>) as functions of $\xx$ alone is typically called a parametrization of the effect of the unobserved $\yy$-variables in the climate sciences <cit.>.

Note that there is considerable freedom in the choosing the noise, since we only require agreement up to the second moment. However, a direct consequence of this weak-coupling parametrization is that realizations of the
noise can be produced by directly integrating the decoupled hidden variables, or by representing it using simple autoregressive models <cit.>. We are assuming here that the uncoupled dynamics leads to a noisy signal; this can be achieved either by the presence of stochastic forcing in the hidden variables <cit.> or by their uncoupled dynamics being chaotic <cit.>.

To summarize, the weak-coupling limit allows one to develop a parametrization of the hidden variables for a system of coupled equations where no separation of time scales is assumed. Moreover, this approach provides explicit approximate expressions for the deterministic, stochastic, and non-Markovian terms in the Mori-Zwanzig formalism of Eq. (<ref>).

There are two sources of error in the parametrization proposed in Eq. (<ref>). First, the truncation performed in the Dyson expansion neglects higher-order effects, which are weighted by the third power of the small coupling parameter. Secondly, averaging over the statistics of the uncoupled dynamics can also introduce errors. Furthermore, the nature of the stochastic correction is not fully determined except for its lagged correlation.

The perturbation operator approach taken here is analogous to that of <cit.>, who only considered the independent or additive-coupling
cases; the latter is expanded upon in Sect. <ref>. Here, though,
we generalize further the parametrization formulas that can be obtained via perturbative expansion of linear operators. In fact, the present approach can also be extended to weakly coupled systems of the form:

\begin{align}
\label{mle12}
\dot{\mathbf{x}}(t) &= \FFF(\mathbf{x}(t)) + \epsilon \CCC^{\xx}\left(\yy(t) \right), \\
\label{mle22}
\dot{\mathbf{y}}(t) &= \GGG(\mathbf{y}(t)) + \epsilon \CCC^{\yy}\left(\xx(t),\yy(t) \right),
\end{align}

where $\CCC^{\yy}$ encodes interactions that need not be separable between the $\xx$- and $\yy$-variables in the hidden layer of the model. Note that the full parametrization of arbitrary couplings was discussed by the two authors of <cit.> in previous work <cit.>, in
which they used a response-theoretic approach.

§.§ The Additive-Coupling Case

The approximate Dyson expansion given in Eq. (<ref>) is exact
in the case of additive coupling. Such systems take the form:

\begin{align}
\label{add_mle12}
\dot{\mathbf{x}}(t) &= \FFF(\mathbf{x}(t)) + \epsilon \CCC^{\xx}\left(\yy(t) \right), \\
\label{add_mle22b}
\dot{\mathbf{y}}(t) &= \GGG(\mathbf{y}(t)) + \epsilon \CCC^{\yy}\left(\xx(t) \right).
\end{align}

Indeed, letting $\CCC ^{\yy}(\xx,\yy)=\CCC ^{\yy}(\xx)$ in Eq. (<ref>) and using Eq. (<ref>) allows us to avoid the truncation of the Dyson expansion and yields the following expression for the memory term:

\begin{align}
\tilde{\mathcal{K}}(t,s,\xx,\yy_0)&=e^{s\LLL}\LLL_1e^{(t-s)\LLL_0}\CCC^{\xx}\left(\yy_0 \right) \\ &= e^{s\LLL}\left(\CCC^{\xx}(\yy_0) \cdot \nabla_{\xx} + \CCC^{\yy}(\xx) \cdot \nabla_{\yy}\right)e^{(t-s)\LLL_0} \CCC^{\xx}\left(\yy_0 \right)\\ &= e^{s\LLL} \left(\CCC^{\yy}(\xx) \right)\cdot \nabla_{\yy}e^{(t-s)\LLL_0} \CCC^{\xx}\left(\yy_0 \right),
\end{align}

which is exact.
Next, taking averages with respect to $\nu$, we obtain:

\begin{align}\label{wl parametrization additive}
\mathcal{K}(t,s,\xx)&=  \int \nu(\dd \yy_0)e^{s\LLL} \CCC^{\yy}(\xx) \cdot \nabla_{\yy}e^{(t-s)\LLL_0}\CCC^{\xx}\left(\yy_0 \right).
\end{align}

Hence, the parametrization in this additive-coupling case is exact, as no
terms proportional to $\epsilon^k$, $k\geq3$ are present. The only assumption made is that the statistics in the $\yy$ variables have reached a steady state according to the unperturbed system. Finally, the full SIDE in
this case takes the form:

\begin{equation}\label{lw parametrization independent}
\dot{\xx}(t)=\FFF(\xx) + \epsilon \eta(t)+\epsilon^2\int_{0}^t\mathcal{K}(t,s,\xx)\dd s,
\end{equation}

where the stochastic process $\eta$ has the mean and correlation properties given by Eqs. (<ref>).
This equation is, thus, exactly the one obtained in <cit.>.

Memory effects represented by integral terms seem unavoidable unless the memory kernels vanish quickly with respect to time. Infinite time scale separation between the two sets of variables leads, though, to the vanishing of the associated integral expressions <cit.>. Here, we are not assuming no such property in the coupled dynamical system under study; see Eqs. (<ref>) and (<ref>). On the other hand, reduced phase spaces can help explain  the statistics of the dynamical system without resorting to delayed effects that entail the integrals in Eqs. (<ref>) and (<ref>). Following <cit.>, we briefly review in Appendix <ref> a criterion based on Koopman operators — and, more generally, Markov operators — that enables one to decide whether memory effects can help explain
the dynamics and statistics in reduced phase spaces.

§.§ Markovian Representation through Leading Koopman Eigenfunctions

In the context of Langevin dynamics, there are known conditions on memory kernels that allow one to recast certain stochastic integro-differential equations into a Markovian SDE by means of extended variables; see <cit.>. The stochastic processes that allow such a procedure are called quasi-Markovian <cit.>.[Such a Markovianization procedure is actually not limited to stochastic processes and it relies on the same type of ideas in other contexts; see refs. <cit.>
and <cit.>.]

This Markovianization theory can be formulated in the setting of near-equilibrium statistical mechanics, where one uses fluctuation-dissipation–like relations that link the decay properties of the memory kernel and the
decorrelation rates of the  fluctuations. Here, we follow the approach in
<cit.> but without making any assumptions on the Hamiltonian behavior of the $\yy$-variables. We need, though, to make assumptions on the spectral properties of the generator of the $\yy$-dynamics, as explained below.

We define the generator $\LLL^{\yy} _0$ of the Koopman semigroup associated with the $\yy$-dynamics by:

\begin{equation}\label{Def_L0}
\LLL_0^{\yy}\Phi(\yy)=\GGG(\yy)\cdot \nabla\Phi(\yy),
\end{equation}

for every real-valued observable $\Phi \in \mathcal{C}^{\infty}\left(\mathbb{R}^{d_2}\right)$ and we denote the associated Koopman operator at time $t$ by $U_t=e^{\LLL^{\yy} _0 t}$; the subscript $\mathbf{y}$ has been
dropped from the $\nabla$ operator herewith, for notational clarity. Recall that the the spectrum of such operators provides useful insights into the statistical properties of the system; this topic is beyond the scope of the present paper but it is treated in detail in <cit.>.

It suffices to show below that the spectrum of $U_t$ allows one to characterize the constitutive ingredients of the WL parameterization (<ref>) and (<ref>), subject to natural assumptions. Even though we have clarified in Eq. (<ref>) how the Koopman operator acts on vector-valued observables of any dimension, we restrict now its action for simplicity to scalar real-valued observables, as in Eq. (<ref>).  In this case, along the lines of the methodology of dynamic mode decomposition <cit.>, we can (formally) decompose the operator as:

\begin{equation}\label{Eq_L0_decomp}
e^{\LLL^{\yy}_0 t} = \sum _{j=1}^{  N} e^{t\lambda _j}\Pi_j + \mathcal{R}(t),
\end{equation}

where $\{\lambda_j\}_{j=1}^{  N}$ are the eigenvalues that form the point spectrum of $\LLL^{\yy} _0$ and $\Pi _j$ is the spectral projector onto the eigenspace spanned by the eigenfunction $\psi _j$. Here, $\mathcal{R}(t)$ is the residual operator associated with the essential spectrum of
$\LLL^{\yy} _0$ and its norm is controlled by a decaying exponential. We assume furthermore that the spectrum of $\LLL^{\yy} _0$ lies in the complex left half-plane, and that in particular $\mathfrak{Re}\lambda _j \leq 0$ for any $j$.

Such a spectral decomposition and its properties can be rigorously justified for a broad class of differential equations perturbed by small noise disturbances; see <cit.>.
Based on these rigorous results, we assume, roughly speaking,  that these properties survive in a certain small-noise limit, and concentrate here
on vector fields $\yy$ given by $\GGG$ in (<ref>) for which a decomposition such as (<ref>) holds and a spectral gap does exist in the appropriate functional space.

In the following lines, we examine the expression of the memory kernel $\mathcal{K}$ appearing in Eq. (<ref>) using the eigendecomposition proposed in Eq. (<ref>). In particular, we study such an integral kernel $\mathcal{K}$ component-wise:

\begin{align}\label{memory kernel spectral decomposition}
[\mathcal{K}(t,s,\xx)]_{i} &= \CCC^{\yy}(\xx(s))\cdot \bigg \langle \nabla\sum_{j=1}^Ne^{\lambda_j (t-s)}\alpha_{i,j} \psi _j (\yy) \bigg \rangle + \mathcal{R}(t-s)[\CCC ^{\xx}]_i \\&\approx \CCC^{\yy}(\xx(s))\cdot
\bigg \langle \nabla\sum_{j=1}^Ne^{\lambda_j (t-s)}\alpha_{i,j} \psi _j
(\yy)  \bigg \rangle \\ &=\CCC^{\yy}(\xx(s))\cdot  \sum_{j=1}^Ne^{\lambda_j (t-s)}\alpha_{i,j}\bigg \langle \nabla\psi _j (\yy)  \bigg \rangle,
\end{align}

for  $i=1,\ldots,d_1$, where
$$\alpha_{i,j}= \left\langle \psi^{\ast}_j, [\CCC ^{\xx}]_i \right\rangle = \int \nu(\dd \yy)\overline{\psi_j^{\ast}(\yy)}[\CCC ^{\xx}(\yy)]_i,$$
and we have neglected the contribution
coming from the essential spectrum. The $(\cdot)^{\ast}$ superscript is used to denote the dual eigenfunction.

This decomposition highlights the fact that the leading eigenvalues of the operator governing the evolution of observables in the uncoupled $\yy$-dynamics set the time scale for the memory kernel. Furthermore, this spectral approximation implies that the correlation functions of the noise have the same decay properties, as will become apparent later in the proof of Theorem <ref>. It follows that the correspondence between the noise and integral time scales allows us to recast the SIDE in the
WL equation (<ref>) into a fully Markovian
version with linearly driven hidden variables that are forced by the observed variables, through a functional dependence that can be nonlinear. More exactly, we have the following theorem.

Consider the system (<ref>) where Eq. (<ref>) is, instead, a scalar equation for a real-valued variable $x(t)$.
Let $\nu$ be the physical invariant measure associated with the equation

= 𝐆(),
i.e., with the flow determined by the vector field $\mathbf{G}$ in system (<ref>), for $\epsilon=0$. Moreover, let  $\LLL_{0}^{\yy}$ be
the (uncoupled) Koopman operator associated with (<ref>) as defined in Eq. (<ref>).

The point spectrum of $\LLL_{0}^{\yy}$ is assumed to be constituted of $N$ simple eigenvalues, whose corresponding eigenpairs
$\left\{ (\lambda_j,\psi_{j}), \; j=1,\ldots,N \right\}$ are ordered as
follows: $0\geq \mathfrak{Re}\lambda_j\geq \mathfrak{Re}\lambda_{j+1}$ and $\lambda_j=\overline{\lambda_{j+1}}$ when $\mathfrak{Im}\lambda_j>0$,
for $j$ in $\{1,\ldots, N\}$.

We assume that  $\CCC^{\xx}$ in (<ref>) lies in the $\mathrm{span}\left\{ \psi_{j}, \; j=1,\ldots,N \right\}$ and has $\nu$-mean zero.

Then, the WL equation (<ref>) associated with system (<ref>) admits a Markovianization of the form:

\begin{align}
\label{remarkovianisation3}	\dot{x}(t) & = \FFF (x(t)) + \epsilon \Lambda\cdot \mathbf{Z}(t), \\
\label{remarkovianisation2} \dd\mathbf{Z}(t) & =
\left(\epsilon \mathbf{R}(x(t))  +
\mathrm{D}\mathbf{Z}(t)\right)\dd t + \Sigma \dd W_t,
\end{align}

where $\Lambda$ and $\mathbf{Z}(t)$ lie in $\mathbb{C}^{N}$ for every $t$,  while the inner product $\Lambda \cdot\mathbf{Z}(t)$ is real.
Here, $\mathbf{R}$ is mapping  $\mathbb{R}$ into $\mathbb{C}^{N}$, $W_t$
is a (real-valued) $N$-dimensional Wiener process with covariance matrix $\Sigma$,  and $\mathrm{D}$ is an $N \times N$ matrix with complex entries, as specified below.

More precisely,

\begin{equation}\label{Eq_lambda}
\Lambda=\left[  \alpha_1^{1/2}\beta_1^{1/2}, \ldots,\alpha_N^{1/2}\beta_N^{1/2} \right]^{\top},
\end{equation}

\begin{align}
\label{eq:alpha}	\alpha_{j}&=\left \langle \psi^{\ast}_j , \CCC^{\xx}
\right \rangle = \int \nu(\dd \yy)\overline{\psi_j^{\ast}(\yy)}\CCC^{\xx}(\yy), \\
\label{eq:beta}	\beta_j&=\left \langle\CCC^{\xx}, \psi_j \right \rangle = \int \nu(\dd \yy)\CCC^{\xx}(\yy)\psi_j(\yy).
\end{align}

The $\mathbb{C}^{N}$-valued mapping $\mathbf{R}$ is defined as

\begin{equation}\label{Eq_Rexpression}
\mathbf{R}\left(x \right)= \left(\CCC^{\yy}(x)\cdot\frac{\alpha_1^{1/2}}{\beta_1^{1/2}}\bigg \langle \nabla\psi _1 (\yy)  \bigg \rangle,\ldots,\CCC^{\yy}(x)\cdot\frac{\alpha_N^{1/2}}{\beta_N^{1/2}}\bigg \langle \nabla\psi _N (\yy)  \bigg \rangle\right),
\end{equation}

where $\CCC^{\yy}$ is defined in  (<ref>), and $\langle \cdot \rangle$ denotes averaging with respect to the invariant measure $\nu$.

Finally, $\mathrm{D}=\mathrm{diag}\left(\lambda_1,\ldots,\lambda_N\right)$ and the covariance matrix $\Sigma$ is given by

\begin{align} \label{eq: def covariance matrix proposition}
\Sigma=\left(-\left(\mathrm{D}+\mathrm{D}^{\ast}\right)\right)^{1/2}\mathrm{H},
\end{align}

where $\mathrm{H}$ is an $N\times N$  matrix whose entries are defined as
follows: If $\lambda_j$ is real, then $\mathrm{H}_{j,j}=1$, and if $\lambda_j=\overline{\lambda_{j+1}}$,

while all other entries are zero.

The full proof appears in Appendix <ref>.

Note that the Koopman operator of interest here is the one associated with the $\yy$-subspace $\mathcal{Y} \subseteq \mathbb{R}^{d_2}$ and not with the entire $(\xx, \yy)$-space $\mathcal{X}\times \mathcal{Y}\subseteq \mathbb{R}^{d_1+d_2}$. Other techniques, like the DMD mentioned in Sect. <ref>, aim at extracting the modes of variability of the full system by means of studying the Koopman operator in the entire phase space through suitably defined observables. To this end, the latter methods employ projections of observables onto the eigenfunctions of the Koopman operator to obtain the so-called Koopman modes, which are susceptible of capturing the underlying dynamics. Notice that in Theorem 2.1, instead, we are using the Koopman eigenfunctions to identify the closure model, while projections only come into play in the definition of the coefficients $\alpha_{j}$ and $\beta_j$; see Eqs. (<ref>) and (<ref>), respectively.

(i) Assumptions on $\FFF$, $\mathbf{R}$ and $\Lambda$ that ensure that (<ref>) possesses a global random attractor — and thus a stable asymptotic behavior in the pullback sense —  appear in <cit.>.

(ii) Note that Theorem <ref> can be viewed as a generalization of other Markovianization results for GLEs that appeared in the literature; see <cit.>. For instance,
the scalar GLE in $\mathbb{R}$ reduces to

ẋ=F (x(t)) -∫_0^t K(t-s) x(s) s +η(t),
where $\lambda$ is in $\mathbb{R}^{n}$, $M$ is a positive definite $n\times n$ matrix and $K(t-s)=\left(e^{M (t-s)} \lambda\right)\cdot \lambda$
determines the autocorrelation of the process $\eta(t)$.  In this setting, Eq. (<ref>) is equivalent to the following SDE:
ẋ=F(x(t))+ λ·z

z  =(xλ-M z) t + ΣW_t,
with $\Sigma \Sigma^{\ast}=M+M^\ast$. Theorem <ref> allows for nonlinear dependence on $x$ in the $z$-equation,
and thus for memory kernels that are more complicated  than in (<ref>). Such a  generalization is of practical importance since the process $z$ can then have a more complex correlation dependence on the observed
variable $x$ than the one afforded by linear memory terms.

When $\LLL_{0}^{\yy}$ is self-adjoint — in a suitable Hilbert space, as outlined in Appendix <ref>, i.e., when $\LLL_{0}^{\yy}=\LLL_{0}^{\yy^{\ast}}$ — the eigenvalues are real and the eigenvectors are mutually orthogonal. Self-adjointness thus  implies that there are no oscillations in the correlation functions of the noise or, equivalently, peaks in their power spectrum.  With respect to Theorem <ref>, the matrix $\mathrm{H}$ in this case would be the identity, since the
eigenvalues and eigenfunctions are real and the Itô solutions of <ref> are, hence, real as well.

The resulting system given by Eq. (<ref>) is now fully Markovian and the only sources of error with respect to the original SIDE
(<ref>) lie (i) in the effects of the essential spectrum, which are neglected herein; and (ii) the assumptions about the coupling terms. Neglecting the essential spectrum is only valid for Koopman operators with a point spectrum capable of capturing the correlations in the
decoupled $\yy$-system; the latter might only hold in the case of Markovian diffusion processes and not for deterministic ones. Also, the assumption that the coupling functions project solely on the point spectrum might
not hold in general.

From a practical perspective, though, a suitable choice of dominant eigendirections can reduce the number of extra dimensions needed to integrate the system. Such a suitable choice boils down to neglecting particular eigendirections and this can be done according to two handy criteria:

(i) The weight determined by the $\alpha_{j}$ and $\beta_j$ coefficients defined in Eqs. (<ref>, <ref>) is small; and

(ii) The eigenvalues $\lambda_j$ of $\LLL_{0}^{\yy}$ satisfy $\mathfrak{Re} \lambda_{j_0} \ll \lambda^{\dagger}$, for some $j_0\in \{1,\ldots,N\}$, in which case $e^{\lambda_{j_0} t}$ decays rapidly as $t$ grows;
here $\lambda^{\dagger} < 0$ and $|\lambda^{\dagger}|$ is some characteristic inverse time for the deterministic system $\FFF$. In addition, if $\alpha_{j_0}>\alpha_{j}$ for $j=1,\ldots,N$ and $j\neq j_0$, both the memory and the noise correlations die out fast. Hence, one can neglect the integral terms and perform a fully Markovian parametrization, which is possible in the presence of white noise.

Theorem <ref> is stated for $x$ scalar for the sake of simplicity, but this result extends to the $d_1$-dimensional  Eq. (<ref>) for the (observed) variables $\xx$. In this Remark, we sketch the main elements that permit such a generalization.

Aside from the obvious generalization of the assumptions in Theorem <ref> to a multidimensional setting, the main  hypothesis consists of assuming that the now vector-valued coupling function $\CCC^{\xx}$ in Eq. (<ref>) has components $\{\left[\CCC^{\xx}\right]_{i}: i=1,\ldots, d_1\}$ that project onto the $N$ simple eigenspaces of the decoupled Koopman operator introduced in Eq. (<ref>). In this case,
the construction of a multilevel Markovianization like Eq. (<ref>) can be done in the following fashion:

\begin{align}
\dot\xx(t) &= \FFF(\xx(t))+\epsilon \Lambda\mathbf{Z}(t), \\
\dd \mathbf{Z}(t) &= \left(\epsilon \mathbf{R}(\xx(t))+\mathcal{D}\mathbf{Z}(t) \right)\dd t + \mathcal{S} \dd W_t.
\end{align}

Here $W_t$ is a $d_1N$-dimensional Wiener process, $\mathbf{Z}(t)$ is a $d_1N$-dimensional vector, $\Lambda$ is a matrix of size $d_1\times d_1N$, $\mathbf{R}:\mathbb{R}^{d_1}\longrightarrow \mathbb{C}^{d_1N}$ and $\mathcal{D}$ and $\mathcal{S}$ are $d_1N\times d_1N$ block-diagonal matrices
given by:

\begin{equation}\label{Eq_multilevel matrices}
\mathcal{D}=\begin{bmatrix}\mathrm{D}_1 &&\\ &\ddots &\\ && \mathrm{D}_N\end{bmatrix} \text{ and } \mathcal{S}=\begin{bmatrix}\Sigma_1 &&\\ &\ddots &\\ && \Sigma_N\end{bmatrix},
\end{equation}

where $\mathrm{D}_j$ and $\Sigma_j$ are $d_1\times d_1$ diagonal matrices with every (non-zero) element being equal to $\lambda_j$ or $\sqrt{-2\mathfrak{Re}\lambda_j}$, respectively. More importantly, the vectors $\mathbf{Z}(t)$ and $\mathbf{R}(\xx)$ are split into $N$ column vectors $\mathbf{z}_j(t)$ and $\mathbf{r}_j(\xx)$ of length $d_1$ with $j$ in $\{1,\ldots, N  \}$. This way, $\mathbf{Z}(t)=\left[ \mathbf{z}^{\top}_1(t),\ldots,\mathbf{z}^{\top}_{N}(t) \right]^{\top}$ and $\mathbf{R}(\xx)=\left[
\mathbf{r}^{\top}_1(t),\ldots,\mathbf{r}^{\top}_{N}(t) \right]^{\top}$.

Therefore Eq. (<ref>) can be written as:

\begin{align}
\dot\xx(t) &= \FFF(\xx(t))+\epsilon \Lambda\mathbf{Z}(t), \\
\dd \mathbf{z}_1(t) &= \left(\epsilon \mathbf{r}_1(\xx(t))+\mathrm{D}_1\mathbf{z}_1(t) \right)\dd t + \Sigma_1
\dd W^{(1)}_t, \\ & \vdots \nonumber \\ \label{eq: multilevel markovianization2 last level}
\dd \mathbf{z}_N(t) &= \left(\epsilon \mathbf{r}_N(\xx(t))+\mathrm{D}_N\mathbf{z}_N(t) \right)\dd t + \Sigma_N \dd W^{(N)}_t,
\end{align}

where $W^{(j)}_t$ is a $d_1$-dimensional Wiener process. The vectors $\mathbf{r}_j(\xx)$ are given by:

\begin{equation}
\mathbf{r}_j(\xx)=\left[ \CCC^{\yy}(\xx)\cdot \gamma_{1,j}\bigg \langle \nabla\psi _1 (\yy)  \bigg \rangle,\ldots, \CCC^{\yy}(\xx)\cdot \gamma_{d_1,j}\bigg \langle \nabla\psi _1 (\yy)  \bigg \rangle \right]^{\top},
\end{equation}

where $\gamma_{i,j}$ are defined in terms of the parameters $(\alpha_{j}, \beta_j)$ introduced in Eqs. (<ref>, <ref>), respectively. Here we don't give the explicit expression of $\Lambda$, but its role is to provide suitable weights, in the spirit of Eq. (<ref>),
to the levels in Eq. (<ref>) so that: (a) the correlation functions match those of the coupling function $\CCC^{\xx}$ in the uncoupled regime; and (b) the resulting term $\Lambda\mathbf{Z}(t)$ is real.

The system Eq. (<ref>) has the general structure one would obtain if the coupling function $\CCC^{\xx}$ projected along all the eigendirections in the point spectrum. This might not be true in general, but the drift matrix $\mathcal{D}$ can be rearranged so that only the relevant modes of variability be modeled — following criteria (i) and (ii), as formulated in Remark <ref> — and still afford a reduction of the number of levels $N$. Notice that the $N$th level variables $\mathbf{z}_{N}$ described by Eq. (<ref>) decorrelate the fastest,
the rest,
since their exponential decorrelation rate is given by $|\mathfrak{Re}\lambda_N |\geq |\mathfrak{Re}\lambda_j|$, for all $j=1,\ldots,N-1$.

The advantages of the Markovian system of Eqs. (<ref>) and (<ref>) over the original WL equation (<ref>) are twofold. First, we identify situations in which the WL equation can be Markovianized by introducing
extended, hidden variables. This idea was already introduced in a preliminary application of the WL parametrization <cit.>, in which the authors resorted to a Markovian system to perform their simulations.  In fact, one of their examples is studied in the present framework; see  Sect. <ref>.

Second, memory equations contain nonlocal terms that are cumbersome and computationally expensive to integrate, as well as requiring much larger storage for the full history of the system's variables. The efficient Markovianization of evolution equations with memory terms is an active field of research in diverse areas of mathematics and the applied sciences; these areas include the study of bifurcations of delay differential equations <cit.>, the reduction of stochastic partial differential equations to stochastic invariant manifolds <cit.>, and material sciences <cit.>, among many others.

§.§ Preliminary Example

As seen earlier in Theorem <ref>, if the coupling function
is resonant with the Koopman operator associated with the $\yy$-dynamics,
one can identify the dominant exponential rates of decay of the memory term and the characteristic decorrelation time of the noise. As a consequence, one can Markovianize the parametrization and greatly facilitate the numerical integrations involved.

To illustrate the above statement, we revisit here the preliminary application of the WL parametrization in the context of multiscale triads <cit.>. In that work, the authors implemented the parametrization for a collection of three-dimensional models that do exhibit time scale separation and compare here the corresponding outputs to those obtained via homogenization. The results are encouraging, since the parametrizations
in <cit.> were obtained only from the decoupled hidden dynamics, in the lines of the present paper as well; see derivation of Eq. (<ref>).

One of the first multiscale triads studied in <cit.> is the following:

\begin{align}
\dot{x}(t) & = \epsilon B^{(0)}y_1y_2,
\label{eq:B0}\\
\dot{y_1}(t) & = \epsilon B^{(1)}xy_2 - \gamma _1y_1 + \sigma _1\dd W^{(1)}_t,  \label{eq:B1}\\
\dot{y_2}(t) & = \epsilon B^{(2)}xy_1 - \gamma _2y_2 + \sigma _2\dd W^{(2)}_t. \label{eq:B2}
\end{align}

Here we require that $\sum_jB^{(j)} = 0$, $\dd W^{(1)}_t$ and $\dd W^{(2)}_t$ are scalar Brownian increments and the parameter $\epsilon$ indicates both the time scale separation and the coupling strength. Notice that
when the system is decoupled, i.e. when $\epsilon=0$, the fast dynamics
evolve according to an Ornstein-Uhlenbeck (OU) process whose steady-state
statistics are governed by Gaussian distributions with explicit mean and variance (see, e.g. <cit.>). Hence, by virtue of the previous formulas or by following <cit.>, the WL parametrization yields the following scalar SIDE:

\begin{equation}\label{eq:scalar_memory}
\dot{x}(t)=\epsilon \eta (t) + \epsilon ^2 \int _0^t\mathcal{K}(s,x(t-s))\dd s;
\end{equation}

\begin{align}
\left\langle \eta(t) \right\rangle & = 0, \\
\left\langle \eta(t+s)\eta(s) \right\rangle & = \left( B^{(0)}\right)^2e^{-(\gamma_1 + \gamma_2)t}\frac{\sigma_1^2}{2\gamma_1}\frac{\sigma_2^2}{2\gamma_2},\\ \label{eq: memory triad 1}
\mathcal{K}(s,x)&= \begin{bmatrix}
x \\ x
\end{bmatrix}\cdot \bigg \langle \begin{bmatrix}
B^{(1)}y_2 \\B^{(2)} y_1
\end{bmatrix}: \nabla_{\yy} y_1(s)y_2(s) \bigg \rangle,
\end{align}

where the angular brackets refer to the ensemble averages according to the already mentioned Gaussian distributions arising from the decoupled model. Expanding these averages Eq. (<ref>) leads to:

\begin{align}
\mathcal{K}(s,x)&= xe^{-(\gamma_1 + \gamma_2)}\bigg \langle B^{(1)}y^2_2 +  B^{(2)}y^2_1\bigg \rangle  \\ & = xB^{(0)}e^{-(\gamma_1 + \gamma_2)s}\left(B^{(1)}  \frac{\sigma_2^2}{2\gamma_2}+ B^{(2)} \frac{\sigma_1^2}{2\gamma_1}\right).
\end{align}

The time scales are indicated by the exponents in the formulas above and they are the same for the noise and the memory kernel. This equality suggests the possibility of Markovianizing the memory equation into the following two-dimensional system:

\begin{align}
\dot{z_1}(t)&=\epsilon B^{(0)}z_2, \\
\dot{z_2}(t)&=-(\gamma_1 + \gamma_2)z_2 + \frac{\sigma_1\sigma_2}{2\gamma_1\gamma_2}
\{2(\gamma_1 + \gamma_2)\}^{1/2}
\dd W_t +\epsilon \left( B^{(1)}\frac{\sigma_2^2}{2\gamma_2} + B^{(2)}\frac{\sigma_1^2}{2\gamma_1} \right)z_1.
\end{align}

Clearly, performing a numerical integration of this system is easier
than for a memory equation like Eq. (<ref>).

The results of Sec. <ref> allow us to carry out the dimension reduction of the multiscale triad by analyzing the spectral properties of the Koopman operator associated with the decoupled $\yy$-dynamics. Since the $\yy$-variables evolve stochastically, the Koopman operator becomes the backward-Kolmogorov equation, which governs the evolution of the expectation values of the observables. Thus, for a generic observable $\Psi$ in the $\yy$ phase space, the evolution of its expectation
value is given by:

\begin{equation}
\partial _t \Psi (y_1,y_2)= \LLL_0^{\yy}\Psi(y_1,y_2) = \begin{bmatrix}
-\gamma_1 y_1 \\ -\gamma_2y_2
\end{bmatrix}\cdot \nabla \Psi (y_1,y_2) + \sigma^2_1\partial^2_{y_1}\Psi(y_1,y_2) + \sigma^2_2\partial^2_{y_2}\Psi(y_1,y_2).
\end{equation}

Now, let $\Psi(y_1,y_2)=y_1y_2$ be the coupling function of the triad system (<ref>), for which we find that:

\begin{equation}
\LLL^{\yy}_{0}\Psi(y_1,y_2) = -(\gamma_1 + \gamma_2)\Psi(y_1,y_2).
\end{equation}

The above equation is an eigenvalue problem, showing that this particular
$\Psi$ is an eigenfunction of the Koopman operator associated with the eigenvalue $e^{-(\gamma_1 + \gamma_2)}$. This is no surprise, since $y_1$ and $y_2$ are respectively the Hermite polynomial eigenfunctions of the backward-Kolmogorov equation of the scalar OU process <cit.>. Hence, the product $y_1y_2$ is also an eigenfunction of the same equation for the joint process. Therefore, we can immediately re-Markovianize the parametrization according to Eqs. (<ref>), where $\mathrm{D}=\gamma_1 + \gamma_2$ and

\begin{equation}
\mathbf{R}(x(t))=\left( B^{(1)}\frac{\sigma_2^2}{2\gamma_2} + B^{(2)}\frac{\sigma_1^2}{2\gamma_1} \right)x(t).
\end{equation}

§ MULTILEVEL STOCHASTIC MODELS AND EMPIRICAL MODEL REDUCTION (EMR)

§.§ Multilevel Stochastic Models (MSMs)

MSMs are a general class of  SDEs that were introduced in <cit.> and are, by their layered structure,  susceptible to provide a
good approximation of the GLE (<ref>) formulated by Mori and Zwanzig when a high-dimensional system is partially observed; see  <cit.>. The MSM framework allows one to provide such approximations that are accompanied by useful dynamical properties, such as the existence of random attractors <cit.>. The conditions on the high-dimensional system's coupling interactions between the resolved and hidden variables are also well understood <cit.>. [Moreover, random attractors with fractal structures that survive highly degenerate noise <cit.>—a situation that might occur when approximating deterministic chaotic dynamics by stochastic pathwise dynamics—can still be present in MSMs; see <cit.>.]

As discussed in <cit.>, MSMs arise in a variety of data-driven protocols for model reduction that typically use successive regressions from partial observations; see Sec. <ref> below. The general form of an MSM is given by <cit.>; we only use herein its most basic version, which has the following structure:

\begin{align}
\dd \xx(t) & = \Big(\FFF (\xx(t)) + \epsilon \Pi \rr(t)\Big) \dd t, \label{msm1} \\
\dd \rr(t) & = \Big(\epsilon \mathrm{C}\xx(t) -\mathrm{D}\rr(t)\Big) \dd t + \Sigma \dd \W_t. \label{msm2}
\end{align}

Here the observed vector variable $\xx(t)$ lies in $\mathbb{R}^{d_1}$ and, for $\epsilon = 0$, the hidden variables $\rr(t) \in \mathbb{R}^{d_2}$ evolve in time independently. Otherwise, the dynamics of the $\xx$-variables is linearly coupled to that of the $\rr$-variables, which act upon (<ref>) as a stochastic forcing, via the canonical projection $\Pi: \mathbb{R}^{d_2} \longrightarrow \mathbb{R}^{d_1}$, while $\W_t$ in (<ref>) is a $d_1$-dimensional Wiener process.

The matrix $\mathrm{C}$ in $\mathbb{R}^{d_2\times d_1}$ models the feedback of the $\xx$-process onto the $\rr$-variables. In the case of $\mathrm{C}\equiv 0$, $\rr$ whould evolve according to an OU process with drift matrix $\mathrm{D}$ and covariance matrix $\Sigma\Sigma ^{\ast}$. For the sake of simplicity, we restrict ourselves to the case $d_1=d_2$ so that
the projection $\Pi$ reduces to the identity.

The more general MSM with nonlinear coupling considered in <cit.> was shown to be equivalent to a SIDE with explicit expressions for the memory kernels and stochastic forcing being obtained; see <cit.>. The noise term there results from successive convolutions of the homogeneous solutions of the lower levels of the system with an OU process. In particular, using the Itô stochastic calculus, one readily obtains a SIDE that is equivalent to an MSM; see
<cit.> and Appendix <ref> below.

We show next that the same SIDE can actually be obtained by using the operator formalism presented in Sec. <ref> above. One might object that an MSM is a stochastic system, due to the
presence of white noise in the hidden layer, whereas the theory presented
above applies to deterministic dynamics. However, as clarified below, the
operator formalism applies equally well to the MSM case.

In fact, given a smooth, $\CC^\infty$ observable
Φ:ℝ^d_1×ℝ^d_2 →ℝ,   (,) ↦Φ(,),
its expected value along a stochastic trajectory $X_t=(\xx(t),\rr(t))^{\top}$ solving  Eq. (<ref>), namely $\mathbb{E}(\Phi(X_t))$, defines a Markov semigroup $P_t$ by
P_t Φ=𝔼(Φ(X_t)),
which solves the backward Kolmogorov equation <cit.> associated with  Eq. (<ref>):

\begin{equation}\label{backward Kolmogorov}
\partial _t \big (P_t\Phi \big) =\begin{bmatrix}\FFF(\xx)+\epsilon \rr\\ \epsilon \mathrm{C}\xx - \mathrm{D}\rr \end{bmatrix}\cdot \nabla P_t\Phi + \frac{1}{2}\begin{bmatrix} 0 \\  \Sigma\Sigma^{\ast}\nabla_{\rr}^2 P_t\Phi  \end{bmatrix};
\end{equation}

the only difference with respect to the transport equation (<ref>) lies in the presence of a second-order differential operator
induced by the white noise.

We introduce the operators:

\begin{align}
&  \LLL_0=\begin{bmatrix} \FFF(\xx)\\ - \mathrm{D}\rr \end{bmatrix}\cdot \nabla +\begin{bmatrix} 0 \\     \Sigma\Sigma^{\ast}\nabla_{\rr}^2 \end{bmatrix}, \label{eq:LLL_0} \\
& \LLL_1=\begin{bmatrix} \rr\\ \mathrm{C}\xx \end{bmatrix}\cdot \nabla,
\label{eq:LLL_1}
\end{align}

which play a role that is analogous to their deterministic relatives in Eq. (<ref>) of the previous section. Again, the operator $\LLL_1$ is viewed as a perturbation to the operator $\LLL_0$ due to the coupling.

If one considers observables $\Phi=\Phi(\xx)$, Eq. (<ref>) becomes at time $t=0$:

\begin{equation}
\partial _t \big (P_t\Phi \big)|_{t=0}=\begin{bmatrix}\FFF(\xx)+\epsilon \rr \end{bmatrix}\cdot \nabla _{\xx} \Phi,
\end{equation}

and we apply now, as in Sec. <ref>, the  Dyson perturbative expansion. By virtue of the formula (<ref>),
the parametrization leads to a reduced equation of the form:

\begin{equation}
\dot{\xx}(t) = \FFF(\xx(t)) + \epsilon \eta (t) + \epsilon^2 \int _0^t\mathcal{K}(s , \xx(t-s))\dd s,
\end{equation}

where the hidden variables in the decoupled regime are governed by an OU process with invariant measure $\mu_{\rr}$. The properties of the stochastic noise $\eta(t)$ are given by:

\begin{align}\label{msm wl}
\left\langle \eta(t)\eta^{\top}(0) \right\rangle  &= \int \dd \mu_{\rr}(\rr_0) e^{t\LLL_0} \rr_0\rr_0^{\top} \\&= \int \dd \mu_{\rr}(\rr_0) \mathbb{E}\left(\rr (t)|\rr_0 \right)\left(\mathbb{E}\left(\rr(0)|\rr_0 \right)\right)^{\top} \\
&= \int \dd \mu_{\rr}(\rr_0) e^{-t\mathrm{D}} \rr_0\rr_0^{\top}
\\
&= \int \dd \mu_{\rr}(\rr_0) e^{-t\mathrm{D}} \rr_0\rr_0^{\top} = e^{-t\mathrm{D}}\Sigma\Sigma^{\ast},
\end{align}

where $\rr$ is a function analogous to the coupling function $\CCC^{\xx}_{\yy}$ in the previous section and the initial condition $\rr_0$ is assumed to be normally distributed with zero mean and variance $\Sigma\Sigma^{\ast}$. The memory kernel is given by:

\begin{align}
\mathcal{K}(s,\xx(t-s))&=\int \dd \mu_{\rr}(\rr_0) \mathrm{C}\xx(t-s)\cdot \nabla_{\rr_0}\mathbb{E}\left( \rr(s)|\rr_0 \right)\\&=\int \dd \mu_{\rr}(\rr_0) \mathrm{C}\xx(t-s)\cdot \nabla_{\rr_0}e^{-s\mathrm{D}}\rr_0
\\&= \mathrm{C}\xx(t-s)\cdot e^{-s\mathrm{D}},
\\&= e^{-s\mathrm{D}}\mathrm{C}\xx(t-s).
\end{align}

Using the intermediate steps above, the explicit parametrization becomes,

\begin{equation}
\dot{\xx}(t) = \FFF(\xx(t)) + \epsilon \eta (t) + \epsilon^2 \int _0^te^{-s\mathrm{D}}\mathrm{C}\xx(t-s)\dd s.
\end{equation}

The integro-differential equation above is the same as Eq. (<ref>)
one obtains using the Itô integration described in Appendix <ref>. This similarity of results occurs because we are considering the case of additive coupling, and the Dyson expansion can be truncated after the memory term proportional to $\epsilon^2$, cf. Sect. <ref>.

§.§ Empirical Model Reduction (EMR)

As discussed in Sects. <ref> and <ref>, and illustrated in Fig. <ref>, the evolution of the resolved variables is forced by fluctuating terms and the effects of the previous state of the system. It is desirable, therefore, to construct a full model of the system even when only capable to partially observe it. The EMR methodology <cit.> aims at achieving this goal; we discuss it below in the broader context of MSMs. Note that EMR provides a solution for the
dynamical closure of partially observed systems and thus it differs from the methodology recently proposed in ref. <cit.>, which requires one to fully observe the system for the data-driven discovery of
its underlying equations to work.

Having a set of reduced $d_1$-dimensional observations $\{\xx_i: i = 1,
\ldots, n\}$ every $\dd t$ time units, one seeks to regress the tendencies $\{\dd \xx_i: i = 1, \ldots, n\}$ of the data onto a quadratic function of the form:

\begin{equation}\label{eq: EMR_vector_field}
\FFF(\mathbf{x})=\mathbf{f} + \mathbf{b}\cdot \mathbf{x} +\mathbf{Q}(\xx),
\end{equation}

where $\mathbf{b}$ in $\mathbb{R}^{d_1\times d_1}$ describes dissipative
processes and $\mathbf{Q}$ is a quadratic form describing self-interaction between the $\xx$ variables. The $i$th component of the quadratic form is given by:

\begin{equation}
\end{equation}

where $\mathbf{A}_i$ in $\mathbb{R}^{d_1 \times d_1}$. The function $\FFF$ is expected to approximate the vector field driving the dynamics in the
absence of hidden external influences. Of course, performing regressions yields an error called residual $\{\rr_i: i = 1, \ldots, n\}$. Hence, the evolution of the $\xx$-variables satisfies the equation:

\begin{equation}
\frac{\dd \xx}{\dd t}  = \FFF(\xx) + \rr.
\end{equation}

At this point, one can study the properties of the residual time series $\{\rr_i\}_{i=1}^{n}$ and construct a model able to reproduce its main statistical features. However, we know that if it is possible to sample all the variables of the dynamical system of interest, one expects that the
residuals are explained by the errors committed exclusively in the regression algorithm. If some sort of subsampling is done, whether spatial or temporal, the residuals are due also to the delayed influence of unresolved processes that are involved in the coupling.

Allowing for the main level variables $\xx$ to be linearly coupled with the residual $\rr$, we are creating a model that is able to incorporate memory effects as well. Hence, for each component $i$ of $\xx$ we have:

\begin{align}
& \frac{\dd  [\mathbf{x}]_i }{\dd t} = [\mathbf{f}]_i + \mathbf{b}_i^{(0)}\cdot \mathbf{x}+ \mathbf{x}^{\top}\mathbf{A}_i\mathbf{x}  + [\rr^{(0)}]_i, \label{EMR_0}\\
& \frac{\dd   [\rr^{(0)}]_i}{\dd t}  = \mathbf{b}_i^{(1)}\cdot [\mathbf{x},\mathbf{r}^{(0)}] + [\rr^{(1)}]_i, \label{EMR_1}\\
& \frac{\dd   [\rr^{(1)}]_i }{\dd t} = \mathbf{b}_i^{(2)}\cdot [\mathbf{x},\mathbf{r}^{(0)},\mathbf{r}^{(1)}] +
[\rr^{(2)}]_i, \label{EMR_2}\\
& \qquad \vdots  \label{EMR_k} \\
& \frac{\dd   [\rr^{(l)}]_i}{\dd t} = \mathbf{b}_i^{(l+1)}\cdot [\mathbf{x},\mathbf{r}^{(0)},\mathbf{r}^{(1)},\ldots, ,\mathbf{r}^{(l)}] 	+ [\rr^{(l+1)}]_i, \label{EMR_l}
\end{align}

where we have introduced new matrices $\mathbf{b}^{(j)}\in \mathbb{R}^{d_1 \times (j+2)d_1}$ that model the linear coupling.  The residual at the last level $[\rr^{(l+1)}]_i$ is assumed to obey a Wiener process for which the correlation matrix is obtained from the last residual time series. The choice of stochastic process in the last step can only be done if the
decorrelation of $[\rr^{(l+1)}]_i$ is sufficiently fast according to the time scale set by $\dd t$. This motivates the problem of choosing the optimal number $l$ of levels.

Several criteria have been established to determine the optimal number of
levels $l$.
The basic idea is that the resulting $(l + 1)$-residual in Eq. (<ref>) should be well approximated by Gaussian white noise <cit.>. One has, therefore, to test whether the residual variables decorrelate at lag $\dd t$ and whether the lag-0 covariance matrix is invariant in the last levels.

Therefore, regression on the tendency of the optimal level $\rr^{(l)}$ should yield:

\begin{equation}
\rr^{(l+1)}-\rr^{(l)} \simeq -\rr^{(l)} + \gamma ^{(l)},
\end{equation}

where $\gamma^{(l)}$ is the residual of the previous regression and is approximately equal to $\rr^{(l+1)}$. Hence, $\gamma^{(l)}$ would become a lagged version of $\rr^{(l+1)}$. Subject to this assumption, it is possible to estimate the optimal value of the coefficient of determination $R^2$:

\begin{equation}\label{emrdiagnostic}
R^2=1-\frac{\sum _{k}\gamma_l^2}{\sum_{k}\left(\rr^{(l+1)} - \rr^{(l)}
\right)^2} \simeq 1 - \frac{\sum _l{\rr^{(l+1)}}^2}{\sum _l {\rr^{(l+1)}}^2 + {\rr^{(l)}}^2} \simeq 0.5.
\end{equation}

This means that, when the amount of unexplained variance of the last regression is $50\%$, one has reached the optimal number of levels.
It is worth stressing that the empirical model (<ref>) has the structure (<ref>) of an MSM, as discussed in <cit.>. It can, therewith, be integrated to transform it into an integro-differential equation with explicit formulas for the fluctuating noise and memory kernel, cf. <cit.>; see also Sec. <ref> for such a transformation from another perspective.

Finally, note that the aforementioned stopping criterion for EMR — namely $R^2\simeq 0.5$, see <cit.>) — is based on decorrelation times and it is also present in the multilevel WL equation (<ref>). We noted, in fact, in Remark <ref> of Sect. <ref> that the last level modeled by Eq. (<ref>) decorrelates the fastest with respect to the rest. Ultimately, making this points amounts to saying that a low number of levels is expected to arise
in the EMR method, provided most of the eigenvalues $\lambda_j$ in Theorem <ref> are located far away from the imaginary axis, except for a very few of them.
Conversely, if the Koopman eigenvalues cluster near the imaginary axis or do not exhibit a spectral gap located at a, suitably defined, small negative real part, many levels are expected to be needed to capture the hidden dynamics; see again Remark <ref>.

§ NUMERICAL EXPERIMENTS

In Sects. <ref> and <ref>, we have shown that both the WL top-down approach and the EMR data-driven method yield a set of multilevel equations for the variables of interest in a multi-scale system.
In particular, both approaches give explicit formulas for the fluctuation
term and memory kernel in the GLE (<ref>) of the Mori-Zwanzig formalism. Furthermore, their Markovian representation share that the
hidden layers are linearly driven with a white noise background (see Eq. <ref> and Eq. (<ref>)). We now compare the two
approaches to model reduction in a simple, conceptual stochastic climate model.

Since the  modeling of geophysical flows is the primary motivation for this research, we consider a set of SDEs proposed in <cit.> as a physically consistent climate “toy” model. In such a model, the main $\xx$-variables are slow and weakly coupled to the fast $\yy$-variables. The latter correspond to weather fluctuations and carry, in
fact, most of the system's variance. The model's governing equations are:

\begin{align}
\dd x_1 & = \left\lbrace -x_2\left( L_{12} + a_1x_1 + a_2x_2\right) -d_1x_1 +F_1 +\epsilon\left(L_{13}y_1+c_{134}y_1y_2\right)\right\rbrace \dd
t, 	\label{stochastic model 2} \\
\dd x_2 & =  \left\lbrace x_1\left( L_{21} + a_1x_1 + a_2x_2 \right) -d_2x_2 +F_2 +
\epsilon L_{24}y_2  \right\rbrace \dd t, \label{stochastic model 3} \\
\dd y_1 & = \left \lbrace \epsilon \left(-L_{13}x_1 + c_{341}y_2x_1\right) + F_3 - \frac{\gamma_1}{h}y_1  \right\rbrace \dd t +\frac{\sigma _1}{\sqrt{h}}\dd W^{(1)}_t, \label{stochastic model 4} \\
\dd y_2 & = \left\lbrace -\epsilon \left(L_{24}x_2 + c_{413}y_1x_2\right) +F_4 - \frac{\gamma_2}{h}y_2 \right \rbrace \dd t + \frac{\sigma_2}{\sqrt{h}}\dd W^{(2)}_t,\label{stochastic model 5}
\end{align}

where $W^{(1)}_t$ and $W^{(2)}_t$ are two independent Wiener processes.

These equations describe the evolution of four real variables $\xx=(x_1,x_2)$ and $\yy=(y_1,y_2)$; their time scale separation is determined by the parameter $h$ and the coupling strength is controlled by $\epsilon$. The parameter values used herein are: $c_{134}=c_{341}=0.25$, $c_{413}=-0.5$, $L_{12}=L_{21}=1$, $L_{24}=-L_{13}=1$, $a_1=-a_2=1$, $d_1=0.2$, $d_2=0.1$, $F_1=-0.25$, $F_2=F_3=F_4=0$, $\gamma _1 = 2 $, $\gamma _2 = 1$ and $\sigma _1 = \sigma _2 = 1$. The time scale separation and the coupling strength are $h=0.1$ and $\epsilon= 0.5$, respectively.

§.§ WL Approximation

Notice that the hidden variables evolve according to a decoupled OU process. Taking advantage of this fact, we calculate the weak-coupling–limit parametrization of the model, according to the formulas presented in Sect. <ref> for the separable coupling functions given by:

\begin{align}
\CCC ^{\xx}(\xx , \yy)=\CCC _{\yy}^{\xx}(\yy) &= \begin{bmatrix}
L_{13}y_1 + c_{134}y_1y_2, \\
\end{bmatrix}, \label{eq:slow} \\
\CCC ^{\yy}(\xx , \yy)=\CCC ^{\yy}_{\xx}(\xx ):\CCC _{\yy}^{\xx}(\yy) &= \begin{bmatrix}
-L_{13}x_1 + c_{341}y_2x_1 \\
-L_{24}x_2 + c_{413}y_1x_2
\end{bmatrix}. \label{eq:fast}
\end{align}

The coupling function $\CCC ^{\xx}$ in the slow equation (<ref>) is independent of the $\xx$-variables, indicating that the noise correction can be additively incorporated and implemented by examining the decoupled hidden process. Note that the functional form of  $\CCC ^{\yy}$ implies that the WL parametrization cannot be exact in $\epsilon$, as noted in
Eqs. (<ref>) and (<ref>). This indicates that the WL reduced model will not only introduce an error in averaging over the decoupled steady state, but also that the Dyson expansion Eq. (<ref>) has to be truncated at $\epsilon ^{3}$, rather than merely at $\epsilon ^{2}$ where no memory effects would be included.

According to the WL parametrization discussed in Sect. <ref>, the fluctuation terms correspond to the decoupled
evolution of the coupling function $\CCC ^{\xx}_{\yy}$, concretely as in Eq. (<ref>). This allows to directly compute the correlation function:

\begin{align}\label{autocorrelation function wl}
\bigg \langle  \CCC ^{\xx}_{\yy}(\yy)\CCC ^{\xx}_{\yy}(\yy (t))^{\top} \bigg \rangle = \begin{bmatrix} L_{13}^2e^{-(\gamma_1/h)t}\frac{\sigma_1^2}{2\gamma_1} + c_{134}^2e^{-(\gamma_1+\gamma_2)t/h}\frac{\sigma_1^2\sigma_2^2}{2\gamma _1 \gamma_2}& 0 \\ 0 & L_{24}^2e^{-(\gamma_2/h)t}\end{bmatrix},
\end{align}

From Eq. (<ref>) we deduce that the noise covariance matrix for the given parameter values is given by:

\begin{equation}\label{covariance matrix wl}
\bigg \langle  \CCC ^{\xx}_{\yy}(\yy)\CCC ^{\xx}_{\yy}(\yy)^{\top} \bigg
\rangle = \begin{bmatrix}
0.2578... & 0 \\ 0 & 0.5
\end{bmatrix}
\end{equation}

The memory kernel $\mathcal{K}$, which is a vector of two components $(\kappa _1, \kappa _2)^{\top}$, is given by:

\begin{align}
\mathcal{K}(s,\xx)=\begin{bmatrix} \kappa _1(s,\xx) \\ \kappa_2(s,\xx) \end{bmatrix} = \bigg \langle  \CCC^{\yy}(\xx , \yy)\cdot \nabla_{\yy} \CCC^{\xx}(\xx (s),\yy (s)) \bigg \rangle;
\end{align}

here the brackets $\left \langle \cdot  \right \rangle$ indicate the averages for the uncoupled equilibrium in the $\yy$-variables, which happen to be a set of independent OU processes. Explicitly,

\begin{align}
\kappa_1(s,\xx) =&  \bigg \langle \left( -L_{13}x_1+c_{341}y_2x_1 \right)\partial _{y_1}\left( L_{13}y_1(s) + c_{134}y_1(s)y_2(s) \right) \bigg
\rangle \\
&-\bigg \langle \left( L_{24}x_2 + c_{413}y_1x_2 \right) \partial _{y_2}\left( L_{13}y_1(s) + c_{134}y_1(s)y_2(s) \right) \bigg  \rangle \\ =& -L_{13}^2e^{-(\gamma _1/h)s}x_1 + c_{341}c_{134}e^{-\left(\gamma_1 + \gamma_2\right)s/h}\frac{\sigma _2^2}{2\gamma_2}x_1+ c_{134}c_{341}e^{-\left(\gamma_1 + \gamma_2\right)s/h}\frac{\sigma _1^2}{2\gamma_1}x_1, \\
%% and \begin{align}
\kappa_2(s,\xx) =&  \bigg \langle \left( -L_{13}x_1+c_{341}y_2x_1 \right)\partial _{y_1}\left( L_{24}y_2(s) \right) \bigg  \rangle
-\bigg \langle L_{24}x_2 \partial _{y_2}\left( L_{24}y_2(s) \right) \bigg
\rangle \\ =& -L_{24}^2e^{-\left(\gamma _2/h\right)s}x_2.
\end{align}

The reduced-order model obtained herewith does give explicit formulas for
the evaluation of the stochastic noise and the memory kernel, independently of the time scale separation $h$, although these formulas are rather complicated. Still, the scheme remains the same when changing parameter values, so it is flexible in studying different scenarios.

§.§ EMR Model and Results

Basic EMR algorithm implementation.
Regarding the data-driven EMR protocol, we integrated the full model with
a time step of $\dd_{\ell} t = 10^{-3}$ time units for a duration of $T_{\ell} = 10^{4}$ time units in order to learn the model parameters. Then, a separate run was performed in order to examine the ability of the inferred model to reproduce the general statistical features. This time, the EMR system was integrated together with the full model using a time step of $\dd_{\tau} t = 10^{-2}$ time units for a total of $T_{\tau} = 10^5$ time units. The equations were solved using a fourth-order Runge-Kutta and a Euler-Maruyama method for the deterministic and stochastic components, respectively.

By sampling every time step, we learned an EMR model whose coefficients were explicitly found. The convergence criterion $R^2 \simeq 0.5$ was attained by adding two extra levels., for a total of three. The convergence was not affected by changes in the time scale separation parameter — namely, $h=0.1$ and $h=1$ in the case at hand. Probably the value of $h$
was not that important here because of the low dimensionality and stochastic nature of the hidden process. However, convergence is likely to be altered in more complicated models, as illustrated in Appendix <ref>.

The climatologies of the slow $\xx$-variables are obtained using data from the full model, the EMR model and the WL parametrization. The two-dimensional probability density functions (PDFs) of the stochastic model (<ref>) in the $(x_1, x_2)$-plane are shown in Fig. <ref>. These PDFs were calculated by employing the Matlab R2019a kernel smoothing function ksdensity. Their respective marginals are shown in Fig. <ref>. The agreement between the two methodologies when approximating the clearly non-Gaussian density arising
from the full model is clearly excellent.

The time scale separation between the $\xx$-variables and the $\yy$-variables is clearly depicted in the left panel of Fig. <ref>, where the fast variables decorrelate almost instantly compared to the slow ones. The approximation of these autocorrelation functions is also obtained using the EMR and WL methods.

Full model

EMR model

WL parametrization

Two-dimensional probability density functions (PDFs) of the stochastic
model (<ref>) in the $(x_1, x_2)$-plane, as obtained with: (a) the full integration; (b) an integration of the EMR model; and (c) the WL parametrization. The time scale separation parameter used is $h = 0.1$. The PDFs shown here and in Fig. <ref>
were obtained by using the Matlab R2019a kernel smoothing function ksdensity.

PDFs of  (a) the $x_1$ variable; and (b) the
$x_2$ variable.
The separation parameter is $h = 0.1$ and colors used for each method are indicated by the legends inside
the panels.

Autocorrelation functions for the four variables
$x_1, x_2, y_1, y_2$ obtained (a) from the full model;
and (b) the comparison of the corresponding results for $x_1, x_2$ with the full model,
the EMR model, and the WL parametrization. See legend for the choice of lines; $h = 0.1$.

In general, cf. <cit.>, the regressions performed in the main level (<ref>) of the EMR model allow one to effectively reconstruct the coefficients of a weakly coupled model; see Appendix <ref>.
The EMR methodology, though, only allows for linear coupling between the slow $\xx$'s and the fast $\yy$'s. The nonlinear coupling between the slow and fast variables in system (<ref>) compromises the
estimation of the main model parameters in Eq. (<ref>), so that we cannot expect to recover the original, full model's behavior given by (<ref>). The EMR model coefficients at the first and second levels are as shown on tables <ref> and <ref>, respectively.

Empirically estimated EMR model coefficients at the first level,
Eq. (<ref>), for $h = 0.1$. First column gives the coefficients for the constant forcing $\mathbf{f}^{(0)}$, the second and third columns indicate the linear component of the vector field $\mathbf{b}^{(0)}$ and the last three columns determine the quadratic form $\mathbf{A}$.

$f$     $x_1$     $x_2$     $x_1^2$     $x_1x_2$     $x_2^2$

-0.31404     -0.50954     -0.065313     0     -1.0092     0.99704

-0.15356     0.12353     0.21979     1.0092     -0.99704     0

Empirically estimated EMR model coefficients at the second
level, Eq. (<ref>), for $h = 0.1$. First column gives the coefficients for the constant forcing $\mathbf{f}^{(1)}$, the next two columns indicate the linear coupling to the main level (i.e. the first two columns of $\mathbf{b}^{(1)}$) and the last two columns determine the linear drift for the second level (i.e. the last two columns of $\mathbf{b}^{(1)}$).

$f^{(1)}$     $x_1$     $x_2$     $r^{(1)}_1$     $r^{(1)}_2$

0     -5.6397e-4     -6.9382e-05     -20.2823     0.016165

0     -1.648e-4     -4.72e-4     0.086621     -10.0426

As discussed in Sect. <ref>, the EMR has the structure of an MSM and it can be recast into an integro-differential equation. If one only considers the first added level, the EMR can be readily integrated giving the following equation for the evolution of the slow variables $\xx = (x_1,x_2)$:

\begin{equation}\label{intdiff equation}
\dot{\xx}(t) = \FFF (\xx(t)) + e^{-\mathrm{D}t}\yy (0) + \int _{0}^te^{-\mathrm{D}(t-s)} \Sigma \dd \mathbf{W}_s + \int _{0}^{t}e^{-\mathrm{D}(t-s)}\mathrm{C}\xx(s)\dd s.
\end{equation}

Here $\mathbf{W}_s$ is an independent two-dimensional Wiener process and

\begin{align}
-19.9982 & 2.1122\cdot 10^{-3}  \\
-0.77528 & -10.116
\end{bmatrix}, \mathrm{C}=\begin{bmatrix}
-5\cdot 10^{-3} & -5\cdot 10^{-4}\\
-1\cdot 10^{-3} & -5\cdot 10^{-3}
\end{bmatrix}, \label{eq:D+C}
\\
&\Sigma =\begin{bmatrix}
0.2626 & -0.0014  \\
-0.0014 & 0.5013
\end{bmatrix}. \label{eq:covar}
\end{align}

First thing to note is that the matrix $\mathrm{C}$ has a small norm and,
by virtue of Eq. (<ref>), it means that memory effects are going to be very small. On the other hand, the eigenvalues of the matrix $\mathrm{D}$ are $\lambda _1 \simeq -20, \lambda _2 \simeq -10$, which
are approximatelly the drift coefficients of the uncoupled OU process driving the $\yy$-variables. This indicates that the exponential kernel is damping the effects of the $\xx$-variables in past times rather quickly.
Moreover, the covariance matrix $\Sigma$ corresponds to that obtained by integrating the $\yy$-dynamics independently, according to Eq. (<ref>).

Regarding the WL approximation, we stress that the $\yy$-variables are no
longer present, after taking the averages in its construction. The memory
kernel $\mathcal{K}$ in this case differs from the matrix $\mathrm{D}$ above, although its dominant terms correspond to its eigenvalues. Note that, if $L_{13} = 0$, the coupling function $\CCC^{\xx}_{\yy}$ would project entirely onto the eigenfunction of the OU process associated with the eigenvalue $-(\gamma _1 + \gamma _2)$. The same statement would hold for $c_{134}=0$, where in this case $\CCC^{\xx}_{\yy}$ projects onto the eigenfunctions associated with the eigenvalue $-\gamma_1$.

Reduced Time Scale Separation. The parameter $h$ controls the
time scale separation in the evolution of the $\xx$- and $\yy$-variables.
Here we set $h=1$ so that this separation is reduced by an order of magnitude as illustrated by the autocorrelation functions in Fig. <ref>. The question to be addressed in this subsection is the effect of such a reduction in the WL and EMR parametrizations and their respective performance.

In the WL parametrization, there is no need to sample the dynamics in order to construct it, since the formulas of Sect. <ref> are explicit and do not depend on $h$. In the case of model (<ref>), the covariance matrix and time correlations of the WL noise correction are thus given by Eq. (<ref>) with no reference to $h$.  The memory term, though, is expected to change as the kernel $\mathcal K$ will decay more slowly by a factor of $10$. Therefore, memory effects are more important, as expected.

The EMR approach, on the other hand, requires a new learning phase for this value of $h=1$. We used the same numerical integration parameters $\dd_{\ell} t = 10^{-3}$ and $T_{\ell} = 10^{4}$ time units as for the previous case. In the first level regression, one observes that the coefficient values listed in Table <ref> are essentially the same from those estimated in the previous case, for $h = 0.1$, and listed in Table <ref>, as well as being even more distant from the original ones in the table's first row. The second level coefficients, for $h=1$, are shown on Table <ref>.

The covariance matrix $\Sigma$ of the noise correction is indicated in
Eq. (<ref>)
Eq. (<ref>) and it agrees fairly well with the previous values, for $h = 0.1$, as given in Eq. (<ref>). The matrix  $\mathrm{C}$ that indicates the strength of the memory effects also has a magnitude that is of the same order as that in the previous case of $h = 0.1$, which is rather suprising, given the factor of $10$ in time scale separation $h$; compare Eqs. (<ref>) and (<ref>).
This observation tells us that the loss of Markovianity might be intrinsic to the nature of the coupling rather than being due to the time scale separation, even though, in the limit case of infinite scale separation, memory effects will dissapear entirely. The memory kernel, as determined by $\mathrm{D}$, scales almost exactly with the time scale separation and it is expected to change depending on how the coupling functions project onto the eigenspaces of the underlying Orstein-Uhlenbeck process, as discussed more generally earlier in Theorem <ref>.

The performance of both parametrization techniques is summarized in Figs. <ref>–<ref>. These figures are the exact counterparts of Figs. <ref>–<ref> for the reduced time scale separation $h = 1$. First of all, we note from Fig. <ref>(a) that, for $h =
1$  there is indeed no strict time scale separation, as indicated by the autocorrelation functions obtained from the full model. Secondly, consideration of Figs. <ref>(a)–(c), <ref>(a,b) and <ref>(b) shows that neither the WL nor the EMR approach seems to be affected by the time scale reduction.

Full model

EMR model

WL parametrization

Two-dimensional smoothed PDFs of the stochastic model (<ref>),  but with a time scale separation of $h = 1$. Panels (a), (b)
and (c) are calculated by integrating the full model, the EMR model and WL approximation, respectively,  as in Fig. <ref>.

PDFs of (a) the $x_1$ variable and (b) the $x_2$ variable, for a time scale separation of $h=1$; compare with Fig. <ref>.

Autocorrelation functions for
the four variables
$x_1, x_2, y_1, y_2$ obtained (a) from the full model;
and (b) the comparison of the corresponding results for $x_1, x_2$ with
the full model,
the EMR model, and the WL parametrization. See legend for the choice of
lines; $h = 1$.

Empirically estimated EMR model coefficients at the first level, for $h = 1.$

$f$     $x_1$     $x_2$     $x_1^2$     $x_1x_2$     $x_2^2$

-0.31181     -0.339     -0.42944     0     -0.93439     0.97958

-0.16925     0.46833     0.17513     0.93439     -0.97958     0

Empirically estimated EMR model coefficients at the second level, for $h = 1$.

$f^{(1)}$     $x_1$     $x_2$     $r^{(1)}_1$     $r^{(1)}_2$

0     -3.9452e-3     -1.5027e-4     -2.1001     0.016509

0     -5.0312e-4     -3.79e-3     0.054861     -1.1214

\begin{align}
-2.1001 & 0.016509  \\
0.054861 & -1.1214
\end{bmatrix} , \mathrm{C}=\begin{bmatrix}
-4\cdot 10^{-3} & -1\cdot 10^{-4} \\
-6\cdot 10^{-4} & -4\cdot 10^{-3}
\end{bmatrix}, \label{cov mat emr epsilon 1} \\
&\Sigma =\begin{bmatrix}
0.2618 & 0.0011  \\
0.0011 & 0.4554
\end{bmatrix}. \label{eq:covar_h2}
\end{align}

§.§ Memory Effects

We would like to end this results section by analyzing the role of the memory effects when performing a reduction of the highly idealized model given by Eqs. (<ref>). For this purpose, we apply the criterion (<ref>) discussed in the corresponding Appendix <ref>. We thus spectrally approximate the autocorrelation functions of the variables $x_1$ and $x_2$ using Eq. (<ref>) with the Koopman operator $\mathcal{T}_{\tau}$ estimated using Ulam's method with a transition time of $\tau=0.5$ time units for the case where $h=0.1$ and $\tau=1$ time unit for $h=1$.

The difference between the two $\tau$-values is due to the fact that we expect the coarse-grained phase space to be sensitive to the system's variability. Hence, if the time scale separation is large, a shorter transition time is requiered in order to capture the influence of the hidden processes. In fact, a range of transition times was tested in the case of $h=1$ to find that the optimal value was $\tau=1$. It follows that the methodology is robust in showing the effects of memory in the projected phase space.

We clearly observe in Figure <ref> that the correlation functions can be accurately reconstructed in the case of large time scale separation $h=0.1$ (figure (a)) but not so for $h=1$ (figure (b)). This indicates, naturally, that memory effects are negligible in the first case and relevant in the second.

Autocorrelation functions for the
two $\xx$-variables obtained from the full model in the spectral reconstruction
using the Kooman operator: (a) $h=0.1$ and (b) $h=1$.

§ CONCLUSIONS

To formulate accurate and efficient parametrizations for multiscale processes is a crucial challenge in many areas of science and technology for one of two reasons: either the numerical simulation of all scales active in a given system is computationally unfeasible; or there is a mismatch between model resolution and the granularity and homogeneity of the observations, as in the case in geophysical flows and in the climate system. Moreover, the construction of parametrizations is instrumental to
help understand the nature of nonlinear fluxes across scales and the physical processes responsible for cascades, instabilities, and feedbacks.

There are two main approaches for constructing parametrizations:  top-down, by deriving the parametrizations directly from the evolution equations
governing the system through the use of suitable approximations; and data-driven, in which the parametrizations are constructed through suitable optimization procedures, which are first tuned in a training phase and then actually used in the prediction phase. Both approaches aim to derive the effective dynamics for the variables of interest: formally, this is achieved by applying the Mori-Zwanzig projection operator <cit.> to the full dynamics. The result of doing so is to describe the impact of the hidden variables by formulating a generalized Langevin equation (GLE) (<ref>) for the variables of interest that includes a deterministic, a stochastic, and a non-Markovian component.

Top-down and data-driven approaches are conceptually complemetary and have different practical advantages and disadvantages. In this paper, we have shown the fundamental equivalence between a top-down and a data-driven approach that have been formulated and applied in the recent literature. This equivalence was illustrated schematically in Fig. <ref>.

We first revisited in Sect. <ref> the
WL parametrization of <cit.>, which relies on an
assumption of weak coupling between the hidden and observed variables, and have extended the previous results by considering more general coupling classes. We have also shown that the perturbative expansion that yields the WL parametrization is exact when the coupling between the hidden
and resolved variables is additive.

The Dyson formalism (<ref>) appears to be essential for computing the effects of the hidden processes on the dynamics of the observed variables, when working at the level of the system's observables. This methodology is explicit, in the sense that no information about the actual coupled process is needed, because the formal computations are performed by considering the limit in which no coupling is present. Other advantages of this approach are that it can be implemented without the need for any hypothesis on the time scale separation between the hidden variables and the observed ones, and that it is also scale adaptive <cit.>.

We addressed systematically the problem of re-Markovianizing the WL memory equation, which was first pointed out in <cit.> and discussed further in Sect. <ref> here. In this example, a system (<ref>) with one observed and two hidden variables that yielded a scalar WL parametrization was re-Markovianized to a Markovian system with just two scalar differential equations. Throughout Sect. <ref>, we provided a broader framework for re-Markovianization. This framework was presented for a scalar equation in Theorem <ref> and described for higher dimensional systems in Remark <ref>. The required assumptions for this treatment boil down to certain
spectral properties of the Koopman operator for the hidden variables.

The multilevel structure of the re-Markovianization obtained in Sect. <ref> motivated the comparison with multilayer stochastic models (MSMs) in Sect. <ref>. Such MSMs arise naturally in data-driven reduction methods and they had been shown in <cit.> to approximate the GLE predicted by Mori <cit.> and Zwanzig <cit.>.

We showed in Sect. <ref> that a seemless application of the WL parametrization solves the MSM of Eq. (<ref>) and coincides with its Itô integration, cf. <cit.>. Note that an MSM can be obtained from partial observations of the coupled system, which amounts to the special case of the data-driven empirical model reduction (EMR) methodology <cit.>.

The EMR methodology was revisited here in Sect. <ref> and it is, in principle, dual to the WL parametrization, in the sense that only partial observations of the coupled system are required, without the need for
knowing the actual equations of motion. Comparing the multilevel structure of Eq. (<ref>) with that of Eq. (<ref>) suggests that the Koopman eigenvalues $\lambda_j$ highlighted in
Theorem <ref> may help provide insights into the number of
levels needed for EMR to converge. This practical role of the $\lambda_j$'s deserves, therewith, a more careful examination
in further work.

Additionally, we considered in Sect. <ref> a conceptual climate model to which we applied both of the methodologies revisited herein. Since both the MSM and the WL parametrization yield a memory equation that involves integrals and stochastic noise, we were able to compare their structure, as well as their statistical outputs. We found that both methodologies produced equivalent numerical results and that the memory kernel and
noise predicted in the WL parametrization agreed with what was found using the data-driven EMR approach.

Concluding, our view point is complementary to the dynamic mode decomposition <cit.> as it uses the basis of eigenvectors of the Koopman operator to construct the projected — in the sense of Mori-Zwanzig — dynamics of the observables of interest, which is then recast in Markovian form using the multilevel Markovian model framework, where the number of levels corresponds to the number of eigenvectors
of the Koopman operators one considers in the reconstruction of the dynamics. In a nutshell, our findings support, on the one hand, the physical basis and robustness of the EMR methodology and, on the other hand, illustrate the practical relevance of the WL perturbative expansion used for deriving the parametrizations.

§ ACKNOWLEDGEMENTS

It is a pleasure to acknowledge useful exchanges with Georg Gottwald, Jeroen Wouters and Gabriele Vissio. VL acknowledges the support received from the European Union's Horizon 2020 research and innovation program through the project CRESCENDO (Grant Agreement No. 641816). This article is TiPES contribution # 56; the TiPES (Tipping Points in the Earth System) project has received funding from the European Union’s Horizon 2020
research and innovation program under Grant Agreement No. 820970. This research was partially supported by the Israeli Council for Higher Education (CHE) via the Weizmann Data Science Research Center, by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 810370), and by the Office of Naval Research (ONR) Multidisciplinary University Research Initiative (MURI) grant N00014-20-1-2023. This paper has also been partially supported by the EIT Climate-KIC; EIT Climate-KIC is supported by the European Institute of Innovation & Technology (EIT), a body of the European Union, under Grant Agreement No. 190733 and by the Russian Science Foundation (Grant No. 20-62-46056).

§ DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author
upon request.

§ PROOF OF THEOREM <REF>

The aim is to show that under the assumptions of this Theorem — which require that the coupling function $\CCC^{\xx}:\mathbb{R}\longrightarrow\mathbb{R}$ projects entirely onto $\mathrm{span}\left\{\psi_{j}, \; j=1,\ldots,N \right\}$ — the memory and noise terms of the WL equation (<ref>)
are obtained from the term $\epsilon \Lambda\cdot \mathbf{Z}(t)$ in Eq. (<ref>), after  integration of Eq. (<ref>).

Step 1. In this step, we expand the memory term and the lag-correlations of the noise in
the WL equation (<ref>)  in terms of the leading eigenelements of the uncoupled Koopman operator $\LLL_{0}^{\yy}$.
These expansions will serve us in Step 2 below to compare the noise and memory terms of the WL-equation with those produced after integration of Eq. (<ref>).

Let the $\lambda_j$ be sorted as in the hypotheses of the theorem. Hence, to distinguish real and complex conjugate eigenvalues and for notational convenience, we introduce the set of indices  $I_r$ and $I_{+}$ defined

\begin{align}
I_r&=\{j \in \{1,\ldots,N\} : \lambda_j \text{ is real} \},\\
I_+&=\{j \in \{1,\ldots,N\} : \mathfrak{Im}\lambda_j>0\}.
\end{align}

An immediate consequence that is used several times throughout the proof
is that the sum of the eigenvalues is real, and may be split  as follows:

\begin{equation}
\sum_{j=1}^{N}\lambda_j=\sum_{j\in I_r}\lambda_j + \sum_{j\in I_+}\lambda_j+ \sum_{j\in I_+}\overline{\lambda_j}=\sum_{j=1}^{N}\mathfrak{Re}\lambda_j \in \mathbb{R}.
\end{equation}

As previously stated, we expand the mean and correlation functions of the scalar noise term $\eta$ in the WL-equation (<ref>) in terms of the eigenpairs. The mean is zero by assumption, but the autocorrelation function can be expanded as follows,
based on Eq. (<ref>),

\begin{align}\label{proof wl noise 3}
\bigg \langle \eta(t)\eta(0) \bigg \rangle = C_{\CCC^{\xx},\CCC^{\xx}}(t)   =\sum_{j=1}^{N}e^{\lambda_jt}\alpha_j\beta_j;
\end{align}

herein, $\alpha_j$ and $\beta_j$ are as defined in (<ref>) and (<ref>), respectively.
The expansion of the correlation function in Eq. (<ref>) is a finite sum by virtue of the assumption that $\CCC^{\xx}$ lies in $\mathrm{span}\left\{\psi_{j}, \; j=1,\ldots,N \right\}$ and therefore there is no contribution from the essential spectrum.

Regarding the complex scalars $\alpha_j$ and $\beta_j$ defined (<ref>) and (<ref>), it follows that for each $j$ such that $\lambda_j=\overline{\lambda_{j+1}}$, we get $\alpha_j=\overline{\alpha_{j+1}}$ and $\beta_j=\overline{\beta_{j+1}}$. Indeed,

\begin{align}
\alpha_j &= \int \nu(\dd \yy)  \overline{\psi_j^{\ast}(\yy)}\CCC^{\xx}(\yy) = \overline{\int \nu(\dd \yy)  \psi_j^{\ast}(\yy)\CCC^{\xx}(\yy)} \\&= \overline{\int \nu(\dd \yy)  \overline{\psi_{j+1}^{\ast}(\yy)}\CCC^{\xx}(\yy)} = \overline{\alpha_{j+1}},
\end{align}

in which we have exploited the fact that $\psi_j=\overline{\psi_{j+1}}$ when $\lambda_{j}$ is complex and $j$ in $\{1,\ldots, N\}$.

The same proof can be repeated for $\beta_j$. Such a conjugacy relation also holds for the gradients of the eigenfunctions $\nabla  \psi_j$, for those $j$ in $\{1\ldots,N\}$ such that $\lambda_j=\overline{\lambda_{j+1}}$. This is observed by the following equality:

\begin{align}\label{Eq_toto1}
\nabla \psi_{j}(\yy)=\nabla \overline{\psi_{j +1}}(\yy)=\overline{\nabla\psi_{j+1}(\yy)},
\end{align}

since $\nabla$ is a differential operator that only involves here differentiation with respect to a real variable.

Exploiting (<ref>), the memory kernel $\mathcal{K}$ then expands as (recalling that $\mathcal{R}(t)=0$):

𝒦(t,s,x)    = ^(x(s))·⟨∇∑_j=1^Ne^λ_j (t-s)α_j ψ_j () ⟩

=^(x(s))·∑_j∈I_re^λ_j (t-s)α_j⟨∇ψ_j ()  ⟩+^(x(s))·∑_j∈I_+e^λ_j (t-s)α_j⟨∇ψ_j ()  ⟩

+^(x(s))·∑_j∈I_+e^λ_j (t-s)α_j⟨∇ψ_j ()  ⟩
which leads to

𝒦(t,s,x)=^(x(s))·∑_j=1^Nℜ𝔢(e^λ_j (t-s)α_j⟨∇ψ_j () ⟩).
Note that to go from (<ref>) to (<ref>),  we have made use of the aforementioned conjugacy
relations, that led to a real-valued  memory kernel at the end. With (<ref>) and (<ref>) at hand, the noise and memory terms in
the WL equation (<ref>) are thus characterized in terms of the leading eigenelements of the uncoupled Koopman operator $\LLL_{0}^{\yy}$.

Step2. The second step consists of analyzing the noise and memory terms produced by integration of  Eq. (<ref>) and to compare these terms with those of the WL-equation.

Performing an Itô integration of Eq. (<ref>) leads

\begin{equation}\label{eq: ito integral proposition}
\mathbf{Z}(t)=e^{\mathrm{D}t}\mathbf{Z}(0)+ \underbrace{ \int _{0}^{t}e^{\mathrm{D}(t-s)}\Sigma \dd W_s}_{\textrm{noise term}} + \epsilon\underbrace{\int _{0}^te^{\mathrm{D}(t-s)}\mathbf{R}(x(s))\dd s}_{\textrm{memory term}},
\end{equation}

where (for simplicity) we have assumed that the initial condition distributes normally with mean zero and variance equal to the identity matrix,
and the function $\mathbf{R}:\mathbb{R}\longrightarrow \mathbb{C}^{N}$ is as defined  in (<ref>).
The noise and memory contributions of  $\mathbf{Z}(t)$ are as indicated by the brackets in (<ref>).

Let us denote the noisy component of Eq. (<ref>) as $\mathbf{q}(t)$ in $\mathbb{R}^N$. Then, it is clear that $\mathbf{q}$ has zero mean and the lag cross-correlations read as:

\begin{align}
\mathbb{E}\left( \mathbf{q}(t)\mathbf{q}^{\top}(0) \right) =e^{\mathrm{D}t}
= \begin{bmatrix}
e^{\lambda_1t} &&\\ &\ddots&\\&&e^{\lambda_Nt}
\end{bmatrix}.
\end{align}

Let $\Lambda$ be defined as in (<ref>) and let us calculate the mean and lag-correlations of the one-dimensional stochastic process $\Lambda\cdot\mathbf{q}$, aimed at approximating $\eta$ in the WL equation.
First, note that
$\Lambda\cdot \mathbf{q}$ is a zero-mean Gaussian process, with lag correlations given by:

\begin{align}\label{eq: proof correlation A Z}
\mathbb{E}\left( \left(\Lambda\cdot \mathbf{q}(t)\right)\left(\Lambda\cdot \mathbf{q}(0)\right)  \right) = \left(e^{\mathrm{D}t}\Lambda\right)\cdot \Lambda.
\end{align}

Now, expanding  $\left(e^{\mathrm{D}t}\Lambda\right)\cdot \Lambda$ in (<ref>) shows that  we recover
the right-hand side (RHS) of (<ref>).

However, the noise term $\eta$ in the WL-equation is a real-valued stochastic process, and we are dealing with complex scalars, so therefore
we still have to show that $\Lambda\cdot \mathbf{q}(t)$ is real for every $t$ in $\mathbb{R}$.  To do so, let us denote by $\mathbf{w}^{\top}=[w_1,\ldots ,w_N]$ any arbitrary row vector in $\mathbb{R}^{N}$.  In other
words,  $\mathbf{w}$ is an arbitrary column vector with real entries.

Consider the following inner product:

\begin{align}
\Lambda\cdot e^{\mathrm{D}t}\Sigma \mathbf{w} &= \Lambda\cdot e^{\mathrm{D}t}\begin{bmatrix}\sqrt{-2\mathfrak{Re}\lambda_1}&& \\ &\ddots & \\
\end{bmatrix}\mathrm{H} \mathbf{w} \\
&=\Lambda\cdot \begin{bmatrix}e^{\lambda_1t}\sqrt{-2\mathfrak{Re}\lambda_1}&& \\ &\ddots  & \\ &&e^{\overline{\lambda_2}t}\sqrt{-2\mathfrak{Re}\lambda_N}
\end{bmatrix}\mathrm{H} \mathbf{w}.
\end{align}

By construction of the matrix $\mathrm{H}$ in Eq. (<ref>)-(<ref>), the product $\mathrm{H}\mathbf{w}$
is given component-wise, for $j=2,\ldots,N$, as:

\begin{equation}
\left[\mathrm{H}\mathbf{w}\right]_j=\begin{cases}w_j, \text{ if } j\in I_r \text{ or, } j\in I_+, \\  w_{j-1}, \text{ if }\lambda_j=\overline{\lambda_{j-1}}, \end{cases}
\end{equation}

while $\left[\mathrm{H}\mathbf{w}\right]_1=w_1$. This implies in particular that $\left[\mathrm{H}\mathbf{w}\right]_j=\left[\mathrm{H}\mathbf{w}\right]_{j+1}$ whenever $\lambda_j=\overline{\lambda_{j+1}}$. As a consequence, we get
Λ·e^DtΣ𝐰=   ∑_j∈I_rα_j^1/2β_j^1/2e^λ_jt√(-2λ_j)w_j+∑_j∈I_+α_j^1/2β_j^1/2e^λ_jt√(-2ℜ𝔢λ_j)w_j

which shows that $\Lambda\cdot e^{\mathrm{D}t}\Sigma$ is a real-valued quantity, and thus  for any realization of the $N$-dimensional Wiener process $W_t$, the product $\Lambda \cdot e^{\mathrm{D}(t-s)}\Sigma \dd W_s$ is real and hence, $\Lambda \cdot \mathbf{q}(t)$ is also real for every $t$.

Finally, we are left with showing that the memory kernel $\mathcal{K}$ of the WL-equation  coincides with that of the memory term in Eq. (<ref>), when multiplied by the vector $\Lambda$. To
do so, we exploit the expansion (<ref>) of
$\mathcal{K}$ for this comparison, namely using the expression of $\mathbf{R}$ in (<ref>), we observe that
Λ·∫_0^te^D(t-s)𝐑(x(s))s    =^(x(s))·∑_j∈I_re^λ_j (t-s)α_j⟨∇ψ_j ()  ⟩

+^(x(s))·∑_j∈I_+e^λ_j (t-s)α_j⟨∇ψ_j ()  ⟩+^(x(s))·∑_j∈I_+e^λ_j (t-s)α_j⟨∇ψ_j ()  ⟩

=    ^(x(s))·∑_j=1^Nℜ𝔢(e^λ_j(t-s)α_j⟨∇ψ_j ()  ⟩),
which indeed coincides with the expression of  $\mathcal{K}$, as desired.  The proof is complete.

§ SEMIGROUP PROPERTY OF THE PROJECTED KOOPMAN OPERATOR FAMILY

It was shown in <cit.> that projection onto a reduced state space is closely related with a coarse graining of the (full) probability transitions on the original system's attractor, while <cit.> dealt recently with the impact of such a projection in terms of reduction of the Koopman semigroup.
In <cit.>, the authors proposed a criterion based on the spectral theory of Markov semigroups to ascertain whether the reduced state
space associated with a given projection can fully explain the statistics
of the desired variables. This approach provides potential insights into the need for modeling non-Markovian effects by inspecting the loss of the
semigroup property, as explained below.

Moreover, it follows from <cit.> that the analysis of correlation functions is not only of physical interest but also of methodological utility. Correlation functions can be defined by means of the Koopman
operator or, dually, by means of the transfer operator's providing the solution of the Liouville equation.

Let $\mu$ denote an ergodic invariant measure of the system and take two observables $\Phi_1,\Phi_2$ in the space $L^2_{\mu}$ of zero-mean functions that are square-integrable with respect to $\mu$. Assume furthermore that the spectrum of the operator $\LLL$ in  $L^2_{\mu}$ is a pure point spectrum, given by the eigenvalues $\{\lambda _j\}_{j=1}^{\infty}$ and their associated eigenfunctions $\{\psi_j\}_{j=1}^{\infty}$, where the eigenvalues are ordered by their decreasing real parts.

Then, the correlation function $C_{\Phi_1,\Phi_2}(t)$ between the functions $\Phi_1$ and $\Phi_2$ is given by:

\begin{equation}\label{spectral correlations0}
C_{\Phi_1,\Phi_2}(t) = \int  \Phi_1\cdot e^{t\LLL }\Phi_2 \dd \mu = \int  e^{t\mathcal{L}^{\ast}}\Phi_1\cdot \Phi_2 \dd \mu
\end{equation}

and it can be expanded, formally, as

C_Φ_1,Φ_2(t) = ∑_j=1^∞e^λ_jt⟨Φ_1,ψ_j ⟩_μ⟨ψ_j^∗,Φ_2 ⟩_μ.
The dual operators in (<ref>) and the adjoint eigenvectors  in (<ref>) are indicated by the superscript $(\cdot)^\ast$, while $\langle \cdot, \cdot \rangle _{\mu}$ denotes the inner product. We refer to <cit.> for a proof of (<ref>) in the context of Markov semigroups. The proof actually applies to the case of the Koopman semigroups considered here as long as the Koopman semigroup $U_t$ defined by (<ref>) is a strongly continuous semigroup in $L^2_\mu$. The  RHS of Eq. (<ref>) consists of a linear combination of exponential
terms whose coefficients are calculated by projecting $\Phi_1$ and $\Phi_2$ onto the corresponding eigenspaces. These coefficients weight each exponential function and they can become exceedingly large if the Koopman operator deviates very much from normality <cit.>.  Note also
that the set of eigenvalues $\lambda_j$'s play a key role in defining the
response of the system to perturbations <cit.>.

The interactions between the resolved and hidden variables that are modeled by the Dyson expansion of the Koopman operator in Sec. <ref> may introduce memory effects into the closed, reduced model for the $\xx$-variables, as given by Eqs. (<ref>)–(<ref>). In certain situations, such memory effects can be neglected, even in the absence of exact slaving relationships between the resolved and hidden variables <cit.>. But the loss of slaving
relationships requires, in general, an explicit representation of memory effects <cit.> to achieve an efficient model reduction.

Furthermore it was shown in <cit.> that the reduction of the Koopman semigroup to observables that act only on the reduced state space leads, in most circumstances, to a family of operators that, while Markovian, no longer satisfy the semigroup property. One might
then ask to which extent this loss of the semigroup property arising from
the reduction, and the related emergence of memory effects, are crucial for providing a faithful reduced model of the observed variables.

When considering reduced state spaces obtained by projection, along with observables $\Phi_1$ and $\Phi_2$ defined on them, Theorem 2 in <cit.> shows the existence of a family of Markov operators $\{\mathcal{T}_t\}_{t\geq 0}$ that satisfies:

\begin{equation}\label{criterion memory}
\int \Phi_1 \cdot \mathcal{T}_t\Phi_2\dd \mu_{\xx} = \int _{0}^t[\Phi_1\circ \pi_{\xx}]\cdot e^{\LLL t}[\Phi_2\circ \pi _{\xx}] \dd \mu = C_{\displaystyle{(\Phi_1\circ \pi_{\xx},\Phi_2 \circ \pi_{\xx}})}(t),
\end{equation}

for every $t\geq0$, where ${\displaystyle{\pi_{\xx}}}$ is the canonical projection onto the reduced subspace and $\mu_{\xx}$ is the disintegrated or sample measure associated with $\pi_{\xx}$; see <cit.>. However, due to the projection, the semigroup property is lost, namely, $\mathcal{T}_{s}\mathcal{T}_t\neq \mathcal{T}_{t+s}$ for some $t,s$.

Following the reasoning above, one can establish a criterion for the need
to model a memory contribution when performing the model reduction. Formally, if there exist $\tau > 0$ and $T\in \mathbb{N}$ such that, for every
$t\in \{k\tau \in \mathbb{R} : 0\leq k \leq T \}$, we have

\begin{equation}\label{criterion memory effects}
C_{\displaystyle{(\Phi_1\circ \pi_{\xx},\Phi_2 \circ \pi_{\xx}})}(t) = \int \Phi_1\cdot \mathcal{T}_{t}\Phi_2\dd \mu = \int \Phi_1\cdot \left(\mathcal{T}_{\tau}\right)^k\Phi_2\dd \mu,
\end{equation}

and one can say that the semigroup is preserved, to some extent, depending on how large $T$ can be in Eq. (<ref>) above. Other such criteria are available in the context of mutually dual Koopman
and transfer operators. Thus, A. Tantet and coauthors <cit.> had already considered empirical ways of quantifying the loss of the semigroup property in reduced dimensions.

The interpretation of $\tau$ comes from the practical implementation of the methodology and it is usually referred to as the transition time. Indeed, numerically, the approximation of such Markov operators is done by Ulam's method, by projecting them onto a finite basis, typically using the characteristic functions of certain domains of phase space. Then, the transitions between domains — after an adequate transition time $\tau$ — are counted to obtain matrix estimates of the operator $\mathcal{T}_{\tau}$ acting on the reduced phase space. Hence, one seeks a suitably
small, or large <cit.>, transition time to obtain the best candidate for applying Eq. (<ref>), see also <cit.>. Thus, Eq. (<ref>) can be implemented in practice this way in order to (potentially)
reconstruct correlation functions on the whole phase space. A very simple
illustration of such a transfer operator calculation is given in <cit.>.

§ ITÔ INTEGRATION OF THE MSM

In the main text, we proposed a solution of the MSM given by Eqs. (<ref>) using the Dyson expansion for the linear operators involved in the backward Kolmogorov equation.
advection acting on functions. Therefore, we substituted nonlinear ordinary differential equations for a partial differential equation, for the sake of having linear operators in hand.
The same solution can be attained by direct integration of the MSM in the form (<ref>). We  convolute, in the Itô sense, Eq. (<ref>) to find an explicit solution for $\rr(t)$ when $d_1 = d_2$:

\begin{equation} \label{eq:Ito}
\rr (t) = e^{-\mathrm{D}t}\rr (0) + \int _{0}^te^{-\mathrm{D}(t-s)} \Sigma \dd W_s + \epsilon \int _{0}^{t}e^{-\mathrm{D}(t-s)}\mathrm{C}\xx(s)\dd s;
\end{equation}

here, $\rr(0)$ indicates an initial state that can be assumed to be distributed in a prescribed way. For the more general, nonlinear MSMs considered there, see <cit.>.

By substituting the expression (<ref>) into Eq. (<ref>), we find an exact expression for the evolution of $\xx (t)$:

\begin{equation}\label{eq:ITO}
\dot{\xx}(t) =\FFF (\xx(t)) + \epsilon e^{-\mathrm{D}t}\rr (0) + \epsilon \int _{0}^te^{-\mathrm{D}(t-s)} \Sigma \dd W_s + \epsilon ^2 \int _{0}^{t}e^{-\mathrm{D}(t-s)}\mathrm{C}\xx(s)\dd s,
\end{equation}

in which the memory effects in the fourth term are of second order in $\epsilon$. Note that the $\epsilon$-order terms arise from a noise realization in the decoupled regime, whereas the memory term is exclusively due to the coupling of the main variables with the hidden ones. Hence, the only degree of freedom left is the distribution of the initial states $\rr(0)$.

§ THE COUPLED L84–L63 SYSTEM

The EMR methodology's ability to capture the statistics of low-dimensional dynamical systems was illustrated in <cit.>, where the authors considered the L63 system <cit.> as a test case in which the phase space can be fully sampled. Moreover, provided that the integration time step is short enough, the parameters of the underlying model can be fully captured with a high degree of confidence.

Here, we repeat the analysis of <cit.> to illustrate the effectiveness of EMR in capturing statistical and dynamical properties in an
extended system. The model we consider is the result of coupling the $\X = (X, Y, Z)$ variables of the L84 system <cit.> with the $\x = (x,y,z)$ variables of the L63 system <cit.>, namely:

\begin{align}
\dot X & = -Y^2-Z^2-aX+a(F_0 + hx),\\
\dot Y & = XY - bXZ - Y +G, \\
\dot Z & = XZ+bXY-Z,\\
\dot x & = \tau s (y-x),\\
\dot y & = \tau (\rho x - y -xz),\\
\dot z & = \tau (xy-\beta z);
\end{align}

the parameter values are: $a=0.25,b=4,$ $F_0=8, G=1,$ and $s=10,$ $\rho = 28, \beta =8/3$, respectively. The parameter $h$ measures the strength of the coupling, while $\tau$ scales the rate of change in the L63 system and, therewith, the time scale ratio between the two subsystems.

This system is a skew-product, in the sense of <cit.>, since the coupling is one-way only, with the L63 system driving the L84 dynamics. Hence, one has — as noted in <cit.> — a fully Markovian parametrization of the L63 variables. Furthermore, the correlation function that defines the stochastic noise $\eta(t)$ can be further expanded
and simplified,with respect to Eq. (<ref>).  One can, in fact,
write explicitly:

\begin{equation}
C(\eta (0),\eta(t))=\left \langle (x(0),0,0)\cdot (x(t),0,0) \right \rangle,
\end{equation}

where the angular brackets $\langle \cdot \rangle$ indicate averages with
respect to the physical measure associated with the L63 system. Since L84
does not feed back into L63, the evolution of $\x(t)$ only obeys the dynamics of L63, and thus the decorrelation of the noise scales with $\tau$.

In most of the numerical experiments, the time scale separation between the two systems is $\tau=5$. The relevance of this time scale parameter was investigated in previous work <cit.>. Here, we focus on the effects of the coupling strength $h$, and we shall study the cases of $h=0.25$ and $0.025$. Partial observations only will be used in these experiments, by sampling the three-dimensional outputs of the L84 system. Then, the observed tendencies are regressed and sequentially layered following the EMR approach, as explained in Sect. <ref>.

§.§ EMR Outputs

The L84–L63 model is integrated for 730 time units that correspond in L84 to 10 natural years, with a time step of $5\cdot 10^{-3}$ time units; two separate runs are made for the coupling strengths $h = 0.25$ and $h = 0.025$. These two full-model runs are used to train the corresponding
EMR model versions, both of which use only the slow $\X$-variables and eliminate the fast $\x$-variables. Then, two separate full-model simulations are run, for testing purposes, over 7 300 time units, and the EMR model's output is compared with it, for the two parameter values. Below we show the main statistical outputs of the EMR methodology compared to the two
reference integrations of the full model. The results for the two separate $h$-values are shown in Figs. <ref>–<ref> and Figs. <ref>–<ref>, respectively.

The region of phase space explored by the EMR model clearly coincides with the one visited by the full model, as seen in Figs. <ref> and <ref>, and the relative occupancies within
this region—as indicated by the smoothed PDFs shown in Figs. <ref> and <ref>, respectively—agree very well.  The time scales are also well captured, as indicated by the good approximation of the autocorrelation functions, cf. Figs. <ref> and <ref>.

Notice that, while the original L84–L63 system is purely deterministic, the EMR model includes white noise acting on the hidden layers of the learned model. This fact could suggest that a smoothing of the invariant measure is inevitable and that the EMR methodology may not be able to capture fractal geometries in phase space, since the EMR model does not satisfy the Hörmander's hypoellipticity condition <cit.>. The numerical evidence in Figs. <ref> and <ref>, however, illustrates a strikingly good approximation of the full model's attractor, including its very fine, and presumably fractal, structure.