3.1 Definition of the Koopman operator and its infinitesimal generator

Let $x_{-Q+1},x_{-Q+2},\ldots ,\ldots ,x_{N-1}$ be a sequence of states $x_{n}=(\boldsymbol{u}_{n},T_{n})\in X$ of the Boussinesq system from § 2, sampled at times $t_{n}=n\unicode[STIX]{x1D70F}$ , where $\unicode[STIX]{x1D70F}$ is a fixed sampling interval and $Q,N$ are positive integers. In what follows, $x_{-Q+1}$ will be the first state sampled after the initial equilibration interval in our simulations. For the DNS data analysis, we will take $N=10\,000$ , $Q=30$ and $\unicode[STIX]{x1D70F}=1$ (the latter in units of the free-fall time $T_{f}$ ). It should be noted that we assign the potentially negative timestamp $t_{-Q+1}$ to the first sample in our dataset in anticipation of the fact that we will be performing delay-coordinate maps with $Q$ delays – this ‘uses up’ the first $Q-1$ samples, so we will be left with $N$ available samples for analysis starting at time $t_{0}=0$ .

Formally, we assume that our Boussinesq system admits a compact invariant set $A\subset X$ (e.g. an attractor) with a compact forward-invariant neighbourhood ${\mathcal{U}}$ in which $x_{0}$ lies. We also assume that $x_{0}$ lies in the basin of a physical ergodic probability measure $\unicode[STIX]{x1D707}$ supported on $A$ (Young Reference Young2002). Under this assumption, given any continuous function $f:X\mapsto \mathbb{C}$ (also referred to as observable), its time average $\bar{f}_{N}$ along the orbit $\{x_{n}\}_{n=0}^{N-1}$ converges to its expectation value with respect to the invariant measure, which translates here to

(3.1) $$\begin{eqnarray}\displaystyle \lim _{N\rightarrow \infty }\bar{f}_{N}=\lim _{N\rightarrow \infty }\frac{1}{N}\mathop{\sum }_{n=0}^{N-1}f(x_{n})=\int _{A}f\,\text{d}\unicode[STIX]{x1D707}. & & \displaystyle\end{eqnarray}$$

These assumptions, which cannot be rigorously justified but are tacitly made in many long-time statistical analyses of such systems, are sufficient to ensure that the results of our data-driven algorithms converge in the limit of large data, $N\rightarrow \infty$ (Das & Giannakis Reference Das and Giannakis2017).

The dynamics of our turbulent convection flow in state space can be characterized through a nonlinear evolution map $\unicode[STIX]{x1D6F7}_{t}:X\mapsto X$ , $t\in \mathbb{R}$ , such that $\unicode[STIX]{x1D6F7}_{t}(x)$ is the state reached at time $t$ under dynamical evolution by the Boussinesq system (2.1) with initial data $x\in X$ . Moreover, associated with $\unicode[STIX]{x1D6F7}_{t}$ is a unitary group of Koopman operators $U_{t}:L^{2}(A,\unicode[STIX]{x1D707})\mapsto L^{2}(A,\unicode[STIX]{x1D707})$ , acting on observables in the Hilbert space $L^{2}(A,\unicode[STIX]{x1D707})$ of square-integrable (with respect to $\unicode[STIX]{x1D707}$ ) complex-valued functions on the attractor (Eisner et al. Reference Eisner, Farkas, Haase and Nagel2015). These linear operators act on $f\in L^{2}(A,\unicode[STIX]{x1D707})$ via composition with the flow map, i.e. $U_{t}\,f=f\circ \unicode[STIX]{x1D6F7}_{t}$ , and for all $t\in \mathbb{R}$ and $f\in L^{2}(A,\unicode[STIX]{x1D707})$ the map $t\mapsto U_{t}\,f$ is continuous. As with all such continuous one-parameter unitary groups, $U_{t}$ is fully determined by its infinitesimal generator, $V:D(V)\mapsto L^{2}(A,\unicode[STIX]{x1D707})$ . This operator is an unbounded skew-adjoint operator with a dense domain $D(V)\subset L^{2}(A,\unicode[STIX]{x1D707})$ , acting on observables $f\in D(V)$ via

(3.2) $$\begin{eqnarray}\displaystyle Vf=\lim _{t\rightarrow 0}\frac{1}{t}(U_{t}\,f-f). & & \displaystyle\end{eqnarray}$$

The Koopman operator at any $t\in \mathbb{R}$ can be recovered via exponentiation of $V$ , namely $U_{t}=\exp (tV)$ . Figure 2 summarizes this connection between the linear Koopman operator and the nonlinear evolution (or flow) map.

Figure 2.

The connection between the nonlinear evolution map and the linear Koopman operator in a commutative diagram.

The action of an ergodic dynamical system on observables can be characterized in terms of the spectral properties of Koopman operators $U_{t}$ (Mezić Reference Mezić2005), or, equivalently, the related transfer operators (Dellnitz & Junge Reference Dellnitz and Junge1999) which are adjoint to $U_{t}$ . Another equivalent characterization, which we will adopt here, is through eigenfunctions of the generator $V$ . In particular, a distinguished class of observables are Koopman eigenfunctions, satisfying the eigenvalue problem $V\unicode[STIX]{x1D713}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D713}$ for some $\unicode[STIX]{x1D713}\in D(V)$ and $\unicode[STIX]{x1D706}\in \mathbb{C}$ . For measure-preserving dynamical systems (including the ergodic systems of interest here), every non-zero eigenvalue $\unicode[STIX]{x1D706}=\text{i}\unicode[STIX]{x1D714}$ of $V$ is purely imaginary and captures an intrinsic frequency $\unicode[STIX]{x1D714}\in \mathbb{R}$ of the system. Moreover, the corresponding eigenfunctions $\unicode[STIX]{x1D713}$ evolve periodically under the dynamics at the frequencies associated with the eigenvalues; that is, $U_{t}\unicode[STIX]{x1D713}=\exp (\text{i}\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D713}$ . Because $Vf^{\ast }=(Vf)^{\ast }$ for any $f\in D(V)$ , the Koopman eigenvalues and eigenfunctions occur as complex-conjugate pairs; i.e. if $\unicode[STIX]{x1D713}$ is an eigenfunction at eigenvalue $\unicode[STIX]{x1D706}$ , then $\unicode[STIX]{x1D713}^{\ast }$ is an eigenfunction at eigenvalue $\unicode[STIX]{x1D706}^{\ast }=-\unicode[STIX]{x1D706}$ . For ergodic systems, all eigenvalues of $V$ are simple. Throughout this paper, we will assume that all Koopman eigenfunctions are continuous (to the best of our knowledge, there are no known continuous-time ergodic dynamical systems with discontinuous Koopman eigenfunctions).

Evaluated along the convection flow trajectory in state space given by $x_{0},x_{1},\ldots ,$ each Koopman eigenfunction induces a time series $\unicode[STIX]{x1D713}(x_{0}),\unicode[STIX]{x1D713}(x_{1}),\ldots$ that can be thought of as an analogue of a temporal (chronos) mode in a POD (Aubry, Guyonnet & Lima Reference Aubry, Guyonnet and Lima1991), but with potentially higher dynamical significance. For our purposes, of particular interest are eigenfunctions at small corresponding frequencies $\unicode[STIX]{x1D714}$ . Such functions will correspond to the LSC and other coherent flow features (Mezić Reference Mezić2013).

When dealing with systems with quasiperiodic or chaotic behaviour exhibiting non-isolated eigenvalues and/or a continuous spectrum, the eigenvalue problem for the raw generator $V$ is numerically ill-conditioned. Instead, one may seek approximate Koopman eigenvalues and eigenfunctions associated with a regularized generator

(3.3) $$\begin{eqnarray}L=V-\unicode[STIX]{x1D701}\unicode[STIX]{x1D6E5},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}>0$ is a small regularization parameter and $\unicode[STIX]{x1D6E5}$ is a positive-semidefinite diffusion operator (Giannakis et al. Reference Giannakis, Slawinska and Zhao2015; Das & Giannakis Reference Das and Giannakis2017; Giannakis Reference Giannakis2017). Under relatively mild assumptions, $L$ has a point spectrum only (Franke et al. Reference Franke, Hwang, Pai and Sheu2010) and its eigenfunctions $\unicode[STIX]{x1D713}$ approximate the eigenfunctions of $V$ as $\unicode[STIX]{x1D701}\rightarrow 0$ . Moreover, the imaginary parts of the eigenvalues $\unicode[STIX]{x1D706}$ of $L$ approximate the Koopman eigenfrequencies $\unicode[STIX]{x1D714}$ . The real part of $\unicode[STIX]{x1D706}$ is non-positive by construction and is related to a Dirichlet energy (roughness) of $\unicode[STIX]{x1D713}$ , as discussed below.

3.2 Data-driven basis from NLSA

In order to solve the eigenvalue problem for $L$ , we first specify an orthonormal basis of $L^{2}(A,\unicode[STIX]{x1D707})$ which can be approximated from data sampled on the orbit $\{x_{n}\}$ . Here, we build our data-driven basis by applying NLSA to the velocity-field snapshots $\boldsymbol{u}_{-Q+1},\ldots ,\boldsymbol{u}_{N-1}$ , given by $\boldsymbol{u}_{n}=F(x_{n})$ , where $F$ is the observation function (projection to velocity-field data) defined in (2.3). Applied to this dataset, NLSA approximates an ergodic compact Markov operator $P:L^{2}(A,\unicode[STIX]{x1D707})\mapsto L^{2}(A,\unicode[STIX]{x1D707})$ with real eigenvalues $\unicode[STIX]{x1D6EC}_{k}$ and a corresponding set of orthonormal eigenfunctions $\unicode[STIX]{x1D719}_{k}$ . This operator is approximated by an $N\times N$ Markov matrix $\unicode[STIX]{x1D64B}=[\unicode[STIX]{x1D617}_{ij}]$ , with $\sum _{j=0}^{N-1}\unicode[STIX]{x1D617}_{ij}=1$ , whose elements are computed through a symmetric kernel function $K:X\times X\mapsto \mathbb{R}_{+}$ operating on delay sequences of velocity-field snapshots. Given an arbitrary pair of states $x,x^{\prime }\in X$ with $x=(\boldsymbol{u},T)$ and $x^{\prime }=(\boldsymbol{u}^{\prime },T^{\prime })$ , we define

(3.4) $$\begin{eqnarray}K(x,x^{\prime })=\text{e}^{-d_{Q}^{2}(x,x^{\prime })/\unicode[STIX]{x1D716}},\quad \text{with }d_{Q}^{2}(x,x^{\prime })=\frac{1}{Q}\mathop{\sum }_{q=0}^{Q-1}\Vert F(\unicode[STIX]{x1D6F7}_{-q\unicode[STIX]{x1D70F}}(x))-F(\unicode[STIX]{x1D6F7}_{-q\unicode[STIX]{x1D70F}}(x^{\prime }))\Vert _{L^{2}}^{2},\end{eqnarray}$$

where $Q$ is a non-negative integer parameter (the number of delays), $\unicode[STIX]{x1D716}$ is a positive kernel bandwidth parameter and $\Vert \cdot \Vert _{L^{2}}$ denotes the norm in velocity-field space $L^{2}(\unicode[STIX]{x1D6FA};\mathbb{R}^{3})$ . Moreover, $d_{Q}$ is a distance function based on data in delay-coordinate space. Intuitively, $K(x,x^{\prime })$ can be thought of as a pairwise measure of similarity between $Q$ -element sequences of velocity fields. In particular, for two states $x_{i}$ and $x_{j}$ in the training data, we have

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D612}_{ij}:=K(x_{i},x_{j})=\exp \left(-\frac{1}{\unicode[STIX]{x1D716}Q}\mathop{\sum }_{q=0}^{Q-1}\Vert \boldsymbol{u}_{i-q}-\boldsymbol{u}_{j-q}\Vert _{L^{2}}^{2}\right). & & \displaystyle\end{eqnarray}$$

In NLSA, the matrix $\unicode[STIX]{x1D64B}$ is formed by performing the normalization procedure introduced in the diffusion maps algorithm on the pairwise kernel values, namely (Coifman & Lafon Reference Coifman and Lafon2006a ; Berry & Sauer Reference Berry and Sauer2016),

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D617}_{ij}=\frac{\unicode[STIX]{x1D612}_{ij}}{\left(\displaystyle \mathop{\sum }_{n=0}^{N-1}\unicode[STIX]{x1D612}_{in}q_{n}^{-1/2}\right)q_{j}^{1/2}},\quad \text{with }q_{i}=\mathop{\sum }_{n=0}^{N-1}\unicode[STIX]{x1D612}_{in}.\end{eqnarray}$$

Let now $f:X\mapsto \mathbb{C}$ be a continuous function. Given the finite dataset $\{x_{n}\}_{n=0}^{N-1}$ , we represent $f$ by the column vector $\boldsymbol{f}=(f(x_{0}),\ldots ,f(x_{N-1}))^{\top }\in \mathbb{C}^{N}$ , whose components are equal to the values of $\boldsymbol{f}$ at the sampled states. It can be shown (e.g. Das & Giannakis Reference Das and Giannakis2017) that for every such function $f$ , the matrix–vector product $\boldsymbol{g}=\unicode[STIX]{x1D64B}\boldsymbol{f}$ is an ergodic average that approximates the action $Pf$ . This, in conjunction with compactness of $P$ , implies that the non-zero eigenvalues and corresponding eigenvectors, $(\unicode[STIX]{x1D6EC}_{k},\unicode[STIX]{x1D753}_{k})$ , of $\unicode[STIX]{x1D64B}$ approximate those of $P$ (von Luxburg, Belkin & Bousquet Reference von Luxburg, Belkin and Bousquet2008). Moreover, as $N\rightarrow \infty$ , the set $\{\unicode[STIX]{x1D753}_{k}\}$ converges in an appropriate sense to a basis for a closed subspace $\tilde{{\mathcal{D}}}$ of $L^{2}(A,\unicode[STIX]{x1D707})$ . One further verifies that the eigenvalues of $\unicode[STIX]{x1D64B}$ are real and admit the ordering $1=\unicode[STIX]{x1D6EC}_{0}>\unicode[STIX]{x1D6EC}_{1}\geqslant \unicode[STIX]{x1D6EC}_{2}\geqslant \cdots \,.$ We will use this ordering throughout the paper. In addition, the eigenvectors $\unicode[STIX]{x1D753}_{k}$ are orthogonal with respect to the weighted inner product

(3.7) $$\begin{eqnarray}\langle \,\boldsymbol{f},\boldsymbol{g}\rangle _{\unicode[STIX]{x1D64B}}=\mathop{\sum }_{i=0}^{N-1}f_{i}^{\ast }g_{i}w_{i},\quad \boldsymbol{f}=(f_{0},\ldots ,f_{N-1})^{\top },\boldsymbol{g}=(g_{0},\ldots ,g_{N-1})^{\top },\end{eqnarray}$$

where the weights $w_{i}$ are equal to the components of the stationary distribution $\boldsymbol{w}=(w_{0},\ldots ,w_{N-1})$ of $\unicode[STIX]{x1D64B}$ , satisfying $\boldsymbol{w}\unicode[STIX]{x1D64B}=\boldsymbol{w}$ and $\sum _{i=0}^{N-1}w_{i}=1$ .

A key property of $\unicode[STIX]{x1D64B}$ constructed via the NLSA kernel in (3.4) is that as the number of delays $Q$ grows, $\tilde{{\mathcal{D}}}$ becomes invariant under the action of the Koopman operator $U_{t}$ , and is also contained in the so-called discrete subspace ${\mathcal{D}}$ of $U_{t}$ spanned by all Koopman eigenfunctions (Das & Giannakis Reference Das and Giannakis2017). We therefore have the invariant inclusions $\tilde{{\mathcal{D}}}\subseteq {\mathcal{D}}\subseteq L^{2}(A,\unicode[STIX]{x1D707})$ . In particular, every Koopman eigenfunction $\unicode[STIX]{x1D713}\in \tilde{{\mathcal{D}}}$ can be approximated by a linear combination $\unicode[STIX]{x1D74D}=\sum _{k=0}^{\ell -1}c_{k}\unicode[STIX]{x1D753}_{k}$ of eigenvectors of $\unicode[STIX]{x1D64B}$ for some spectral truncation parameter $\ell \leqslant N-1$ , and in a suitable limit $\ell ,Q,N\rightarrow \infty$ , $\unicode[STIX]{x1D74D}$ converges to $\unicode[STIX]{x1D713}$ .

Another important property of $\unicode[STIX]{x1D64B}$ is that it induces a measure of roughness of observables through its eigenvalues. Specifically, to each $\unicode[STIX]{x1D753}_{k}$ we assign the quantity $\unicode[STIX]{x1D702}_{k}$ , with

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{k}=\frac{1}{\unicode[STIX]{x1D716}}(1-\unicode[STIX]{x1D6EC}_{k}^{-1}),\end{eqnarray}$$

and $\unicode[STIX]{x1D702}_{0},\unicode[STIX]{x1D702}_{1},\ldots$ forms an increasing sequence with $\unicode[STIX]{x1D702}_{0}=0$ . Here, the parameter $\unicode[STIX]{x1D716}$ stands again for the kernel bandwidth (see (3.4)). If the invariant set $A$ has a Riemannian manifold structure, the $\unicode[STIX]{x1D702}_{k}$ converge in a suitable asymptotic limit $\unicode[STIX]{x1D716}\rightarrow 0$ and $N\rightarrow \infty$ , and, up to a proportionality constant, to the values $E_{k}={\mathcal{E}}(\unicode[STIX]{x1D719}_{k})$ of the Dirichlet energy functional ${\mathcal{E}}(\unicode[STIX]{x1D719}_{k})=\int _{A}\Vert \text{grad}\,\unicode[STIX]{x1D719}_{k}\Vert ^{2}\,\text{d}\unicode[STIX]{x1D707}$ , where $\unicode[STIX]{x1D719}_{k}$ is again the eigenfunction of $P$ to which $\unicode[STIX]{x1D753}_{k}$ converges (Coifman & Lafon Reference Coifman and Lafon2006a ). In this limit, the $E_{k}$ form together with the $\unicode[STIX]{x1D719}_{k}$ eigenvalue–eigenfunction pairs of the Laplacian $\unicode[STIX]{x1D6E5}=-\text{div}\,\text{grad}$ on the manifold. Given our convention of ordering the $\unicode[STIX]{x1D702}_{k}$ in increasing order, for every non-negative integer $\ell$ , $\{\unicode[STIX]{x1D719}_{0},\ldots ,\unicode[STIX]{x1D719}_{\ell -1}\}$ corresponds to an $\ell$ -element orthonormal set of functions in $L^{2}(A,\unicode[STIX]{x1D707})$ with the least Dirichlet energy.

Intuitively, the Dirichlet energy can be thought of as a measure of ‘roughness’, or spatial variability, of functions on a manifold $A$ . In particular, a highly oscillatory function $f:A\mapsto \mathbb{R}$ has, on average, large values of its gradient norm $\Vert \text{grad}\,f\Vert$ , and thus large Dirichlet energy ${\mathcal{E}}(f)$ , whereas ${\mathcal{E}}(f)$ will be small if $f$ is spatially smooth. It should be noted that the Dirichlet energy should not be understood as a kinetic or thermal energy. Rather, if one were to think of the function $f$ as describing the displacement of a linearly elastic membrane on $A$ , ${\mathcal{E}}(f)$ would measure the elastic potential energy stored in the membrane.

While the invariant set $A$ of our turbulent convection flow is not a smooth manifold, and we only have access to finitely many snapshots, one can still think of the $\unicode[STIX]{x1D702}_{k}$ from (3.8) as Dirichlet energies, measuring the roughness of the corresponding eigenvectors $\unicode[STIX]{x1D753}_{k}$ . In particular, we observe that because $\unicode[STIX]{x1D64B}$ is Markov, $\unicode[STIX]{x1D753}_{0}=(1,1,\ldots ,1)^{\top }\in \mathbb{C}^{N}$ is the constant eigenvector that gets assigned zero Dirichlet energy, and $\unicode[STIX]{x1D753}_{1},\unicode[STIX]{x1D753}_{2},\ldots$ are mutually orthogonal eigenvectors with respect to the inner product in (3.7) of increasingly higher energy. Moreover, despite $A$ lacking a manifold structure, it is still possible to construct a self-adjoint positive-semidefinite unbounded operator $\unicode[STIX]{x1D6E5}:D(\unicode[STIX]{x1D6E5})\mapsto L^{2}(A,\unicode[STIX]{x1D707})$ , analogous to the Laplacian on Riemannian manifolds. This operator has a purely discrete spectrum and acts on functions $f=\sum _{k=0}^{\infty }c_{k}\unicode[STIX]{x1D719}_{k}$ through the formula

(3.9) $$\begin{eqnarray}\unicode[STIX]{x0394}f=\mathop{\sum }_{k=0}^{\infty }E_{k}c_{k}\unicode[STIX]{x1D719}_{k},\end{eqnarray}$$

where the $E_{k}$ are quantities to which the $\unicode[STIX]{x1D702}_{k}$ converge as $N\rightarrow \infty$ (Das & Giannakis Reference Das and Giannakis2017). The data-driven analogue $\unicode[STIX]{x1D63F}$ of this operator is then defined through

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D63F}\boldsymbol{f}=\mathop{\sum }_{k=0}^{N-1}\unicode[STIX]{x1D702}_{k}c_{k}\unicode[STIX]{x1D753}_{k},\end{eqnarray}$$

where $\boldsymbol{f}=\sum _{k=0}^{N-1}c_{k}\unicode[STIX]{x1D753}_{k}$ .

From a numerical standpoint, it should be noted that the $(\unicode[STIX]{x1D6EC}_{k},\unicode[STIX]{x1D753}_{k})$ can be computed stably and efficiently by solving the eigenvalue problem for the $N\times N$ Markov matrix $\unicode[STIX]{x1D64B}$ . In particular, at large sample numbers $N$ , that matrix can be made sparse by zeroing all but the largest $k_{nn}\ll N$ elements in each row of the unnormalized kernel matrix $\unicode[STIX]{x1D646}=[\unicode[STIX]{x1D612}_{ij}]$ , and symmetrizing the resulting sparse matrix. Moreover, the computation to form $\unicode[STIX]{x1D64B}$ scales linearly with the data space dimension $d$ (in this case, the number of velocity-field gridpoint values) and is trivially parallelizable. Indeed, once $\unicode[STIX]{x1D646}$ has been computed, the cost of forming $\unicode[STIX]{x1D64B}$ and computing $(\unicode[STIX]{x1D6EC}_{k},\unicode[STIX]{x1D753}_{k})$ , as well as the cost of the Koopman eigenvalue problem described in § 3.3 ahead, is independent of the ambient space dimension. In the problem studied here, $N=10\,000$ is small enough so as not to require sparsification of $\unicode[STIX]{x1D64B}$ , but this would become important at larger sample numbers (e.g. for the analysis of convection at higher Rayleigh numbers). However, the number of degrees of freedom of the velocity field in the DNS model has the value $d=3\times 128^{3}\approx 6.3\times 10^{6}$ , which is more than 600 times larger than $N$ , and the favourable computational cost of our approach with respect to $d$ becomes crucial.

3.3 Galerkin method for the Koopman eigenvalue problem

We compute approximate Koopman eigenfunctions through a Petrov–Galerkin method for the eigenvalue problem of the regularized generator $L$ in (3.3), defined using the diffusion operator $\unicode[STIX]{x1D6E5}$ in (3.9). This scheme will be formulated in a data-driven basis consisting of eigenfunctions of $\unicode[STIX]{x1D64B}$ . First, for any $p\in \mathbb{N}$ , we introduce the Sobolev spaces

(3.11) $$\begin{eqnarray}\displaystyle H^{p}=\left\{f=\mathop{\sum }_{k=1}^{\infty }c_{k}\unicode[STIX]{x1D719}_{k}\in L^{2}(A,\unicode[STIX]{x1D707}):\mathop{\sum }_{k=1}^{\infty }E_{k}^{p}|c_{k}^{2}|<\infty \right\} & & \displaystyle\end{eqnarray}$$

and equip these spaces with the inner products

(3.12) $$\begin{eqnarray}\displaystyle \langle \,f,g\rangle _{H^{p}}=\mathop{\sum }_{k=1}^{\infty }(1+E_{k}+\cdots +E_{k}^{p})c_{k}^{\ast }d_{k} & & \displaystyle\end{eqnarray}$$

and thus the norms $\Vert f\Vert _{H^{p}}=\langle \,f,f\rangle _{H^{p}}^{1/2}$ , where

(3.13a,b ) $$\begin{eqnarray}\displaystyle f=\mathop{\sum }_{k=1}^{\infty }c_{k}\unicode[STIX]{x1D719}_{k}\quad \text{and}\quad g=\mathop{\sum }_{k=1}^{\infty }d_{k}\unicode[STIX]{x1D719}_{k}. & & \displaystyle\end{eqnarray}$$

It should be noted that the infinite-dimensional spaces $H^{p}$ contain functions of higher regularity than $L^{2}$ by imposing a condition on the rate of decay of the expansion coefficients $c_{k}$ , which becomes stronger as $p$ grows. Moreover, all $H^{p}$ spaces contain zero-mean functions since they are $L^{2}$ -orthogonal to the constant function $\unicode[STIX]{x1D719}_{0}$ . In what follows, we will only require the spaces with $p\leqslant 2$ .

With these definitions in place, the eigenvalue problem for $L$ is approximated by the following well-posed regularized Koopman eigenvalue problem in weak form:

(3.14) $$\begin{eqnarray}{\mathcal{L}}(f,\unicode[STIX]{x1D713}):=\langle \,f,V\unicode[STIX]{x1D713}\rangle _{H^{0}}-\unicode[STIX]{x1D701}\langle \,f,\unicode[STIX]{x0394}\unicode[STIX]{x1D713}\rangle _{H^{0}}=\unicode[STIX]{x1D706}\langle \,f,\unicode[STIX]{x1D713}\rangle _{H^{0}},\quad \forall f\in H^{0}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D706}\in \mathbb{C}$ and $\unicode[STIX]{x1D713}\in H^{2}$ are weak eigenvalues and eigenfunctions of $L$ . Furthermore, $f$ is a test function and ${\mathcal{L}}:H^{0}\times H^{2}\mapsto \mathbb{C}$ is bounded.

Two steps are involved in passing to a data-driven approximation of the eigenvalue problem, namely (i) approximation of the trial and test spaces by appropriate ‘data-driven’ Sobolev spaces and (ii) approximation of the action of the generator $V$ on functions by finite differences or matrix logarithms.

Our data-driven Sobolev spaces are defined in direct analogy to their infinite-dimensional counterparts as $H_{N}^{p}:=\text{span}\{\unicode[STIX]{x1D753}_{1},\ldots ,\unicode[STIX]{x1D753}_{N-1}\}$ , and are equipped with the inner products and norms

(3.15a,b ) $$\begin{eqnarray}\displaystyle \langle \,\boldsymbol{f},\boldsymbol{g}\rangle _{H_{N}^{p}}=\mathop{\sum }_{k=1}^{N-1}(1+\unicode[STIX]{x1D702}_{k}+\cdots +\unicode[STIX]{x1D702}_{k}^{p})c_{k}^{\ast }d_{k},\quad \Vert \boldsymbol{f}\Vert _{H_{N}^{p}}=\langle \,\boldsymbol{f},\boldsymbol{f}\rangle _{H_{N}^{p}}^{1/2} & & \displaystyle\end{eqnarray}$$

respectively, where

(3.16a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{f}=\mathop{\sum }_{k=1}^{N-1}c_{k}\unicode[STIX]{x1D753}_{k}\quad \text{and}\quad \boldsymbol{g}=\mathop{\sum }_{k=1}^{N-1}d_{k}\unicode[STIX]{x1D753}_{k}. & & \displaystyle\end{eqnarray}$$

It should be noted that $\dim (H_{N}^{p})=N-1$ and that $\unicode[STIX]{x1D753}_{0}$ (which is not included in our definition of $H_{N}^{p}$ ) has constant elements. Moreover, the inner product of $H_{N}^{0}$ is equivalent to the weighted inner product in (3.7). We equip the data-driven Sobolev spaces $H_{N}^{p}$ with the normalized basis functions $\unicode[STIX]{x1D753}_{k}^{(p)}=\unicode[STIX]{x1D753}_{k}/\unicode[STIX]{x1D702}_{k}^{p/2}$ for $k\in \{1,2,\ldots ,N-1\}$ . The normalization factors $\unicode[STIX]{x1D702}_{k}$ are given by (3.8). These basis functions have the important property that for any $\boldsymbol{f}=\sum _{k=1}^{N-1}c_{k}\unicode[STIX]{x1D753}_{k}^{(p)}$ , converging as $N\rightarrow \infty$ to $f\in H^{p}$ , the $H_{N}^{p}$ norm can be bounded by $C\Vert \boldsymbol{c}\Vert$ , where $C$ is a constant independent of $N$ and $\Vert \boldsymbol{c}\Vert$ is the 2-norm of the vector of expansion coefficients $\boldsymbol{c}=(c_{1},\ldots ,c_{N-1})\in \mathbb{C}^{N-1}$ . This means that as $N\rightarrow \infty$ , $\{\unicode[STIX]{x1D753}_{k}^{(p)}\}_{k=1}^{N-1}$ converges to a basis of $H^{p}$ , whereas the unnormalized basis functions $\unicode[STIX]{x1D753}_{k}$ do not.

To approximate the action of the generator on functions, we recall from § 3.1 that $Vf$ corresponds to a ‘time derivative’ of observable $f$ along the dynamical flow and that $V$ is related to the Koopman operator $U_{t}$ via exponentiation, $U_{t}=\exp (tV)$ . For a sufficiently smooth function $f:X\mapsto \mathbb{C}$ , these facts respectively suggest that we can approximate $Vf$ through finite differences in time or via matrix logarithms. As a concrete example of the former approach, the quantity

(3.17) $$\begin{eqnarray}\displaystyle V_{\unicode[STIX]{x1D70F}}\,f=\frac{U_{\unicode[STIX]{x1D70F}}f-U_{-\unicode[STIX]{x1D70F}}f}{2\unicode[STIX]{x1D70F}} & & \displaystyle\end{eqnarray}$$

is a second-order central finite-difference approximation of $Vf$ for the sampling interval $\unicode[STIX]{x1D70F}$ . On evaluating this expression at a state $x_{n}$ of the training data, we have

(3.18) $$\begin{eqnarray}\displaystyle V_{\unicode[STIX]{x1D70F}}\,f(x_{n})=\frac{f(x_{n+1})-f(x_{n-1})}{2\unicode[STIX]{x1D70F}}. & & \displaystyle\end{eqnarray}$$

Moreover, since the components $\unicode[STIX]{x1D719}_{nk}$ of the basis vectors $\unicode[STIX]{x1D753}_{k}=(\unicode[STIX]{x1D719}_{0k},\ldots ,\unicode[STIX]{x1D719}_{N-1,k})$ correspond to the values of a function at $x_{n}\in X$ (see § 3.2), we can define an approximate generator $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}:H_{N}^{0}\mapsto H_{N}^{0}$ based on the action of $V_{\unicode[STIX]{x1D70F}}$ on the $\unicode[STIX]{x1D753}_{k}$ basis elements. Specifically, by introducing the vectors $\unicode[STIX]{x1D753}_{k}^{\prime }=(\unicode[STIX]{x1D719}_{0k}^{\prime },\ldots ,\unicode[STIX]{x1D719}_{N-1,k}^{\prime })$ with elements

(3.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{nk}^{\prime }=\frac{\unicode[STIX]{x1D719}_{n+1,k}-\unicode[STIX]{x1D719}_{n-1,k}}{2\unicode[STIX]{x1D70F}},\quad 1\leqslant n\leqslant N-2, & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D719}_{0k}^{\prime }=\unicode[STIX]{x1D719}_{N-1,k}^{\prime }=0$ , the action $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}\boldsymbol{f}=\boldsymbol{g}$ on an observable $\boldsymbol{f}=\sum _{k=1}^{N-1}c_{k}\unicode[STIX]{x1D753}_{k}\in H_{N}^{0}$ is given by $\boldsymbol{g}=\sum _{j=1}^{N-1}d_{j}\unicode[STIX]{x1D753}_{j}$ , where

(3.20) $$\begin{eqnarray}d_{j}=\mathop{\sum }_{k=1}^{N-1}\unicode[STIX]{x1D61D}_{\unicode[STIX]{x1D70F},jk}c_{k},\quad \unicode[STIX]{x1D61D}_{\unicode[STIX]{x1D70F},jk}=\langle \unicode[STIX]{x1D753}_{j},\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D753}_{k}\rangle _{H_{N}^{0}}=\langle \unicode[STIX]{x1D753}_{j},\unicode[STIX]{x1D753}_{k}^{\prime }\rangle _{H_{N}^{0}}.\end{eqnarray}$$

In the above, the quantities $\unicode[STIX]{x1D61D}_{\unicode[STIX]{x1D70F},jk}$ are the matrix elements of $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}$ in the $\{\unicode[STIX]{x1D753}_{k}\}_{k=1}^{N-1}$ basis of $H_{N}^{0}$ . Alternatively, to construct an approximation of $V$ based on matrix logarithms, we proceed as above, replacing $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}$ in (3.20) with

(3.21a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D70F}^{-1}\log \unicode[STIX]{x1D650}_{\unicode[STIX]{x1D70F}},\quad \unicode[STIX]{x1D650}_{\unicode[STIX]{x1D70F}}=[U_{\unicode[STIX]{x1D70F},jk}],\quad U_{\unicode[STIX]{x1D70F},jk}=\langle \unicode[STIX]{x1D753}_{j},\unicode[STIX]{x1D753}_{k}^{+}\rangle _{H_{N}^{0}},\end{eqnarray}$$

where $\unicode[STIX]{x1D753}_{k}^{+}=(\unicode[STIX]{x1D719}_{0k}^{+},\ldots ,\unicode[STIX]{x1D719}_{N-1,k}^{+})$ is the vector with elements

(3.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{nk}^{+}=\unicode[STIX]{x1D719}_{n+1,k},\quad 1\leqslant n\leqslant N-1, & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D719}_{0,k}^{+}=0$ .

With these definitions, we pose the following variational eigenvalue problem, which is a data-driven analogue of the regularized Koopman eigenvalue problem in (3.14). Fixing a spectral truncation parameter $\ell \leqslant N-1$ , the problem consists of finding eigenvalues $\unicode[STIX]{x1D706}\in \mathbb{C}$ and eigenfunctions $\unicode[STIX]{x1D74D}\in \text{span}\{\unicode[STIX]{x1D753}_{1}^{(2)},\ldots ,\unicode[STIX]{x1D753}_{\ell }^{(2)}\}$ such that for all $\boldsymbol{f}\in \text{span}\{\unicode[STIX]{x1D753}_{1},\ldots ,\unicode[STIX]{x1D753}_{\ell }\}$ ,

(3.23) $$\begin{eqnarray}\hat{{\mathcal{L}}}(\boldsymbol{f},\unicode[STIX]{x1D74D})=\unicode[STIX]{x1D706}\langle \,\boldsymbol{f},\unicode[STIX]{x1D74D}\rangle _{H_{N}^{0}}\end{eqnarray}$$

holds, where $\hat{{\mathcal{L}}}:H_{N}^{0}\times H_{N}^{2}\mapsto \mathbb{C}$ is the sesquilinear form defined as

(3.24) $$\begin{eqnarray}\displaystyle \hat{{\mathcal{L}}}(\boldsymbol{f},\unicode[STIX]{x1D74D})=\langle \,\boldsymbol{f},\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D74D}\rangle _{H_{N}^{0}}-\unicode[STIX]{x1D701}\langle \,\boldsymbol{f},\unicode[STIX]{x1D63F}\unicode[STIX]{x1D74D}\rangle _{H_{N}^{0}}, & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D70F}}$ is an approximation of the generator, obtained via either the finite-difference or the matrix logarithm approach in (3.20) and (3.21) respectively. Numerically, this variational problem is equivalent to solving a matrix generalized eigenvalue problem

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D647}\boldsymbol{c}_{k}=\unicode[STIX]{x1D706}_{k}\unicode[STIX]{x1D63D}\boldsymbol{c}_{k},\quad \unicode[STIX]{x1D706}_{k}\in \mathbb{C},\end{eqnarray}$$

where $\unicode[STIX]{x1D647}$ and $\unicode[STIX]{x1D63D}$ are $\ell \times \ell$ matrices with elements

(3.26a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D613}_{ij}=\hat{{\mathcal{L}}}(\unicode[STIX]{x1D753}_{i},\unicode[STIX]{x1D753}_{j}^{(2)})=\unicode[STIX]{x1D61D}_{\unicode[STIX]{x1D70F},ij}\unicode[STIX]{x1D702}_{j}^{-1}-\unicode[STIX]{x1D701}\unicode[STIX]{x1D6FF}_{ij},\quad \unicode[STIX]{x1D609}_{ij}=\langle \unicode[STIX]{x1D753}_{i},\unicode[STIX]{x1D753}_{j}^{(2)}\rangle _{H_{N}^{0}}=\unicode[STIX]{x1D702}_{i}^{-1}\unicode[STIX]{x1D6FF}_{ij}\end{eqnarray}$$

respectively and $\boldsymbol{c}_{k}=(c_{1k},\ldots ,c_{\ell k})^{\top }$ is a column vector in $\mathbb{C}^{\ell }$ containing the expansion coefficients of the solution $\unicode[STIX]{x1D74D}_{k}$ in the $\{\unicode[STIX]{x1D753}_{j}^{(2)}\}$ basis of $H_{N}^{2}$ , namely $\unicode[STIX]{x1D74D}_{k}=\sum _{j=1}^{\ell }c_{jk}\unicode[STIX]{x1D753}_{j}^{(2)}$ . It should be noted that matrix eigenvalue problems of this class are encountered frequently in finite-element methods for elliptic partial differential equations (Babuška & Osborn Reference Babuška, Osborn, Ciarlet and Lions1991), where the matrices $\unicode[STIX]{x1D647}$ and $\unicode[STIX]{x1D63D}$ are referred to as stiffness and mass matrices respectively.

It can be shown that under mild assumptions on the convergence of the $\unicode[STIX]{x1D702}_{k}$ to the Dirichlet energies $E_{k}$ , for any $\ell$ , the solutions of the eigenvalue problem in (3.23) converge as $N\rightarrow \infty$ to the solutions of the problem in (3.14), restricted to the $\ell$ -dimensional subspace of $H^{2}$ spanned by $\{\unicode[STIX]{x1D719}_{1},\ldots ,\unicode[STIX]{x1D719}_{\ell }\}$ . It should be noted again that it is important to work with the normalized basis $\{\unicode[STIX]{x1D753}_{k}^{(2)}\}$ of $H_{N}^{2}$ as opposed to the unnormalized basis $\{\unicode[STIX]{x1D753}_{k}\}$ in order to ensure good conditioning of our scheme at large $\ell$ . This is reflected from the fact that $\unicode[STIX]{x1D63F}$ is represented in this basis by the identity matrix, whose condition number is 1 for all $\ell$ . In practice, one typically works with $\ell \ll N-1$ , so that (3.25) can be feasibly solved at large sample numbers. As with the NLSA basis functions from § 3.2, we order our solutions $(\unicode[STIX]{x1D706}_{k},\unicode[STIX]{x1D74D}_{k})$ , $k\in \{1,\ldots ,\ell \}$ , in order of increasing Dirichlet energy,

(3.27) $$\begin{eqnarray}\hat{{\mathcal{E}}}(\unicode[STIX]{x1D74D}_{k}):=\frac{\langle \unicode[STIX]{x1D74D}_{k},\unicode[STIX]{x1D63F}\unicode[STIX]{x1D74D}_{k}\rangle _{H_{N}^{0}}}{\Vert \unicode[STIX]{x1D74D}_{k}\Vert _{H_{N}^{0}}^{2}}=\frac{\displaystyle \mathop{\sum }_{i=1}^{\ell }|c_{ik}|^{2}/\unicode[STIX]{x1D702}_{i}}{\displaystyle \mathop{\sum }_{j=1}^{\ell }|c_{jk}|^{2}/\unicode[STIX]{x1D702}_{j}^{2}}.\end{eqnarray}$$

As stated in § 3.1, the Koopman eigenfunctions corresponding to non-zero eigenvalues form complex-conjugate pairs (with equal Dirichlet energies), and this is usually the case for the data-driven $(\unicode[STIX]{x1D706}_{k},\unicode[STIX]{x1D74D}_{k})$ . However, due to the presence of diffusion, it is possible that $\unicode[STIX]{x1D647}$ has purely real negative eigenvalues with corresponding real eigenfunctions. We will, in fact, encounter such a solution in § 4 ahead. It should be noted that the Koopman eigenvalue problem also has the trivial (constant) solution $\unicode[STIX]{x1D706}_{0}=0$ and $\unicode[STIX]{x1D74D}_{0}=\unicode[STIX]{x1D753}_{0}$ , with vanishing Dirichlet energy, but our numerical eigenvalue problem does not yield this solution since our trial space is orthogonal to $\unicode[STIX]{x1D74D}_{0}$ .

Table 1.

Summary of the several levels of description towards the data-driven and Galerkin-approximated eigenvalue problem for the regularized Koopman generator in weak form. The step from ${\mathcal{L}}$ to $\hat{{\mathcal{L}}}$ requires the generation of a data-driven basis that consists of the eigenfunctions of the Markov operator $\unicode[STIX]{x1D64B}$ obtained by NLSA.

Table 1 summarizes the different levels of description, starting from the original Koopman operator eigenvalue problem to the data-driven and Galerkin-approximated eigenvalue problem for the regularized Koopman generator.

3.4 Spatiotemporal reconstruction

Solution of the eigenvalue problem in (3.25) results in the data-driven complex Koopman eigenvalues $\unicode[STIX]{x1D706}_{k}\in \mathbb{C}$ and the corresponding eigenfunctions $\unicode[STIX]{x1D74D}_{k}\in \mathbb{C}^{N}$ , with $1\leqslant k\leqslant \ell \ll N-1$ , and $N$ being the number of turbulent convection snapshots. The spatiotemporal patterns corresponding to the eigenfunctions are computed through the standard approach employed in singular spectrum analysis algorithms (the analogue of POD in delay-coordinate space) by Ghil et al. (Reference Ghil, Allen, Dettinger, Ide, Kondrashov, Mann, Robertson, Saunders, Tian, Varadi and Yiou2002), which is also employed in NLSA. Here, we describe the reconstruction procedure for the velocity field itself, but note that the same approach can be used to reconstruct other observables such as the temperature field or combined quantities of both fields.

First, given a lead/lag $t^{\prime }\in \mathbb{R}$ , we compute the projections $\boldsymbol{A}_{k}(t^{\prime })$ of the observable $F$ onto the Koopman eigenfunctions $\unicode[STIX]{x1D713}_{k}$ ,

(3.28) $$\begin{eqnarray}\boldsymbol{A}_{k}(t^{\prime })=\langle \unicode[STIX]{x1D713}_{k},F\circ \unicode[STIX]{x1D6F7}_{t^{\prime }}\rangle _{H^{0}}=\int _{A}\unicode[STIX]{x1D713}_{k}^{\ast }(x)F(\unicode[STIX]{x1D6F7}_{t^{\prime }}(x))\,\text{d}\unicode[STIX]{x1D707}(x).\end{eqnarray}$$

It should be noted that for each time $t^{\prime }$ , $\boldsymbol{A}_{k}(t^{\prime })$ is a velocity-field snapshot in $L^{2}(\unicode[STIX]{x1D6FA};\mathbb{R}^{3})$ . In applications, we evaluate this quantity at the discrete times $t_{q}^{\prime }=q\unicode[STIX]{x1D70F}$ , with $q\in \mathbb{Z}$ , and approximate the integral with respect to the invariant measure $\unicode[STIX]{x1D707}$ in (3.28) by time averages. This leads to a complex-valued vector field given by

(3.29) $$\begin{eqnarray}\hat{\boldsymbol{A}}_{k}(q\unicode[STIX]{x1D70F})=\frac{1}{N^{\prime }}\mathop{\sum }_{n=0}^{N^{\prime }-1}\unicode[STIX]{x1D713}_{nk}^{\ast }\boldsymbol{u}_{n+q},\quad \text{with }N^{\prime }=\min \{N,N-q\}.\end{eqnarray}$$

In (3.29), we made use of the fact that the observable $F$ evaluated on a state $(\boldsymbol{u},T)$ of the Boussinesq system is simply equal to the velocity field $\boldsymbol{u}$ ; see § 2. In particular, $\boldsymbol{u}_{n+q}$ is the velocity field of the convection flow sampled at time $t=(n+q)\unicode[STIX]{x1D70F}$ , and, moreover, $\unicode[STIX]{x1D713}_{nk}$ is the $n$ th component of the Koopman eigenvector $\unicode[STIX]{x1D74D}_{k}=(\unicode[STIX]{x1D713}_{0k},\ldots ,\unicode[STIX]{x1D713}_{N-1,k})^{\top }$ . In the Koopman operator literature, $\hat{\boldsymbol{A}}_{k}(0)$ is referred to as the Koopman mode associated with the observable $F$ (Mezić Reference Mezić2005; Budisić, Mohr & Mezić Reference Budisić, Mohr and Mezić2012).

The corresponding data-derived sequence of (complex-valued) spatiotemporal velocity-field patterns associated with the pair $(\unicode[STIX]{x1D74D}_{k},\hat{\boldsymbol{A}}_{k})$ is given by

(3.30) $$\begin{eqnarray}\hat{\boldsymbol{u}}_{n}^{(k)}=\frac{1}{Q^{\prime }}\mathop{\sum }_{q=0}^{Q^{\prime }-1}\hat{\boldsymbol{A}}_{k}(-q\unicode[STIX]{x1D70F})\unicode[STIX]{x1D713}_{n+q,k},\quad \text{with }Q^{\prime }=\min \{Q,N-1-n\},0\leqslant n\leqslant N-1.\end{eqnarray}$$

Since the data-driven Koopman eigenfunctions usually form complex-conjugate pairs, $(\unicode[STIX]{x1D74D}_{k},\unicode[STIX]{x1D74D}_{k^{\prime }})$ , we compute the sum $\hat{\boldsymbol{u}}_{k}+\hat{\boldsymbol{u}}_{k^{\prime }}$ to obtain a real-valued physical velocity field. As is evident from (3.30), a reconstruction at time $n\unicode[STIX]{x1D70F}$ is a superposition of $Q$ snapshots. This is a result of the time-lagged embedding and different from a standard POD analysis.

4.1 Time-averaging analysis at probe points

Before coming to the Koopman analysis, we show in figure 3 time series of the vertical velocity, measured at eight different simulation gridpoints in the midplane in their statistically steady state. This analysis is similar to the one conducted by Foroozani et al. (Reference Foroozani, Niemela, Armenio and Sreenivasan2017). Four of these points are in the four corner regions, slightly away from the faces, and the other four points are located halfway between the corner points. Specifically, point #1 in the $\unicode[STIX]{x1D703}=5\unicode[STIX]{x03C0}/4$ corner is found at $(x_{1},y_{1})=(0.078,-0.422)$ , point #2 at $(0.5,-0.422)$ and point #3 in the $\unicode[STIX]{x1D703}=7\unicode[STIX]{x03C0}/4$ corner at $(0.922,-0.422)$ . The other five points $(x_{4},y_{4}),\ldots ,(x_{8},y_{8})$ follow correspondingly by symmetry in counterclockwise orientation. In addition, we apply a temporal and spatial smoothing of the strongly fluctuating signals to highlight the large-scale motion and to reduce the amount of available information. In detail, we substitute the vertical velocity $u_{z}(x_{i},y_{i},0,t)=w_{x,y}(t)$ at the monitoring points $(x_{1},y_{1}),\ldots (x_{8},y_{8})$ by

(4.1) $$\begin{eqnarray}\tilde{w}_{x,y}(t)=\frac{1}{4}\left[\frac{w_{x+h,y}(t)+w_{x-h,y}(t)+w_{x,y+h}(t)+w_{x,y-h}(t)+4w_{x,y}(t)}{2}\right],\end{eqnarray}$$

which corresponds to a local nearest-neighbour averaging on the simulation mesh with mesh width $h$ . In addition, functions are smoothed with Bézier splines using a $10T_{f}$ temporal window. The results confirm that all four diagonal states are visited over our analysis time series. In what follows, we will describe how Koopman eigenfunction analysis yields a representation of these states without requiring spatial or temporal smoothing, while also providing access to LSC patterns on shorter time scales (the secondary eigenfunctions), which would be difficult to extract using averaging approaches.

4.2 Primary eigenfunctions

The primary eigenfunction family consists of the leading three data-driven Koopman eigenfunctions, $\unicode[STIX]{x1D74D}_{1},\unicode[STIX]{x1D74D}_{2},\unicode[STIX]{x1D74D}_{3}\in \mathbb{C}^{N}$ , ordered with respect to increasing Dirichlet energy, $\hat{{\mathcal{E}}}(\unicode[STIX]{x1D74D}_{1})=\hat{{\mathcal{E}}}(\unicode[STIX]{x1D74D}_{2})<\hat{{\mathcal{E}}}(\unicode[STIX]{x1D74D}_{3})$ . This family includes a complex-conjugate pair, $(\unicode[STIX]{x1D74D}_{1},\unicode[STIX]{x1D74D}_{2})=\text{Re}(\unicode[STIX]{x1D74D}_{1})\pm \text{i}\text{Im}(\unicode[STIX]{x1D74D}_{1})$ , and a purely real mode, $\unicode[STIX]{x1D74D}_{3}=\text{Re}(\unicode[STIX]{x1D74D}_{3})$ . Thus, by defining $y_{1}(n)=\text{Re}(\unicode[STIX]{x1D713}_{n1})$ , $y_{2}(n)=\text{Im}(\unicode[STIX]{x1D713}_{n1})$ and $y_{3}(n)=\unicode[STIX]{x1D713}_{n3}$ , with $n\in \{0,\ldots ,N-1\}$ , the evolution of the LSC configuration reconstructed from the primary eigenfunctions can be represented via a trajectory in a three-dimensional phase space with coordinates $(y_{1},y_{2},y_{3})$ , as shown in figure 4. As stated in § 3.1, the trajectory is sampled at the $N=10\,000$ temporally consecutive states $x_{n}$ at times $t_{n}=n\unicode[STIX]{x1D70F}$ , which are separated by one free-fall time unit, $\unicode[STIX]{x1D70F}=1$ , from one another. Visually, it can readily be seen that the trajectory yields four clustered regions. These regions can also be identified quantitatively via $K$ -means clustering performed on the $(y_{1},y_{2},y_{3})$ coordinates. Residence times in these clusters can be well in excess of 1000 free-fall time units.

Figure 4.

The dynamical trajectory in the three-dimensional phase space with coordinates $(y_{1}(n),y_{2}(n),y_{3}(n))$ , constructed from the primary Koopman eigenfunctions. The trajectory is sampled at the $N=10\,000$ analysis states at times $t_{n}=n\unicode[STIX]{x1D70F}$ , where the sampling interval $\unicode[STIX]{x1D70F}$ is equal to one free-fall time unit. Four clusters can be assigned from a hierarchical $K$ -means clustering performed on this point cloud. One can notice the larger number of transitions between clusters A and B, or C and D, compared with B and C, or D and A.

In order to gain insight into the properties of the LSC states captured by the primary eigenfunctions, it is useful to examine their corresponding projected velocity-field patterns $\hat{\boldsymbol{A}}_{k}(q\unicode[STIX]{x1D70F})$ from (3.29). Figure 5 displays these patterns for $q=0$ for horizontal cross-sections near the top plate ( $z=0.45$ ) and at the midplane ( $z=0$ ), as well as for a vertical cross-section at $y=0$ . We also show the time-averaged state, which corresponds to the $\hat{\boldsymbol{A}}_{0}(0)$ projection pattern associated with the trivial (constant) Koopman eigenfunction $\unicode[STIX]{x1D74D}_{0}=(1,\ldots ,1)^{\top }\in \mathbb{C}^{N}$ . It is evident from the results in figure 5 that the real and imaginary parts of $\hat{\boldsymbol{A}}_{1}(0)$ capture diagonal LSC patterns, with the circulation locked along one of the two diagonals of the square cross-section of the flow domain. On the other hand, at least when examined on the $z=0$ plane, the LSC patterns in $\hat{\boldsymbol{A}}_{3}(0)$ and the time-averaged state appear to be oriented along one of the two sidewall pairs. By examining $\hat{\boldsymbol{A}}_{k}(q\unicode[STIX]{x1D70F})$ at different values of $q$ (not shown here), it can be verified that the patterns in figure 5 remain largely constant over the $Q=30$ delay-embedding window. Due to this, in the ensuing discussion, we will drop the $q\unicode[STIX]{x1D70F}$ argument from $\hat{\boldsymbol{A}}_{k}(q\unicode[STIX]{x1D70F})$ in the interest of a simple notation.

Figure 5.

Projected velocity-field snapshots $\hat{\boldsymbol{A}}_{k}(q\unicode[STIX]{x1D70F})$ at lag $q=0$ for the primary Koopman eigenfunctions and the trivial (constant) Koopman eigenfunction, corresponding to a time average. The colour plots in (a–d) and (e–h) show the vertical velocity-field components at horizontal cross-sections near the top ( $z=0.45$ ) and in the middle ( $z=0$ ) of the domain. The vector plots in (i–l) show the $(\boldsymbol{e}_{x},\boldsymbol{e}_{z})$ velocity-field components at a vertical cross-section perpendicular to $\boldsymbol{e}_{y}$ through the middle ( $y=0$ ) of the domain. Panels (a,e,i), (b,f,j), (c,g,k) and (d,h,l) show $\text{Re}\,\hat{\boldsymbol{A}}_{1}$ , $\text{Im}\,\hat{\boldsymbol{A}}_{1}$ , $\hat{\boldsymbol{A}}_{3}$ and $\hat{\boldsymbol{A}}_{0}$ respectively.

Upon initial inspection, the vertical velocity-field configuration in $\text{Re}\,\hat{\boldsymbol{A}}_{1}$ and $\text{Im}\,\hat{\boldsymbol{A}}_{1}$ suggests that these patterns are related by $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ (anticlockwise) rotations. From a dynamical symmetry standpoint, this is consistent with the fact that the Boussinesq system in a cubical domain is equivariant under the $\mathbb{Z}_{4}$ cyclic group associated with rotations by integer multiples of $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ . That is, if $(\boldsymbol{u}(\boldsymbol{x},t),T(\boldsymbol{x},t))$ is a solution, then $(\unicode[STIX]{x1D64D}_{g}\boldsymbol{u}(\unicode[STIX]{x1D64D}_{g}^{-1}\boldsymbol{x},t),T(\unicode[STIX]{x1D64D}_{g}^{-1}\boldsymbol{x},t))$ is also a solution, where $\unicode[STIX]{x1D64D}_{g}$ is a $3\times 3$ rotation matrix representing group element $g\in \mathbb{Z}_{4}$ . One can verify that if the system samples such solutions with equal probability in the course of dynamical evolution (more specifically, if the symmetry group action preserves the ergodic invariant measure of the dynamics), then the real and imaginary parts of the Koopman projected patterns $\hat{\boldsymbol{A}}_{k}$ should span two-dimensional unitary representation spaces of the symmetry group, which implies in turn that $\text{Re}\,\hat{\boldsymbol{A}}_{k}$ and $\text{Im}\,\hat{\boldsymbol{A}}_{k}$ should be relatable by symmetry operations (for more details, see, e.g., Giannakis et al. Reference Giannakis, Ourmazd, Slawinska and Zhao2017). The fact that $\text{Re}\,\hat{\boldsymbol{A}}_{1}$ and $\text{Im}\,\hat{\boldsymbol{A}}_{1}$ are at least qualitatively relatable by a $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ rotation is consistent with this picture. What is unexpected, however, is a difference in the velocity-field magnitudes exhibited by $\text{Re}\,\hat{\boldsymbol{A}}_{1}$ and $\text{Im}\,\hat{\boldsymbol{A}}_{1}$ . In particular, as can be seen in figure 5, typical vertical velocity-field magnitudes are approximately a factor of two smaller in $\text{Im}\,\hat{\boldsymbol{A}}_{1}$ than $\text{Re}\,\hat{\boldsymbol{A}}_{1}$ . This difference can be attributed to the fact that $|\text{Im}\,\unicode[STIX]{x1D74D}_{1}|$ is typically approximately a factor of two smaller than $|\text{Re}\,\unicode[STIX]{x1D74D}_{1}|$ , but, aside from our analysis not having sampled all LSC states in equal proportion, there is no obvious mechanism underlying such an asymmetry. We will return to the repercussions of this asymmetry on the properties of the reconstructed spatiotemporal patterns associated with the primary Koopman eigenfunctions below.

Next, turning to the $\hat{\boldsymbol{A}}_{3}$ and $\hat{\boldsymbol{A}}_{0}$ patterns, we first note that despite the $\mathbb{Z}_{4}$ equivariance of the Boussinesq system, the time-averaged velocity field is not fully invariant under this symmetry group. This is manifested by the vertical velocity-field pattern in $\hat{\boldsymbol{A}}_{0}$ at the $z=0$ horizontal cross-section. This pattern exhibits a clear bipolar structure oriented in the $\unicode[STIX]{x1D703}=0$ direction, indicating that the system does not visit all $\mathbb{Z}_{4}$ symmetry states with equal probability in our analysis interval. Intriguingly, at this cross-section, the $\hat{\boldsymbol{A}}_{3}$ and $\hat{\boldsymbol{A}}_{0}$ patterns can be related via a $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ rotation. However, elsewhere in the domain, these patterns have pronounced differences, which cannot be accounted for by the $\mathbb{Z}_{4}$ symmetry group. For instance, in the $z=0.45$ horizontal cross-sections in figure 5, the vertical velocity field associated with $\hat{\boldsymbol{A}}_{3}$ exhibits a quadrupolar pattern, whereas that associated with $\hat{\boldsymbol{A}}_{0}$ is unipolar.

We now examine the spatiotemporal reconstructions associated with the primary Koopman eigenfunctions. Due to the non-trivial structure of the time-averaged state with respect to $\mathbb{Z}_{4}$ rotations, and its relationship with $\hat{\boldsymbol{A}}_{3}$ described above, in what follows, we will include the time-averaged pattern in our reconstructions. That is, we consider the reconstructed velocity field

(4.2) $$\begin{eqnarray}\hat{\boldsymbol{u}}_{n}^{primary}:=\hat{\boldsymbol{u}}_{n}^{(0)}+\hat{\boldsymbol{u}}_{n}^{(1)}+\hat{\boldsymbol{u}}_{n}^{(2)}+\hat{\boldsymbol{u}}_{n}^{(3)},\quad 0\leqslant n\leqslant N-1,\end{eqnarray}$$

with $\hat{\boldsymbol{u}}_{n}^{(1)}$ , $\hat{\boldsymbol{u}}_{n}^{(2)}$ and $\hat{\boldsymbol{u}}_{n}^{(3)}$ given by (3.30) and $\hat{\boldsymbol{u}}_{n}^{(0)}=\hat{\boldsymbol{A}}_{0}(0)$ . It should be noted that the resulting velocity-field projection is real-valued since the Koopman eigenfunctions $\unicode[STIX]{x1D74D}_{0}$ and $\unicode[STIX]{x1D74D}_{3}$ are real-valued, and $\unicode[STIX]{x1D74D}_{1}$ and $\unicode[STIX]{x1D74D}_{2}$ form a complex-conjugate pair. We found that inclusion of $\hat{\boldsymbol{u}}_{0}^{(n)}$ aids the physical interpretation of our results, as, without this pattern, our reconstructions would describe anomalies relative to the time-averaged state, which are somewhat cumbersome to visually interpret against a non-trivial time-averaged background. With this convention, due to the fact that the $\hat{\boldsymbol{A}}_{k}(q\unicode[STIX]{x1D70F})$ patterns are nearly constant over the delay-embedding window for the primary Koopman eigenfunctions, our reconstructed velocity-field patterns are, to a good approximation, given by the time-averaged state $\hat{\boldsymbol{A}}_{0}(0)$ plus linear combinations of the $\text{Re}\,\hat{\boldsymbol{A}}_{1}(0)$ , $\text{Im}\,\hat{\boldsymbol{A}}_{1}(0)$ and $\hat{\boldsymbol{A}}_{3}(0)$ patterns, with coefficients determined from the time averages of the $y_{1}$ , $y_{2}$ and $y_{3}$ coordinates respectively over the time interval $[n\unicode[STIX]{x1D70F},n\unicode[STIX]{x1D70F}+(Q-1)\unicode[STIX]{x1D70F}]$ . The nature of the resulting spatiotemporal patterns will therefore depend on the temporal relationships of the $y_{k}$ , and particularly their clustering behaviour, as we now discuss.

Figure 6 shows representative velocity-field snapshots $\hat{\boldsymbol{u}}_{n}^{primary}$ , visualized at the same horizontal ( $z=0.45$ and $z=0$ ) and vertical ( $y=0$ ) planes as in figure 5, for each of the four clusters in the $(y_{1},y_{2},y_{3})$ phase space. As stated above, the system can spend periods exceeding 1000 free-fall time units in each of these clusters, and displays infrequent rapid transitions between them. Based on the velocity-field patterns at $z=0.45$ , it is evident that clusters A, B, C and D are associated with diagonal LSC states with impingement points at the $\unicode[STIX]{x1D703}=5\unicode[STIX]{x03C0}/4$ , $7\unicode[STIX]{x03C0}/4$ , $\unicode[STIX]{x03C0}/4$ and $3\unicode[STIX]{x03C0}/4$ corners respectively. Due to the asymmetry between the real and imaginary parts of the eigenfunctions $\hat{\unicode[STIX]{x1D74D}}_{1,2}$ and the corresponding projected patterns, the $\unicode[STIX]{x1D703}=3\unicode[STIX]{x03C0}/4$ and $7\unicode[STIX]{x03C0}/4$ impingement points are somewhat less sharply resolved than their $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/4$ and $5\unicode[STIX]{x03C0}/4$ counterparts. Intriguingly, at the $z=0$ cross-section, the vertical velocity-field components of the four clusters feature stronger activity adjacent to the sidewalls. In line with this observation, clusters A and B (C and D) appear to form ‘superclusters’, whose vertical velocity-field patterns are approximately relatable by $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}$ rotations.

Figure 6.

Reconstructed velocity-field snapshots $\hat{\boldsymbol{u}}_{n}^{primary}$ associated with the primary Koopman eigenfunctions (a–d, e–h and i–l) and the corresponding phase space coordinates $(y_{1},y_{2},y_{3})$ (m–p). The reconstructed velocity fields were computed via (3.30) and (4.2), and are displayed for four time instances $n\unicode[STIX]{x1D70F}$ , each associated with residence of the system state in one of the four phase space clusters. The cross-sections in (a–l) are as in figure 5. The dark solid lines in (m–p) indicate the phase space coordinates over the time interval $[n\unicode[STIX]{x1D70F},n\unicode[STIX]{x1D70F}+(Q-1)\unicode[STIX]{x1D70F}]$ employed for reconstruction. The snapshots at $n\unicode[STIX]{x1D70F}=200,2397,7595$ and 8661 exhibit the four possible diagonal LSC states associated with clusters A, B, C and D respectively (see figure 3).

Large-scale circulation states with the qualitative features described above have previously been found in laboratory experiments (Bai et al. Reference Bai, Ji and Brown2016), as well as in large-eddy simulations (Foroozani et al. Reference Foroozani, Niemela, Armenio and Sreenivasan2014, Reference Foroozani, Niemela, Armenio and Sreenivasan2017), of turbulent Rayleigh–Bénard convection in a cubical domain, both at higher Rayleigh numbers. In addition, the experiments were run over longer time intervals than our DNS study. Those studies found the flow to be trapped for long periods of the order of a few hundred free-fall times in one of the four diagonal flow states associated with the $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/4$ , $3\unicode[STIX]{x03C0}/4$ , $5\unicode[STIX]{x03C0}/4$ and $7\unicode[STIX]{x03C0}/4$ corners, and to subsequently switch to a rotationally adjacent macrostate. This behaviour was monitored in several different ways, e.g. by temperature probes at $z=0$ in the laboratory experiments (Bai et al. Reference Bai, Ji and Brown2016) or by the vertical velocity measured as in our figure 3 (Foroozani et al. Reference Foroozani, Niemela, Armenio and Sreenivasan2017). The latter study also found that during transitions between the diagonal macrostates, the LSC exhibits wall-aligned states. Such states were called transition states. In our reconstructions, we have also observed such wall-aligned states when the system transitions between clusters A and B, or C and D.

With regards to time scales, the angular frequency results in table 2 indicate a time scale of $T_{1,2}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{1,2}\approx 2\times 10^{4}$ for the expected time taken for the system to undergo a full cycle through all diagonal and wall-aligned LSC states, which is actually larger than our total observation interval. Even though this fact indicates that this $T_{1,2}$ value should be interpreted with caution (and should be thought of as an order of magnitude estimate of the time scale of the LSC cycling process), it suggests that a sampling of all LSC states would require analysis of an even longer DNS, which was not feasible with our available resources.

We also found that switching among the four stable diagonal LSC states is not equally probable in our analysis dataset. In particular, as can be seen in figure 3, the system exhibits periods where the LSC oscillates preferentially either between clusters A and B, or between clusters C and D, while transitions between A and D, or B and C, are significantly less frequent. This is noteworthy since A $\leftrightarrow$ B, B $\leftrightarrow$ C, C $\leftrightarrow$ D and D $\leftrightarrow$ A should be equiprobable transitions due to the square symmetry of horizontal cross-sections of the flow domain. The system thus appears to display a kind of broken ergodicity over time.

4.3 Statistical ensemble analysis

In order to analyse the asymmetry of residencies and transitions between clusters A, B, C and D more closely, we performed an additional series of 10 long-term simulations, each evolving for a time period of 5000 free-fall times in the statistically steady regime. Snapshots were written out again for each free-fall time unit. Each of the 10 cases started with a unique distribution of initial random perturbations of the equilibrium state, with the linear temperature profile and the fluid being at rest, and was given enough time to relax into a statistically steady state. The simulations were conducted at a smaller grid resolution of $64^{3}$ points.

Figure 7.

Time series of the midplane $(z=0)$ filtered vertical velocity $\tilde{w}_{x,y}$ from (4.1), sampled at four points near the corners of this cross-section (labelled #1, #3, #5, #7) and at four points near the midpoints of the sides (labelled #2, #4, #6, #8) for 10 different ensemble runs (shown in separate panels). Within each panel, pairs of time series are shifted vertically to improve visibility. Dotted lines mark $\tilde{w}_{x,y}=0$ axes. The numbering of the eight probe points is as in figure 3. The magnitude is again $|\tilde{w}_{x,y}|\leqslant 0.5$ . The time interval shown here is $[0,5000]$ in units of the free-fall time, which corresponds to $n\unicode[STIX]{x1D70F}\in \{0,\ldots ,5000\}$ .

Figure 7 displays again time series of the vertical velocity resulting from an analysis similar to figure 3. This simulation ensemble exhibits a variety of different flow regimes. In several of the 10 cases, longer time intervals are detected with dominating diagonal flow regimes. For example, in simulations 1, 2, 6, 7 and 9, one observes an LSC roll locked in the diagonal that connects points #1 and #5. In cases 1, 2, 3, 4, 8 and 9, the roll is locked in the other diagonal connecting points #3 and #7. Longer periods of low-amplitude oscillations at the diagonal corner points #1, #3, #5 and #7, as well as faster switches by $90^{\circ }$ from one diagonal state to another, are also observed in some of the cases. To summarize this analysis: we confirm that all possible LSC states can again be found; in addition, it is even possible to detect candidates for LSC reversals. During such events, all eight time series cross the zero value and the diagonal roll changes its spin, staying in the same pair of corners. This behaviour can be observed, for example, in case 3 at $t\approx 2900\unicode[STIX]{x1D70F}$ for corner points #1 and #5, or in case 4 at $t\approx 2500\unicode[STIX]{x1D70F}$ for corner points #3 and #7.

This ensemble analysis demonstrates that the long-term evolution may exhibit various different scenarios and switching patterns. The evolution gets locked in some cases in a particular regime for a longer time period. It is thus likely that our original long-term trajectory is not long enough to adequately sample all such regimes and transition states. This geometrically constrained convection flow thus obeys a kind of broken ergodicity, as is the case in other systems such as spin ensembles in condensed matter physics (Palmer Reference Palmer1982). It is beyond the scope of this work to address this problem in full completeness. We also wish to point once more to the experiments of Bai et al. (Reference Bai, Ji and Brown2016) and stress that they were actually run for 10 times longer than our DNS, namely 500 000 s, which corresponds to approximately 100 000 $T_{f}$ , all this at a Rayleigh number larger than ours by a factor of 50.

4.4 Secondary eigenfunctions

The three primary eigenfunctions are followed by a group of eight secondary eigenfunctions, $\unicode[STIX]{x1D74D}_{4},\ldots ,\unicode[STIX]{x1D74D}_{11}$ , which form four complex-conjugate pairs, all of which exhibit nearly equal Koopman frequencies $\unicode[STIX]{x1D714}_{k}$ , as shown in table 2. Figure 8 displays representative eigenfunction time series from each of the four pairs. These time series are compared with the time evolution of the primary eigenfunction $\unicode[STIX]{x1D74D}_{1}$ . All four secondary eigenfunction time series display time intervals with enhanced and strongly fluctuating amplitudes which are interrupted by periods of ambient low-amplitude variations.

The comparison with the primary mode indicates clearly that high-fluctuation intervals of each of the four pairs are synchronized with one of the four diagonal macrostates (see figure 3). For example, when the LSC flow is in cluster A (see also figure 4), the complex-conjugate pair $(\unicode[STIX]{x1D74D}_{10},\unicode[STIX]{x1D74D}_{11})$ fluctuates strongly in comparison with the other three pairs of secondary eigenfunctions. Similarly, the switch into macrostates D, C and B is in line with a switch to enhanced fluctuations of $(\unicode[STIX]{x1D74D}_{4},\unicode[STIX]{x1D74D}_{5})$ , $(\unicode[STIX]{x1D74D}_{6},\unicode[STIX]{x1D74D}_{7})$ and $(\unicode[STIX]{x1D74D}_{8},\unicode[STIX]{x1D74D}_{9})$ respectively. This suggests that the four pairs of secondary eigenfunctions could play an important role in exchanging kinetic energy between the four diagonal LSC configurations and degrees of freedom at smaller scales. If so, these patterns would be important building blocks for low-dimensional dynamical models of the LSC of this systems. The construction of such models is beyond the scope of this work.

Figure 8.

Modulating relationships between primary and secondary Koopman eigenfunctions. We display real parts only. (a) The first primary eigenfunction. Dashed horizontal lines indicate the four clusters (diagonal macrostates) which are also seen in the central trajectory plot in figure 4. Different diagonal macrostates in correspondence with figure 3 are indicated by different background colours (A $=$ green, B $=$ blue, C $=$ red, D $=$ yellow). Panels (b–e) display one eigenfunction from each of the four pairs of secondary eigenfunctions. In each of panels (b–e), we mark with coloured backgrounds intervals of residence in one of the macrostates that correspond to periods of enhanced fluctuations of the secondary eigenfunctions.

The typical time scale $T_{k}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{k}$ for the secondary eigenfunctions determined from their corresponding eigenfrequencies is approximately 40 free-fall time units. These quantities correspond to the time scale on which the secondary modes undergo oscillations in their active state. As with the primary eigenfunctions, we should interpret this time scale as an order of magnitude estimate, although it should be noted that in this case $T_{k}$ is significantly shorter than our analysis timespan and thus has higher robustness. The synchronization to the LSC states displayed in figure 8 suggests that $T_{k}$ is a mean turnover time of a fluid parcel in the LSC roll. It is this time scale that drives corner vortices and horizontal swirls which in turn feed the LSC or cause the crossover to another macrostate. The order of magnitude of $T_{k}$ agrees well with turnover times that have been found in a Lagrangian analysis of a turbulent convection flow in a cylindrical cell of aspect ratio 1 at the same Rayleigh number (Emran & Schumacher Reference Emran and Schumacher2010).

Figure 9 displays velocity-field reconstructions based on the secondary eigenfunctions at the same time instants as the data shown in figure 6. As in the case of the primary eigenfunctions, we display the vertical velocity component at two horizontal cross-sections and the $(x,z)$ components in a vertical cross-section. At each time instance, we show the reconstruction from those eigenfunctions that fluctuate most strongly when the corresponding macrostate is established (see also figure 8 for comparison). For example, figure 6(b,f,j,n) exhibits macrostate B, and this should be considered together with figure 9(b,f,j,n), showing reconstructions based on $\unicode[STIX]{x1D74D}_{8}$ and $\unicode[STIX]{x1D74D}_{9}$ . The phase portraits in figure 9(m–p) confirm this relationship. That is, the eigenfunction values participating in each reconstruction (indicated by solid dark lines) are always found at larger radii where the activated secondary modes oscillate. The phase portraits in figure 9(m–p) also demonstrate that the real and imaginary parts of the secondary eigenfunctions are to a good approximation in a $90^{\circ }$ phase difference, representing a coherent amplitude-modulated oscillator.

Figure 9.

Velocity-field reconstructions based on secondary Koopman eigenfunctions. The time instances and displayed fields in panels (a–l) are as in figure 6. In each column, we only use those eigenfunctions that fluctuate most strongly for the corresponding diagonal state. From left to right these are $(\unicode[STIX]{x1D74D}_{10},\unicode[STIX]{x1D74D}_{11})$ , in connection with macrostate A, $(\unicode[STIX]{x1D74D}_{8},\unicode[STIX]{x1D74D}_{9})$ , in connection with B, $(\unicode[STIX]{x1D74D}_{6},\unicode[STIX]{x1D74D}_{7})$ , in connection with C, and $(\unicode[STIX]{x1D74D}_{4},\unicode[STIX]{x1D74D}_{5})$ , in connection with D. Panels (m–p) show phase portraits of the secondary eigenfunctions in each case, with each time instance coloured in accordance with the primary cluster affiliation as in figure 3. The dark solid line corresponds to the eigenfunction values used for reconstruction via (3.30), indicating the strong activity of the secondary eigenfunction pair used for reconstruction in each case.

As shown in figure 9, the secondary eigenfunctions correspond to secondary structures in the convection flow which can appear in multiple configurations, including vortices that co- or counter-rotate with the LSC, corner vortices and prominent up- or downdraughts in the midplane or close to the plate. The vortices in the bulk and the corners can be interpreted as the flow structures that generate the oscillatory patterns in the time series of figure 8. They are potentially also responsible for the eventual switch of the LSC roll out of one of the four stable macrostates into another one. It should be noted also that the amplitudes of the primary eigenfunctions are not significantly larger (i.e. larger by more than an order of magnitude) than those of the secondary eigenfunctions. This confirms that it requires a strong perturbation to drive the LSC roll out of the quasistable equilibrium configuration into a different one.

4.5 Sensitivity analysis

In separate calculations, we have verified that the properties of the eigenfunctions listed in table 2 are robust under changes of NLSA and Koopman parameters. Among these parameters, a particularly important one is the number of delays $Q$ , as it controls the ability of the NLSA eigenfunctions to span the discrete subspace ${\mathcal{D}}$ of the Koopman operator (see § 3.2). Here, we have computed NLSA and Koopman eigenfunctions using $Q=15$ and $60$ delays, and found that the recovered primary and secondary eigenfunction time series remain qualitatively unchanged. The frequencies $\unicode[STIX]{x1D714}_{k}$ of the secondary eigenfunctions typically changed by less than 5 % in these calculations. As expected from the fact that they lie close to zero (in comparison with the frequency resolution $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}$ afforded by our dataset), the frequencies of the primary eigenfunctions exhibited significantly larger, 50 %–75 %, changes, but still captured an $O(10^{-4})$ frequency scale. The primary frequencies were also more sensitive to the diffusion regularization parameter $\unicode[STIX]{x1D701}$ (probably for the same reasons). In particular, we found that an increase of $\unicode[STIX]{x1D701}$ to $\simeq 10^{-3}$ caused the primary eigenfunctions to become purely diffusive, with $\unicode[STIX]{x1D714}_{k}$ collapsing to zero. While the LSC states discussed in § 4.2 are still captured by these diffusive modes, the ability to separate the two pairs of diagonal states into distinct modes is somewhat degraded (i.e. the diagonal states in figure 5 become mixed). Similarly, we found that the primary eigenfunctions are more sensitive to the effects of diffusion if the finite-difference scheme is used instead of the logarithm scheme to approximate the Koopman generator $V$ (see § 3.3). As a final sensitivity test, we have examined the dependence of our results on the dimension $\ell$ of the Galerkin approximation space, and found that satisfactory primary eigenfunctions can be obtained with $\ell$ as low as 100. Secondary eigenfunctions capturing similar time scales to those listed in table 2 can also be obtained using $\ell \simeq 100$ , but the modulating relationships between primary and secondary eigenfunctions discussed in § 4.4 were found to require larger $\ell$ .

Overall, the observations described above highlight the fact that successful data-driven recovery of LSC patterns in turbulent convective flows via Koopman operator techniques depends strongly on (i) the basis functions $\unicode[STIX]{x1D753}_{k}$ used for the Galerkin approximation (here, depending on the number of delays $Q$ ), (ii) the properties of diffusion regularization (quantified by $\unicode[STIX]{x1D701}$ ) and (iii) the approximation scheme for the generator $V$ .

It should be noted that, in addition to the three points listed above, the approximation accuracy of our scheme also depends on the sampling interval $\unicode[STIX]{x1D70F}$ . Due to the high computational cost of the NLSA basis calculation (which would have to be repeated if $\unicode[STIX]{x1D70F}$ were to be changed), we have not performed a detailed sensitivity analysis for this parameter. In general, for a fixed number of analysis samples $N$ , which may be limited by the available computational resources, one seeks a compromise between approximation accuracy for the generator (which increases with decreasing $\unicode[STIX]{x1D70F}$ ) and accuracy of ergodic time averaging (which increases with increasing $N\unicode[STIX]{x1D70F}$ ). Our choice $\unicode[STIX]{x1D70F}=1T_{f}$ reflects such a compromise for patterns evolving on time scales between $O(10T_{f})$ (secondary modes) and $O(10^{3}T_{f})$ (primary modes). An alternative numerical approach, lying beyond the scope of the present work, could be to maintain a $\unicode[STIX]{x1D70F}=O(1)$ sampling interval in the computation of the NLSA basis, but use a smaller interval $\unicode[STIX]{x1D70F}^{\prime }$ (perhaps as small as a simulation time step) to locally approximate the action of the generator on the NLSA basis functions. This would require the evaluation of NLSA eigenfunctions at states $x_{n}^{\prime }$ sampled at times $t_{n}^{\prime }=n(\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D70F}^{\prime })$ , which can be performed efficiently using out-of-sample extension techniques for Laplacian eigenfunctions (Coifman & Lafon Reference Coifman and Lafon2006b ).